Filter Library

Camera

Interface Physics

ScanFilterFFT.cpp

/*
 * ScanFilterFFT.cpp - forward FFT.
 *
 * This file is part of the Camera Filter Library.
 * Computer Aided Measurement Environment for Realtime Atomic imaging (Camera)
 *
 * Copyright (C) 2004-2005, Leiden Probe Microscopy.
 * Copyright (C) 2004-2005, Universiteit Leiden.
 *
 * Authors: Martin J. Moene
 *
 * $Id: ScanFilterFFT.cpp 400 2006-09-28 13:34:02Z moene $
 */

/*
 * Note: adding anything before stdafx has no effect>
 */
#include "stdafx.h"                     // for common (pre-compiled) headers
//#include <cfl/stdafx.h>                 // for common (pre-compiled) headers
#include <cfl/resource.h>               // for messages and other resources

#include <Camera/InterfaceDll.h>        // Camera--Filter Library interface

/*
 * configuration defines:
 */
// #define ScanFilter_Use_FFTW          // use FFTW routines; NOTE: this is not yet in order (MO 20060719)
// #define USE_SEPARATE_AMPLITUDE_SCALE // use separate scaling factor for each power spectrum (L2R, R2L)

/*
 * end of configuration defines.
 */

#include <cfl/ScanFilterFFT.h>          // this filter's header
#include <cfl/ToolCilImage.h>           // for Image type and ToImage conversion shim
#include <cil/Math.h>                   // for real(), //min(), max(), sqrt(), sin(), Pi

#if defined( ScanFilter_Use_FFTW )
#   include <cil/FFTWImageFilter.h>             // FFTW image filter
#   include <cil/FFTWInverseImageFilter.h>      // FFTW inverse image filter
#else
#   include <cil/FFTImageFilter.h>              // FFT image filter
#   include <cil/FFTInverseImageFilter.h>       // FFT inverse image filter
#endif // defined( ScanFilter_Use_FFTW )

#include <cil/FFTPowerSpectrumImageFilter.h>    // FFT power spectrum image filter
#include <cil/FFTSwapQuadrantsImageFilter.h>    // FFT swap quadrants image filter
#include <cil/ShiftScaleImageFilter.h>          // FFT scale intensity image filter

// #pragma warning( disable : 4663 )       // template<> Class<type>
// #include <limits>                       // for numeric_limits<>
// #pragma warning( default : 4663 )

#undef min                              // use cil::min instead of macro
#undef max                              // use cil::max instead of macro

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

#define CFL_SCHEMAVERSION_FILTERFFT 0   // permanent storage schema version
                                        // this filter's name
LPCTSTR CScanFilterFFT::m_lpcsFilterName = _T("Fast-Fourier-Transform");
                                        // this filter's short name
LPCTSTR CScanFilterFFT::m_lpcsShortFilterName = _T( "FFT" );
                                        // the maximum amplitude to scale power spectrum to
const CScanFilterFFT::ValueType CScanFilterFFT::maximumAmplitude = 20000;
                                        // serialization
IMPLEMENT_SERIAL(CScanFilterFFT, CScanFilter, CFL_SCHEMAVERSION_FILTERFFT )

/**
 * constructor.
 */
CScanFilterFFT::CScanFilterFFT()
   : m_progress( 0 )
{
   m_csFilterName = m_lpcsFilterName;
}

/**
 * destructor.
 */
CScanFilterFFT::~CScanFilterFFT()
{
   ; // do nothing
}

/**
 * no dialog.
 */
BOOL CScanFilterFFT::CanDoDialogEntry() const
{
   return FALSE;
}

/**
 * store or retrieve the object's settings;
 * see also CScanFilterNull::Serialize().
 */
void CScanFilterFFT::Serialize(CArchive& ar)
{
   CScanFilter::Serialize( ar );

   if ( ar.IsStoring() )
   {
      ; // do nothing
   }
   else
   {
      ; // do nothing
   }
};

/**
 * no parameters.
 */
#pragma warning( push )
#pragma warning( disable : 4100 )       // unreferenced formal parameter

BOOL CScanFilterFFT::SetParameters(LPCTSTR lpParameters)
{
   return TRUE;
}

#pragma warning( pop )

/*
 * no parameters.
 */
LPCTSTR CScanFilterFFT::GetParameters() const
{
   return NULL;
}

/**
 * check parameters, take care of a properly sized output buffer,
 * set its name and copy filter parameters and process;
 * see also CScanFilterNull::ApplyCore().
 */
BOOL CScanFilterFFT::Apply()
{
//   Q_LOG( _T("CScanFilterCrossCorrelate::Apply()") );

   if ( Q_INVALID( NULL == m_lpsbIn && "ensure input buffer") )
      return FALSE;

   if ( Q_INVALID( NULL == m_lpsbOut && "ensure output buffer" ) )
      return FALSE;

   if ( Q_INVALID( !m_lpsbIn->IsCompleteFrame() && "ensure complete frame" ) )
      return FALSE;

   /*
    * copy the filter settings (not the data):
    */
   Q_RETURN( m_lpsbOut->CreateOutputBufferFor( *m_lpsbIn ) );

   m_lpsbOut->m_csBufferName = m_lpsbIn->m_csBufferName + _T(" - ") + m_lpcsShortFilterName;
   m_lpsbOut->m_csDataName.Format( _T( "10Log %s(%s)" ), m_lpcsShortFilterName, m_lpsbIn->m_csDataName );

   /*
    * create a new suspended thread to do the computation;
    * the progressdialog resumes the thread and displays the progress:
    */
   CProgressDlg dlg( AfxGetMainWnd() ); // FFT has no dialog

   ThreadArgs threadArgs( &dlg, this );

   CWinThread *hThread = AfxBeginThread(
      ComputationThreadFunction,        // the controlling function for the worker thread
      &threadArgs,                      // the arguments and return value
      THREAD_PRIORITY_BELOW_NORMAL,     // the priority
      0,                                // 0: the stack size defaults to the same size
                                        //    stack as the creating thread.
      CREATE_SUSPENDED,                 // create it in suspended state
      NULL                              // security attributes structure
   );

   /*
    * show progress dialog, try to start thread:
    */
   if ( dlg.DoModal( hThread ) <= 0 )
   {
      return FALSE;
   };

   /*
    * return thread result:
    */
   return threadArgs.bReturn;
}

/**
 * create a thread that computes the result.
 */
UINT CScanFilterFFT::ComputationThreadFunction( LPVOID lpContext )
{
//   Q_LOG( _T("CScanFilterFFT::ComputationThreadFunction()") );

   CScanFilterFFT::ThreadArgs *lpTreadArgs = static_cast< CScanFilterFFT::ThreadArgs* >( lpContext );

   lpTreadArgs->bReturn = lpTreadArgs->pThis->Compute( lpTreadArgs->pDlg );

   lpTreadArgs->pDlg->Terminate();

   return lpTreadArgs->bReturn;
}

/**
 * compute the desired result.
 */
BOOL CScanFilterFFT::Compute( CProgressDlg *pDlg )
{
   ASSERT( m_lpsbIn );
   ASSERT( m_lpsbOut );

   /*
    * prepare some local varables:
    */
   BOOL bL2R = m_lpsbIn->HasLeftToRightScan();
   BOOL bR2L = m_lpsbIn->HasRightToLeftScan();

   /*
    * update progress dlg; don't touch dialog from other thread!
    */
   pDlg->SetText ( _T("Computing 2D fast fourier transformation" ) );
   pDlg->SetRange( 0, bR2L && bL2R ? 8 : 4 );

#if defined( ScanFilter_Use_FFTW )
   int  nrow = cil::GetSizeForFFT( m_lpsbIn->Rows()    );
   int  ncol = cil::GetSizeForFFT( m_lpsbIn->Columns() );
#else
   int  nrow = m_lpsbIn->GetRowsForFFT();
   int  ncol = m_lpsbIn->GetColumnsForFFT();
#endif // defined( ScanFilter_Use_FFTW )

   /*
    * (re)allocate buffers:
    */
   Q_RETURN( m_lpsbOut->ResizeDataBuffer( CSize( ncol, nrow ) ) );
   Q_RETURN( m_lpsbOut->CreateFourierBuffer()                   );

   /*
    * perform the left-right and right-left transformations:
    */
   m_progress = 0;

#ifdef USE_SEPARATE_AMPLITUDE_SCALE

   if ( bL2R )
   {
      ComputePart(
         pDlg,
         ToImage( m_lpsbIn , false ),
         ToImage( m_lpsbOut, false ),
         m_lpsbOut->FourierDataL2R()
      );
   }

   if ( bR2L )
   {
      ComputePart(
         pDlg,
         ToImage( m_lpsbIn , true ),
         ToImage( m_lpsbOut, true ),
         m_lpsbOut->FourierDataR2L()
      );

#else // ! USE_SEPARATE_AMPLITUDE_SCALE

   using cil::Size2D;

   /*
    * maximum amplitudes found for buffers:
    */
   RealType maxAmpL2R = 0;
   RealType maxAmpR2L = 0;

   /*
    * buffers for power spectrum:
    */
   Size2D zeroSize  (    0, 0    );
   Size2D fftSize   ( ncol, nrow );

   Size2D fftSizeL2R( bL2R ? fftSize : zeroSize );
   Size2D fftSizeR2L( bR2L ? fftSize : zeroSize );

   PowerImageType powerImageL2R( fftSizeL2R );
   PowerImageType powerImageR2L( fftSizeR2L );

   if ( bL2R )
   {
      maxAmpL2R = ComputePower(
         pDlg,
         ToImage( m_lpsbIn , false ),
         m_lpsbOut->FourierDataL2R(),
         powerImageL2R
      );
   }

   if ( bR2L )
   {
      maxAmpR2L = ComputePower(
         pDlg,
         ToImage( m_lpsbIn , true ),
         m_lpsbOut->FourierDataR2L(),
         powerImageR2L
      );
   }

   /*
    * select maximum amplitude to use:
    */
   enum ScaleSelector { e_Same, e_L2R, e_R2L };
   ScaleSelector s = e_Same;

   switch ( s )
   {
      case e_Same:
         maxAmpL2R = maxAmpR2L = lpm::absmax( maxAmpL2R, maxAmpR2L );
         break;

      case e_L2R:
         maxAmpR2L = maxAmpL2R ;
         break;

      case e_R2L:
         maxAmpL2R = maxAmpR2L;
         break;
   }

   /*
    * compute output:
    */
   if ( bL2R )
   {
      ComputeOutput(
         pDlg,
         maxAmpL2R,
         powerImageL2R,
         ToImage( m_lpsbOut, false )
      );
   }

   if ( bR2L )
   {
      ComputeOutput(
         pDlg,
         maxAmpR2L,
         powerImageR2L,
         ToImage( m_lpsbOut, true )
      );
   }
#endif // USE_SEPARATE_AMPLITUDE_SCALE

   /*
    * check if a stop is requested from the progress dialog:
    */
   if ( pDlg->StopRequested() )
   {
      return FALSE;
   }
   else
   {
      return TRUE;
   }
}

#ifdef USE_SEPARATE_AMPLITUDE_SCALE

/**
 * compute a frame.
 */
void CScanFilterFFT::ComputePart( CProgressDlg *pDlg, const Image& src, Image& dst, FourierElementType** pFFT )
{
   using   cil::Size2D;

   typedef std::complex< FourierElementType > Complex;

   typedef      Image             InputImageType;
   typedef      Image             OutputImageType;

   typedef cil::Image< RealType > PowerImageType;
   typedef cil::Image< Complex  > FFTImageType;

#if defined( ScanFilter_Use_FFTW )
   typedef cil::FFTWImageFilter             < InputImageType , FFTImageType    > FFTImageFilter;
#else
   typedef cil::FFTImageFilter              < InputImageType , FFTImageType    > FFTImageFilter;
#endif // defined( ScanFilter_Use_FFTW )
   typedef cil::FFTPowerSpectrumImageFilter < FFTImageType   , PowerImageType  > FFTPowerSpectrumImageFilter;
   typedef cil::FFTSwapQuadrantsImageFilter < OutputImageType                  > FFTSwapQuadrantsImageFilter;
   typedef cil::ShiftScaleImageFilter       < PowerImageType , OutputImageType > ScaleImageFilter;

   /*
    * the FFT Image view on the CScanBaseBuffer FFT databuffer:
    *
    * the FFT Image is a 1-dimensional array of complex<FourierElementType>
    * the CScanBaseBuffer FFT databuffer is a 2-dimensional array of FourierElementType
    *
    *            fftImage.Data () ----------->   cmplx
    * m_lpsbOut->FourierDataL2R() => [row0] -> [ re,im . . . ]  <- important: contiguous block
    *                                   .      [   .   . . . ]
    *                                   .      [   .   . . . ]
    *                                [rowN] -> [   .   . . . ]
    *
    */
   Size2D fftSize( Size2D( dst.Width(), dst.Height() ) );
   FFTImageType fftImage( reinterpret_cast<FFTImageType::PixelType*>( cil::To1DView( pFFT ) ), fftSize );

#if defined( USE_FFT_FFTI )     // for test, perform fft-fft inverse
   /*
    * filters:
    */
   typedef cil::FFTWInverseImageFilter< FFTImageType, OutputImageType > FFTInverseImageFilter;

   FFTImageFilter        fftFilter;
   FFTInverseImageFilter fftiFilter;

   /*
    * perform transformations:
    */
    fftFilter.Apply( src     , fftImage );
   fftiFilter.Apply( fftImage, dst      );
#else // !USE_FFT_FFTI
   /*
    * buffer for power spectrum:
    */
   PowerImageType powerImage( fftSize );  // this is temporary

   /*
    * filters:
    */
   FFTImageFilter                 fftFilter;
   FFTPowerSpectrumImageFilter    powerSpectrumFilter;
   FFTSwapQuadrantsImageFilter    swapQuadrantsFilter;
   ScaleImageFilter               scaleFilter;

   /*
    * perform transformations:
    */
             fftFilter.Apply( src       , fftImage   ); pDlg->SetPos( ++m_progress );
   powerSpectrumFilter.Apply( fftImage  , powerImage ); pDlg->SetPos( ++m_progress );

        scaleFilter.SetScale( maximumAmplitude / powerSpectrumFilter.GetAbsoluteMaximumAmplitude() );

           scaleFilter.Apply( powerImage, dst        ); pDlg->SetPos( ++m_progress );
   swapQuadrantsFilter.Apply( dst                    ); pDlg->SetPos( ++m_progress );
#endif
}

/*
 * stubs
 */
#pragma warning( disable : 4100 )       // unreferenced formal parameter
CScanFilterFFT::RealType CScanFilterFFT::ComputePower( CProgressDlg *pDlg, const Image& src, FourierElementType** pFFT, PowerImageType& pwr )
{
}

void CScanFilterFFT::ComputeOutput( CProgressDlg *pDlg, RealType maxAmpIn, const PowerImageType& pwr, Image& dst )
{
}
#pragma warning( default : 4100 )

#else // ! USE_SEPARATE_AMPLITUDE_SCALE

/*
 * stub
 */
#pragma warning( disable : 4100 )       // unreferenced formal parameter
void CScanFilterFFT::ComputePart( CProgressDlg *pDlg, const Image& src, Image& dst, FourierElementType** pFFT )
{
}
#pragma warning( default : 4100 )

#if defined( _CIL_MSC_6 )
#pragma optimize( "atpw", on )
#endif
/**
 * compute power for a frame.
 */
CScanFilterFFT::RealType CScanFilterFFT::ComputePower( CProgressDlg *pDlg, const Image& src, FourierElementType** pFFT, PowerImageType& pwr )
{
   using   cil::Size2D;

   typedef std::complex< FourierElementType > Complex;

   typedef      Image             InputImageType;
   typedef      Image             OutputImageType;

   typedef cil::Image< Complex  > FFTImageType;

#if defined( ScanFilter_Use_FFTW )
   typedef cil::FFTWImageFilter             < InputImageType , FFTImageType    > FFTImageFilter;
#else
   typedef cil::FFTImageFilter              < InputImageType , FFTImageType    > FFTImageFilter;
#endif // defined( ScanFilter_Use_FFTW )
   typedef cil::FFTPowerSpectrumImageFilter < FFTImageType   , PowerImageType  > FFTPowerSpectrumImageFilter;

   /*
    * the FFT Image view on the CScanBaseBuffer FFT databuffer:
    *
    * the FFT Image is a 1-dimensional array of complex<FourierElementType>
    * the CScanBaseBuffer FFT databuffer is a 2-dimensional array of FourierElementType
    *
    *            fftImage.Data () ----------->   cmplx
    * m_lpsbOut->FourierDataL2R() => [row0] -> [ re,im . . . ]  <- important: contiguous block
    *                                   .      [   .   . . . ]
    *                                   .      [   .   . . . ]
    *                                [rowN] -> [   .   . . . ]
    *
    */
   Size2D fftSize( Size2D( pwr.Width(), pwr.Height() ) );
   FFTImageType fftImage( reinterpret_cast<FFTImageType::PixelType*>( cil::To1DView( pFFT ) ), fftSize );

   /*
    * filters:
    */
   FFTImageFilter                 fftFilter;
   FFTPowerSpectrumImageFilter    powerSpectrumFilter;

   /*
    * perform transformations:
    */
             fftFilter.Apply( src       , fftImage ); pDlg->SetPos( ++m_progress );
   powerSpectrumFilter.Apply( fftImage  , pwr      ); pDlg->SetPos( ++m_progress );

   return powerSpectrumFilter.GetAbsoluteMaximumAmplitude();
}

/**
 * compute output for a frame.
 */
void CScanFilterFFT::ComputeOutput( CProgressDlg *pDlg, RealType maxAmpIn, const PowerImageType& pwr, Image& dst )
{
   using   cil::Size2D;

   typedef Image OutputImageType;

   typedef cil::ShiftScaleImageFilter       < PowerImageType , OutputImageType > ScaleImageFilter;
   typedef cil::FFTSwapQuadrantsImageFilter < OutputImageType                  > FFTSwapQuadrantsImageFilter;

   /*
    * filters:
    */
   ScaleImageFilter            scaleFilter;
   FFTSwapQuadrantsImageFilter swapQuadrantsFilter;

   /*
    * perform transformations:
    */
        scaleFilter.SetScale( maximumAmplitude / maxAmpIn );

           scaleFilter.Apply( pwr, dst ); pDlg->SetPos( ++m_progress );
   swapQuadrantsFilter.Apply( dst      ); pDlg->SetPos( ++m_progress );
}
#pragma optimize( "", on )      // restore

#endif // #ifdef USE_SEPARATE_AMPLITUDE_SCALE

/*
 * end of file
 */

---