xsmmath.C

Go to the documentation of this file.

/* Gxsm - Gnome X Scanning Microscopy
 * universal STM/AFM/SARLS/SPALEED/... controlling and
 * data analysis software
 * 
 * Copyright (C) 1999,2000,2001,2002,2003 Percy Zahl
 *
 * Authors: Percy Zahl <zahl@users.sf.net>
 * additional features: Andreas Klust <klust@users.sf.net>
 * WWW Home: http://gxsm.sf.net
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
 */

/* -*- Mode: C++; indent-tabs-mode: nil; c-basic-offset: 8 c-style: "K&R" -*- */

#include <math.h>

#include "xsmmath.h"
#include "xsmtypes.h"
#include "glbvars.h"

#include "bench.h"
#include "regress.h"


char *MathErrString[] = {
  "Math OK",
  "Size of sources is not the same !",
  "Size of selection is invalid !",
  "Division by zero occured !",
  "undefined error"
};


void MkMausSelect(Point2D *List, MOUSERECT *msel, int mx, int my){
  msel->xLeft   = MAX(0, MIN( MIN(List[0].x, List[1].x), mx));
  msel->xRight  = MAX(0, MIN( MAX(List[0].x, List[1].x), mx));
  msel->xSize   = msel->xRight - msel->xLeft;
  msel->xRatio  = (double)msel->xSize / (double)mx;

  msel->yTop    = MAX(0, MIN( MIN(List[0].y, List[1].y), my));
  msel->yBottom = MAX(0, MIN( MAX(List[0].y, List[1].y), my));
  msel->ySize   = msel->yBottom - msel->yTop;
  msel->yRatio  = (double)msel->ySize / (double)my;

  if(msel->ySize > 0)
    msel->Aspect  = (double)msel->xSize / (double)msel->ySize;
  else
    msel->Aspect  = 0.;

  msel->Area    = msel->xSize * msel->ySize;
  msel->Radius2 = (msel->xSize)*(msel->xSize) + (msel->ySize)*(msel->ySize);
  msel->xCenter = List[0].x;
  msel->yCenter = List[0].y;

  XSM_DEBUG (DBG_L3, "MkMausSelect: (" << msel->xLeft << ", " << msel->yBottom << ")-(" << msel->xRight  << ", " << msel->yTop << ")"); 
  XSM_DEBUG (DBG_L3, "MkMausSelect Size: (" << msel->xSize << ", " << msel->ySize << ")");
}

//
// als Vorlage für Implementation neuer Filter
//
// Muster Filter "NULL EFFEKT"
//==============================
// MATHOPPARAMS:  
// Scan *Src:  Datenquelle
// Scan *Dest: Datenziel
//
// ===============================
// *Dest ist ein neuer, leerer Scan mit identischen Parametern wie *Src !
// !!!!!!  *Src sollte NIE UND NIMMER hier verändert werden !!!!!!
//
// Größe:
// ===============================
// nx = Dest->mem2d->GetNx(), ... GetNy()
// Dest->Resize(Dest->data.s.nx, Dest->data.s.ny) :  Größe von ZielScan ändern
//
// Datenzugriff über Mem2d Objekt:
// ===============================
// * Kopieren eines rechteckigen Bereichs (Typ Src == Typ Dest):
// Dest->mem2d->CopyFrom(Src->mem2d, int x0, int y0, int tox, int toy, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
// * Kopieren eines rechteckigen Bereichs mit ggf. Typenkonversion:
// Dest->mem2d->ConvertFrom(Src->mem2d, int x0, int y0, int tox, int toy, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
//
//
// [[ Src->mem2d->GetDataLine(int line, ZTYPE *buffer)   :    ganze Zeile ]] please do not use !!
// double value = Src->mem2d->GetDataPkt(int x, int y) :    einen Punkt
// ... oder auch (etwas schneller, da direkter)
// double value = Src->mem2d->data->Z(int x, int y) :    einen Punkt
// value = Src->mem2d->GetDataPkt(int x, int y) :    einen Punkt
// value = Src->mem2d->GetDataPktLineReg(int x, int y) :    einen Punkt, mit Line Regession !!
// [[ Dest->mem2d->PutDataLine(int line, ZTYPE *buffer)  :    ganze Zeile ]] please do not use !!
// Dest->mem2d->PutDataPkt(ZTYPE value, int x, int y) :    einen Punkt
// ... oder auch (etwas schneller, da ohne invalidate lineregress parameters)
// Src->mem2d->data->Z(double value, int x, int y) :    einen Punkt
//
// for large Data Transfer use:
// ============================================================
// for same types:
// inline int CopyFrom(Mem2d *src, int x, int y, int tox, int toy, int nx, int ny=1)
//
// for different types:
// inline int ConvertFrom(Mem2d *src, int x, int y, int tox, int toy, int nx, int ny=1)
//
// for quick linear data access of elements in one line use:
// ============================================================
// void   mem2d->data->SetPtr(int x, int y) to set internal pointer to x,y  -- no range check !!!
// double mem2d->data->GetNext() to access Element x,y and point to next one in line y -- no range check !!!
// double mem2d->data->SetNext(double z) to access Element x,y and point to next one in line y -- no range check !!!
//
/* For Example use constructs like:
  ZData  *SrcZ, *DestZ;

  SrcZ  =  Src->mem2d->data;
  DestZ = Dest->mem2d->data;

  for ( line=0; line < Dest->data.s.ny; line++){ 
    DestZ->SetPtr(0, line);
    SrcZ ->SetPtr(0, line);
    DestZ->SetNext(SrcZ->GetNext() + ....));
  }

  or a advanced example:  => F2D_LineInterpol(MATHOPPARAMS) 

 */
// Parameter:
// ===============================
// Dest->data.s.xxx:  Scanparameterstruktur
// Src->Pkt2d[0..15]: Maus-Punkte

// uses new Object methods now !
gboolean CopyScan(MATHOPPARAMS){
  XSM_DEBUG (DBG_L3, "Copy Scan");
  Dest->mem2d->Resize(Src->mem2d->GetNx(), Src->mem2d->GetNy());

  BenchStart(mbs,"Copy Scan",Dest->mem2d->GetNx()*Dest->mem2d->GetNy());

  Dest->mem2d->ConvertFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());

  BenchStop(mbs);

  return MATH_OK;
}

// Ausschnitt herausschneiden
gboolean CropScan(MATHOPPARAMS){
  MOUSERECT msr;
  XSM_DEBUG (DBG_L3, "Crop Scan");

  MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());
    
  if( msr.xSize  < 1 || msr.ySize < 1){
    XSM_DEBUG (DBG_L3, "Crop:" << msr.xSize << " " << msr.ySize);
    return MATH_SELECTIONERR;
  }

  if(gapp->xsm->MausMode() == MRECT){
    Dest->data.s.nx = msr.xSize;
    Dest->data.s.ny = msr.ySize;
    Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
    Dest->data.s.x0 += (msr.xLeft + msr.xSize/2 - Src->data.s.nx/2)*Src->data.s.dx;
    Dest->data.s.y0 -= msr.yTop*Src->data.s.dy;

    // Anpassen an nächstmöglichen Wert  
    Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
    Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;
    
    Dest->mem2d->CopyFrom(Src->mem2d, msr.xLeft,msr.yTop, 0,0, Dest->data.s.nx, Dest->data.s.ny);

    return MATH_OK;
  }
  else if(gapp->xsm->MausMode() == MCIRC){
    int Radius, xx,x,y,r2,y2;
    int col;
    const SHT bgval=(SHT)gapp->xsm->mradius;    // Get Radius Field as BgVal
    Dest->data.s.nx = 2*(Radius = (int)sqrt((double)msr.Radius2));
    Dest->data.s.ny = 2*Radius;
    Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
    r2 = Radius*Radius;

    // Anpassen an nächstmöglichen Wert  
    Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
    Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;
    
    for(int line = 0; line < Dest->data.s.ny; line++){
      for(col = 0; col < Dest->data.s.nx; col++)
        Dest->mem2d->PutDataPkt(bgval, col, line);

      y=(int)(line+msr.yCenter-Radius);
      y2  = line-Radius;
      y2 *= y2;
      if(y>0 && y<Src->mem2d->GetNy())
        for(col = 0; col < Dest->data.s.nx; col++){
          x   = (int)(col+msr.xCenter-Radius);
          xx  = col-Radius;
          xx *= xx;
          if(xx < (r2-y2) && x >= 0 && x < Src->mem2d->GetNx())
            Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt(x,y), col, line);
        }
    }
    return MATH_OK;
  }
  return MATH_SELECTIONERR;
}

// Größe halbieren, 2x2 Mittelung
// Essential for Gxsm core!!!
gboolean TR_QuenchScan(MATHOPPARAMS){
  long col, line;
  ZData  *SrcZ, *SrcZn, *DestZ;
  double val;
  XSM_DEBUG (DBG_L3, "Quench Scan");

  SrcZ  =  Src->mem2d->data;

  DestZ = Dest->mem2d->data;

  Dest->data.s.nx = (Src->data.s.nx/2) & ~1;
  Dest->data.s.dx = Src->data.s.dx*2.;
  Dest->data.s.ny = (Src->data.s.ny/2) & ~1;
  Dest->data.s.dy = Src->data.s.dy*2.;
  Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);

  Mem2d NextLine(Src->mem2d->GetNx(), 1, Src->mem2d->GetTyp());
  SrcZn = NextLine.data;
  for(line = 0; line < Dest->data.s.ny; line++){
    NextLine.CopyFrom(Src->mem2d, 0,2*line+1, 0,0, Src->mem2d->GetNx());
    DestZ->SetPtr(0, line);
    SrcZ ->SetPtr(0, 2*line);
    SrcZn->SetPtr(0, 0);
    for(col = 0; col < Dest->data.s.nx; col++){
      val  = SrcZ ->GetNext();
      val += SrcZ ->GetNext();
      val += SrcZn->GetNext();
      val += SrcZn->GetNext();
      val /= 4.;
      DestZ->SetNext(val);
    }
  }

  return MATH_OK;
}

// In Ausschnitt hineinzoomen 
gboolean ZoomInScan(MATHOPPARAMS){
  long col, line;
  MOUSERECT msr;
  double FXY, facx, facy;

  MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1 || msr.ySize < 1)
    return MATH_SELECTIONERR;
   
  // Seitenverhältnis und Fläche beibehalten:
  FXY = (double)(Src->data.s.nx*Src->data.s.ny);
  Dest->data.s.nx = (int)(sqrt(FXY * msr.Aspect));
  Dest->data.s.ny = (int)(sqrt(FXY / msr.Aspect));

  Dest->data.s.x0 += (msr.xLeft + msr.xSize/2 - Src->data.s.nx/2)*Src->data.s.dx;
  Dest->data.s.y0 -= msr.yTop*Src->data.s.dy;

  facx = (double)Dest->data.s.nx/msr.xSize;
  facy = (double)Dest->data.s.ny/msr.ySize;
  Dest->data.s.dx /= facx;
  Dest->data.s.dy /= facy;

  /*
  // Anpassen an nächstmöglichen Wert  ... erst bevor Scanning
  Dest->data.s.dx = R2INT(Dest->data.s.dx/Inst->XResolution())*Inst->XResolution();
  Dest->data.s.dy = R2INT(Dest->data.s.dy/Inst->YResolution())*Inst->YResolution();
  */

  Dest->data.s.rx = (Dest->data.s.nx-1)*Dest->data.s.dx;
  Dest->data.s.ry = (Dest->data.s.ny-1)*Dest->data.s.dy;

  Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);

  for(line=0; line<Dest->data.s.ny; line++)
    for(col=0; col<Dest->data.s.nx; col++)
      Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt((int)(msr.xLeft+col/facx),
                                                     (int)(msr.yTop+line/facy)),
                              col, line);
  return MATH_OK;
}

// auf Ausschnitt herauszoomen (*1/3)
gboolean ZoomOutScan(MATHOPPARAMS){
  int col, line;
  int x0,y0;
  XSM_DEBUG(DBG_L2, "ZoomOut" );

  // Seitenverhältnis und Fläche beibehalten:
  Dest->data.s.nx *= 3;
  Dest->data.s.ny *= 3;
  Dest->mem2d->Resize(Dest->data.s.nx, Dest->data.s.ny);
  //  Dest->data.s.x0 += 0;
  Dest->data.s.y0 += Src->data.s.ny*Src->data.s.dy;

  // Anpassen an nächstmöglichen Wert  
  Dest->data.s.rx = Dest->data.s.nx*Dest->data.s.dx;
  Dest->data.s.ry = Dest->data.s.ny*Dest->data.s.dy;

  x0=Src->data.s.nx;
  y0=Src->data.s.ny;
  for(line = 0; line < Src->mem2d->GetNy(); line++)
    for(col = 0; col < Src->mem2d->GetNx(); col++)
      Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPktLineReg(col, line), col+x0, line+y0);
  return MATH_OK;
}

// Copy Scan, do Lin Regress
gboolean BgLin1DScan(MATHOPPARAMS){
  int line, col;
  XSM_DEBUG (DBG_L3, "BgLin1D Scan");

  for(line=0; line<Dest->data.s.ny; line++)
    for(col=0; col<Dest->data.s.nx; col++)
      Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPktLineReg(col, line), col, line);
  return MATH_OK;
}

// Do E Regress
gboolean BgERegress(MATHOPPARAMS){
        double sumi, sumiy, sumyq, sumiq, sumy, val;
        double ysumi, ysumiy, ysumyq, ysumiq, ysumy;
        double a, b, n, imean, ymean,ax,bx,ay,by;
        long line, i;
        MOUSERECT msr;

        XSM_DEBUG (DBG_L3, "BgERegress Scan");

        MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

        if( msr.xSize  < 1 || msr.ySize < 1)
                return MATH_SELECTIONERR;

        sumi = sumiq = sumy = sumyq = sumiy = ax = ay = bx = by = 0.;
        ysumi = ysumiq = ysumy = ysumyq = ysumiy = 0.;
        for (val=msr.xLeft; val < msr.xRight; val += 1.){
                sumiq += val*val;
                sumi  += val;
        }
        n = msr.xSize;
        imean = sumi / n;
        for (line = msr.yTop; line<msr.yBottom; ++line) {
                sumy = sumyq = sumiy = 0.;
                for (i=msr.xLeft; i < msr.xRight; ++i) {
                        val = Src->mem2d->GetDataPkt(i, line);
                        sumy   += val;
                        sumyq  += val*val;
                        sumiy  += val*i;
                }
                ymean = sumy / n;
                a = (sumiy- n * imean * ymean ) / (sumiq - n * imean * imean);
                b = ymean - a * imean;
                ax += a;
                bx += b;
                /* Werte fuer y-linregress */
                val = line;
                ysumy  += b;
                ysumyq += b*b;
                ysumiy += b*val;
                ysumiq += val*val;
                ysumi  += val;
        }
        n = msr.ySize;
        ax = ax/n;
        bx = bx/n;
        imean = ysumi / n;
        ymean = ysumy/  n;
        ay = (ysumiy- n * imean * ymean ) / (ysumiq - n * imean * imean);
        by = ymean - ay * imean;
        sumy = sumi = 0.;

        for (line = 0; line<Dest->data.s.ny; ++line, sumy += ay) {
                sumi = by + sumy + .5;
                for (i=0; i < Dest->data.s.nx; ++i, sumi += ax)
                        Dest->mem2d->PutDataPkt (Src->mem2d->GetDataPkt (i, line) - sumi,
                                                 i, line);
        }
        return MATH_OK;
}

// Do Parabol Regress
gboolean BgParabolRegress(MATHOPPARAMS){
        double sumi, sumiy, sumyq, sumiq, sumy, val;
        double ysumi, ysumiy, ysumyq, ysumiq, ysumy;
        double a, b, n, imean, ymean,ax,bx,ay,by;
        int line, i;
        MOUSERECT msr;

        XSM_DEBUG (DBG_L3, "BgParabolRegress Scan");

        MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

        if( msr.xSize  < 1 || msr.ySize < 1)
                return MATH_SELECTIONERR;


        // plane regress first

        sumi = sumiq = sumy = sumyq = sumiy = ax = ay = bx = by = 0.;
        ysumi = ysumiq = ysumy = ysumyq = ysumiy = 0.;

        for (val=msr.xLeft; val < msr.xRight; val += 1.){
                sumiq += val*val;
                sumi  += val;
        }
        n = msr.xSize;
        imean = sumi / n;
        for (line = msr.yTop; line<msr.yBottom; ++line) {
                sumy = sumyq = sumiy = 0.;
                for (i=msr.xLeft; i < msr.xRight; ++i) {
                        val = Src->mem2d->GetDataPkt(i, line);
                        sumy   += val;
                        sumyq  += val*val;
                        sumiy  += val*i;
                }
                ymean = sumy / n;
                a = (sumiy- n * imean * ymean ) / (sumiq - n * imean * imean);
                b = ymean - a * imean;
                ax += a;
                bx += b;
                /* Werte fuer y-linregress */
                val = line;
                ysumy  += b;
                ysumyq += b*b;
                ysumiy += b*val;
                ysumiq += val*val;
                ysumi  += val;
        }
        n = msr.ySize;
        ax = ax/n;
        bx = bx/n;
        imean = ysumi / n;
        ymean = ysumy/  n;
        ay = (ysumiy- n * imean * ymean ) / (ysumiq - n * imean * imean);
        by = ymean - ay * imean;
        sumy = sumi = 0.;

        for (line = 0; line<Dest->data.s.ny; ++line, sumy += ay) {
                sumi = by + sumy + .5;
                for (i=0; i < Dest->data.s.nx; ++i, sumi += ax)
                        Dest->mem2d->PutDataPkt (Src->mem2d->GetDataPkt (i, line) + sumi,
                                                 i, line);
        }

        // now evaluate 2nd order terms
        double M0, M1, M2, li;
        double D, N;
        double qx0, qx, qxx, qy0, qy, qyy, lb;

        // avg. 2nd order terms in X
        N = (double)Src->mem2d->GetNx();
        D = (N*N-1.0)*(N*N-4.0);
        qx0 = qx = qxx = 0.;
        for (int line=0; line < Src->mem2d->GetNy(); line++) {
                M0 = 0.0; M1 = 0.0; M2 = 0.0;

                for (int col = 0; col < Src->mem2d->GetNx(); col++) {
                        li = (double)col;
                        M0 += lb = Src->mem2d->GetDataPkt (col, line);
                        M1 += lb * li;
                        M2 += lb * li*li;
                }

                M0 /= N;
                M1 /= N;
                M2 /= N;

                qx0 += (3.*(N+1.)*(N+2.)/D)*((3.*N*N+3.*N+2.)*M0-6*(2*N+1.)*M1+10.*M2);
                qx  += -6./D*(3*(N+1.)*(N+2.)*(2.*N+1.)*M0-2.*(2.*N+1.)*(8.*N+11.)*M1+30.*(N+1.)*M2);
                qxx += 30./D*((N+1.)*(N+2.)*M0-6.*(N+1.)*M1+6.*M2);
        }
        qx0 /= (double) Src->mem2d->GetNy();
        qx  /= (double) Src->mem2d->GetNy();
        qxx /= (double) Src->mem2d->GetNy();
         
        // avg. 2nd order terms in Y
        // swaped line/col here
        N = (double)Src->mem2d->GetNy();
        D = (N*N-1.0)*(N*N-4.0);
        qy0 = qy = qyy = 0.;
        for (int line=0; line < Src->mem2d->GetNx(); line++) {
                M0 = 0.0; M1 = 0.0; M2 = 0.0;

                for (int col = 0; col < Src->mem2d->GetNy(); col++) {
                        li = (double)col;
                        M0 += lb = Src->mem2d->GetDataPkt (line, col);
                        M1 += lb * li;
                        M2 += lb * li*li;
                }

                M0 /= N;
                M1 /= N;
                M2 /= N;

                qy0 += (3.*(N+1.)*(N+2.)/D)*((3.*N*N+3.*N+2.)*M0-6*(2*N+1.)*M1+10.*M2);
                qy  += -6./D*(3*(N+1.)*(N+2.)*(2.*N+1.)*M0-2.*(2.*N+1.)*(8.*N+11.)*M1+30.*(N+1.)*M2);
                qyy += 30./D*((N+1.)*(N+2.)*M0-6.*(N+1.)*M1+6.*M2);
        }
        qy0 /= (double) Src->mem2d->GetNx();
        qy  /= (double) Src->mem2d->GetNx();
        qyy /= (double) Src->mem2d->GetNx();
        

        for (line = 0; line<Dest->data.s.ny; ++line)
                for (i=0; i < Dest->data.s.nx; ++i, sumi += ax){
                        sumi = qx0 + qy0 + qx*i + qy*line + qxx*i*i + qyy*line*line;
                        Dest->mem2d->PutDataPkt (Dest->mem2d->GetDataPkt (i, line) - sumi, i, line);
        }

        return MATH_OK;
}

//======================================================================================
//
// 1D Filterfunktionen 
//
//======================================================================================


gboolean F1D_Despike(MATHOPPARAMS){
  int i;
  int n=Dest->mem2d->GetNx();
  int m=Dest->mem2d->GetNy();

  XSM_DEBUG (DBG_L3, "F1D Despike");

  Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());
  for(i=1;i<m;i++){
    int q, j;
    double avg, avgdev;
    
    {
      double sum=0.;
      double absdevsum=0.;
      double diff;
      for(j=0;j<n;j++){
        diff = Src->mem2d->GetDataPkt(j,i) - Src->mem2d->GetDataPkt(j,i-1);
        sum += diff;
        absdevsum += fabs(diff);
      }
      avg = sum/n;
      avgdev = absdevsum/(n-1);
    }
    for(q=0;q<4;q++)    {
      double sum, dsum, na;
      for(na=dsum=sum=0.0,j=0;j<n;j++)  {
        const double sdev = (Src->mem2d->GetDataPkt(j,i) - Src->mem2d->GetDataPkt(j,i-1) - avg)/avgdev;
        const double k = 1/(3+sdev*sdev);
        sum += k*sdev;
        dsum += k*sdev*sdev;
        na += k;
      }
      avg += avgdev*sum/na;
      avgdev = 0.2*avgdev + 0.8*sqrt(dsum/(na-1));
    }
    
    {
      SHT iavg = (SHT)(rint(avg));
      for(j=0;j<n;j++)
        Dest->mem2d->PutDataPkt(Src->mem2d->GetDataPkt(j,i) - iavg, j,i);
    }
  }
  return MATH_OK;
}


// compute 1D power spectrum (row by row)
gboolean F1D_LogPowerSpec(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F1D LogPowerSpec");

  ZData  *SrcZ, *DestZ;

  SrcZ  =  Src->mem2d->data;
  DestZ = Dest->mem2d->data;

  XSM_DEBUG (DBG_L3, "F1D LogPowerSpec: using libfftw");

  // get memory for complex data
  double *in  = new double [Src->mem2d->GetNx()];
  fftw_complex *out = new fftw_complex [Src->mem2d->GetNx()];

  // create plan for fft
  fftw_plan plan = fftw_plan_dft_r2c_1d (Src->mem2d->GetNx(), in, out, FFTW_ESTIMATE);
  if (plan == NULL) {
    XSM_DEBUG (DBG_L3, "F1D LogPowerSpec: libfftw: Error creating plan");
    return MATH_LIB_ERR;
  }

  Dest->mem2d->Resize (Src->mem2d->GetNx()/2, Src->mem2d->GetNy(), ZD_DOUBLE);
  Dest->data.s.nx = Dest->mem2d->GetNx();
  Dest->data.s.x0 = 0.;
  Dest->data.s.y0 = 0.;
  Dest->data.s.rx = Src->data.s.rx/2.;
  Dest->data.s.ry = Src->data.s.ry;
  Dest->data.s.dl = Src->data.s.rx/2.;
  Dest->data.s.dz = Src->data.s.dz;
  Dest->data.s.dy = Src->data.s.dy;
 
  // compute 1D fourier transform for every row
  for (int line = 0; line < Src->mem2d->GetNy(); line++) {
    SrcZ ->SetPtr(0, line);
    
    // prepare image row data for fftw
    for (int col = 0; col < Src->mem2d->GetNx(); col++)
       in[col] = SrcZ->GetNext();

    // compute transform
    fftw_execute (plan);

    // convert complex data to image row
    // scale data to 16Bit integer
    // reorder data that freq 0 is in the middle and positive freq go to the right

    DestZ->SetPtr(0, line);

    for (int col = Src->mem2d->GetNx()/2-1; col >= 0; --col)
            DestZ->SetNext ( ZEROVALUE + ( c_re(out[col])*c_re(out[col]) + 
                                           c_im(out[col])*c_im(out[col]) ));
  }

  int n = Dest->mem2d->GetNx();
  double s = Src->data.s.rx/n;
  for (int i=0; i<n; ++i)
    Dest->mem2d->data->SetXLookup(i, (n-i)*s);
    
  // delete plan:
  fftw_destroy_plan (plan);

  // free temp data memory
  delete in;
  delete out;

  return MATH_OK;
}

// 2D Filterfunktionen

// remove bad area
gboolean F2D_RemoveRect(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D Remove Rect");

  // get selection
  MOUSERECT msr;
  MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  // check for valid selection
  if( msr.xSize  < 1 || msr.yBottom > Src->mem2d->GetNy() || msr.yTop < 0)
    return MATH_SELECTIONERR;

  // resize Dest
  Dest->mem2d->Resize(Src->mem2d->GetNx(), Src->mem2d->GetNy()-msr.ySize);

  // copy data except selection
  for (int line = 0; line < msr.yTop; line++)
    Dest->mem2d->CopyFrom(Src->mem2d, 0,line, 0, line, Dest->mem2d->GetNx());

  for (int line = msr.yBottom+1; line < Src->mem2d->GetNy(); line++)
    Dest->mem2d->CopyFrom(Src->mem2d, 0,line, 0, line-msr.ySize-1, Dest->mem2d->GetNx());

  return MATH_OK;
}


// Do 2D DeSpike
gboolean F2D_Despike(MATHOPPARAMS){
  int i,j,k,l,za=0,nx,ny;
  int anz=1;
  int num=1;
  double reihe1[20],reihe2[20],mi;

  XSM_DEBUG (DBG_L3, "F2D Despike");

  BenchStart(mbs,"2D Despike",Dest->mem2d->GetNx()*Dest->mem2d->GetNy());
  
  Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, 
                        nx=Dest->mem2d->GetNx(),ny=Dest->mem2d->GetNy());

  for(j=anz; j<ny-anz; ++j)  {
    for(i=anz; i<nx-anz; ++i) {
      for(k=j-anz, l=0; k<=j+anz; k++,l++)
        reihe1[l] = Src->mem2d->data->Z(i,k);
      for(k=0; k<2*anz+1;k++)  {
        mi = 32650;
        for(l=0; l<2*anz+1;l++) {
          if(reihe1[l]<mi) {  mi=reihe1[l]; reihe2[k]=mi; za=l; }
          
        }
        reihe1[za]=32650;
      }
      Dest->mem2d->data->Z(reihe2[num], i,j);
    }  /*i*/
  } /*j*/
  BenchStop(mbs);

  return MATH_OK;
}



// remove line shifts
// uses new Object methods now !
extern gboolean F2D_LineShifts(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D LineShifts");
  double avg;
  int nx,ny;
  int x1,x2,y1,y2;
  ZData  *SrcZ, *DestZ;
  
  nx=Dest->mem2d->GetNx();
  ny=Dest->mem2d->GetNy();
  BenchStart(mbs,"F2D LineShifts",nx*ny);

  x1=0; x2=nx; y1=y2=0;

  if (Src->number_of_object () > 0){
          XSM_DEBUG(DBG_L2, "found object(s)" );
          if( Src->get_object_data (0) -> get_num_points () == 2){
                  double x=0.;
                  double y=0.;
                  Src->get_object_data (0) -> get_xy (0, x, y);
                  Src->World2Pixel (x,y, x1,y1);
                  Src->get_object_data (0) -> get_xy (1, x, y);
                  Src->World2Pixel (x,y, x2,y2);
                  XSM_DEBUG(DBG_L2, "got x range " << x1 << " .. " << x2 );
          }

          if (x1 > x2){
                  int tmp = x1;
                  x1 = x2;
                  x2 = tmp;
          }
  }

  Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, nx,ny);
  SrcZ  =  Src->mem2d->data;
  DestZ = Dest->mem2d->data;
  for (int line = 0; line < ny; line++) {
    avg = 0.;
    SrcZ->SetPtr(0,line);

    for (int col = x1; col < x2; col++)
      avg += SrcZ->GetNext();

    avg /= (x2-x1);
    SrcZ->SetPtr(0,line);
    DestZ->SetPtr(0,line);

    for (int col = 0; col < nx; col++)
      DestZ->SetNext(SrcZ->GetNext() - avg);
  }
  BenchStop(mbs);

  return MATH_OK;
}



// "fix" bad line
// using an interpolation between the neighbour lines
// uses new Object methods now !
extern gboolean F2D_LineInterpol(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D LineInterpol");

  MOUSERECT msr;

  MkMausSelect(Src->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1)
    return MATH_SELECTIONERR;

  // copy scan data:
  Dest->mem2d->CopyFrom(Src->mem2d, 0,0, 0,0, Dest->mem2d->GetNx(),Dest->mem2d->GetNy());

  // interpolation if in between first and last line
  if(msr.yTop > 0){
    if(msr.yTop < (Src->mem2d->GetNy()-1)){
      // create help Mem2d object (only one line)
      Mem2d PreLine(Src->mem2d->GetNx(), 1, Src->mem2d->GetTyp());
      // with "PreLine" and set internal Ptr to start
      PreLine.CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,0, Src->mem2d->GetNx());
      PreLine.data->SetPtr(0,0);
      // set Src Ptr to "PostLine"
      Src->mem2d->data->SetPtr(0,msr.yTop+1);
      // set Dest Ptr to "Line"
      Dest->mem2d->data->SetPtr(0,msr.yTop);
      for(int col = 0; col < Src->mem2d->GetNx(); col++)
        Dest->mem2d->data->SetNext((PreLine.data->GetNext()+Src->mem2d->data->GetNext())/2.);
    }
    else  // last line: copy previous line
      Dest->mem2d->CopyFrom(Src->mem2d, 0,msr.yTop-1, 0,msr.yTop, Src->mem2d->GetNx());
  }
  else  // first line: copy 2nd line
    Dest->mem2d->CopyFrom(Src->mem2d, 0,msr.yTop+1, 0,msr.yTop, Src->mem2d->GetNx());

  return MATH_OK;
}




// 2D power spectrum in logarithmic scale
gboolean F2D_LogPowerSpec(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D FFT");
  ZData  *SrcZ;

  SrcZ  =  Src->mem2d->data;
  Dest->mem2d->Resize (Dest->data.s.nx, Dest->data.s.ny, ZD_COMPLEX);
  // set "Complex" layer param defaults
  Dest->data.s.nvalues=3;
  Dest->data.s.ntimes=1;
  Dest->mem2d->data->SetVLookup(0, 0);
  Dest->mem2d->data->SetVLookup(1, 1);

  XSM_DEBUG (DBG_L3, "F2D FFT: using libfftw");

  // allocate memory for real and inplace complex data
  int xlen=2*(Src->mem2d->GetNx ()/2+1);
  double *in  = new double [ Src->mem2d->GetNx () * Src->mem2d->GetNy () ];
  fftw_complex *dat  = new fftw_complex [ xlen*Src->mem2d->GetNy () ];

  if (dat == NULL) {
    XSM_DEBUG (DBG_L3, "F2D FFT: libfftw: Error allocating memory for complex data");
    return MATH_NOMEM;
  }

  memset(in, 0, sizeof(in));

  // create plan for in-place transform
  fftw_plan plan    = fftw_plan_dft_r2c_2d (Src->mem2d->GetNy (), Src->mem2d->GetNx (), 
                                            in, dat, FFTW_ESTIMATE);

  if (plan == NULL) {
    XSM_DEBUG (DBG_L3, "F2D FFT: libfftw: Error creating plan");
    return MATH_LIB_ERR;
  }

  // convert image data to fftw_real
  for (int line=0; line < Src->mem2d->GetNy(); line++) {
          SrcZ ->SetPtr(0, line);
          for (int col=0; col < Src->mem2d->GetNx(); col++)
                  in[line*Src->mem2d->GetNx() + col] = SrcZ->GetNext();
  }

  // compute 2D transform using in-place fourier transform
  fftw_execute (plan);

  // convert complex data to image
  // and flip image that spatial freq. (0,0) is in the middle
  // (quadrant swap using the QSWP macro defined in xsmmath.h)
  
  XSM_DEBUG (DBG_L3, "F2D FFT done, converting data to ABS/RE/IM");

  xlen /= 2; // jetzt für complex

  double I;
  int m,n;
 
  for (int line=0; line < Src->mem2d->GetNy(); line++) {
          for (int i,col=0; col < Src->mem2d->GetNx()/2; col++) {
                  i=line*xlen + col;
                  m = QSWP(col, Src->mem2d->GetNx());
                  n = QSWP(line, Src->mem2d->GetNy());
                  I = ZEROVALUE + sqrt (c_re(dat[i]) * c_re(dat[i]) +  c_im(dat[i]) * c_im(dat[i]));
                  Dest->mem2d->PutDataPkt(I, m, n, 0);
                  Dest->mem2d->PutDataPkt(c_re(dat[i]), m, n, 1);
                  Dest->mem2d->PutDataPkt(c_im(dat[i]), m, n, 2);
                  if(line != (Src->mem2d->GetNy()/2)){
                          m = Src->mem2d->GetNx() - m;
                          n = Src->mem2d->GetNy() - n;
                          Dest->mem2d->PutDataPkt(I, m, n, 0);
                          Dest->mem2d->PutDataPkt(c_re(dat[i]), m, n, 1); 
                          Dest->mem2d->PutDataPkt(c_im(dat[i]), m, n, 2);  
                  }
          }
  }
  
  // destroy plan
  fftw_destroy_plan(plan); 

  // free real/complex data memory
  delete in;
  delete dat;

  return MATH_OK;
}


// Echte 1D Fourier Filter :=)

// FFTW-Kernel, applys some function to the spectra
// Forward FT => do something with spectra => backward FT
gboolean F1D_ift_ft(MATH2OPPARAMS,
                    gboolean (*spkfkt)(MATH2OPPARAMS,
                                       fftw_complex *dat, 
                                       int line)
        ){
  XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) )");
  gboolean ret=MATH_OK;
  ZData  *SrcZ, *DestZ;

  SrcZ  = Src1->mem2d->data;
  DestZ = Dest->mem2d->data;

  // error reporting for invalid selections:
  // Src1 and Src2 scans must have the same dimensions
  // if Src2 == NULL filter spkfkt needs no selection, see e.g. SpkAutoCorr
  if(Src2 != NULL && (Src1->data.s.nx != Src2->data.s.nx || Src1->data.s.ny != Src2->data.s.ny))
    return MATH_SELECTIONERR;

  // ======== 1. fourier transform Src1 => complex data spectrum =======

  int line, col;

  // allocate memory for complex data
  double  *in  = new double [ Src1->mem2d->GetNx() ];
  fftw_complex *dat = new fftw_complex [ Src1->mem2d->GetNx() ];
  double scale;

  fftw_plan plan    = fftw_plan_dft_r2c_1d (Src1->mem2d->GetNx(), in, dat, FFTW_ESTIMATE);
  fftw_plan planinv = fftw_plan_dft_c2r_1d (Src1->mem2d->GetNx(), dat, in, FFTW_ESTIMATE);

  if (!plan || !planinv) {
    XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error creating plans");
    if(plan)
      fftw_destroy_plan (plan); 
    if(planinv)
      fftw_destroy_plan (planinv); 
    return MATH_LIB_ERR;
  }



  if (!in || !dat) {
    XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error allocating memory for complex data");
    if(in)  delete in;
    if(dat) delete dat;
    return MATH_NOMEM;
  }

  memset(in, 0, sizeof(in));

  // convert image data to fftw_real
  for (line=0; line < Src1->mem2d->GetNy(); line++) {
          SrcZ ->SetPtr(0, line);
          for (col=0; col < Src1->mem2d->GetNx(); col++)
                  in[col] = SrcZ->GetNext();
    
          // 1. compute 1D transform using in-place fourier transform
          fftw_execute (plan);
          
          // now the complex spectrum is in dat too !
          
          // 2. apply something to the spektra
          if (spkfkt (MATH2OPVARS, dat, line) == MATH_OK){
      
                  // 3. inverse fourier transform modified spectrum
      
                  fftw_execute (planinv);
      
                  // 4. take the real part as data and apply the fft normalization factor scale
      
                  scale = 1. / Src1->mem2d->GetNx();
                  DestZ->SetPtr(0, line);
                  for(col=0; col<Dest->mem2d->GetNx(); col++)
                          DestZ->SetNext(scale*(in[col]));
          }
          else
                  ret=MATH_SELECTIONERR;
  }

  // free [complex] data memory
  delete in;

  // destroy plan
  fftw_destroy_plan (plan); 
  fftw_destroy_plan (planinv); 
  
  return ret;
}


// Filter Kernel 1D "Window"
gboolean SpkWindow1D(MATH2OPPARAMS, fftw_complex *dat, int line){
  int col;
  // ====== 2. multiply complex data spectrum with Src2, eg. windowing =====
  // and take into account that spatial freq. (0,0) is in the middle

  // apply Window

  int n    = Src1->mem2d->GetNx();
  for (col=0; col < n/2; col++)
    if (Src2->mem2d->GetDataPkt(n/2-col-1, line) == ZEROVALUE){
            c_re (dat[col]) = 0.;
            c_im (dat[col]) = 0.;
            c_re (dat[n-col-1]) = 0.;
            c_im (dat[n-col-1]) = 0.;
    }
  return MATH_OK;
}


// Filter Kernel "Gauss Stopp"
gboolean SpkGaussStop1D(MATH2OPPARAMS, fftw_complex *dat, int line){
  MOUSERECT msr;
  MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1 || msr.ySize < 1)
    return MATH_SELECTIONERR;

  // apply gaussian window to complex data
  // scale = 1 - exp(-r2/msr.Radius2)
  double r, scale;
  int n = Src1->mem2d->GetNx();
  for (int col = 0; col < n/2; col++) {
    r = (n/2-col-1) - msr.xCenter; r *= r;
    scale = 1 - exp(-r / msr.Radius2);
    c_re (dat[col])     *= scale;
    c_re (dat[n-col-1]) *= scale;
  }
  return MATH_OK;
}

// Filter Kernel "Gauss Pass"
gboolean SpkGaussPass1D(MATH2OPPARAMS, fftw_complex *dat, int line){
  MOUSERECT msr;
  MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1 || msr.ySize < 1)
    return MATH_SELECTIONERR;

  // apply gaussian window to complex data
  // scale = exp(-r2/msr.Radius2)
  double r, scale;
  int n = Src1->mem2d->GetNx();
  for (int col = 0; col < n/2; col++) {
    r = (n/2-col-1) - msr.xCenter; r *= r;
    scale = exp(-r / msr.Radius2);
    c_re (dat[col])   *= scale;
    c_re (dat[n-col-1]) *= scale;
  }
  return MATH_OK;
}


// Fourierfilter:
// Pass frequency window with gaussian profile
gboolean F1D_FT_Window(MATH2OPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F1D FT Window");
  return(F1D_ift_ft(MATH2OPVARS, SpkWindow1D));
  return MATH_OK;
}

gboolean F1D_FT_GaussStop(MATH2OPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F1D FT GaussStop");
  return(F1D_ift_ft(MATH2OPVARS, SpkGaussStop1D));
  return MATH_OK;
}

gboolean F1D_FT_GaussPass(MATH2OPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F1D FT GaussPass");
  return(F1D_ift_ft(MATH2OPVARS, SpkGaussPass1D));
  return MATH_OK;
}


// Echte 2D Fourier Filter :=)

// FFTW-Kernel, applys some function to the spectra
// Forward FT => do something with spectra => backward FT
gboolean F2D_ift_ft(MATH2OPPARAMS, gboolean (*spkfkt)(MATH2OPPARAMS, fftw_complex *dat)){
  XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) )");
  gboolean ret=MATH_OK;
  ZData  *SrcZ, *DestZ;

  SrcZ  = Src1->mem2d->data;
  DestZ = Dest->mem2d->data;

  // error reporting for invalid selections:
  // Src1 and Src2 scans must have the same dimensions
  // if Src2 == NULL filter spkfkt needs no selection, see e.g. SpkAutoCorr
  if(Src2 != NULL && (Src1->data.s.nx != Src2->data.s.nx || Src1->data.s.ny != Src2->data.s.ny))
    return MATH_SELECTIONERR;

  // ======== 1. fourier transform Src1 => complex data spectrum =======

  int line, col;

  // allocate memory for complex data
  int xlen=2*(Src1->mem2d->GetNx()/2+1);
  double    *in  = new double [ xlen*Src1->mem2d->GetNy() ];
  fftw_complex *dat = (fftw_complex*)(in); // use same memory - saves a lot space...
  double scale;

  if (!in || !dat) {
    XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error allocating memory for complex data");
    if(in)  delete in;
    if(dat) delete dat;
    return MATH_NOMEM;
  }

  // create plan for in-place transform
  fftw_plan plan    = fftw_plan_dft_r2c_2d (Src1->mem2d->GetNy(), Src1->mem2d->GetNx(), in, dat, FFTW_ESTIMATE);
  fftw_plan planinv = fftw_plan_dft_c2r_2d (Src1->mem2d->GetNy(), Src1->mem2d->GetNx(), dat, in, FFTW_ESTIMATE);

  if (!plan || !planinv) {
          XSM_DEBUG (DBG_L3, "iftXft: libfftw: Error creating plans");
          if(plan)
                  fftw_destroy_plan (plan); 
          if(planinv)
                  fftw_destroy_plan (planinv); 
          return MATH_LIB_ERR;
  }

  memset(in, 0, sizeof(in));

  // convert image data to fftw_real
  for (line=0; line < Src1->mem2d->GetNy(); line++) {
          SrcZ ->SetPtr(0, line);
          for (col=0; col < Src1->mem2d->GetNx(); col++)
                  in[line*xlen + col] = SrcZ->GetNext();
  }

  // 1. compute 2D transform using in-place fourier transform

  fftw_execute (plan);
  // now the complex spectrum is in dat too !

  // 2. apply something to the spektra
  if (spkfkt (MATH2OPVARS, dat) == MATH_OK){

    // 3. inverse fourier transform modified spectrum

    fftw_execute (planinv);

    // 4. take the real part as data and apply the fft normalization factor scale

    scale = 1. / Src1->mem2d->GetNx() / Src1->mem2d->GetNy();
    for(line=0; line<Dest->mem2d->GetNy(); line++){
            DestZ->SetPtr(0, line);
            for(col=0; col<Dest->mem2d->GetNx(); col++)
                    DestZ->SetNext(scale*(in[line*xlen + col]));
    }
  }
  else
    ret=MATH_SELECTIONERR;

  // free [complex] data memory
  delete in;

  // destroy plan
  fftw_destroy_plan (plan); 
  fftw_destroy_plan (planinv); 

  return ret;
}

// Filter Kernel "Window"
gboolean SpkWindow(MATH2OPPARAMS, fftw_complex *dat){
  int line, col;
  // ====== 2. multiply complex data spectrum with Src2, eg. windowing =====
  // and take into account that spatial freq. (0,0) is in the middle

  // apply Window

  int xlen = Src1->mem2d->GetNx()/2+1;
  for (line=0; line < Src1->mem2d->GetNy(); line++)
    for (col=0; col < xlen; col++)
      if(Src2->mem2d->GetDataPkt(QSWP(col,Src1->mem2d->GetNx()), 
                                 QSWP(line,Src1->mem2d->GetNy()))
         == ZEROVALUE){
        c_re(dat[line*xlen + col]) = 0.;
        c_im(dat[line*xlen + col]) = 0.;
      }
  return MATH_OK;
}


// Src1 : Daten Quelle
// Src2 : Fensterfunktion
//        kann ggf. obiges PowerSpektrum sein, wobei
//        folgendes gilt:
//        Dest = Re(IFT( FT(Src1) * (Src2 == ZEROVALUE ? ZEROVALUE : 1) ))
//        d.h. alle Spektralanteile, die in Src2 Null sind werden im Complexen Spektrum ausgeblendet !
//        Dazu muessen natuerlich vorher im Power-Spektrum entspechende Bereiche auf Null gesetzt werden :=)
//        --- that´s it --- PZ
gboolean F2D_iftXft(MATH2OPPARAMS){
  XSM_DEBUG (DBG_L3, "Doing iftXft: RE( ift(X * FT(Active)) )  -- Window(s)");
  return(F2D_ift_ft(MATH2OPVARS, SpkWindow));
  return MATH_OK;
}


// Filter Kernel "Gauss Stopp"
gboolean SpkGaussStop(MATH2OPPARAMS, fftw_complex *dat){
  MOUSERECT msr;
  MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1 || msr.ySize < 1)
    return MATH_SELECTIONERR;

  // apply gaussian window to complex data
  // scale = 1 - exp(-r2/msr.Radius2)
  double r2, scale;
  int xlen = Src1->mem2d->GetNx()/2+1;
  for (int line = 0; line < Src1->mem2d->GetNy(); line++)
    for (int col = 0; col < xlen; col++) {
      r2 = (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) * (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) + 
        (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter) * (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter);
      scale = 1 - exp(-r2 / msr.Radius2);
      c_re(dat[line*xlen + col]) *= scale;
      c_im(dat[line*xlen + col]) *= scale;
      }
  return MATH_OK;
}

// Fourierfilter:
// Pass frequency window with gaussian profile
gboolean F2D_FT_GaussStop(MATH2OPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D FT GaussStop");
  return(F2D_ift_ft(MATH2OPVARS, SpkGaussStop));
  return MATH_OK;
}


// Filter Kernel "Gauss Pass"
gboolean SpkGaussPass(MATH2OPPARAMS, fftw_complex *dat){
  MOUSERECT msr;
  MkMausSelect(Src2->Pkt2d, &msr, Dest->mem2d->GetNx(), Dest->mem2d->GetNy());

  if( msr.xSize  < 1 || msr.ySize < 1)
    return MATH_SELECTIONERR;

  // apply gaussian window to complex data
  // scale = exp(-r2/msr.Radius2)
  double r2, scale;
  int xlen = Src1->mem2d->GetNx()/2+1;
  for (int line = 0; line < Src1->mem2d->GetNy(); line++)
    for (int col = 0; col < xlen; col++) {
      r2 = (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) * (QSWP(line,Src1->mem2d->GetNy()) - msr.yCenter) + 
        (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter) * (QSWP(col,Src1->mem2d->GetNx()) - msr.xCenter);
      scale = exp(-r2 / msr.Radius2);
      c_re(dat[line*xlen + col]) *= scale;
      c_im(dat[line*xlen + col]) *= scale;
      }
  return MATH_OK;
}


// Fourierfilter:
// Pass frequency window with gaussian profile
gboolean F2D_FT_GaussPass(MATH2OPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D FT GaussPass");
  return(F2D_ift_ft(MATH2OPVARS, SpkGaussPass));
}


// Filter Kernel "Abs. Square"
// (for auto correlation)
gboolean SpkAutoCorr(MATH2OPPARAMS, fftw_complex *dat)
{
  // take abs. value of complex data
  int xlen = Src1->mem2d->GetNx()/2+1;
  for (int line = 0; line < Src1->mem2d->GetNy(); line++)
    for (int col = 0; col < xlen; col++) {
      c_re(dat[line*xlen + col]) = sqrt( c_re(dat[line*xlen + col])*c_re(dat[line*xlen + col]) +
                                         c_im(dat[line*xlen + col])*c_im(dat[line*xlen + col]) );
      c_im(dat[line*xlen + col]) = 0.0;
      }


  return MATH_OK;
}


// auto correlation
gboolean F2D_AutoCorr(MATHOPPARAMS)
{
  XSM_DEBUG (DBG_L3, "F2D AutoCorr");
  return(F2D_ift_ft(Src, NULL, Dest, SpkAutoCorr));
}

---