📄 ximadsp.cpp
字号:
for (j=0;j<h;j++) {
for (k=0;k<w;k++) {
real[k] = grid[k][j].x;
imag[k] = grid[k][j].y;
}
if (bXpow2) FFT(direction,m,real,imag);
else DFT(direction,w,real,imag,real2,imag2);
for (k=0;k<w;k++) {
grid[k][j].x = real[k];
grid[k][j].y = imag[k];
}
}
free(real);
free(imag);
/* Transform the columns */
real = (double *)malloc(h * sizeof(double));
imag = (double *)malloc(h * sizeof(double));
m=0;
while((1<<m)<h) m++;
for (k=0;k<w;k++) {
for (j=0;j<h;j++) {
real[j] = grid[k][j].x;
imag[j] = grid[k][j].y;
}
if (bYpow2) FFT(direction,m,real,imag);
else DFT(direction,h,real,imag,real2,imag2);
for (j=0;j<h;j++) {
grid[k][j].x = real[j];
grid[k][j].y = imag[j];
}
}
free(real);
free(imag);
free(real2);
free(imag2);
/* converting from double to byte, there is a HUGE loss in the dynamics
"nn" tries to keep an acceptable SNR, but 8bit=48dB: don't ask more */
double nn=pow((double)2,(double)log((double)max(w,h))/(double)log((double)2)-4);
//reversed gain for reversed transform
if (direction==-1) nn=1/nn;
//bMagnitude : just to see it on the screen
if (bMagnitude) nn*=4;
for (j=0;j<h;j++) {
for (k=0;k<w;k++) {
if (bMagnitude){
tmpReal->SetPixelIndex(k,j,(BYTE)max(0,min(255,(nn*(3+log(_cabs(grid[k][j])))))));
if (grid[k][j].x==0){
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128+(atan(grid[k][j].y/0.0000000001)*nn)))));
} else {
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128+(atan(grid[k][j].y/grid[k][j].x)*nn)))));
}
} else {
tmpReal->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128 + grid[k][j].x*nn))));
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128 + grid[k][j].y*nn))));
}
}
}
for (k=0;k<w;k++) free (grid[k]);
free (grid);
if (srcReal==0 && dstReal==0) delete tmpReal;
if (srcImag==0 && dstImag==0) delete tmpImag;
return true;
}
////////////////////////////////////////////////////////////////////////////////
bool CxImage::IsPowerof2(long x)
{
long i=0;
while ((1<<i)<x) i++;
if (x==(1<<i)) return true;
return false;
}
////////////////////////////////////////////////////////////////////////////////
/**
This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays of n=2^m points.
o(n)=n*log2(n)
dir = 1 gives forward transform
dir = -1 gives reverse transform
Written by Paul Bourke, July 1998
FFT algorithm by Cooley and Tukey, 1965
*/
bool CxImage::FFT(int dir,int m,double *x,double *y)
{
long nn,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;
/* Calculate the number of points */
nn = 1<<m;
/* Do the bit reversal */
i2 = nn >> 1;
j = 0;
for (i=0;i<nn-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<nn;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<nn;i++) {
x[i] /= (double)nn;
y[i] /= (double)nn;
}
}
return true;
}
////////////////////////////////////////////////////////////////////////////////
/**
Direct fourier transform o(n)=n^2
Written by Paul Bourke, July 1998
*/
bool CxImage::DFT(int dir,long m,double *x1,double *y1,double *x2,double *y2)
{
long i,k;
double arg;
double cosarg,sinarg;
for (i=0;i<m;i++) {
x2[i] = 0;
y2[i] = 0;
arg = - dir * 2.0 * PI * i / (double)m;
for (k=0;k<m;k++) {
cosarg = cos(k * arg);
sinarg = sin(k * arg);
x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
}
}
/* Copy the data back */
if (dir == 1) {
for (i=0;i<m;i++) {
x1[i] = x2[i] / m;
y1[i] = y2[i] / m;
}
} else {
for (i=0;i<m;i++) {
x1[i] = x2[i];
y1[i] = y2[i];
}
}
return true;
}
////////////////////////////////////////////////////////////////////////////////
/**
* Combines different color components into a single image
* \param r,g,b: color channels
* \param a: alpha layer, can be NULL
* \param colorspace: 0 = RGB, 1 = HSL, 2 = YUV, 3 = YIQ, 4 = XYZ
* \return true if everything is ok
*/
bool CxImage::Combine(CxImage* r,CxImage* g,CxImage* b,CxImage* a, long colorspace)
{
if (r==0 || g==0 || b==0) return false;
long w = r->GetWidth();
long h = r->GetHeight();
Create(w,h,24);
g->Resample(w,h);
b->Resample(w,h);
if (a) {
a->Resample(w,h);
#if CXIMAGE_SUPPORT_ALPHA
AlphaCreate();
#endif //CXIMAGE_SUPPORT_ALPHA
}
RGBQUAD c;
for (long y=0;y<h;y++){
info.nProgress = (long)(100*y/h); //<Anatoly Ivasyuk>
for (long x=0;x<w;x++){
c.rgbRed=r->GetPixelIndex(x,y);
c.rgbGreen=g->GetPixelIndex(x,y);
c.rgbBlue=b->GetPixelIndex(x,y);
switch (colorspace){
case 1:
SetPixelColor(x,y,HSLtoRGB(c));
break;
case 2:
SetPixelColor(x,y,YUVtoRGB(c));
break;
case 3:
SetPixelColor(x,y,YIQtoRGB(c));
break;
case 4:
SetPixelColor(x,y,XYZtoRGB(c));
break;
default:
SetPixelColor(x,y,c);
}
#if CXIMAGE_SUPPORT_ALPHA
if (a) AlphaSet(x,y,a->GetPixelIndex(x,y));
#endif //CXIMAGE_SUPPORT_ALPHA
}
}
return true;
}
////////////////////////////////////////////////////////////////////////////////
/**
* Smart blurring to remove small defects, dithering or artifacts.
* \param radius: normally between 0.01 and 0.5
* \param niterations: should be trimmed with radius, to avoid blurring should be (radius*niterations)<1
* \param colorspace: 0 = RGB, 1 = HSL, 2 = YUV, 3 = YIQ, 4 = XYZ
* \return true if everything is ok
*/
bool CxImage::Repair(float radius, long niterations, long colorspace)
{
if (!IsValid()) return false;
long w = GetWidth();
long h = GetHeight();
CxImage r,g,b;
r.Create(w,h,8);
g.Create(w,h,8);
b.Create(w,h,8);
switch (colorspace){
case 1:
SplitHSL(&r,&g,&b);
break;
case 2:
SplitYUV(&r,&g,&b);
break;
case 3:
SplitYIQ(&r,&g,&b);
break;
case 4:
SplitXYZ(&r,&g,&b);
break;
default:
SplitRGB(&r,&g,&b);
}
for (int i=0; i<niterations; i++){
RepairChannel(&r,radius);
RepairChannel(&g,radius);
RepairChannel(&b,radius);
}
CxImage* a=NULL;
#if CXIMAGE_SUPPORT_ALPHA
if (AlphaIsValid()){
a = new CxImage();
AlphaSplit(a);
}
#endif
Combine(&r,&g,&b,a,colorspace);
delete a;
return true;
}
////////////////////////////////////////////////////////////////////////////////
bool CxImage::RepairChannel(CxImage *ch, float radius)
{
if (ch==NULL) return false;
CxImage tmp(*ch);
if (!tmp.IsValid()) return false;
long w = ch->GetWidth()-1;
long h = ch->GetHeight()-1;
double correction,ix,iy,ixx,ixy,iyy,den,num;
int x,y,xy0,xp1,xm1,yp1,ym1;
for(x=1; x<w; x++){
for(y=1; y<h; y++){
xy0 = ch->GetPixelIndex(x,y);
xm1 = ch->GetPixelIndex(x-1,y);
xp1 = ch->GetPixelIndex(x+1,y);
ym1 = ch->GetPixelIndex(x,y-1);
yp1 = ch->GetPixelIndex(x,y+1);
ix= (xp1-xm1)/2.0;
iy= (yp1-ym1)/2.0;
ixx= xp1 - 2.0 * xy0 + xm1;
iyy= yp1 - 2.0 * xy0 + ym1;
ixy=(ch->GetPixelIndex(x+1,y+1)+ch->GetPixelIndex(x-1,y-1)-
ch->GetPixelIndex(x-1,y+1)-ch->GetPixelIndex(x+1,y-1))/4.0;
num= (1.0+iy*iy)*ixx - ix*iy*ixy + (1.0+ix*ix)*iyy;
den= 1.0+ix*ix+iy*iy;
correction = num/den;
tmp.SetPixelIndex(x,y,(BYTE)min(255,max(0,(xy0 + radius * correction))));
}
}
for (x=0;x<=w;x++){
tmp.SetPixelIndex(x,0,ch->GetPixelIndex(x,0));
tmp.SetPixelIndex(x,h,ch->GetPixelIndex(x,h));
}
for (y=0;y<=h;y++){
tmp.SetPixelIndex(0,y,ch->GetPixelIndex(0,y));
tmp.SetPixelIndex(w,y,ch->GetPixelIndex(w,y));
}
ch->Transfer(tmp);
return true;
}
////////////////////////////////////////////////////////////////////////////////
/**
* Enhance the variations between adjacent pixels.
* Similar results can be achieved using Filter(),
* but the algorithms are different both in Edge() and in Contour().
* \return true if everything is ok
*/
bool CxImage::Contour()
{
if (!pDib) return false;
long Ksize = 3;
long k2 = Ksize/2;
long kmax= Ksize-k2;
long i,j,k;
BYTE maxr,maxg,maxb;
RGBQUAD pix1,pix2;
CxImage tmp(*this,pSelection!=0,true,true);
if (!tmp.IsValid()) return false;
long xmin,xmax,ymin,ymax;
if (pSelection){
xmin = info.rSelectionBox.left; xmax = info.rSelectionBox.right;
ymin = info.rSelectionBox.bottom; ymax = info.rSelectionBox.top;
} else {
xmin = ymin = 0;
xmax = head.biWidth; ymax=head.biHeight;
}
for(long y=ymin; y<ymax; y++){
info.nProgress = (long)(100*y/head.biHeight);
if (info.nEscape) break;
for(long x=xmin; x<xmax; x++){
#if CXIMAGE_SUPPORT_SELECTION
if (SelectionIsInside(x,y))
#endif //CXIMAGE_SUPPORT_SELECTION
{
pix1 = GetPixelColor(x,y);
maxr=maxg=maxb=0;
for(j=-k2, i=0;j<kmax;j++){
for(k=-k2;k<kmax;k++, i++){
pix2=GetPixelColor(x+j,y+k);
if ((pix2.rgbBlue-pix1.rgbBlue)>maxb) maxb = pix2.rgbBlue;
if ((pix2.rgbGreen-pix1.rgbGreen)>maxg) maxg = pix2.rgbGreen;
if ((pix2.rgbRed-pix1.rgbRed)>maxr) maxr = pix2.rgbRed;
}
}
pix1.rgbBlue=(BYTE)(255-maxb);
pix1.rgbGreen=(BYTE)(255-maxg);
pix1.rgbRed=(BYTE)(255-maxr);
tmp.SetPixelColor(x,y,pix1);
}
}
}
Transfer(tmp);
return true;
}
////////////////////////////////////////////////////////////////////////////////
/**
* Adds a random offset to each pixel in the image
* \param radius: maximum pixel displacement
* \return true if everything is ok
*/
bool CxImage::Jitter(long radius)
{
if (!pDib) return false;
long nx,ny;
CxImage tmp(*this,pSelection!=0,true,true);
if (!tmp.IsValid()) return false;
long xmin,xmax,ymin,ymax;
if (pSelection){
xmin = info.rSelectionBox.left; xmax = info.rSelectionBox.right;
ymin = info.rSelectionBox.bottom; ymax = info.rSelectionBox.top;
} else {
xmin = ymin = 0;
xmax = head.biWidth; ymax=head.biHeight;
}
for(long y=ymin; y<ymax; y++){
info.nProgress = (long)(100*y/head.biHeight);
if (info.nEscape) break;
for(long x=xmin; x<xmax; x++){
#if CXIMAGE_SUPPORT_SELECTION
if (SelectionIsInside(x,y))
#endif //CXIMAGE_SUPPORT_SELECTION
{
nx=x+(long)((rand()/(float)RAND_MAX - 0.5)*(radius*2));
ny=y+(long)((rand()/(float)RAND_MAX - 0.5)*(radius*2));
if (!IsInside(nx,ny)) {
nx=x;
ny=y;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -