二维变换算法
⑴ Radon变换及逆变换公式的来历及证明过程,要详细的,不懂得不要回答
radon变换 两维情况下radon变换大致可以这样理解:一个平面内沿不同的直线(直线与原点的距离为d,方向角为alfa)对f(x,y)做线积分,得到的像F(d,alfa)就是函数f的Radon变换。也就是说,平面(d,alfa)的每个点的像函数值对应了原始函数的某个线积分值。一个更直观的理解是,假设你的手指被一个很强的平行光源透射,你迎着光源看到的手指图像就是手指的光衰减系塌尘数的三维Radon变换(小小的推广)在给定方向(两个角坐标)的时候的值, 一个最简单而直接的应用就是拿来检测图像里面含有的直线成分,很显然地,任何直线都会导致Randon像在该直线对应(d,alfa)处的极值。 具体的CT断层团逗禅影像重建算法当中其实没怎么用到Radon变换,或者说Radon变换仅仅只有一点点理论上的意义。原因是:CT机做扫描:球管发出X-ray,经过人体,被吸收一部分,进入检测器队列(球管是旋转的,检测器呈扇形分布,很老的和很新的除外,老式的ct做平行扫描,效率低,很新式的什么多层螺旋扫描,我也不知道咋回事)显然检测器读数就是人体的x-ray吸收系数(是空间的函数)对相应路径的线积分,所以这样转一圈下来再把所有的检测器读数值按照(d,alfa)的方式排列一下就算完成了某个被检测截面的Radon变换了,这个过程是人体和X-rayscaner一起完成的,显然不干软件什么事。接下来,照理说是要靠计算机把获得的数据做一个逆Radon变换,就能得到被检测截面的X-ray吸收系数的分布图像了。CT的图像其实就是一个吸收系数的图,类似的B超或者声纳之类的图像是大致是一个弹性模量的图(反射声波)... 但是接下来这里有一个问题就是Radon变换是不是可逆,google了下好像是可逆的,我的理解: 1)有另外一种求逆方法,就是解代数方程,简化地说可以大致设想把整个截面离散网格化,每个格子对应一个吸收系数,把每根扫描积分路径经过的格子按照权重(显然透心凉和擦点皮对吸收的贡献不同)作累加,令他们等于相应Radon值(积分变成了加权累加而已)显然设计好的话,这个方程组肯定是有解的(不过运算量会很庞大,比如一个512X512的网格...) 2)工程师不问这么无聊不切实际的问题,所以以前学的时候就压根没想到。 3)最重要的原因,是下面要说的求逆问题,竟然变成了二维的fourier逆变换。所以忘掉Radon变换吧。 有这样一个事实:把某个角度坐标alfa对应的一“条”Radon值(一系列检测器的读数,也实际上就是原始截面吸收系数在方向为alfa+-Pi/2直线簇上的线积分值)作一个fourier变换,得到的就是整个原始被检测截面(吸收系数)的二维fourier像在某条直线上的值(这条直线经过频域的原点并且方向为alfa)如果把所有角度的Radon值作一维Fourier变换,然后按照合适的角度(alfa)经过原点把这些一维fourier像值放在频域平面上,就得到了整个二维fourier像!!!这个其实直观上很容易想象其合理性,还是以手指头为例(不考虑它指向的方向)对着光源看,从左至右,透光率不同产生明暗的变化,亮暗本身是沿前后方向的积分结果决定,但是相邻的亮暗变化却反应了整个手指截面的从左至右这个方向上的频域信息,看到的细节越多,频域的高频分量越多(与前后方向的细节毫无关系,因为被radon积分掉了)。 以上关于CT其实是过分简化的描述,只能提供一个大致的原理。实际情况会有些不同,首先检测器读数是有限空间的,这就是相当于理想的投影函数乘了一个窗函数(某段区间内为一,其他地方为零的函数),在频域内窗函数会“扩散”所以他们频域做卷积的结果是频域的扩展。也可以说成是,对于非周期函数(包括周期不为无穷大)的fourier级数在边界的间断处只能是平均收敛,“平均”的结果就是在光滑的地方拟合的很好,在间断点处发生振荡。工程中管这个叫做吉布斯(Gibbs)效应,它告诉我们:用有限项指消级数的和去表示一个函数,随着项数的增加,振荡发生的位置会越来越接近间断点,但是它的摆幅不变(写到这忽然觉得它的名字似乎翻译成“挤不死”更贴切)另外,检测器只能读出空间上分立的数值,所谓的取样过程就是投影函数乘一个迪拉克函数组成的序列(假设周期为L)而迪拉克序列变换到频域仍然是一个迪拉克序列,只是周期变成了1/L。投影和取样序列相乘在频域就是卷积,出来的结果就是具有了周期频谱,显然可用的只能是原点(DC)所在的一个周期内的数据。当L越来越小的时候,频谱周期越来越大,空间分辨率越来越高。当L为有限的时候,分辨率如果用频率来表示的话,从原点(“直流”分量)开始算,由于周期性缘故显然最高到1/2L处。 设想一间黑屋子,唯一的光源是一个可调节频率的频闪光源,一台电风扇。假定光源闪烁频率为w,显然理论上能够检测到的风扇转速u将允许加上任意整数个w。比方说,每秒亮一下,你看到了风扇转动了1/4圈,那么你可以认为风扇每秒转动1/4圈,但也可以是5/4圈(多转了一圈,有何不可?),9/4圈...也可以是(反着转)-3/4圈,-7/4圈...原因就是前面说的,用一个脉冲序列(光源频闪)去做取样,必然会得到周期性的频谱。接下来,当光源的闪烁频率和风扇的转速(用转/秒来表述)相等的时候,你将看到风扇是停止的,当光源频率高于风扇转速的两倍时,你才能看到风扇正常的转动,如果光源频率介于风扇转速一倍和两倍之间,那你会看到风扇倒着转了。这里的情况被称为频谱混叠。此类现象生活中常遇到。另外,函数变换本身还带来了坐标平移一类的问题。实际当中用的最多的是一种叫做滤波反投影的算法来实现断层重建,说穿了关键就针对以上一些问题设计合理的滤波器。 另外值得一提的是,这里用到的数学大概一百年前就有了,但是随着计算机技术的进步,具体实现的时候,出现过不同的版本。譬如说,70年代的商业运行的CT(256X256),带一台长得像电冰箱般的“卷积器(convolver)”,顾名思义,它的滤波器实现大概是用DSP做卷积的(离散的卷积就是一系列的移位连乘连加)。而现代,随着硬件技能的突飞猛进,FFT不成问题了,这个交给CPU在频域内作乘法就能搞定。退一步说,我甚至怀疑,那个形体巨大的Convolver做卷积的性能恐怕未必能赶上我正在码字的电脑。此刻,它正在运行音乐播放软件foorbar,同时一起运行的还有一堆插件(也可看作卷积器),比如老式电子管音色,教堂混响,耳机模拟现场音效之类的... 以上这些基本上是相关领域的abc,没有深入,基本凭借记忆,说法可能和专门的教材不完全一样,而且很多地方一知半解,肯定会有谬误,大家随便看看不可当真,当然欢迎拍砖。
⑵ 图像算法中的二维傅立叶变换(DCT)及蝶形算法(butterfly algorithm)
本文讨论图像算法中的核心概念,即二维傅立叶变换(DCT)与蝶形算法(butterfly algorithm)。首先,我们概述了离散2维Fourier的基础知识,虽然这里不进行深入讲解,但有需要时可参考fourier.eng.hmc.e/e10...的资料。
在图像处理和压缩中,如JPEG等,傅立叶变换发挥着关键作用。我们通常不会对整个图像进行傅立叶变换,而是将其分割成小块,如4x4或8x8像素块,以达到计算效率与质量平衡的目的。
接下来,以4x4像素块为例,解析DCT变换及蝶形算法的原理。首先,我们定义了4x4矩阵的DCT变换矩阵。通过推导,我们发现变换矩阵为[公式]。在实际应用中,文献通常采用简化形式,其矩阵为[公式],其中系数C影响直流系数的大小,对图像亮度处理尤为重要。尽管C的引入改变了变换矩阵的形式,但它不会对图像处理结果造成实质性的影响。
基于此,我们直接使用简化后的DCT变换矩阵,即[公式],进行图像变换。随后,我们引入蝶形算法,通过合并重复的计算组合,大大降低运算量。以4x4矩阵为例,蝶形算法图示展示了如何合并计算路径,形成依赖关系,从而减少重复计算,提高计算效率。
总结,二维傅立叶变换(DCT)与蝶形算法(butterfly algorithm)在图像算法中扮演关键角色,通过分割图像、变换矩阵的简化以及蝶形算法的应用,实现高效且精确的图像处理与压缩。
⑶ 如何用findhomography计算二维仿射变换矩阵
1.计算方法不同:通过跟踪源码,发现getPerspectiveTransform用的是SVD分解,findHomography看不出是用什么方法(没注释枣迅,一堆等式)。但两者计算结果是一样的。
2.输入参数不同:getPerspectiveTransform只会拿前4个点去计算,findHomography则会拿一堆点(>=4)去计算(其是不断从一堆点中重复拿出4个点去计算出一个结果,再采用一些优化算法凳搏此RANSAC/银神LMEDS去筛选出最优解)。
⑷ 怎样用C语言实现FFT算法啊
1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下:
先对各行逐一进行一维FFT,然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示:
for (int i=0; i<M; i++)
FFT_1D(ROW[i],N);
for (int j=0; j<N; j++)
FFT_1D(COL[j],M);
其中,ROW[i]表示矩阵的第i行。注意这只是一个简单的记法,并不能完全照抄。还需要通过一些语句来生成各行的数据。同理,COL[i]是对矩阵的第i列的一种简单表示方法。
所以,关键是一维FFT算法的实现。
2、例程:
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineN1000
/*定义复数类型*/
typedefstruct{
doublereal;
doubleimg;
}complex;
complexx[N],*W;/*输入序列,变换核*/
intsize_x=0;/*输入序列的大小,在本程序中仅限2的次幂*/
doublePI;/*圆周率*/
voidfft();/*快速傅里叶变换*/
voidinitW();/*初始化变换核*/
voidchange();/*变址*/
voidadd(complex,complex,complex*);/*复数加法*/
voidmul(complex,complex,complex*);/*复数乘法*/
voidsub(complex,complex,complex*);/*复数减法*/
voidoutput();
intmain(){
inti;/*输出结果*/
system("cls");
PI=atan(1)*4;
printf("Pleaseinputthesizeofx: ");
scanf("%d",&size_x);
printf("Pleaseinputthedatainx[N]: ");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return0;
}
/*快速傅里叶变换*/
voidfft(){
inti=0,j=0,k=0,l=0;
complexup,down,proct;
change();
for(i=0;i<log(size_x)/log(2);i++){/*一级蝶形运算*/
l=1<<i;
for(j=0;j<size_x;j+=2*l){/*一组蝶形运算*/
for(k=0;k<l;k++){/*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],&proct);
add(x[j+k],proct,&up);
sub(x[j+k],proct,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化变换核*/
voidinitW(){
inti;
W=(complex*)malloc(sizeof(complex)*size_x);
for(i=0;i<size_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*变址计算,将x(n)码位倒置*/
voidchange(){
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for(i=0;i<size_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while((t--)>0){
j=j<<1;
j|=(k&1);
k=k>>1;
}
if(j>i){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*输出傅里叶变换的结果*/
voidoutput(){
inti;
printf("Theresultareasfollows ");
for(i=0;i<size_x;i++){
printf("%.4f",x[i].real);
if(x[i].img>=0.0001)printf("+%.4fj ",x[i].img);
elseif(fabs(x[i].img)<0.0001)printf(" ");
elseprintf("%.4fj ",x[i].img);
}
}
voidadd(complexa,complexb,complex*c){
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul(complexa,complexb,complex*c){
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub(complexa,complexb,complex*c){
c->real=a.real-b.real;
c->img=a.img-b.img;
}
⑸ 计算机图形学直线段裁剪算法或二维基本变换算法,能复制运行的来啊,谢谢大侠们了 啊
#include"graphics.h"
#include"math.h"
typedef struct Point /* 点 */
{
int x;
int y;
}Point;
/* 点的平移变换 */
void PinYi(int * x,int * y,int tx,int ty)
{
*x = *x + tx;
*y = *y + ty;
}
/* 点的旋转变换 */
void XuanZhuan(int * x,int * y,float q)
{
int m ;
int n;
float Q;
Q = (3.1415926/180)*q;
m = (*x);
n = (*y);
*x = m * cos(Q) - n * sin(Q);
*y = m * sin(Q) + n * cos(Q);
}
void XuanZhuan_RY(int m,int n,int * x,int * y,float q)/* 绕(m,n)的旋转*/
{
PinYi(x,y,-m,-n);
XuanZhuan(x,y,q);
PinYi(x,y,m,n);
}
void line_dda(int x1,int y1,int x2,int y2 ,COLORREF color)
{
float x = 0.0;
float y = 0.0;
float x3 = ( float )x1;
float y3 = ( float )y1;
float n;
n=( float )( abs( x2-x1 ) >= abs( y2-y1 ) ? abs( x2-x1 ) : abs( y2-y1 ) );
if( n != 0.0 )
{
x=( ( float )( x2-x1 ) ) / n;
y=( ( float )( y2-y1 ) ) / n;
}
while( n >= 0 )
{
putpixel( ( int )x3 , ( int )y3 , color );
x3 += x ;y3 += y ;n -= 1.0;
}
}
void line(Point i,Point j,COLORREF color){
line_dda(i.x,i.y,j.x,j.y,color);
}
Point p[4]={-50,50,50,50,-50,-50,50,-50};
int PY[2]={0,0},XZ=0;
void main()
{
int width=600,height=480;
int zbx=200,zby=200;
COLORREF color=0x00ff00;
char ch;
initgraph(width, height);
ch=getch();
for(int o=0;o<4;o++)
PinYi(&p[o].x,&p[o].y,zbx,zby);
while(ch=getch())
{
PY[0]=0;PY[1]=0;
XZ=0;
switch(ch){
case 'q': XZ--;break;
case 'e': XZ++;break;
case 'w': PY[1]--;break;
case 's': PY[1]++;break;
case 'a': PY[0]--;break;
case 'd': PY[0]++;break;
//default :return 0;
}
for(int o=0;o<4;o++)
{
PinYi(&p[o].x,&p[o].y,PY[0],PY[1]);
PinYi(&zbx,&zby,PY[0],PY[1]);
XuanZhuan_RY(zbx,zby,&p[o].x,&p[o].y,XZ);
}
cleardevice();
line(p[0],p[1],color);
line(p[1],p[2],color);
line(p[2],p[3],color);
line(p[3],p[0],color);
}
}
有个问题是,,旋转的时候因为π取的3.1415926七位,所以越转越小!!!
需要改动!!就不弄了,要睡了!!
⑹ 试讨论一维DCT、二维DCT变换不采用快速算法时所需的加法、乘法次数。
【答案】:(1)一维DCT:
Y=AX
对每一个变换系数Y(K),需作N次乘法,N-1次加法。共有N个变换系数,需N2次乘法,N(N-1)次加法。
(2)二维DCT:
F=[DCT],[DCT]T
对每一个变换系数Y(u,v),需作2次矩阵的相乘。每次相乘需N次乘法,N-1次加法,所以要得到Y(u,v),需2N次乘法,2(N-1)次加法。共有N2个变换系数(N×N像素块),因此共需2N3次乘法,2N2(N-1)次加法。