imagingUtils.hpp

Go to the documentation of this file.

/** \file imagingUtils.hpp
 * \brief Utilities for modeling image formation
 * \ingroup imaging
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 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 imagingUtils_hpp
#define imagingUtils_hpp
 
#include <cmath>
 
#include "../mxlib.hpp"
#include "../improc/eigenImage.hpp"
#include "../math/constants.hpp"
 
#include "imagingArray.hpp"
 
namespace mx
{
 
namespace wfp
{
 
/// Calculate the angular plate scale (radians per pixel) of an image after propagation by FFT.
/**
 *
 * \returns the platescale of the wavefront after propagation by FFT, in radians per pixel.
 *
 * \tparam realT a real floating point type
 *
 * \ingroup imaging
 */
template <typename realT>
realT fftPlateScale( size_t pixels,        ///< [in] the linear dimension of the FFT (including 0 pad, etc.)
                     realT metersPerPixel, ///< [in] the scale of the input wavefront [m/pix]
                     realT lambda          ///< [in] the wavelength of the wavefront [m]
)
{
    return ( lambda / metersPerPixel ) * ( 1. / pixels );
}
 
/// Fill in an Eigen-like array with a circular pupil mask.
/** Sets any pixel which is at rad <= r < rad+(1.0/overscan) pixels to rho = 1,
 *
 *
 * \retval 0 on success
 * \retval -1 on error
 *
 * \ingroup imaging
 */
template <class arrayT>
int circularPupil( arrayT &m, ///< [in.out] is the allocated Array.  Dimensions are used to create the pupil.
                   typename arrayT::Scalar eps = 0, ///< [in] [optional] is the central obscuration.  0-1, default is 0.
                   typename arrayT::Scalar rad = 0, ///< [in] [optional] is the desired radius. If rad <= 0, then the
                                                    ///< maximum radius based on dimensions of m is used. Default is 0.
                   typename arrayT::Scalar overscan = 0 ///< [in] [optional] overscan in fractional pixels, to include
                                                        ///< partial pixels on the edge. Default is 0.
)
{
 
    if( eps < 0 )
    {
        internal::mxlib_error_report( error_t::invalidarg, "Central obscuration can not be < 0." );
        return -1;
    }
 
    if( eps > 1 )
    {
        internal::mxlib_error_report( error_t::invalidarg, "Central obscuration can not be > 1." );
        return -1;
    }
 
    size_t l0 = m.rows();
    size_t l1 = m.cols();
 
    typename arrayT::Scalar r;
    typename arrayT::Scalar xc = 0.5 * ( l0 - 1 );
    typename arrayT::Scalar yc = 0.5 * ( l1 - 1 );
 
    if( rad <= 0 )
        rad = 0.5 * std::min( l0 - 1, l1 - 1 );
 
    for( size_t i = 0; i < l0; i++ )
    {
        for( size_t j = 0; j < l1; j++ )
        {
            r = std::sqrt( std::pow( i - xc, 2 ) + std::pow( j - yc, 2 ) );
 
            if( r <= rad + ( overscan ) && r >= eps * rad )
                m( i, j ) = 1;
            else
                m( i, j ) = 0;
        }
    }
 
    return 0;
}
 
/// Draw a line in an image
/** \todo should handle width much more intelligently, this only works for ~45 degree lines.
 *
 * \tparam arrayT is an Eigen-like array with public typedef Scalar
 */
template <class arrayT>
void drawLine( arrayT &im,                       ///< [in.out] The input image, modified.
               typename arrayT::Scalar x0,       ///< [in] the x value, relative to image center, of the starting point
               typename arrayT::Scalar y0,       ///< [in] the y value, relative to image center, of the starting point
               typename arrayT::Scalar x1,       ///< [in] the x value, relative to image center, of the end point
               typename arrayT::Scalar y1,       ///< [in] the y value, relative to image center, of the end point
               typename arrayT::Scalar halfWidth ///< [in] the half-width of the line.
)
{
 
    int d1 = im.rows();
    int d2 = im.cols();
 
    typename arrayT::Scalar xc, yc;
 
    xc = 0.5 * ( im.rows() - 1 );
    yc = 0.5 * ( im.cols() - 1 );
 
    typename arrayT::Scalar m = ( y1 - y0 ) / ( x1 - x0 );
    int y;
 
    if( x1 > x0 )
    {
        for( int x = x0; x <= x1; ++x )
        {
            y = y0 + ( x - x0 ) * m;
 
            for( int i = 0; i <= halfWidth; ++i )
            {
                if( x + xc >= 0 && x + xc < d1 && y + yc + i >= 0 && y + yc + i < d2 )
                    im( x + xc, y + yc + i ) = 0;
                if( x + xc >= 0 && x + xc < d1 && y + yc - i >= 0 && y + yc - i < d2 )
                    im( x + xc, y + yc - i ) = 0;
            }
        }
    }
    else
    {
        for( int x = x1; x <= x0; ++x )
        {
            y = y0 + ( x - x0 ) * m;
 
            for( int i = 0; i <= halfWidth; ++i )
            {
                if( x + xc >= 0 && x + xc < d1 && y + yc + i >= 0 && y + yc + i < d2 )
                    im( x + xc, y + yc + i ) = 0;
                if( x + xc >= 0 && x + xc < d1 && y + yc - i >= 0 && y + yc - i < d2 )
                    im( x + xc, y + yc - i ) = 0;
            }
        }
    }
}
 
/// Create a complex pupil plane wavefront from a real amplitude mask.
/** The input real amplitude mask is placed in the center of a 0-padded complex array.
 *
 * \param [out] complexPupil the complex pupil plane wavefront
 * \param [in] realPupil a real amplitude mask.
 * \param [in] wavefrontSizePixels the desired size of the ouput wavefront, should be at least as big as realPupil
 *
 * \ingroup imaging
 */
template <typename arrayOutT, typename arrayInT>
void makeComplexPupil( arrayOutT &complexPupil, const arrayInT &realPupil, int wavefrontSizePixels )
{
 
    complexPupil.resize( wavefrontSizePixels, wavefrontSizePixels );
    complexPupil.setZero();
 
    // Lower-left corner of insertion region
    int bl = 0.5 * ( complexPupil.rows() - 1 ) - 0.5 * ( realPupil.rows() - 1. );
 
    for( int i = 0; i < realPupil.cols(); ++i )
    {
        for( int j = 0; j < realPupil.rows(); ++j )
        {
            complexPupil( bl + i, bl + j ) =
                typename arrayOutT::Scalar( realPupil( i, j ), 0 ); //*exp( typename arrayOutT::Scalar(0,1));
        }
    }
    // complexPupil.block(bl, bl, realPupil.rows(), realPupil.rows()) = realPupil*std::complex<realT>(1,0);
}
 
/// Create a complex wavefront from a real amplitude and a real phase.
/** The wavefront is placed in the center of a 0-padded complex array.
 *
 * \param [out] complexWavefront the complex pupil plane wavefront
 * \param [in] realAmplitude is the real-valued amplitude.
 * \param [in] realPhase is the real-valued phase in radians, same size as realAmplitude
 * \param [in] wavefrontSizePixels the desired size of the ouput wavefront, should be at least as big as the real arrays
 *
 * \ingroup imaging
 */
template <typename arrayOutT, typename arrayInT>
void makeComplexPupil( arrayOutT &complexWavefront,
                       const arrayInT &realAmplitude,
                       const arrayInT &realPhase,
                       int wavefrontSizePixels )
{
 
    complexWavefront.resize( wavefrontSizePixels, wavefrontSizePixels );
    complexWavefront.setConstant( typename arrayOutT::Scalar( 0, 0 ) );
 
    // Lower-left corner of insertion region
    int bl = 0.5 * ( complexWavefront.rows() - 1 ) - 0.5 * ( realAmplitude.rows() - 1. );
 
    for( int i = 0; i < realAmplitude.cols(); ++i )
    {
        for( int j = 0; j < realAmplitude.rows(); ++j )
        {
            complexWavefront( bl + j, bl + i ) =
                realAmplitude( j, i ) * exp( ( typename arrayOutT::Scalar( 0, 1 ) ) * realPhase( j, i ) );
        }
    }
}
 
/// Apply a tilt to a wavefront
/**
 * \param complexWavefront [in.out] the complex wavefront to tilt, will be modified on output
 * \param xTilt [input] the amount of tilt in the x direction, in pixels
 * \param yTilt [input] the amount of tilt in the y direction, in pixels
 *
 * \ingroup imaging
 */
template <typename wavefrontT>
void tiltWavefront( wavefrontT &complexWavefront,
                    typename wavefrontT::Scalar::value_type xTilt,
                    typename wavefrontT::Scalar::value_type yTilt )
{
    typedef typename wavefrontT::Scalar complexT;
    typedef typename wavefrontT::Scalar::value_type realT;
 
    realT pi = math::pi<realT>();
 
    int wfsSizeX = complexWavefront.cols();
    int wfsSizeY = complexWavefront.rows();
 
    realT xCen = 0.5 * ( wfsSizeX - 1 );
    realT yCen = 0.5 * ( wfsSizeY - 1 );
 
    realT argX = 2.0 * pi / ( wfsSizeX - 1.0 );
    realT argY = 2.0 * pi / ( wfsSizeY - 1.0 );
 
    for( int ii = 0; ii < wfsSizeX; ++ii )
    {
        for( int jj = 0; jj < wfsSizeY; ++jj )
        {
            complexWavefront( ii, jj ) =
                complexWavefront( ii, jj ) *
                exp( complexT( (realT)0., argX * xTilt * ( ii - xCen ) + argY * yTilt * ( jj - yCen ) ) );
        }
    }
}
 
/// Extract a block from one image and insert it into a second
template <typename imageT1,
          typename imageT2>
void extractBlock( imageT1 &dest, ///< [in/out] the image in which to place the extracted block.  Must be pre-allocated.
                   int imX0,      ///< [in] the x/row-coord of the lower left pixel of the destination region.
                   int imXsz,     ///< [in] the x size (number of rows) of the block
                   int imY0,      ///< [in] the y/col-coord of the lower left pixel of the desination region
                   int imYsz,     ///< [in] the y=size (number of cols) ov the block
                   imageT2 &src,  ///< [in] the source of the block.  Must be large enough.
                   int wfX0,      ///< [in] the x/row-coord of the lower left pixel of the source region
                   int wfY0 )     ///< [in] the y/col-coord of the lower left pixel of the source region
{
    int dest_cols = dest.cols();
 
    int src_cols = src.cols();
 
    typedef typename imageT1::Scalar dataT;
 
    dataT *dest_data;
    dataT *src_data;
 
    for( int j = 0; j < imYsz; ++j )
    {
        dest_data = &dest.data()[imX0 + ( imY0 + j ) * dest_cols];
        src_data = &src.data()[wfX0 + ( wfY0 + j ) * src_cols];
 
        memcpy( dest_data, src_data, sizeof( dataT ) * imXsz );
    }
}
 
/// Extract a pixels from one image and insert them into a second based on a mask
/** Only pixels with a non-zero value in mask are changed in dest to have the value in src.  Other pixels are not
 * modified.
 */
template <typename imageT1,
          typename imageT2,
          typename imageT3>
void extractMaskedPixels( imageT1 &dest, ///< [in/out] the image in which to place the extracted pixels.  Must be the
                                         ///< same size as src and mask.
                          const imageT2 &src, ///< [in] the source of the pixels.  Must be the same size as mask.
                          const imageT3 &mask ///< [in] the mask image, where any value other than 0 indicates a pixel
                                              ///< to extract.  Must be the same size as src.
)
{
    if( dest.rows() != src.rows() )
    {
        throw mx::exception( error_t::sizeerr, "dest and src do not have same size (rows)" );
    }
 
    if( dest.cols() != src.cols() )
    {
        throw mx::exception( error_t::sizeerr, "dest and src do not have same size (cols)" );
    }
 
    if( src.rows() != mask.rows() )
    {
        throw mx::exception( error_t::sizeerr, "src and mask do not have same size (rows)" );
    }
 
    if( src.cols() != mask.cols() )
    {
        throw mx::exception( error_t::sizeerr, "src and mask do not have same size (cols)" );
    }
 
    for( int cc = 0; cc < dest.cols(); ++cc )
    {
        for( int rr = 0; rr < dest.rows(); ++rr )
        {
            if( mask( rr, cc ) != 0 )
            {
                dest( rr, cc ) = src( rr, cc );
            }
        }
    }
}
 
template <typename realImageT, typename complexImageT>
void extractIntensityImage(
    realImageT &im, int imX0, int imXsz, int imY0, int imYsz, complexImageT &wf, int wfX0, int wfY0 )
{
    int im_rows = im.cols();
 
    int wf_rows = wf.cols();
 
    typename realImageT::Scalar *im_data;
    typename complexImageT::Scalar *wf_data;
 
    for( int j = 0; j < imXsz; ++j )
    {
        im_data = &im.data()[imX0 + ( imY0 + j ) * im_rows];
        wf_data = &wf.data()[wfX0 + ( wfY0 + j ) * wf_rows];
        for( int i = 0; i < imYsz; ++i )
        {
            im_data[i] = norm( wf_data[i] );
        }
    }
}
 
/// Extract the intensity image from a complex wavefront and accumulate the result
/** Can work on a subset of either or both of the image and the wavefront.  This way only
 * a small part of the image can be extracted for an oversampled wavefront.
 */
template <typename realImageT, typename complexImageT, int cudaGPU = 0>
void extractIntensityImageAccum( realImageT &im,    /**<[out] The real image to populate with intensity */
                                 int imX0,          /**<[in] The x-coord of the starting image pixel*/
                                 int imXsz,         /**<[in] The number of x image pixels*/
                                 int imY0,          /**<[in] The y-coord of the starting image pixel */
                                 int imYsz,         /**<[in] The number of y image pixels */
                                 complexImageT &wf, /**<[in] The complex wavefront*/
                                 int wfX0,          /**<[in] The x-coord of the starting wavefront pixel*/
                                 int wfY0           /**<[in] The y-coord of the starting wavefront pixel*/
);
 
 
 
// CPU versions
template <>
void extractIntensityImageAccum<improc::eigenImage<float>, improc::eigenImage<std::complex<float>>, 0>(
    improc::eigenImage<float> &im, int imX0, int imXsz, int imY0, int imYsz, improc::eigenImage<std::complex<float>> &wf, int wfX0, int wfY0 );
 
template <>
void extractIntensityImageAccum<improc::eigenImage<double>, improc::eigenImage<std::complex<double>>, 0>(
    improc::eigenImage<double> &im, int imX0, int imXsz, int imY0, int imYsz, improc::eigenImage<std::complex<double>> &wf, int wfX0, int wfY0 );
 
 
#ifdef MXLIB_CUDA
 
 
// GPU versions
template <>
void extractIntensityImageAccum<cuda::cudaPtr<float>, cuda::cudaPtr<std::complex<float>>, 1>(
    cuda::cudaPtr<float> &im, int imX0, int imXsz, int imY0, int imYsz, cuda::cudaPtr<std::complex<float>> &wf, int wfX0, int wfY0 );
 
template <>
void extractIntensityImageAccum<cuda::cudaPtr<double>, cuda::cudaPtr<std::complex<double>>, 1>(
    cuda::cudaPtr<double> &im, int imX0, int imXsz, int imY0, int imYsz, cuda::cudaPtr<std::complex<double>> &wf, int wfX0, int wfY0 );
 
#endif
 
} // namespace wfp
} // namespace mx
 
#endif //__imagingUtils_hpp__