fraunhoferPropagator.hpp

Go to the documentation of this file.

/** \file fraunhoferPropagator.hpp
 * \brief Declares and defines a class for Fraunhofer propagation of optical wavefronts
 * \ingroup imaging
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2015-2020 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib 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 3 of the License, or
// (at your option) any later version.
//
// mxlib 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 mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef wfp_fraunhoferPropagator_hpp
#define wfp_fraunhoferPropagator_hpp
 
#include "../math/constants.hpp"
#include "imagingArray.hpp"
#include "imagingUtils.hpp"
 
#include "../math/fft/fft.hpp"
 
namespace mx
{
 
namespace wfp
{
 
/// Class to perform Fraunhofer propagation between pupil and focal planes
/** This class uses the FFT to propagate between planes, and normalizes so that flux
 * is conserved.  For propagation from pupil to focal plane, the pupil wavefront is tilted so that
 * the focal-plane image is centered at the geometric center of the array.  After propagation from
 * focal plane to pupil plane, the pupil plane wavefront is un-tilted to restore the
 * pupil to its original position.
 *
 * \tparam _wavefrontT is an Eigen::Array-like type, with std::complex values.
 *
 * \todo check if we should reverse the FFT orders to be proper
 *
 * \ingroup imaging
 */
template <typename _wavefrontT>
class fraunhoferPropagator
{
 
  public:
    /// The wavefront data type
    typedef _wavefrontT wavefrontT;
 
    /// The complex data type
    typedef typename wavefrontT::Scalar complexT;
 
    /// The real data type
    typedef typename wavefrontT::Scalar::value_type realT;
 
  protected:
    /// The size of the wavefront in pixels
    int m_wavefrontSizePixels{ 0 };
 
    realT m_xcen{ 0 }; ///< x-coordinate of focal plane center, in pixels
    realT m_ycen{ 0 }; ///< x-coordinate of focal plane center, in pixels
 
    /// Determines how the image is centered.
    /** If 0 (default) it is at 0.5*(wfSz-1), if true it is shifted by 0.5*m_wholePixel in each axis.
     */
    realT m_wholePixel{ 0 };
 
    /// Phase screen for tilting the pupil plane so that the focal plane image is centered.
    wavefrontT m_centerFocal;
 
    /// Phase screen for un-tilting the pupil plane after propagating from a centered focal plane.
    wavefrontT m_centerPupil;
 
    /// FFT object for forward FFTs
    math::fft::fftT<complexT, complexT, 2, 0> m_fft_fwd;
 
    /// FFT object for backward FFTs
    math::fft::fftT<complexT, complexT, 2, 0> m_fft_back;
 
    /// Initialize members
    void initialize();
 
  public:
    /// Constructor
    fraunhoferPropagator();
 
    /// Destructor
    ~fraunhoferPropagator();
 
    /// Get the value of the wholePixel parameter
    /** The wholePixel parameter determines how the image is centered.  If 0 (default) it is at 0.5*(wfSz-1), if true it
     * is shifted by 0.5*m_wholePixel in each axis.
     *
     * \returns the value of m_wholePixel
     */
    int wholePixel();
 
    /// Set the value of the wholePixel parameter
    /** The wholePixel parameter determines how the image is centered.  If 0 (default) it is at 0.5*(wfSz-1), if true it
     * is shifted by 0.5*m_wholePixel in each axis.
     */
    void wholePixel( realT wp /**< [in] the new wholePixel value */ );
 
    /// Apply the shift to a pupil wavefront which will center the resultant focal plane image, and apply the
    /// normalization.
    /** You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
     * setWavefrontSizePixels.
     */
    void shiftPupil( wavefrontT &complexPupil /**< [in.out] the complex pupil plane wavefront to shift*/ );
 
    /// Apply the shift to a pupil wavefront which will restore it to a centered pupil image, with correct flux.
    /** You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
     * setWavefrontSizePixels.
     */
    void unshiftPupil( wavefrontT &complexPupil /**< [in.out] the complex pupil plane wavefront to shift*/ );
 
    /// Propagate the wavefront from the pupil plane to the focal plane
    /** The pupil plane wavefront (complexPupil) is multiplied by a tilt to place the
     * image in the geometric center of the focal plane.  This can be prevented by
     * setting doCenter to false.
     *
     */
    void propagatePupilToFocal(
        wavefrontT
            &complexFocal, ///< [out] the focal plane wavefront.  Must be pre-allocated to same size as complexPupil.
        wavefrontT &complexPupil, ///< [in] the pupil plane wavefront. Modified due to application of centering tilt.
        bool doCenter = true      ///< [in] [optional] set to false to not apply the centering shift
    );
 
    /// Propagate the wavefront from Focal plane to Pupil plane
    /** After the fourier transform, the output pupil plane wavefront is de-tilted, restoring it
     * to the state prior to calling \ref propagatePupilToFocal.  This can be prevented by
     * setting doCenter to false.
     *
     */
    void propagateFocalToPupil( wavefrontT &complexPupil, ///< [out] the pupil plane wavefront. Must be pre-allocated to
                                                          ///< same size as complexFocal.
                                wavefrontT &complexFocal, ///< [in] the focal plane wavefront.
                                bool doCenter = true ///< [in] [optional] set to false to not apply the centering shift
    );
 
    /// Set the size of the wavefront, in pixels
    /** Checks if the size changes, does nothing if no change.  Otherwise,  calls
     * \ref makeShiftPhase to pre-calculate the tilt arrays and plans the FFTs.
     *
     */
    void setWavefrontSizePixels( int wfsPix /**< [in] the desired new size of the wavefront */ );
 
  protected:
    /// Calculate the complex tilt arrays for centering and normalizing the wavefronts
    /**
     */
    void makeShiftPhase();
 
}; // class fraunhoferPropagator
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::initialize()
{
    m_wavefrontSizePixels = 0;
 
    m_xcen = 0;
    m_ycen = 0;
}
 
template <typename wavefrontT>
fraunhoferPropagator<wavefrontT>::fraunhoferPropagator()
{
}
 
template <typename wavefrontT>
fraunhoferPropagator<wavefrontT>::~fraunhoferPropagator()
{
}
 
template <typename wavefrontT>
int fraunhoferPropagator<wavefrontT>::wholePixel()
{
    return m_wholePixel;
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::wholePixel( realT wp )
{
    m_wholePixel = wp;
 
    // Re-make the shift phase if size already set.
    if( m_wavefrontSizePixels > 0 )
        makeShiftPhase();
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::shiftPupil( wavefrontT &complexPupil )
{
    complexPupil *= m_centerFocal;
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::unshiftPupil( wavefrontT &complexPupil )
{
    complexPupil *= m_centerPupil;
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::propagatePupilToFocal( wavefrontT &complexFocal,
                                                              wavefrontT &complexPupil,
                                                              bool doCenter /*default = true*/
)
{
    // First setup the tilt screens (does nothing if there's no change in size)
    setWavefrontSizePixels( complexPupil.rows() );
 
    // Apply the centering shift -- this adjusts by 0.5 pixels and normalizes
    if( doCenter )
        shiftPupil( complexPupil );
 
    // fft_fwd.fft(complexPupil.data(), complexFocal.data() );
    m_fft_fwd( complexFocal.data(), complexPupil.data() );
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::propagateFocalToPupil( wavefrontT &complexPupil,
                                                              wavefrontT &complexFocal,
                                                              bool doCenter /*default = true*/
)
{
    // First setup the tilt screens (does nothing if there's no change in size)
    setWavefrontSizePixels( complexPupil.rows() );
 
    // fft_back.fft( complexFocal.data(), complexPupil.data());
    m_fft_back( complexPupil.data(), complexFocal.data() );
 
    // Unshift the wavefront and normalize
    if( doCenter )
        unshiftPupil( complexPupil );
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::setWavefrontSizePixels( int wfsPix )
{
    // If no change in size, do nothing
    if( wfsPix == m_centerFocal.rows() )
        return;
 
    m_wavefrontSizePixels = wfsPix;
 
    m_xcen = 0.5 * ( wfsPix - 1.0 );
    m_ycen = 0.5 * ( wfsPix - 1.0 );
 
    makeShiftPhase();
 
    m_fft_fwd.plan( wfsPix, wfsPix );
 
    m_fft_back.plan( wfsPix, wfsPix, MXFFT_BACKWARD );
}
 
template <typename wavefrontT>
void fraunhoferPropagator<wavefrontT>::makeShiftPhase()
{
    constexpr realT pi = math::pi<realT>();
 
    // The normalization is included in the tilt.
    realT norm = 1. / ( m_wavefrontSizePixels * sqrt( 2 ) );
    complexT cnorm = complexT( norm, norm );
 
    // Resize the center phases
    m_centerFocal.resize( m_wavefrontSizePixels, m_wavefrontSizePixels );
    m_centerPupil.resize( m_wavefrontSizePixels, m_wavefrontSizePixels );
 
    // Shift by 0.5 pixels
    realT arg = -2.0 * pi * 0.5 * ( m_wavefrontSizePixels - m_wholePixel ) / ( m_wavefrontSizePixels - 1 );
 
    for( int ii = 0; ii < m_wavefrontSizePixels; ++ii )
    {
        for( int jj = 0; jj < m_wavefrontSizePixels; ++jj )
        {
            m_centerFocal( ii, jj ) = cnorm * exp( complexT( 0., arg * ( ( ii - m_xcen ) + ( jj - m_ycen ) ) ) );
            m_centerPupil( ii, jj ) =
                cnorm * exp( complexT( 0., 0.5 * pi - arg * ( ( ii - m_xcen ) + ( jj - m_ycen ) ) ) );
        }
    }
}
 
} // namespace wfp
} // namespace mx
 
#endif // wfp_fraunhoferPropagator_hpp