---

pass_band_filter.h

Go to the documentation of this file.

#ifndef __PASS_BAND_FILTER__
#define __PASS_BAND_FILTER__

#include <complex>

#include "function_2D.h"

#include "fftw3.h"

#include "fft/support.h"


namespace FFT{


  template<typename Function>
  inline void difference_of_gaussians(const Function& f,
                                      const float     small_wave_length,
                                      const float     large_wave_length,
                                      Function* const g);



  template<typename Function>
  inline void  difference_of_gaussians(const Function& f,
                                       const float     small_wave_length,
                                       const float     large_wave_length,
                                       Function* const g,
                                       Support&        support);

  template<typename Function>
  inline void difference_of_gaussians(const Function& f,
                                      const float     center_wave_length,
                                      Function* const g);



  template<typename Function>
  inline void  difference_of_gaussians(const Function& f,
                                       const float     center_wave_length,
                                       Function* const g,
                                       Support&        support);


  template<typename Function>
  inline void difference_of_gaussians(const Support::buffer_type& F,
                                      const float     center_wave_length,
                                      Function* const g,
                                      Support&        support);

  
  
  template<typename Function>
  void difference_of_gaussians(const Support::buffer_type& F,
                               const float     small_wave_length,
                               const float     large_wave_length,
                               Function* const g,
                               Support&        support);




  namespace Fourier_domain{ 
  
    inline void
    difference_of_gaussians(const Support::buffer_type& F,
                            const float                 center_wave_length,
                            Support::buffer_type* const G,
                            Support&                    support);
  
    inline void
    difference_of_gaussians(const Support::buffer_type& F,
                            const float                 small_wave_length,
                            const float                 large_wave_length,
                            Support::buffer_type* const G,
                            Support&                    support);
  }


  

/*
  ###########################################################
  ###########################################################
  ###########################################################
  ##############                               ##############
  ##############  I M P L E M E N T A T I O N  ##############
  ##############                               ##############
  ###########################################################
  ###########################################################
  ###########################################################
*/

  


  template<typename Function>
  void difference_of_gaussians(const Function& f,
                               const float     center_wave_length,
                               Function* const g){

    difference_of_gaussians(f,
                            center_wave_length*M_SQRT1_2*M_SQRT1_2,
                            center_wave_length*M_SQRT2*M_SQRT1_2,
                            g);

  }



  template<typename Function>
  void  difference_of_gaussians(const Function& f,
                                const float     center_wave_length,
                                Function* const g,
                                Support&        support){

    difference_of_gaussians(f,
                            center_wave_length*M_SQRT1_2*M_SQRT1_2,
                            center_wave_length*M_SQRT2*M_SQRT1_2,
                            g,
                            support);
  }


  
  template<typename Function>
  void difference_of_gaussians(const Function& f,
                               const float     small_wave_length,
                               const float     large_wave_length,
                               Function* const g,
                               Support&        support){

    // We need a temporary buffer.
    Support::buffer_type buffer = support.create_buffer();

    // Store f data in the fftw structures.
    support.load_space_data(f);
    
    // Compute the FFT.
    support.space_to_frequency();
    
    // Copy the result in a safe place.
    support.save_frequency_data_into_buffer(buffer);

    difference_of_gaussians(buffer,
                            small_wave_length,
                            large_wave_length,
                            g,
                            support);
  }


  
  template<typename Function>
  void difference_of_gaussians(const Function& f,
                               const float     small_wave_length,
                               const float     large_wave_length,
                               Function* const g){

    Support s(f.x_size(),f.y_size());

    difference_of_gaussians(f,small_wave_length,large_wave_length,g,s);
  }



  
  template<typename Function>
  void difference_of_gaussians(const Support::buffer_type& F,
                               const float                 center_wave_length,
                               Function* const             g,
                               Support&                    support){

    difference_of_gaussians(F,
                            center_wave_length*M_SQRT1_2,
                            center_wave_length*M_SQRT2,
                            g,
                            support);
  }

  
  

  
  
  template<typename Function>
  void difference_of_gaussians(const Support::buffer_type& F,
                               const float                 small_wave_length,
                               const float                 large_wave_length,
                               Function* const             g,
                               Support&                    support){

    using namespace Function_2D;


    // Type definition
    typedef float                              real_type;
    typedef std::complex<float>                complex_type;
    typedef FFT_normalized_gaussian            gaussian_type;
    typedef Minus<gaussian_type,gaussian_type> dog_type;

    // Creation of the data structures
    Array_2D<complex_type> dog_kernel(support.x_frequency_size(),
                                      support.y_frequency_size(),
                                      complex_type());

    const real_type width  = static_cast<real_type>(support.x_space_size());
    const real_type height = static_cast<real_type>(support.y_space_size());
    
    const real_type small_sigma_x = 2 * width / small_wave_length;
    const real_type small_sigma_y = 2 * height / small_wave_length;

    const real_type large_sigma_x = 2 * width / large_wave_length;
    const real_type large_sigma_y = 2 * height / large_wave_length;

    gaussian_type small_gaussian(0,0,small_sigma_x,small_sigma_y,0);
    gaussian_type large_gaussian(0,0,large_sigma_x,large_sigma_y,0);

    const dog_type dog(small_gaussian,large_gaussian);

    // Creation of DoG kernel
    FFT::fill_real_part(dog,&dog_kernel);

    // Load the data in the FFT engine
    support.load_frequency_data(dog_kernel);

    // Multiplication "in place" in the frequency space.    
    support.multiply_frequency_data_by_buffer(F,true);

    // Compute inverse FFT transform.
    support.frequency_to_space();

    // Return data into the original structures    
    support.save_space_data(g);
  }



  namespace Fourier_domain{ 

    void difference_of_gaussians(const Support::buffer_type& F,
                                 const float                 center_wave_length,
                                 Support::buffer_type* const G,
                                 Support&                    support){

      FFT::Fourier_domain::difference_of_gaussians(F,
                                                   center_wave_length*M_SQRT1_2*M_SQRT1_2,
                                                   center_wave_length*M_SQRT2*M_SQRT1_2,
                                                   G,
                                                   support);
    }




    void difference_of_gaussians(const Support::buffer_type& F,
                                 const float                 small_wave_length,
                                 const float                 large_wave_length,
                                 Support::buffer_type* const G,
                                 Support&                    support){

      using namespace Function_2D;


      // Type definition
      typedef float                              real_type;
      typedef std::complex<float>                complex_type;
      typedef FFT_normalized_gaussian            gaussian_type;
      typedef Minus<gaussian_type,gaussian_type> dog_type;

      // Creation of the data structures
      Array_2D<complex_type> dog_kernel(support.x_frequency_size(),
                                        support.y_frequency_size(),
                                        complex_type());

      const real_type width  = static_cast<real_type>(support.x_space_size());
      const real_type height = static_cast<real_type>(support.y_space_size());
    
      const real_type small_sigma_x = 2 * width / small_wave_length;
      const real_type small_sigma_y = 2 * height / small_wave_length;

      const real_type large_sigma_x = 2 * width / large_wave_length;
      const real_type large_sigma_y = 2 * height / large_wave_length;

      gaussian_type small_gaussian(0,0,small_sigma_x,small_sigma_y,0);
      gaussian_type large_gaussian(0,0,large_sigma_x,large_sigma_y,0);

      const dog_type dog(small_gaussian,large_gaussian);

      // Creation of DoG kernel
      FFT::fill_real_part(dog,&dog_kernel);

      // Load the data in the FFT engine
      support.load_frequency_data(dog_kernel);

      // Multiplication "in place" in the frequency space.    
      support.multiply_frequency_data_by_buffer(F,false);

      // Return data into a buffer.
      support.save_frequency_data_into_buffer(*G);
    }

  }
  
}

#endif

---