Filter Library

Camera

Interface Physics

ScanFilterFFTInverse.cpp

/*
 * ScanFilterFFTInverse.cpp - inverse 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: M. Seynen (original), Martin J. Moene
 *
 * $Id: ScanFilterFFTInverse.cpp 511 2007-04-03 10:19:11Z moene $
 */

#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

#include <cfl/ScanFilterFFTInverse.h>   // this filter's header
#include <cfl/ScanFilterClip.h>         // clipping (cropping) filter
#include <cfl/ToolCilImage.h>           // for Image type and ToImage conversion shim

#include <cil/Math.h>                   // for real(), //min(), max(), sqrt(), sin(), Pi
#include <cil/Traits.h>                 // for NumericTraits<>::Min()
#include <cil/FFTInverseImageFilter.h>          // FFT inverse image filter
#include <cil/RescaleIntensityImageFilter.h>    // to rescale inverse FFT image

#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

/*
 * configuration defines:
 */

// currently none
// #define USE_FULL_SCALE               // rescale output
// #define USE_FULL_IFFT                // do not crop image

/*
 * end of configuration defines.
 */

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

#define CFL_SCHEMAVERSION_FILTERIFFT 0  // permanent storage schema version
                                        // this filter's name
LPCTSTR CScanFilterFFTInverse::m_lpcsFilterName = _T("Inverse-Fast-Fourier-Transform");
                                        // this filter's short name
LPCTSTR CScanFilterFFTInverse::m_lpcsShortFilterName = _T( "FFT-1" );
                                        // serialization
IMPLEMENT_SERIAL( CScanFilterFFTInverse, CScanFilter, CFL_SCHEMAVERSION_FILTERIFFT)

/**
 * constructor.
 */
CScanFilterFFTInverse::CScanFilterFFTInverse()
   : m_lpsbTmp( new CScanBaseBuffer )
{
//   Q_LOG( _T("CScanFilterBgsDiff::CScanFilterBgsDiff()") );

   m_csFilterName = m_lpcsFilterName;
}

/**
 * destructor.
 */
CScanFilterFFTInverse::~CScanFilterFFTInverse()
{
//   Q_LOG( _T("CScanFilterFFTInverse::~CScanFilterFFTInverse()") );

   delete m_lpsbTmp;
   ; // do nothing
}

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

/**
 * store or retrieve the object's settings;
 * see also CScanFilterNull::Serialize().
 */
void CScanFilterFFTInverse::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 CScanFilterFFTInverse::SetParameters(LPCTSTR lpParameters)
{
   return TRUE;
}

#pragma warning( pop )

/*
 * no parameters.
 */
LPCTSTR CScanFilterFFTInverse::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 CScanFilterFFTInverse::Apply()
{
//   Q_LOG( _T("CScanFilterCrossCorrelate::Apply()") );

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

   if ( Q_INVALID( NULL == m_lpsbTmp ) && "ensure valid temporary 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_lpsbTmp->CreateOutputBufferFor( *m_lpsbIn ) );

   m_lpsbTmp->m_csBufferName = m_lpsbIn->m_csBufferName + _T(" - ") + m_lpcsShortFilterName;
   m_lpsbTmp->m_csDataName.Format(_T("%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() ); // IFFT 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 CScanFilterFFTInverse::ComputationThreadFunction( LPVOID lpContext )
{
//   Q_LOG( _T("CScanFilterFFTInverse::ComputationThreadFunction()") );

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

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

   lpTreadArgs->pDlg->Terminate();

   return lpTreadArgs->bReturn;
}

/**
 * compute the desired result.
 */
BOOL CScanFilterFFTInverse::Compute( CProgressDlg *pDlg )
{
//   Q_LOG( _T("CScanFilterFFTInverse::Compute()") );

   ASSERT( m_lpsbIn  );
   ASSERT( m_lpsbTmp );
   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 inverse fast fourier transformation" ) );
   pDlg->SetRange( 0, bR2L && bL2R ? 2: 1 );

   /*
    * check for FFT input to filter:
    */
   if ( bL2R && 0 == m_lpsbIn->m_pdFourierL2R ||
        bR2L && 0 == m_lpsbIn->m_pdFourierR2L   )
   {
      MessageBox(
         IsInteractive() ? m_pDlg->m_hWnd : AfxGetMainWnd()->m_hWnd,
         _T( "This scanbuffer does not contain a fourier spectrum." ),
         _T( "Cannot transform" ),
         MB_ICONEXCLAMATION | MB_OK
      );

      return FALSE;
   }

   /*
    * allocate fourier buffer:
    *
    * Note if we come this far, we know that the temporary buffer is created
    * from an input buffer that contains an FFT, so the temporary buffer's size
    * is already appropriate for handling FFT.
    */
   Q_RETURN( m_lpsbTmp->CopyFourierBufferFrom( *m_lpsbIn ) );

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

   if ( bL2R )
   {
      ComputePart( m_lpsbTmp->FourierDataL2R(), ToImage( m_lpsbTmp, false ) );

      pDlg->SetPos( ++progress );
   }

   if ( bR2L )
   {
      ComputePart( m_lpsbTmp->FourierDataR2L(), ToImage( m_lpsbTmp, true ) );

      pDlg->SetPos( ++progress );
   }

   /*
    * release the allocated fourier memory, because it does not contain an FFT
    * spectrum anymore:
    */
   m_lpsbTmp->DeleteFourierBuffer();

   /*
    * crop the image to the original size, unless USE_FULL_IFFT is defined;
    */
   Q_RETURN( CropImage() );

   /*
    * set the output buffer name (prevent Clip from showing):
    */
   m_lpsbOut->m_csBufferName.Format( _T("%s - %s"), m_lpsbIn->m_csBufferName, m_lpcsShortFilterName );

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

/**
 * compute a frame.
 */
void CScanFilterFFTInverse::ComputePart( FourierElementType** pFFT, Image& dst )
{
   using   cil::Size2D;

   typedef std::complex< FourierElementType > Complex;

   typedef cil::Image< Complex              > FFTImageType;
   typedef cil::Image< FourierElementType   > IFFTImageType;
   typedef      Image                         OutputImageType;

   typedef cil::FFTInverseImageFilter       < FFTImageType , IFFTImageType   > FFTInverseImageFilter;
   typedef cil::RescaleIntensityImageFilter < IFFTImageType, OutputImageType > RescaleIntensityImageFilter;

   /*
    * 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 );

   IFFTImageType ifftImage;

   /*
    * do the inverse FFT transformation:
    */
   FFTInverseImageFilter fftiFilter;
   fftiFilter.Apply( fftImage, ifftImage );

   /*
    * scale to result to the (integral) destination image:
    */
   RescaleIntensityImageFilter scaleFilter;
   scaleFilter.SetOutputMinimum( lpm::NumericTraits< OutputImageType::ValueType >::Min() );
   scaleFilter.SetOutputMaximum( lpm::NumericTraits< OutputImageType::ValueType >::Max() );
   scaleFilter.Apply( ifftImage, dst );
}

/**
 * crop image to original size (remove border with zero data).
 */
BOOL CScanFilterFFTInverse::CropImage()
{
   typedef CScanFilterClip CropFilter;

   CropFilter cropFilter;

   cropFilter.SetDragInput   ( false, GetUsedRectangle() );
   cropFilter.SetInputBuffer ( 0, m_lpsbTmp );
   cropFilter.SetOutputBuffer( 0, m_lpsbOut );

   return cropFilter.Apply();
}

/**
 * true if any edge contains non-zero data.
 */
bool CScanFilterFFTInverse::HasNonzeroEdge() const
{
   Pointer  data = m_lpsbTmp->Data();
   SizeType ncol = m_lpsbTmp->Columns();
   SizeType size = m_lpsbTmp->GetFrameSize();

   return 0 != data[  0              ]  // top-left
       || 0 != data[        ncol - 1 ]  // top-right
       || 0 != data[ size - ncol - 1 ]  // bottom-left
       || 0 != data[ size        - 1 ]; // bottom-right
}

/**
 * the rectangle that contains non-zero data;
 * from the outside to the inside we search for non-zero data;
 * if USE_FULL_IFFT is defined, rectange is full IFFT result.
 */
CRect CScanFilterFFTInverse::GetUsedRectangle() const
{
   typedef Pointer Iterator;

#ifdef _CIL_MSC_6
   typedef std::reverse_iterator<Pointer, ValueType> ReverseIterator;
#else
   typedef std::reverse_iterator<Pointer> ReverseIterator;
#endif

   SizeType nrow   = m_lpsbTmp->Rows();
   SizeType ncol   = m_lpsbTmp->Columns();
   SizeType size   = m_lpsbTmp->GetFrameSize();

   /*
    * set default cropping coordinates:
    */
   SizeType left   = 0;
   SizeType right  = ncol - 0;
   SizeType top    = 0;
   SizeType bottom = nrow - 0;

#ifdef USE_FULL_IFFT
   return CRect( left, top, right, bottom );
#endif

   /*
    * if any edge contains non-zero data, return full rectangle:
    */
   if ( HasNonzeroEdge() )
   {
      return CRect( left, top, right, bottom );
   }

   /*
    * determine used rectangle:
    */
   Iterator itopleft  = m_lpsbTmp->Data();
   Iterator imidleft  = itopleft + size / 2;
   Iterator icenter   = imidleft + ncol / 2;
   Iterator ibotright = itopleft + size;

   /*
    * downward search each line for non-zero data:
    */
   Iterator pos = std::find_if(
      itopleft,                                         // start
      icenter ,                                         // end
      std::bind2nd( std::not_equal_to<ValueType>(), 0 ) // criterion
   );

   if ( pos != icenter )
   {
      top = 0 + ( pos - itopleft ) / ncol;
   }

   /*
    * upward search each line for non-zero data:
    */
   ReverseIterator rpos = std::find_if(
      ReverseIterator( ibotright ),                     // start
      ReverseIterator( icenter   ),                     // end
      std::bind2nd( std::not_equal_to<ValueType>(), 0 ) // criterion
   );

   if ( rpos != ReverseIterator( icenter ) )
   {
      bottom = nrow - ( rpos - ReverseIterator( ibotright ) ) / ncol;
   }

   /*
    * rightward search each column for non-zero data:
    */

   for ( Iterator icol = itopleft; icol < itopleft + ncol / 2; ++icol )
   {
      for ( Iterator it = icol; it < icenter; it += ncol )
      {
         if ( 0 != *it )
         {
            left = 0 + ( it - itopleft ) % ncol;
            goto foundLeft;
         }
      }
   }
foundLeft: ;

   /*
    * leftward search each column for non-zero data:
    */

   for ( ReverseIterator ricol = ReverseIterator( ibotright );
                         ricol < ReverseIterator( itopleft + ncol / 2 ); ++ricol )
   {
      for ( ReverseIterator it = ricol; it < ReverseIterator( icenter ); it += ncol )
      {
         if ( 0 != *it )
         {
            right = ncol - ( it - ReverseIterator( ibotright ) ) % ncol;
            goto foundRight;
         }
      }
   }
foundRight: ;

   return CRect( left, top, right, bottom );
}

/*
 * end of file
 */

---