xsmmath.C

Go to the documentation of this file.
00001 /* Gxsm - Gnome X Scanning Microscopy
00002  * universal STM/AFM/SARLS/SPALEED/... controlling and
00003  * data analysis software
00004  * 
00005  * Copyright (C) 1999,2000,2001,2002,2003 Percy Zahl
00006  *
00007  * Authors: Percy Zahl <zahl@users.sf.net>
00008  * additional features: Andreas Klust <klust@users.sf.net>
00009  * WWW Home: http://gxsm.sf.net
00010  *
00011  * This program is free software; you can redistribute it and/or modify
00012  * it under the terms of the GNU General Public License as published by
00013  * the Free Software Foundation; either version 2 of the License, or
00014  * (at your option) any later version.
00015  *
00016  * This program is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU General Public License
00022  * along with this program; if not, write to the Free Software
00023  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
00024  */
00025 
00026 /* -*- Mode: C++; indent-tabs-mode: nil; c-basic-offset: 8 c-style: "K&R" -*- */
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 // als Vorlage für Implementation neuer Filter
00074 //
00075 // Muster Filter "NULL EFFEKT"
00076 //==============================
00077 // MATHOPPARAMS:  
00078 // Scan *Src:  Datenquelle
00079 // Scan *Dest: Datenziel
00080 //
00081 // ===============================
00082 // *Dest ist ein neuer, leerer Scan mit identischen Parametern wie *Src !
00083 // !!!!!!  *Src sollte NIE UND NIMMER hier verändert werden !!!!!!
00084 //
00085 // Größe:
00086 // ===============================
00087 // nx = Dest->mem2d->GetNx(), ... GetNy()
00088 // Dest->Resize(Dest->data.s.nx, Dest->data.s.ny) :  Größe von ZielScan ändern
00089 //
00090 // Datenzugriff über Mem2d Objekt:
00091 // ===============================
00092 // * Kopieren eines rechteckigen Bereichs (Typ Src == Typ Dest):
00093 // Dest->mem2d->CopyFrom(Src->mem2d, int x0, int y0, int tox, int toy, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00094 // * Kopieren eines rechteckigen Bereichs mit ggf. Typenkonversion:
00095 // Dest->mem2d->ConvertFrom(Src->mem2d, int x0, int y0, int tox, int toy, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00096 //
00097 //
00098 // [[ Src->mem2d->GetDataLine(int line, ZTYPE *buffer)   :    ganze Zeile ]] please do not use !!
00099 // double value = Src->mem2d->GetDataPkt(int x, int y) :    einen Punkt
00100 // ... oder auch (etwas schneller, da direkter)
00101 // double value = Src->mem2d->data->Z(int x, int y) :    einen Punkt
00102 // value = Src->mem2d->GetDataPkt(int x, int y) :    einen Punkt
00103 // value = Src->mem2d->GetDataPktLineReg(int x, int y) :    einen Punkt, mit Line Regession !!
00104 // [[ Dest->mem2d->PutDataLine(int line, ZTYPE *buffer)  :    ganze Zeile ]] please do not use !!
00105 // Dest->mem2d->PutDataPkt(ZTYPE value, int x, int y) :    einen Punkt
00106 // ... oder auch (etwas schneller, da ohne invalidate lineregress parameters)
00107 // Src->mem2d->data->Z(double value, int x, int y) :    einen Punkt
00108 //
00109 // for large Data Transfer use:
00110 // ============================================================
00111 // for same types:
00112 // inline int CopyFrom(Mem2d *src, int x, int y, int tox, int toy, int nx, int ny=1)
00113 //
00114 // for different types:
00115 // inline int ConvertFrom(Mem2d *src, int x, int y, int tox, int toy, int nx, int ny=1)
00116 //
00117 // for quick linear data access of elements in one line use:
00118 // ============================================================
00119 // void   mem2d->data->SetPtr(int x, int y) to set internal pointer to x,y  -- no range check !!!
00120 // double mem2d->data->GetNext() to access Element x,y and point to next one in line y -- no range check !!!
00121 // double mem2d->data->SetNext(double z) to access Element x,y and point to next one in line y -- no range check !!!
00122 //
00123 /* For Example use constructs like:
00124   ZData  *SrcZ, *DestZ;
00125 
00126   SrcZ  =  Src->mem2d->data;
00127   DestZ = Dest->mem2d->data;
00128 
00129   for ( line=0; line < Dest->data.s.ny; line++){ 
00130     DestZ->SetPtr(0, line);
00131     SrcZ ->SetPtr(0, line);
00132     DestZ->SetNext(SrcZ->GetNext() + ....));
00133   }
00134 
00135   or a advanced example:  => F2D_LineInterpol(MATHOPPARAMS) 
00136 
00137  */
00138 // Parameter:
00139 // ===============================
00140 // Dest->data.s.xxx:  Scanparameterstruktur
00141 // Src->Pkt2d[0..15]: Maus-Punkte
00142 
00143 // uses new Object methods now !
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 // Ausschnitt herausschneiden
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     // Anpassen an nächstmöglichen Wert  
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;    // Get Radius Field as BgVal
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     // Anpassen an nächstmöglichen Wert  
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 // Größe halbieren, 2x2 Mittelung
00219 // Essential for Gxsm core!!!
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 // In Ausschnitt hineinzoomen 
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   // Seitenverhältnis und Fläche beibehalten:
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   // Anpassen an nächstmöglichen Wert  ... erst bevor Scanning
00282   Dest->data.s.dx = R2INT(Dest->data.s.dx/Inst->XResolution())*Inst->XResolution();
00283   Dest->data.s.dy = R2INT(Dest->data.s.dy/Inst->YResolution())*Inst->YResolution();
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 // auf Ausschnitt herauszoomen (*1/3)
00300 gboolean ZoomOutScan(MATHOPPARAMS){
00301   int col, line;
00302   int x0,y0;
00303   XSM_DEBUG(DBG_L2, "ZoomOut" );
00304 
00305   // Seitenverhältnis und Fläche beibehalten:
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   //  Dest->data.s.x0 += 0;
00310   Dest->data.s.y0 += Src->data.s.ny*Src->data.s.dy;
00311 
00312   // Anpassen an nächstmöglichen Wert  
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 // Copy Scan, do Lin Regress
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 // Do E Regress
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                 /* Werte fuer y-linregress */
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 // Do Parabol Regress
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         // plane regress first
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                 /* Werte fuer y-linregress */
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         // now evaluate 2nd order terms
00462         double M0, M1, M2, li;
00463         double D, N;
00464         double qx0, qx, qxx, qy0, qy, qyy, lb;
00465 
00466         // avg. 2nd order terms in X
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         // avg. 2nd order terms in Y
00493         // swaped line/col here
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 // 1D Filterfunktionen 
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 // compute 1D power spectrum (row by row)
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   // get memory for complex data
00596   double *in  = new double [Src->mem2d->GetNx()];
00597   fftw_complex *out = new fftw_complex [Src->mem2d->GetNx()];
00598 
00599   // create plan for fft
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   // compute 1D fourier transform for every row
00617   for (int line = 0; line < Src->mem2d->GetNy(); line++) {
00618     SrcZ ->SetPtr(0, line);
00619     
00620     // prepare image row data for fftw
00621     for (int col = 0; col < Src->mem2d->GetNx(); col++)
00622        in[col] = SrcZ->GetNext();
00623 
00624     // compute transform
00625     fftw_execute (plan);
00626 
00627     // convert complex data to image row
00628     // scale data to 16Bit integer
00629     // reorder data that freq 0 is in the middle and positive freq go to the right
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   // delete plan:
00644   fftw_destroy_plan (plan);
00645 
00646   // free temp data memory
00647   delete in;
00648   delete out;
00649 
00650   return MATH_OK;
00651 }
00652 
00653 // 2D Filterfunktionen
00654 
00655 // remove bad area
00656 gboolean F2D_RemoveRect(MATHOPPARAMS)
00657 {
00658   XSM_DEBUG (DBG_L3, "F2D Remove Rect");
00659 
00660   // get selection
00661   MOUSERECT msr;
00662   MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
00663 
00664   // check for valid selection
00665   if( msr.xSize  < 1 || msr.yBottom > Src->mem2d->GetNy() || msr.yTop < 0)
00666     return MATH_SELECTIONERR;
00667 
00668   // resize Dest
00669   Dest->mem2d->Resize(Src->mem2d->GetNx(), Src->mem2d->GetNy()-msr.ySize);
00670 
00671   // copy data except selection
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 // Do 2D DeSpike
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     }  /*i*/
00710   } /*j*/
00711   BenchStop(mbs);
00712 
00713   return MATH_OK;
00714 }
00715 
00716 
00717 
00718 // remove line shifts
00719 // uses new Object methods now !
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 // "fix" bad line
00778 // using an interpolation between the neighbour lines
00779 // uses new Object methods now !
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   // copy scan data:
00792   Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
00793 
00794   // interpolation if in between first and last line
00795   if(msr.yTop > 0){
00796     if(msr.yTop < (Src->mem2d->GetNy()-1)){
00797       // create help Mem2d object (only one line)
00798       Mem2d PreLine(Src->mem2d->GetNx(), 1, Src->mem2d->GetTyp());
00799       // with "PreLine" and set internal Ptr to start
00800       PreLine.CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,0, Src->mem2d->GetNx());
00801       PreLine.data->SetPtr(0,0);
00802       // set Src Ptr to "PostLine"
00803       Src->mem2d->data->SetPtr(0,msr.yTop+1);
00804       // set Dest Ptr to "Line"
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  // last line: copy previous line
00810       Dest->mem2d->CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,msr.yTop, Src->mem2d->GetNx());
00811   }
00812   else  // first line: copy 2nd line
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 // 2D power spectrum in logarithmic scale
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   // set "Complex" layer param defaults
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   // allocate memory for real and inplace complex data
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   // create plan for in-place transform
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   // convert image data to fftw_real
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   // compute 2D transform using in-place fourier transform
00866   fftw_execute (plan);
00867 
00868   // convert complex data to image
00869   // and flip image that spatial freq. (0,0) is in the middle
00870   // (quadrant swap using the QSWP macro defined in xsmmath.h)
00871   
00872   XSM_DEBUG (DBG_L3, "F2D FFT done, converting data to ABS/RE/IM");
00873 
00874   xlen /= 2; // jetzt für complex
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   // destroy plan
00899   fftw_destroy_plan(plan); 
00900 
00901   // free real/complex data memory
00902   delete in;
00903   delete dat;
00904 
00905   return MATH_OK;
00906 }
00907 
00908 
00909 // Echte 1D Fourier Filter :=)
00910 
00911 // FFTW-Kernel, applys some function to the spectra
00912 // Forward FT => do something with spectra => backward FT
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   // error reporting for invalid selections:
00926   // Src1 and Src2 scans must have the same dimensions
00927   // if Src2 == NULL filter spkfkt needs no selection, see e.g. SpkAutoCorr
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   // ======== 1. fourier transform Src1 => complex data spectrum =======
00932 
00933   int line, col;
00934 
00935   // allocate memory for complex data
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   // convert image data to fftw_real
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           // 1. compute 1D transform using in-place fourier transform
00970           fftw_execute (plan);
00971           
00972           // now the complex spectrum is in dat too !
00973           
00974           // 2. apply something to the spektra
00975           if (spkfkt (MATH2OPVARS, dat, line) == MATH_OK){
00976       
00977                   // 3. inverse fourier transform modified spectrum
00978       
00979                   fftw_execute (planinv);
00980       
00981                   // 4. take the real part as data and apply the fft normalization factor scale
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   // free [complex] data memory
00993   delete in;
00994 
00995   // destroy plan
00996   fftw_destroy_plan (plan); 
00997   fftw_destroy_plan (planinv); 
00998   
00999   return ret;
01000 }
01001 
01002 
01003 // Filter Kernel 1D "Window"
01004 gboolean SpkWindow1D(MATH2OPPARAMS, fftw_complex *dat, int line){
01005   int col;
01006   // ====== 2. multiply complex data spectrum with Src2, eg. windowing =====
01007   // and take into account that spatial freq. (0,0) is in the middle
01008 
01009   // apply Window
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 // Filter Kernel "Gauss Stopp"
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   // apply gaussian window to complex data
01032   // scale = 1 - exp(-r2/msr.Radius2)
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 // Filter Kernel "Gauss Pass"
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   // apply gaussian window to complex data
01053   // scale = exp(-r2/msr.Radius2)
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 // Fourierfilter:
01067 // Pass frequency window with gaussian profile
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 // Echte 2D Fourier Filter :=)
01091 
01092 // FFTW-Kernel, applys some function to the spectra
01093 // Forward FT => do something with spectra => backward FT
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   // error reporting for invalid selections:
01103   // Src1 and Src2 scans must have the same dimensions
01104   // if Src2 == NULL filter spkfkt needs no selection, see e.g. SpkAutoCorr
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   // ======== 1. fourier transform Src1 => complex data spectrum =======
01109 
01110   int line, col;
01111 
01112   // allocate memory for complex data
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); // use same memory - saves a lot space...
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   // create plan for in-place transform
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   // convert image data to fftw_real
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   // 1. compute 2D transform using in-place fourier transform
01148 
01149   fftw_execute (plan);
01150   // now the complex spectrum is in dat too !
01151 
01152   // 2. apply something to the spektra
01153   if (spkfkt (MATH2OPVARS, dat) == MATH_OK){
01154 
01155     // 3. inverse fourier transform modified spectrum
01156 
01157     fftw_execute (planinv);
01158 
01159     // 4. take the real part as data and apply the fft normalization factor scale
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   // free [complex] data memory
01172   delete in;
01173 
01174   // destroy plan
01175   fftw_destroy_plan (plan); 
01176   fftw_destroy_plan (planinv); 
01177 
01178   return ret;
01179 }
01180 
01181 // Filter Kernel "Window"
01182 gboolean SpkWindow(MATH2OPPARAMS, fftw_complex *dat){
01183   int line, col;
01184   // ====== 2. multiply complex data spectrum with Src2, eg. windowing =====
01185   // and take into account that spatial freq. (0,0) is in the middle
01186 
01187   // apply Window
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 // Src1 : Daten Quelle
01203 // Src2 : Fensterfunktion
01204 //        kann ggf. obiges PowerSpektrum sein, wobei
01205 //        folgendes gilt:
01206 //        Dest = Re(IFT( FT(Src1) * (Src2 == ZEROVALUE ? ZEROVALUE : 1) ))
01207 //        d.h. alle Spektralanteile, die in Src2 Null sind werden im Complexen Spektrum ausgeblendet !
01208 //        Dazu muessen natuerlich vorher im Power-Spektrum entspechende Bereiche auf Null gesetzt werden :=)
01209 //        --- that´s it --- PZ
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 // Filter Kernel "Gauss Stopp"
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   // apply gaussian window to complex data
01226   // scale = 1 - exp(-r2/msr.Radius2)
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 // Fourierfilter:
01241 // Pass frequency window with gaussian profile
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 // Filter Kernel "Gauss Pass"
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   // apply gaussian window to complex data
01259   // scale = exp(-r2/msr.Radius2)
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 // Fourierfilter:
01275 // Pass frequency window with gaussian profile
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 // Filter Kernel "Abs. Square"
01284 // (for auto correlation)
01285 gboolean SpkAutoCorr(MATH2OPPARAMS, fftw_complex *dat)
01286 {
01287   // take abs. value of complex data
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 // auto correlation
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 

Generated on Sat Apr 1 09:04:08 2006 for GXSM by  doxygen 1.4.6