00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include <math.h>
00029
00030 #include "xsmmath.h"
00031 #include "xsmtypes.h"
00032 #include "glbvars.h"
00033
00034 #include "bench.h"
00035 #include "regress.h"
00036
00037
00038 char *MathErrString[] = {
00039 "Math OK",
00040 "Size of sources is not the same !",
00041 "Size of selection is invalid !",
00042 "Division by zero occured !",
00043 "undefined error"
00044 };
00045
00046
00047 void MkMausSelect(Point2D *List, MOUSERECT *msel, int mx, int my){
00048 msel->xLeft = MAX(0, MIN( MIN(List[0].x, List[1].x), mx));
00049 msel->xRight = MAX(0, MIN( MAX(List[0].x, List[1].x), mx));
00050 msel->xSize = msel->xRight - msel->xLeft;
00051 msel->xRatio = (double)msel->xSize / (double)mx;
00052
00053 msel->yTop = MAX(0, MIN( MIN(List[0].y, List[1].y), my));
00054 msel->yBottom = MAX(0, MIN( MAX(List[0].y, List[1].y), my));
00055 msel->ySize = msel->yBottom - msel->yTop;
00056 msel->yRatio = (double)msel->ySize / (double)my;
00057
00058 if(msel->ySize > 0)
00059 msel->Aspect = (double)msel->xSize / (double)msel->ySize;
00060 else
00061 msel->Aspect = 0.;
00062
00063 msel->Area = msel->xSize * msel->ySize;
00064 msel->Radius2 = (msel->xSize)*(msel->xSize) + (msel->ySize)*(msel->ySize);
00065 msel->xCenter = List[0].x;
00066 msel->yCenter = List[0].y;
00067
00068 XSM_DEBUG (DBG_L3, "MkMausSelect: (" << msel->xLeft << ", " << msel->yBottom << ")-(" << msel->xRight << ", " << msel->yTop << ")");
00069 XSM_DEBUG (DBG_L3, "MkMausSelect Size: (" << msel->xSize << ", " << msel->ySize << ")");
00070 }
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 gboolean CopyScan(MATHOPPARAMS){
00145 XSM_DEBUG (DBG_L3, "Copy Scan");
00146 Dest->mem2d->Resize(Src->mem2d->GetNx(), Src->mem2d->GetNy());
00147
00148 BenchStart(mbs,"Copy Scan",Dest->mem2d->GetNx()*Dest->mem2d->GetNy());
00149
00150 Dest->mem2d->ConvertFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00151
00152 BenchStop(mbs);
00153
00154 return MATH_OK;
00155 }
00156
00157
00158 gboolean CropScan(MATHOPPARAMS){
00159 MOUSERECT msr;
00160 XSM_DEBUG (DBG_L3, "Crop Scan");
00161
00162 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00163
00164 if( msr.xSize < 1 || msr.ySize < 1){
00165 XSM_DEBUG (DBG_L3, "Crop:" << msr.xSize << " " << msr.ySize);
00166 return MATH_SELECTIONERR;
00167 }
00168
00169 if(gapp->xsm->MausMode() == MRECT){
00170 Dest->data.s.nx = msr.xSize;
00171 Dest->data.s.ny = msr.ySize;
00172 Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
00173 Dest->data.s.x0 += (msr.xLeft + msr.xSize/2 - Src->data.s.nx/2)*Src->data.s.dx;
00174 Dest->data.s.y0 -= msr.yTop*Src->data.s.dy;
00175
00176
00177 Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
00178 Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;
00179
00180 Dest->mem2d->CopyFrom(Src->mem2d, msr.xLeft,msr.yTop, 0,0, Dest->data.s.nx, Dest->data.s.ny);
00181
00182 return MATH_OK;
00183 }
00184 else if(gapp->xsm->MausMode() == MCIRC){
00185 int Radius, xx,x,y,r2,y2;
00186 int col;
00187 const SHT bgval=(SHT)gapp->xsm->mradius;
00188 Dest->data.s.nx = 2*(Radius = (int)sqrt((double)msr.Radius2));
00189 Dest->data.s.ny = 2*Radius;
00190 Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
00191 r2 = Radius*Radius;
00192
00193
00194 Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
00195 Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;
00196
00197 for(int line = 0; line < Dest->data.s.ny; line++){
00198 for(col = 0; col < Dest->data.s.nx; col++)
00199 Dest->mem2d->PutDataPkt(bgval, col, line);
00200
00201 y=(int)(line+msr.yCenter-Radius);
00202 y2 = line-Radius;
00203 y2 *= y2;
00204 if(y>0 && y<Src->mem2d->GetNy())
00205 for(col = 0; col < Dest->data.s.nx; col++){
00206 x = (int)(col+msr.xCenter-Radius);
00207 xx = col-Radius;
00208 xx *= xx;
00209 if(xx < (r2-y2) && x >= 0 && x < Src->mem2d->GetNx())
00210 Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt(x,y), col, line);
00211 }
00212 }
00213 return MATH_OK;
00214 }
00215 return MATH_SELECTIONERR;
00216 }
00217
00218
00219
00220 gboolean TR_QuenchScan(MATHOPPARAMS){
00221 long col, line;
00222 ZData *SrcZ, *SrcZn, *DestZ;
00223 double val;
00224 XSM_DEBUG (DBG_L3, "Quench Scan");
00225
00226 SrcZ = Src->mem2d->data;
00227
00228 DestZ = Dest->mem2d->data;
00229
00230 Dest->data.s.nx = (Src->data.s.nx/2) & ~1;
00231 Dest->data.s.dx = Src->data.s.dx*2.;
00232 Dest->data.s.ny = (Src->data.s.ny/2) & ~1;
00233 Dest->data.s.dy = Src->data.s.dy*2.;
00234 Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
00235
00236 Mem2d NextLine(Src->mem2d->GetNx(), 1, Src->mem2d->GetTyp());
00237 SrcZn = NextLine.data;
00238 for(line = 0; line < Dest->data.s.ny; line++){
00239 NextLine.CopyFrom(Src->mem2d, 0,2*line+1, 0,0, Src->mem2d->GetNx());
00240 DestZ->SetPtr(0, line);
00241 SrcZ ->SetPtr(0, 2*line);
00242 SrcZn->SetPtr(0, 0);
00243 for(col = 0; col < Dest->data.s.nx; col++){
00244 val = SrcZ ->GetNext();
00245 val += SrcZ ->GetNext();
00246 val += SrcZn->GetNext();
00247 val += SrcZn->GetNext();
00248 val /= 4.;
00249 DestZ->SetNext(val);
00250 }
00251 }
00252
00253 return MATH_OK;
00254 }
00255
00256
00257 gboolean ZoomInScan(MATHOPPARAMS){
00258 long col, line;
00259 MOUSERECT msr;
00260 double FXY, facx, facy;
00261
00262 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00263
00264 if( msr.xSize < 1 || msr.ySize < 1)
00265 return MATH_SELECTIONERR;
00266
00267
00268 FXY = (double)(Src->data.s.nx*Src->data.s.ny);
00269 Dest->data.s.nx = (int)(sqrt(FXY * msr.Aspect));
00270 Dest->data.s.ny = (int)(sqrt(FXY / msr.Aspect));
00271
00272 Dest->data.s.x0 += (msr.xLeft + msr.xSize/2 - Src->data.s.nx/2)*Src->data.s.dx;
00273 Dest->data.s.y0 -= msr.yTop*Src->data.s.dy;
00274
00275 facx = (double)Dest->data.s.nx/msr.xSize;
00276 facy = (double)Dest->data.s.ny/msr.ySize;
00277 Dest->data.s.dx /= facx;
00278 Dest->data.s.dy /= facy;
00279
00280
00281
00282
00283
00284
00285
00286 Dest->data.s.rx = (Dest->data.s.nx-1)*Dest->data.s.dx;
00287 Dest->data.s.ry = (Dest->data.s.ny-1)*Dest->data.s.dy;
00288
00289 Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
00290
00291 for(line=0; line<Dest->data.s.ny; line++)
00292 for(col=0; col<Dest->data.s.nx; col++)
00293 Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt((int)(msr.xLeft+col/facx),
00294 (int)(msr.yTop+line/facy)),
00295 col, line);
00296 return MATH_OK;
00297 }
00298
00299
00300 gboolean ZoomOutScan(MATHOPPARAMS){
00301 int col, line;
00302 int x0,y0;
00303 XSM_DEBUG(DBG_L2, "ZoomOut" );
00304
00305
00306 Dest->data.s.nx *= 3;
00307 Dest->data.s.ny *= 3;
00308 Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
00309
00310 Dest->data.s.y0 += Src->data.s.ny*Src->data.s.dy;
00311
00312
00313 Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
00314 Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;
00315
00316 x0=Src->data.s.nx;
00317 y0=Src->data.s.ny;
00318 for(line = 0; line < Src->mem2d->GetNy(); line++)
00319 for(col = 0; col < Src->mem2d->GetNx(); col++)
00320 Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPktLineReg(col, line), col+x0, line+y0);
00321 return MATH_OK;
00322 }
00323
00324
00325 gboolean BgLin1DScan(MATHOPPARAMS){
00326 int line, col;
00327 XSM_DEBUG (DBG_L3, "BgLin1D Scan");
00328
00329 for(line=0; line<Dest->data.s.ny; line++)
00330 for(col=0; col<Dest->data.s.nx; col++)
00331 Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPktLineReg(col, line), col, line);
00332 return MATH_OK;
00333 }
00334
00335
00336 gboolean BgERegress(MATHOPPARAMS){
00337 double sumi, sumiy, sumyq, sumiq, sumy, val;
00338 double ysumi, ysumiy, ysumyq, ysumiq, ysumy;
00339 double a, b, n, imean, ymean,ax,bx,ay,by;
00340 long line, i;
00341 MOUSERECT msr;
00342
00343 XSM_DEBUG (DBG_L3, "BgERegress Scan");
00344
00345 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00346
00347 if( msr.xSize < 1 || msr.ySize < 1)
00348 return MATH_SELECTIONERR;
00349
00350 sumi = sumiq = sumy = sumyq = sumiy = ax = ay = bx = by = 0.;
00351 ysumi = ysumiq = ysumy = ysumyq = ysumiy = 0.;
00352 for (val=msr.xLeft; val < msr.xRight; val += 1.){
00353 sumiq += val*val;
00354 sumi += val;
00355 }
00356 n = msr.xSize;
00357 imean = sumi / n;
00358 for (line = msr.yTop; line<msr.yBottom; ++line) {
00359 sumy = sumyq = sumiy = 0.;
00360 for (i=msr.xLeft; i < msr.xRight; ++i) {
00361 val = Src->mem2d->GetDataPkt(i, line);
00362 sumy += val;
00363 sumyq += val*val;
00364 sumiy += val*i;
00365 }
00366 ymean = sumy / n;
00367 a = (sumiy- n * imean * ymean ) / (sumiq - n * imean * imean);
00368 b = ymean - a * imean;
00369 ax += a;
00370 bx += b;
00371
00372 val = line;
00373 ysumy += b;
00374 ysumyq += b*b;
00375 ysumiy += b*val;
00376 ysumiq += val*val;
00377 ysumi += val;
00378 }
00379 n = msr.ySize;
00380 ax = ax/n;
00381 bx = bx/n;
00382 imean = ysumi / n;
00383 ymean = ysumy/ n;
00384 ay = (ysumiy- n * imean * ymean ) / (ysumiq - n * imean * imean);
00385 by = ymean - ay * imean;
00386 sumy = sumi = 0.;
00387
00388 for (line = 0; line<Dest->data.s.ny; ++line, sumy += ay) {
00389 sumi = by + sumy + .5;
00390 for (i=0; i < Dest->data.s.nx; ++i, sumi += ax)
00391 Dest->mem2d->PutDataPkt (Src->mem2d->GetDataPkt (i, line) - sumi,
00392 i, line);
00393 }
00394 return MATH_OK;
00395 }
00396
00397
00398 gboolean BgParabolRegress(MATHOPPARAMS){
00399 double sumi, sumiy, sumyq, sumiq, sumy, val;
00400 double ysumi, ysumiy, ysumyq, ysumiq, ysumy;
00401 double a, b, n, imean, ymean,ax,bx,ay,by;
00402 int line, i;
00403 MOUSERECT msr;
00404
00405 XSM_DEBUG (DBG_L3, "BgParabolRegress Scan");
00406
00407 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00408
00409 if( msr.xSize < 1 || msr.ySize < 1)
00410 return MATH_SELECTIONERR;
00411
00412
00413
00414
00415 sumi = sumiq = sumy = sumyq = sumiy = ax = ay = bx = by = 0.;
00416 ysumi = ysumiq = ysumy = ysumyq = ysumiy = 0.;
00417
00418 for (val=msr.xLeft; val < msr.xRight; val += 1.){
00419 sumiq += val*val;
00420 sumi += val;
00421 }
00422 n = msr.xSize;
00423 imean = sumi / n;
00424 for (line = msr.yTop; line<msr.yBottom; ++line) {
00425 sumy = sumyq = sumiy = 0.;
00426 for (i=msr.xLeft; i < msr.xRight; ++i) {
00427 val = Src->mem2d->GetDataPkt(i, line);
00428 sumy += val;
00429 sumyq += val*val;
00430 sumiy += val*i;
00431 }
00432 ymean = sumy / n;
00433 a = (sumiy- n * imean * ymean ) / (sumiq - n * imean * imean);
00434 b = ymean - a * imean;
00435 ax += a;
00436 bx += b;
00437
00438 val = line;
00439 ysumy += b;
00440 ysumyq += b*b;
00441 ysumiy += b*val;
00442 ysumiq += val*val;
00443 ysumi += val;
00444 }
00445 n = msr.ySize;
00446 ax = ax/n;
00447 bx = bx/n;
00448 imean = ysumi / n;
00449 ymean = ysumy/ n;
00450 ay = (ysumiy- n * imean * ymean ) / (ysumiq - n * imean * imean);
00451 by = ymean - ay * imean;
00452 sumy = sumi = 0.;
00453
00454 for (line = 0; line<Dest->data.s.ny; ++line, sumy += ay) {
00455 sumi = by + sumy + .5;
00456 for (i=0; i < Dest->data.s.nx; ++i, sumi += ax)
00457 Dest->mem2d->PutDataPkt (Src->mem2d->GetDataPkt (i, line) + sumi,
00458 i, line);
00459 }
00460
00461
00462 double M0, M1, M2, li;
00463 double D, N;
00464 double qx0, qx, qxx, qy0, qy, qyy, lb;
00465
00466
00467 N = (double)Src->mem2d->GetNx();
00468 D = (N*N-1.0)*(N*N-4.0);
00469 qx0 = qx = qxx = 0.;
00470 for (int line=0; line < Src->mem2d->GetNy(); line++) {
00471 M0 = 0.0; M1 = 0.0; M2 = 0.0;
00472
00473 for (int col = 0; col < Src->mem2d->GetNx(); col++) {
00474 li = (double)col;
00475 M0 += lb = Src->mem2d->GetDataPkt (col, line);
00476 M1 += lb * li;
00477 M2 += lb * li*li;
00478 }
00479
00480 M0 /= N;
00481 M1 /= N;
00482 M2 /= N;
00483
00484 qx0 += (3.*(N+1.)*(N+2.)/D)*((3.*N*N+3.*N+2.)*M0-6*(2*N+1.)*M1+10.*M2);
00485 qx += -6./D*(3*(N+1.)*(N+2.)*(2.*N+1.)*M0-2.*(2.*N+1.)*(8.*N+11.)*M1+30.*(N+1.)*M2);
00486 qxx += 30./D*((N+1.)*(N+2.)*M0-6.*(N+1.)*M1+6.*M2);
00487 }
00488 qx0 /= (double) Src->mem2d->GetNy();
00489 qx /= (double) Src->mem2d->GetNy();
00490 qxx /= (double) Src->mem2d->GetNy();
00491
00492
00493
00494 N = (double)Src->mem2d->GetNy();
00495 D = (N*N-1.0)*(N*N-4.0);
00496 qy0 = qy = qyy = 0.;
00497 for (int line=0; line < Src->mem2d->GetNx(); line++) {
00498 M0 = 0.0; M1 = 0.0; M2 = 0.0;
00499
00500 for (int col = 0; col < Src->mem2d->GetNy(); col++) {
00501 li = (double)col;
00502 M0 += lb = Src->mem2d->GetDataPkt (line, col);
00503 M1 += lb * li;
00504 M2 += lb * li*li;
00505 }
00506
00507 M0 /= N;
00508 M1 /= N;
00509 M2 /= N;
00510
00511 qy0 += (3.*(N+1.)*(N+2.)/D)*((3.*N*N+3.*N+2.)*M0-6*(2*N+1.)*M1+10.*M2);
00512 qy += -6./D*(3*(N+1.)*(N+2.)*(2.*N+1.)*M0-2.*(2.*N+1.)*(8.*N+11.)*M1+30.*(N+1.)*M2);
00513 qyy += 30./D*((N+1.)*(N+2.)*M0-6.*(N+1.)*M1+6.*M2);
00514 }
00515 qy0 /= (double) Src->mem2d->GetNx();
00516 qy /= (double) Src->mem2d->GetNx();
00517 qyy /= (double) Src->mem2d->GetNx();
00518
00519
00520 for (line = 0; line<Dest->data.s.ny; ++line)
00521 for (i=0; i < Dest->data.s.nx; ++i, sumi += ax){
00522 sumi = qx0 + qy0 + qx*i + qy*line + qxx*i*i + qyy*line*line;
00523 Dest->mem2d->PutDataPkt (Dest->mem2d->GetDataPkt (i, line) - sumi, i, line);
00524 }
00525
00526 return MATH_OK;
00527 }
00528
00529
00530
00531
00532
00533
00534
00535
00536 gboolean F1D_Despike(MATHOPPARAMS){
00537 int i;
00538 int n=Dest->mem2d->GetNx();
00539 int m=Dest->mem2d->GetNy();
00540
00541 XSM_DEBUG (DBG_L3, "F1D Despike");
00542
00543 Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00544 for(i=1;i<m;i++){
00545 int q, j;
00546 double avg, avgdev;
00547
00548 {
00549 double sum=0.;
00550 double absdevsum=0.;
00551 double diff;
00552 for(j=0;j<n;j++){
00553 diff = Src->mem2d->GetDataPkt(j,i) - Src->mem2d->GetDataPkt(j,i-1);
00554 sum += diff;
00555 absdevsum += fabs(diff);
00556 }
00557 avg = sum/n;
00558 avgdev = absdevsum/(n-1);
00559 }
00560 for(q=0;q<4;q++) {
00561 double sum, dsum, na;
00562 for(na=dsum=sum=0.0,j=0;j<n;j++) {
00563 const double sdev = (Src->mem2d->GetDataPkt(j,i) - Src->mem2d->GetDataPkt(j,i-1) - avg)/avgdev;
00564 const double k = 1/(3+sdev*sdev);
00565 sum += k*sdev;
00566 dsum += k*sdev*sdev;
00567 na += k;
00568 }
00569 avg += avgdev*sum/na;
00570 avgdev = 0.2*avgdev + 0.8*sqrt(dsum/(na-1));
00571 }
00572
00573 {
00574 SHT iavg = (SHT)(rint(avg));
00575 for(j=0;j<n;j++)
00576 Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt(j,i) - iavg, j,i);
00577 }
00578 }
00579 return MATH_OK;
00580 }
00581
00582
00583
00584 gboolean F1D_LogPowerSpec(MATHOPPARAMS)
00585 {
00586 XSM_DEBUG (DBG_L3, "F1D LogPowerSpec");
00587
00588 ZData *SrcZ, *DestZ;
00589
00590 SrcZ = Src->mem2d->data;
00591 DestZ = Dest->mem2d->data;
00592
00593 XSM_DEBUG (DBG_L3, "F1D LogPowerSpec: using libfftw");
00594
00595
00596 double *in = new double [Src->mem2d->GetNx()];
00597 fftw_complex *out = new fftw_complex [Src->mem2d->GetNx()];
00598
00599
00600 fftw_plan plan = fftw_plan_dft_r2c_1d (Src->mem2d->GetNx(), in, out, FFTW_ESTIMATE);
00601 if (plan == NULL) {
00602 XSM_DEBUG (DBG_L3, "F1D LogPowerSpec: libfftw: Error creating plan");
00603 return MATH_LIB_ERR;
00604 }
00605
00606 Dest->mem2d->Resize (Src->mem2d->GetNx()/2, Src->mem2d->GetNy(), ZD_DOUBLE);
00607 Dest->data.s.nx = Dest->mem2d->GetNx();
00608 Dest->data.s.x0 = 0.;
00609 Dest->data.s.y0 = 0.;
00610 Dest->data.s.rx = Src->data.s.rx/2.;
00611 Dest->data.s.ry = Src->data.s.ry;
00612 Dest->data.s.dl = Src->data.s.rx/2.;
00613 Dest->data.s.dz = Src->data.s.dz;
00614 Dest->data.s.dy = Src->data.s.dy;
00615
00616
00617 for (int line = 0; line < Src->mem2d->GetNy(); line++) {
00618 SrcZ ->SetPtr(0, line);
00619
00620
00621 for (int col = 0; col < Src->mem2d->GetNx(); col++)
00622 in[col] = SrcZ->GetNext();
00623
00624
00625 fftw_execute (plan);
00626
00627
00628
00629
00630
00631 DestZ->SetPtr(0, line);
00632
00633 for (int col = Src->mem2d->GetNx()/2-1; col >= 0; --col)
00634 DestZ->SetNext ( ZEROVALUE + ( c_re(out[col])*c_re(out[col]) +
00635 c_im(out[col])*c_im(out[col]) ));
00636 }
00637
00638 int n = Dest->mem2d->GetNx();
00639 double s = Src->data.s.rx/n;
00640 for (int i=0; i<n; ++i)
00641 Dest->mem2d->data->SetXLookup(i, (n-i)*s);
00642
00643
00644 fftw_destroy_plan (plan);
00645
00646
00647 delete in;
00648 delete out;
00649
00650 return MATH_OK;
00651 }
00652
00653
00654
00655
00656 gboolean F2D_RemoveRect(MATHOPPARAMS)
00657 {
00658 XSM_DEBUG (DBG_L3, "F2D Remove Rect");
00659
00660
00661 MOUSERECT msr;
00662 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00663
00664
00665 if( msr.xSize < 1 || msr.yBottom > Src->mem2d->GetNy() || msr.yTop < 0)
00666 return MATH_SELECTIONERR;
00667
00668
00669 Dest->mem2d->Resize(Src->mem2d->GetNx(), Src->mem2d->GetNy()-msr.ySize);
00670
00671
00672 for (int line = 0; line < msr.yTop; line++)
00673 Dest->mem2d->CopyFrom(Src->mem2d, 0,line, 0, line, Dest->mem2d->GetNx());
00674
00675 for (int line = msr.yBottom+1; line < Src->mem2d->GetNy(); line++)
00676 Dest->mem2d->CopyFrom(Src->mem2d, 0,line, 0, line-msr.ySize-1, Dest->mem2d->GetNx());
00677
00678 return MATH_OK;
00679 }
00680
00681
00682
00683 gboolean F2D_Despike(MATHOPPARAMS){
00684 int i,j,k,l,za=0,nx,ny;
00685 int anz=1;
00686 int num=1;
00687 double reihe1[20],reihe2[20],mi;
00688
00689 XSM_DEBUG (DBG_L3, "F2D Despike");
00690
00691 BenchStart(mbs,"2D Despike",Dest->mem2d->GetNx()*Dest->mem2d->GetNy());
00692
00693 Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0,
00694 nx=Dest->mem2d->GetNx(),ny=Dest->mem2d->GetNy());
00695
00696 for(j=anz; j<ny-anz; ++j) {
00697 for(i=anz; i<nx-anz; ++i) {
00698 for(k=j-anz, l=0; k<=j+anz; k++,l++)
00699 reihe1[l] = Src->mem2d->data->Z(i,k);
00700 for(k=0; k<2*anz+1;k++) {
00701 mi = 32650;
00702 for(l=0; l<2*anz+1;l++) {
00703 if(reihe1[l]<mi) { mi=reihe1[l]; reihe2[k]=mi; za=l; }
00704
00705 }
00706 reihe1[za]=32650;
00707 }
00708 Dest->mem2d->data->Z(reihe2[num], i,j);
00709 }
00710 }
00711 BenchStop(mbs);
00712
00713 return MATH_OK;
00714 }
00715
00716
00717
00718
00719
00720 extern gboolean F2D_LineShifts(MATHOPPARAMS)
00721 {
00722 XSM_DEBUG (DBG_L3, "F2D LineShifts");
00723 double avg;
00724 int nx,ny;
00725 int x1,x2,y1,y2;
00726 ZData *SrcZ, *DestZ;
00727
00728 nx=Dest->mem2d->GetNx();
00729 ny=Dest->mem2d->GetNy();
00730 BenchStart(mbs,"F2D LineShifts",nx*ny);
00731
00732 x1=0; x2=nx; y1=y2=0;
00733
00734 if (Src->number_of_object () > 0){
00735 XSM_DEBUG(DBG_L2, "found object(s)" );
00736 if( Src->get_object_data (0) -> get_num_points () == 2){
00737 double x=0.;
00738 double y=0.;
00739 Src->get_object_data (0) -> get_xy (0, x, y);
00740 Src->World2Pixel (x,y, x1,y1);
00741 Src->get_object_data (0) -> get_xy (1, x, y);
00742 Src->World2Pixel (x,y, x2,y2);
00743 XSM_DEBUG(DBG_L2, "got x range " << x1 << " .. " << x2 );
00744 }
00745
00746 if (x1 > x2){
00747 int tmp = x1;
00748 x1 = x2;
00749 x2 = tmp;
00750 }
00751 }
00752
00753 Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, nx,ny);
00754 SrcZ = Src->mem2d->data;
00755 DestZ = Dest->mem2d->data;
00756 for (int line = 0; line < ny; line++) {
00757 avg = 0.;
00758 SrcZ->SetPtr(0,line);
00759
00760 for (int col = x1; col < x2; col++)
00761 avg += SrcZ->GetNext();
00762
00763 avg /= (x2-x1);
00764 SrcZ->SetPtr(0,line);
00765 DestZ->SetPtr(0,line);
00766
00767 for (int col = 0; col < nx; col++)
00768 DestZ->SetNext(SrcZ->GetNext() - avg);
00769 }
00770 BenchStop(mbs);
00771
00772 return MATH_OK;
00773 }
00774
00775
00776
00777
00778
00779
00780 extern gboolean F2D_LineInterpol(MATHOPPARAMS)
00781 {
00782 XSM_DEBUG (DBG_L3, "F2D LineInterpol");
00783
00784 MOUSERECT msr;
00785
00786 MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00787
00788 if( msr.xSize < 1)
00789 return MATH_SELECTIONERR;
00790
00791
00792 Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00793
00794
00795 if(msr.yTop > 0){
00796 if(msr.yTop < (Src->mem2d->GetNy()-1)){
00797
00798 Mem2d PreLine(Src->mem2d->GetNx(), 1, Src->mem2d->GetTyp());
00799
00800 PreLine.CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,0, Src->mem2d->GetNx());
00801 PreLine.data->SetPtr(0,0);
00802
00803 Src->mem2d->data->SetPtr(0,msr.yTop+1);
00804
00805 Dest->mem2d->data->SetPtr(0,msr.yTop);
00806 for(int col = 0; col < Src->mem2d->GetNx(); col++)
00807 Dest->mem2d->data->SetNext((PreLine.data->GetNext()+Src->mem2d->data->GetNext())/2.);
00808 }
00809 else
00810 Dest->mem2d->CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,msr.yTop, Src->mem2d->GetNx());
00811 }
00812 else
00813 Dest->mem2d->CopyFrom(Src->mem2d, 0,msr.yTop+1, 0,msr.yTop, Src->mem2d->GetNx());
00814
00815 return MATH_OK;
00816 }
00817
00818
00819
00820
00821
00822 gboolean F2D_LogPowerSpec(MATHOPPARAMS)
00823 {
00824 XSM_DEBUG (DBG_L3, "F2D FFT");
00825 ZData *SrcZ;
00826
00827 SrcZ = Src->mem2d->data;
00828 Dest->mem2d->Resize (Dest->data.s.nx, Dest->data.s.ny, ZD_COMPLEX);
00829
00830 Dest->data.s.nvalues=3;
00831 Dest->data.s.ntimes=1;
00832 Dest->mem2d->data->SetVLookup(0, 0);
00833 Dest->mem2d->data->SetVLookup(1, 1);
00834
00835 XSM_DEBUG (DBG_L3, "F2D FFT: using libfftw");
00836
00837
00838 int xlen=2*(Src->mem2d->GetNx ()/2+1);
00839 double *in = new double [ Src->mem2d->GetNx () * Src->mem2d->GetNy () ];
00840 fftw_complex *dat = new fftw_complex [ xlen*Src->mem2d->GetNy () ];
00841
00842 if (dat == NULL) {
00843 XSM_DEBUG (DBG_L3, "F2D FFT: libfftw: Error allocating memory for complex data");
00844 return MATH_NOMEM;
00845 }
00846
00847 memset(in, 0, sizeof(in));
00848
00849
00850 fftw_plan plan = fftw_plan_dft_r2c_2d (Src->mem2d->GetNy (), Src->mem2d->GetNx (),
00851 in, dat, FFTW_ESTIMATE);
00852
00853 if (plan == NULL) {
00854 XSM_DEBUG (DBG_L3, "F2D FFT: libfftw: Error creating plan");
00855 return MATH_LIB_ERR;
00856 }
00857
00858
00859 for (int line=0; line < Src->mem2d->GetNy(); line++) {
00860 SrcZ ->SetPtr(0, line);
00861 for (int col=0; col < Src->mem2d->GetNx(); col++)
00862 in[line*Src->mem2d->GetNx() + col] = SrcZ->GetNext();
00863 }
00864
00865
00866 fftw_execute (plan);
00867
00868
00869
00870
00871
00872 XSM_DEBUG (DBG_L3, "F2D FFT done, converting data to ABS/RE/IM");
00873
00874 xlen /= 2;
00875
00876 double I;
00877 int m,n;
00878
00879 for (int line=0; line < Src->mem2d->GetNy(); line++) {
00880 for (int i,col=0; col < Src->mem2d->GetNx()/2; col++) {
00881 i=line*xlen + col;
00882 m = QSWP(col, Src->mem2d->GetNx());
00883 n = QSWP(line, Src->mem2d->GetNy());
00884 I = ZEROVALUE + sqrt (c_re(dat[i]) * c_re(dat[i]) + c_im(dat[i]) * c_im(dat[i]));
00885 Dest->mem2d->PutDataPkt(I, m, n, 0);
00886 Dest->mem2d->PutDataPkt(c_re(dat[i]), m, n, 1);
00887 Dest->mem2d->PutDataPkt(c_im(dat[i]), m, n, 2);
00888 if(line != (Src->mem2d->GetNy()/2)){
00889 m = Src->mem2d->GetNx() - m;
00890 n = Src->mem2d->GetNy() - n;
00891 Dest->mem2d->PutDataPkt(I, m, n, 0);
00892 Dest->mem2d->PutDataPkt(c_re(dat[i]), m, n, 1);
00893 Dest->mem2d->PutDataPkt(c_im(dat[i]), m, n, 2);
00894 }
00895 }
00896 }
00897
00898
00899 fftw_destroy_plan(plan);
00900
00901
00902 delete in;
00903 delete dat;
00904
00905 return MATH_OK;
00906 }
00907
00908
00909
00910
00911
00912
00913 gboolean F1D_ift_ft(MATH2OPPARAMS,
00914 gboolean (*spkfkt)(MATH2OPPARAMS,
00915 fftw_complex *dat,
00916 int line)
00917 ){
00918 XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) )");
00919 gboolean ret=MATH_OK;
00920 ZData *SrcZ, *DestZ;
00921
00922 SrcZ = Src1->mem2d->data;
00923 DestZ = Dest->mem2d->data;
00924
00925
00926
00927
00928 if(Src2 != NULL && (Src1->data.s.nx != Src2->data.s.nx || Src1->data.s.ny != Src2->data.s.ny))
00929 return MATH_SELECTIONERR;
00930
00931
00932
00933 int line, col;
00934
00935
00936 double *in = new double [ Src1->mem2d->GetNx() ];
00937 fftw_complex *dat = new fftw_complex [ Src1->mem2d->GetNx() ];
00938 double scale;
00939
00940 fftw_plan plan = fftw_plan_dft_r2c_1d (Src1->mem2d->GetNx(), in, dat, FFTW_ESTIMATE);
00941 fftw_plan planinv = fftw_plan_dft_c2r_1d (Src1->mem2d->GetNx(), dat, in, FFTW_ESTIMATE);
00942
00943 if (!plan || !planinv) {
00944 XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error creating plans");
00945 if(plan)
00946 fftw_destroy_plan (plan);
00947 if(planinv)
00948 fftw_destroy_plan (planinv);
00949 return MATH_LIB_ERR;
00950 }
00951
00952
00953
00954 if (!in || !dat) {
00955 XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error allocating memory for complex data");
00956 if(in) delete in;
00957 if(dat) delete dat;
00958 return MATH_NOMEM;
00959 }
00960
00961 memset(in, 0, sizeof(in));
00962
00963
00964 for (line=0; line < Src1->mem2d->GetNy(); line++) {
00965 SrcZ ->SetPtr(0, line);
00966 for (col=0; col < Src1->mem2d->GetNx(); col++)
00967 in[col] = SrcZ->GetNext();
00968
00969
00970 fftw_execute (plan);
00971
00972
00973
00974
00975 if (spkfkt (MATH2OPVARS, dat, line) == MATH_OK){
00976
00977
00978
00979 fftw_execute (planinv);
00980
00981
00982
00983 scale = 1. / Src1->mem2d->GetNx();
00984 DestZ->SetPtr(0, line);
00985 for(col=0; col<Dest->mem2d->GetNx(); col++)
00986 DestZ->SetNext(scale*(in[col]));
00987 }
00988 else
00989 ret=MATH_SELECTIONERR;
00990 }
00991
00992
00993 delete in;
00994
00995
00996 fftw_destroy_plan (plan);
00997 fftw_destroy_plan (planinv);
00998
00999 return ret;
01000 }
01001
01002
01003
01004 gboolean SpkWindow1D(MATH2OPPARAMS, fftw_complex *dat, int line){
01005 int col;
01006
01007
01008
01009
01010
01011 int n = Src1->mem2d->GetNx();
01012 for (col=0; col < n/2; col++)
01013 if (Src2->mem2d->GetDataPkt(n/2-col-1, line) == ZEROVALUE){
01014 c_re (dat[col]) = 0.;
01015 c_im (dat[col]) = 0.;
01016 c_re (dat[n-col-1]) = 0.;
01017 c_im (dat[n-col-1]) = 0.;
01018 }
01019 return MATH_OK;
01020 }
01021
01022
01023
01024 gboolean SpkGaussStop1D(MATH2OPPARAMS, fftw_complex *dat, int line){
01025 MOUSERECT msr;
01026 MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
01027
01028 if( msr.xSize < 1 || msr.ySize < 1)
01029 return MATH_SELECTIONERR;
01030
01031
01032
01033 double r, scale;
01034 int n = Src1->mem2d->GetNx();
01035 for (int col = 0; col < n/2; col++) {
01036 r = (n/2-col-1) - msr.xCenter; r *= r;
01037 scale = 1 - exp(-r / msr.Radius2);
01038 c_re (dat[col]) *= scale;
01039 c_re (dat[n-col-1]) *= scale;
01040 }
01041 return MATH_OK;
01042 }
01043
01044
01045 gboolean SpkGaussPass1D(MATH2OPPARAMS, fftw_complex *dat, int line){
01046 MOUSERECT msr;
01047 MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
01048
01049 if( msr.xSize < 1 || msr.ySize < 1)
01050 return MATH_SELECTIONERR;
01051
01052
01053
01054 double r, scale;
01055 int n = Src1->mem2d->GetNx();
01056 for (int col = 0; col < n/2; col++) {
01057 r = (n/2-col-1) - msr.xCenter; r *= r;
01058 scale = exp(-r / msr.Radius2);
01059 c_re (dat[col]) *= scale;
01060 c_re (dat[n-col-1]) *= scale;
01061 }
01062 return MATH_OK;
01063 }
01064
01065
01066
01067
01068 gboolean F1D_FT_Window(MATH2OPPARAMS)
01069 {
01070 XSM_DEBUG (DBG_L3, "F1D FT Window");
01071 return(F1D_ift_ft(MATH2OPVARS, SpkWindow1D));
01072 return MATH_OK;
01073 }
01074
01075 gboolean F1D_FT_GaussStop(MATH2OPPARAMS)
01076 {
01077 XSM_DEBUG (DBG_L3, "F1D FT GaussStop");
01078 return(F1D_ift_ft(MATH2OPVARS, SpkGaussStop1D));
01079 return MATH_OK;
01080 }
01081
01082 gboolean F1D_FT_GaussPass(MATH2OPPARAMS)
01083 {
01084 XSM_DEBUG (DBG_L3, "F1D FT GaussPass");
01085 return(F1D_ift_ft(MATH2OPVARS, SpkGaussPass1D));
01086 return MATH_OK;
01087 }
01088
01089
01090
01091
01092
01093
01094 gboolean F2D_ift_ft(MATH2OPPARAMS, gboolean (*spkfkt)(MATH2OPPARAMS, fftw_complex *dat)){
01095 XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) )");
01096 gboolean ret=MATH_OK;
01097 ZData *SrcZ, *DestZ;
01098
01099 SrcZ = Src1->mem2d->data;
01100 DestZ = Dest->mem2d->data;
01101
01102
01103
01104
01105 if(Src2 != NULL && (Src1->data.s.nx != Src2->data.s.nx || Src1->data.s.ny != Src2->data.s.ny))
01106 return MATH_SELECTIONERR;
01107
01108
01109
01110 int line, col;
01111
01112
01113 int xlen=2*(Src1->mem2d->GetNx()/2+1);
01114 double *in = new double [ xlen*Src1->mem2d->GetNy() ];
01115 fftw_complex *dat = (fftw_complex*)(in);
01116 double scale;
01117
01118 if (!in || !dat) {
01119 XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error allocating memory for complex data");
01120 if(in) delete in;
01121 if(dat) delete dat;
01122 return MATH_NOMEM;
01123 }
01124
01125
01126 fftw_plan plan = fftw_plan_dft_r2c_2d (Src1->mem2d->GetNy(), Src1->mem2d->GetNx(), in, dat, FFTW_ESTIMATE);
01127 fftw_plan planinv = fftw_plan_dft_c2r_2d (Src1->mem2d->GetNy(), Src1->mem2d->GetNx(), dat, in, FFTW_ESTIMATE);
01128
01129 if (!plan || !planinv) {
01130 XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error creating plans");
01131 if(plan)
01132 fftw_destroy_plan (plan);
01133 if(planinv)
01134 fftw_destroy_plan (planinv);
01135 return MATH_LIB_ERR;
01136 }
01137
01138 memset(in, 0, sizeof(in));
01139
01140
01141 for (line=0; line < Src1->mem2d->GetNy(); line++) {
01142 SrcZ ->SetPtr(0, line);
01143 for (col=0; col < Src1->mem2d->GetNx(); col++)
01144 in[line*xlen + col] = SrcZ->GetNext();
01145 }
01146
01147
01148
01149 fftw_execute (plan);
01150
01151
01152
01153 if (spkfkt (MATH2OPVARS, dat) == MATH_OK){
01154
01155
01156
01157 fftw_execute (planinv);
01158
01159
01160
01161 scale = 1. / Src1->mem2d->GetNx() / Src1->mem2d->GetNy();
01162 for(line=0; line<Dest->mem2d->GetNy(); line++){
01163 DestZ->SetPtr(0, line);
01164 for(col=0; col<Dest->mem2d->GetNx(); col++)
01165 DestZ->SetNext(scale*(in[line*xlen + col]));
01166 }
01167 }
01168 else
01169 ret=MATH_SELECTIONERR;
01170
01171
01172 delete in;
01173
01174
01175 fftw_destroy_plan (plan);
01176 fftw_destroy_plan (planinv);
01177
01178 return ret;
01179 }
01180
01181
01182 gboolean SpkWindow(MATH2OPPARAMS, fftw_complex *dat){
01183 int line, col;
01184
01185
01186
01187
01188
01189 int xlen = Src1->mem2d->GetNx()/2+1;
01190 for (line=0; line < Src1->mem2d->GetNy(); line++)
01191 for (col=0; col < xlen; col++)
01192 if(Src2->mem2d->GetDataPkt(QSWP(col,Src1->mem2d->GetNx()),
01193 QSWP(line,Src1->mem2d->GetNy()))
01194 == ZEROVALUE){
01195 c_re(dat[line*xlen + col]) = 0.;
01196 c_im(dat[line*xlen + col]) = 0.;
01197 }
01198 return MATH_OK;
01199 }
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210 gboolean F2D_iftXft(MATH2OPPARAMS){
01211 XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) ) -- Window(s)");
01212 return(F2D_ift_ft(MATH2OPVARS, SpkWindow));
01213 return MATH_OK;
01214 }
01215
01216
01217
01218 gboolean SpkGaussStop(MATH2OPPARAMS, fftw_complex *dat){
01219 MOUSERECT msr;
01220 MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
01221
01222 if( msr.xSize < 1 || msr.ySize < 1)
01223 return MATH_SELECTIONERR;
01224
01225
01226
01227 double r2, scale;
01228 int xlen = Src1->mem2d->GetNx()/2+1;
01229 for (int line = 0; line < Src1->mem2d->GetNy(); line++)
01230 for (int col = 0; col < xlen; col++) {
01231 r2 = (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) * (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) +
01232 (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter) * (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter);
01233 scale = 1 - exp(-r2 / msr.Radius2);
01234 c_re(dat[line*xlen + col]) *= scale;
01235 c_im(dat[line*xlen + col]) *= scale;
01236 }
01237 return MATH_OK;
01238 }
01239
01240
01241
01242 gboolean F2D_FT_GaussStop(MATH2OPPARAMS)
01243 {
01244 XSM_DEBUG (DBG_L3, "F2D FT GaussStop");
01245 return(F2D_ift_ft(MATH2OPVARS, SpkGaussStop));
01246 return MATH_OK;
01247 }
01248
01249
01250
01251 gboolean SpkGaussPass(MATH2OPPARAMS, fftw_complex *dat){
01252 MOUSERECT msr;
01253 MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
01254
01255 if( msr.xSize < 1 || msr.ySize < 1)
01256 return MATH_SELECTIONERR;
01257
01258
01259
01260 double r2, scale;
01261 int xlen = Src1->mem2d->GetNx()/2+1;
01262 for (int line = 0; line < Src1->mem2d->GetNy(); line++)
01263 for (int col = 0; col < xlen; col++) {
01264 r2 = (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) * (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) +
01265 (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter) * (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter);
01266 scale = exp(-r2 / msr.Radius2);
01267 c_re(dat[line*xlen + col]) *= scale;
01268 c_im(dat[line*xlen + col]) *= scale;
01269 }
01270 return MATH_OK;
01271 }
01272
01273
01274
01275
01276 gboolean F2D_FT_GaussPass(MATH2OPPARAMS)
01277 {
01278 XSM_DEBUG (DBG_L3, "F2D FT GaussPass");
01279 return(F2D_ift_ft(MATH2OPVARS, SpkGaussPass));
01280 }
01281
01282
01283
01284
01285 gboolean SpkAutoCorr(MATH2OPPARAMS, fftw_complex *dat)
01286 {
01287
01288 int xlen = Src1->mem2d->GetNx()/2+1;
01289 for (int line = 0; line < Src1->mem2d->GetNy(); line++)
01290 for (int col = 0; col < xlen; col++) {
01291 c_re(dat[line*xlen + col]) = sqrt( c_re(dat[line*xlen + col])*c_re(dat[line*xlen + col]) +
01292 c_im(dat[line*xlen + col])*c_im(dat[line*xlen + col]) );
01293 c_im(dat[line*xlen + col]) = 0.0;
01294 }
01295
01296
01297 return MATH_OK;
01298 }
01299
01300
01301
01302 gboolean F2D_AutoCorr(MATHOPPARAMS)
01303 {
01304 XSM_DEBUG (DBG_L3, "F2D AutoCorr");
01305 return(F2D_ift_ft(Src, NULL, Dest, SpkAutoCorr));
01306 }
01307
01308