ADIobservation.hpp

Go to the documentation of this file.

/** \file ADIobservation.hpp
 * \author Jared R. Males
 * \brief Defines the ADI high contrast imaging data type.
 * \ingroup hc_imaging_files
 * \ingroup image_processing_files
 *
 */
 
#ifndef __ADIobservation_hpp__
#define __ADIobservation_hpp__
#include "../ioutils/fileUtils.hpp"
 
#include "HCIobservation.hpp"
#include "../ioutils/fits/fitsHeader.hpp"
 
#include "imagePads.hpp"
 
namespace mx
{
namespace improc
{
namespace HCI
{
/// Fake injection PSF file specification methods
/** \ingroup hc_imaging_enums
 */
enum fakeMethods
{
    single, ///< A single PSF is used
    list    ///< A list of PSF files, one per input image, is used.
};
 
/// Get the string name of a fake injection method
/**
 * \returns the string name corresponding to the fake injection method
 */
std::string fakeMethodsStr( int method /**< [in] the fake injection method */ );
 
/// Get the fake injection method from its string name
/**
 * \returns the corresponding member of the fakeMethods enum
 */
int fakeMethodFmStr( const std::string &method /**< [in] the fake injection method name*/ );
 
} // namespace HCI
 
/// Process an angular differential imaging (ADI) observation
/** Angular differential imaging (ADI) uses sky rotation to differentiate real objects from
 * speckles.
 *
 * \tparam realT is the floating point type in which to do calculations
 *
 * \tparam _derotFunctObj
 * \parblock
 * is the derotation object with the following minimum interface:
 * \code
 * template<typename _realT>
 * struct derotF
 * {
 *    typedef  _realT realT;
 *
 *    //Vector of keywords to extract from the fits headers
 *    std::vector<std::string> keywords;
 *
 *    //Vector(s) to hold the keyword values
 *    std::vector<realT> keyValue1;
 *
 *    ///To allow ADIobservation to check for errors.
 *    bool isSetup()
 *    {
 *      if( <any condition indicating not set up>) return false;
 *      return true;
 *    }
 *
 *    //Method called by HCIobservation to get keyword-values
 *    void extractKeywords(vector<fitsHeader> & heads)
 *    {
 *       keyValue1 = headersToValues<float>(heads, "KEYWORD1");
 *    }
 *
 *    //Calculate the derotation angle for a given image number
 *    realT derotAngle(size_t imno) const
 *    {
 *       return some_function_of(keyValue1[imno]); //This function uses keyValue1[imno] to produce the derotation angle
 * in radians.
 *    }
 * };
 * \endcode
 * \endparblock
 *
 * \ingroup hc_imaging
 */
template <typename _realT, class _derotFunctObj>
struct ADIobservation : public HCIobservation<_realT>
{
    typedef _realT realT;
    typedef _derotFunctObj derotFunctObj;
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> eigenImageT;
 
    derotFunctObj m_derotF;
 
    derotFunctObj m_RDIderotF;
 
    bool m_doDerotate{ true };
 
    bool m_postMedSub{ false };
 
    ADIobservation();
 
    ADIobservation( const std::string &dir,    ///< [in] the directory to search.
                    const std::string &prefix, ///< [in] the initial part of the file name.  Can be empty "".
                    const std::string &ext     ///< [in] the extension to append to the file name, must include the '.'.
    );
 
    explicit ADIobservation( const std::string &fileListFile /**< [in] a file name path to read.*/ );
 
    ADIobservation(
        const std::string &dir,    ///< [in] the directory to search.
        const std::string &prefix, ///< [in] the initial part of the file name.  Can be empty "".
        const std::string &ext,    ///< [in] the extension to append to the file name, must include the '.'.
        const std::string &RDIdir, ///< [in] the directory to search for the reference files.
        const std::string
            &RDIprefix, ///< [in] the initial part of the file name for the reference files.  Can be empty "".
        const std::string &RDIext = "" ///< [in] [optional] the extension to append to the RDI file name, must include
                                       ///< the '.'.  If empty "" then same extension as target files is used.
    );
 
    ADIobservation( const std::string &fileListFile,   ///< [in] a file name path to read for the target file names.
                    const std::string &RDIfileListFile ///< [in] a file name path to read for the reference file names.
    );
 
    /// Read in the target files
    /** First sets up the keywords, then calls HCIobservation readFiles
     */
    int readFiles();
 
    /// Actions to take after the files are first read in by HCIobservation
    /** Exracts the angle keywords from the FITS headers, and checks that a valid value is found for each file.
     * Performs fake planet injectio if configured.
     *
     * \throws
     * \returns -1 if an error occurs.
     */
    virtual int postReadFiles();
 
    /// Post target coadd actions.
    /** Here updates derotation for new average values.
     */
    virtual int postCoadd();
 
    /// Read in the RDI files
    /** First sets up the keywords, then calls HCIobservation readRDIFiles
     */
    int readRDIFiles();
 
    /// Post reference read actions, including fake injection
    virtual int postRDIReadFiles();
 
    /// Post reference coadd actions.
    /** Here updates derotation for new average values.
     */
    virtual int postRDICoadd();
 
    /// Read in already PSF-subtracted files
    /** Used to take up final processing after applying some non-klipReduce processing steps to
     * PSF-subtracted images.
     */
    int readPSFSub( const std::string &dir, const std::string &prefix, const std::string &ext, size_t nReductions );
 
    /** \name Fake Planets
     * @{
     */
    int m_fakeMethod{ HCI::single };   ///< Method for reading fake files, either HCI::single or HCI::list.
 
    std::string m_fakeFileName;        ///< FITS file containing the fake planet PSF to inject or a list of fake images
 
    std::string m_fakeScaleFileName;   ///< One-column text file containing a scale factor for each point in time.
 
    std::vector<realT> m_fakeSep;      ///< Separation(s) of the fake planet(s)
    std::vector<realT> m_fakePA;       ///< Position angles(s) of the fake planet(s)
    std::vector<realT> m_fakeContrast; ///< Contrast(s) of the fake planet(s)
 
    realT m_RDIFluxScale{ 1 };         /**< Flux scaling to apply to fake planets injected in RDI.
                                            Would depend on the assumed spectrum in SDI.*/
    realT m_RDISepScale{ 1 };          /**< Scaling to apply to fake planet separation in RDI.
                                            Would be ratio of wavelengths for SDI.*/
 
    /// Inkect the fake plants
    /**
     *  \todo should pad the fake before calling the single image version
     *  \todo throw exceptions for all errors, and switch to void
     */
    int injectFake( eigenCube<realT> &ims,              ///< [in.out] the image cube in which to inject the fakes.
                    std::vector<std::string> &fileList, ///< [in] a list of file paths used for per-image fake PSFs.  If
                                                        ///< empty, then m_fakeFileName is used.
                    derotFunctObj &derotF,
                    realT RDIfluxScale, ///< [in] the flux scaling for RDI.  In SDI, this is from the planet spectrum.
                    realT RDISepScale   ///< [in] the separation scale for RDI.  In SDI, this is the ratio of wavlengths
                                        ///< after lambda/D scaling.
    );
 
    /// Inject the fake planet into a single image
    /**
     *  \todo should pad the fake before this point
     *  \todo throw exceptions for all errors, and switch to void
     */
    int injectFake( eigenImageT &fakePSF,
                    eigenCube<realT> &ims,
                    int image_i,
                    realT derotAngle,
                    realT PA,
                    realT sep,
                    realT contrast,
                    realT scale,
                    realT RDIfluxScale,
                    realT RDISepScale );
 
    /// @}
 
    void stdFitsHeader( fits::fitsHeader *head );
 
    virtual void makeMaskCube();
 
    /// De-rotate the PSF subtracted images
    void derotate();
 
    double t_fake_begin{ 0 };
    double t_fake_end{ 0 };
 
    double t_derotate_begin{ 0 };
    double t_derotate_end{ 0 };
};
 
template <typename _realT, class _derotFunctObj>
ADIobservation<_realT, _derotFunctObj>::ADIobservation()
{
}
 
template <typename _realT, class _derotFunctObj>
ADIobservation<_realT, _derotFunctObj>::ADIobservation( const std::string &dir,
                                                        const std::string &prefix,
                                                        const std::string &ext )
    : HCIobservation<realT>( dir, prefix, ext )
{
}
 
template <typename _realT, class _derotFunctObj>
ADIobservation<_realT, _derotFunctObj>::ADIobservation( const std::string &fileListFile )
    : HCIobservation<realT>( fileListFile )
{
}
 
template <typename _realT, class _derotFunctObj>
ADIobservation<_realT, _derotFunctObj>::ADIobservation( const std::string &dir,
                                                        const std::string &prefix,
                                                        const std::string &ext,
                                                        const std::string &RDIdir,
                                                        const std::string &RDIprefix,
                                                        const std::string &RDIext )
    : HCIobservation<realT>( dir, prefix, ext, RDIdir, RDIprefix, RDIext )
{
}
 
template <typename _realT, class _derotFunctObj>
ADIobservation<_realT, _derotFunctObj>::ADIobservation( const std::string &fileListFile,
                                                        const std::string &RDIfileListFile )
    : HCIobservation<realT>( fileListFile, RDIfileListFile )
{
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::readFiles()
{
    this->m_keywords.clear();
 
    if( !m_derotF.isSetup() )
    {
        mxError( "ADIobservation::readFiles", MXE_PARAMNOTSET, "Derotator is not configured." );
        return -1;
    }
 
    /*----- Append the ADI keywords to propagate them if needed -----*/
 
    for( size_t i = 0; i < m_derotF.m_keywords.size(); ++i )
    {
        this->m_keywords.push_back( m_derotF.m_keywords[i] );
    }
 
    if( HCIobservation<realT>::readFiles() < 0 )
    {
        return -1;
    }
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::postReadFiles()
{
    std::optional<std::vector<size_t>> bad = m_derotF.extractKeywords( this->m_heads );
 
    if( bad )
    {
        for( size_t n = 0; n < bad->size(); ++n )
        {
            std::cerr << this->m_fileList[( *bad )[n]] << " conversion failed for " << m_derotF.m_angleKeyword << "\n";
        }
 
        mxThrowException( mx::err::invalidarg,
                          "ADIobservation::postReadFiles",
                          "bad derotation angles in FITS header" );
    }
 
    if( m_fakeFileName != "" && !this->m_skipPreProcess )
    {
        std::cerr << "Injecting fakes in target images...\n";
        if( injectFake( this->m_tgtIms, this->m_fileList, m_derotF, 1, 1 ) < 0 )
        {
            return -1;
        }
    }
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::postCoadd()
{
    m_derotF.extractKeywords( this->m_heads );
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::readRDIFiles()
{
    this->m_RDIkeywords.clear();
 
    if( !m_RDIderotF.isSetup() )
    {
        mxError( "ADIobservation::readRDIFiles", MXE_PARAMNOTSET, "Derotator is not configured." );
        return -1;
    }
 
    /*----- Append the ADI keywords to propagate them if needed -----*/
 
    for( size_t i = 0; i < m_RDIderotF.m_keywords.size(); ++i )
    {
        this->m_RDIkeywords.push_back( m_RDIderotF.m_keywords[i] );
    }
 
    if( HCIobservation<realT>::readRDIFiles() < 0 )
        return -1;
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::postRDIReadFiles()
{
    m_RDIderotF.extractKeywords( this->m_RDIheads );
 
    /*if(m_fakeFileName != ""  && !this->m_skipPreProcess)
    {
       if( injectFake(this->m_refIms, this->m_RDIfileList, m_RDIderotF, m_RDIFluxScale, m_RDISepScale) < 0) return -1;
    }*/
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::postRDICoadd()
{
    m_RDIderotF.extractKeywords( this->m_RDIheads );
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::readPSFSub( const std::string &dir,
                                                        const std::string &prefix,
                                                        const std::string &ext,
                                                        size_t nReductions )
{
    // Load first file to condigure based on its header.
    std::vector<std::string> flist = ioutils::getFileNames( dir, prefix, "000", ext );
 
    fits::fitsHeader fh;
    eigenImage<realT> im;
    fits::fitsFile<realT> ff;
 
    ff.read( im, fh, flist[0] );
 
    if( !m_derotF.isSetup() )
    {
        mxError( "ADIobservation::readFiles", MXE_PARAMNOTSET, "Derotator is not configured." );
        return -1;
    }
 
    if( fh.count( "POSTMEDS" ) != 0 )
    {
        m_postMedSub = fh["POSTMEDS"].Int();
        std::cerr << "postMedSub: " << m_postMedSub << "\n";
    }
 
    if( fh.count( "FAKEFILE" ) != 0 )
    {
        m_fakeFileName = fh["FAKEFILE"].String();
        std::cerr << "fakeFileName: " << m_fakeFileName << "\n";
    }
 
    if( fh.count( "FAKESCFL" ) != 0 )
    {
        m_fakeScaleFileName = fh["FAKESCFL"].String();
        std::cerr << "fakeScaleFileName: " << m_fakeScaleFileName << "\n";
    }
 
    if( fh.count( "FAKESEP" ) != 0 )
    {
        ioutils::parseStringVector( m_fakeSep, fh["FAKESEP"].String(), "," );
 
        if( m_fakeSep.size() == 0 )
        {
            mxError( "KLIPReduction", MXE_PARSEERR, "FAKESEP vector did not parse correctly." );
            return -1;
        }
        std::cerr << "fakeSep: " << fh["FAKESEP"].String() << "\n";
    }
 
    if( fh.count( "FAKEPA" ) != 0 )
    {
        ioutils::parseStringVector( m_fakePA, fh["FAKEPA"].String(), "," );
 
        if( m_fakePA.size() == 0 )
        {
            mxError( "KLIPReduction", MXE_PARSEERR, "FAKEPA vector did not parse correctly." );
            return -1;
        }
        std::cerr << "fakePA: " << fh["FAKEPA"].String() << "\n";
    }
 
    if( fh.count( "FAKECONT" ) != 0 )
    {
        ioutils::parseStringVector( m_fakeContrast, fh["FAKECONT"].String(), "," );
 
        if( m_fakeContrast.size() == 0 )
        {
            mxError( "KLIPReduction", MXE_PARSEERR, "FAKECONT vector did not parse correctly." );
            return -1;
        }
        std::cerr << "fakeContrast: " << fh["FAKECONT"].String() << "\n";
    }
 
    /*----- Append the ADI keywords to propagate them if needed -----*/
 
    for( size_t i = 0; i < m_derotF.m_keywords.size(); ++i )
    {
        this->m_keywords.push_back( m_derotF.m_keywords[i] );
    }
 
    if( HCIobservation<realT>::readPSFSub( dir, prefix, ext, nReductions ) < 0 )
        return -1;
 
    m_derotF.extractKeywords( this->m_heads );
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::injectFake( eigenCube<realT> &ims,
                                                        std::vector<std::string> &fileList,
                                                        derotFunctObj &derotF,
                                                        realT RDIFluxScale,
                                                        realT RDISepScale )
{
    t_fake_begin = sys::get_curr_time();
 
    // typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imT;
    eigenImageT fakePSF;
    std::vector<std::string> fakeFiles; // used if m_fakeMethod == HCI::list
 
    fits::fitsFile<realT> ff;
    std::ifstream scaleFin; // for reading the scale file.
 
    // Fake Scale -- default to 1, read from file otherwise
    std::vector<realT> fakeScale( ims.planes(), 1.0 );
    if( m_fakeScaleFileName != "" )
    {
        std::vector<std::string> sfileNames;
        std::vector<realT> imS;
 
        // Read the scale file and load it into a map
        if( ioutils::readColumns( m_fakeScaleFileName, sfileNames, imS ) < 0 )
            return -1;
 
        if( sfileNames.size() != imS.size() )
        {
            std::cerr << "fake scale file must be two columns of: fileName scale\n";
            exit( -1 );
        }
        std::map<std::string, realT> scales;
        for( size_t i = 0; i < sfileNames.size(); ++i )
            scales[ioutils::pathFilename( sfileNames[i].c_str() )] = imS[i];
 
        for( size_t i = 0; i < fileList.size(); ++i )
        {
            if( scales.count( ioutils::pathFilename( fileList[i].c_str() ) ) > 0 )
            {
                fakeScale[i] = scales[ioutils::pathFilename( fileList[i].c_str() )];
            }
            else
            {
                std::cerr << "File name not found in fakeScaleFile:\n";
                std::cerr << ioutils::pathFilename( fileList[i].c_str() ) << "\n";
                exit( -1 );
            }
        }
    } // if(fakeScaleFileName != "")
 
    if( m_fakeMethod == HCI::single )
    {
        if( ff.read( fakePSF, m_fakeFileName ) < 0 )
            return -1;
    }
 
    if( m_fakeMethod == HCI::list )
    {
        if( ioutils::readColumns( m_fakeFileName, fakeFiles ) < 0 )
            return -1;
    }
 
    for( int i = 0; i < ims.planes(); ++i )
    {
        if( m_fakeMethod == HCI::list )
        {
            ff.read( fakePSF, fakeFiles[i] );
        }
 
        for( size_t j = 0; j < m_fakeSep.size(); ++j )
        {
            if( injectFake( fakePSF,
                            ims,
                            i,
                            derotF.derotAngle( i ),
                            m_fakePA[j],
                            m_fakeSep[j],
                            m_fakeContrast[j],
                            fakeScale[j],
                            RDIFluxScale,
                            RDISepScale ) < 0 )
                return -1;
        }
    }
 
    t_fake_end = sys::get_curr_time();
 
    return 0;
}
 
template <typename _realT, class _derotFunctObj>
int ADIobservation<_realT, _derotFunctObj>::injectFake( eigenImageT &fakePSF,
                                                        eigenCube<realT> &ims,
                                                        int image_i,
                                                        realT derotAngle,
                                                        realT PA,
                                                        realT sep,
                                                        realT contrast,
                                                        realT scale,
                                                        realT RDIFluxScale,
                                                        realT RDISepScale )
{
    // Check for correct sizing
    if( ( fakePSF.rows() < ims.rows() && fakePSF.cols() >= ims.cols() ) ||
        ( fakePSF.rows() >= ims.rows() && fakePSF.cols() < ims.cols() ) )
    {
        mxThrowException( err::sizeerr,
                          "ADIobservation::injectFake",
                          "fake PSF has different dimensions and can't be sized properly" );
    }
 
    // Check if fake needs to be padded out
    if( fakePSF.rows() < ims.rows() && fakePSF.cols() < ims.cols() )
    {
        eigenImageT pfake( ims.rows(), ims.cols() );
        padImage( pfake, fakePSF, 0.5 * ( ims.rows() - fakePSF.rows() ), 0 );
        fakePSF = pfake;
    }
 
    // Check if fake needs to be cut down
    if( fakePSF.rows() > ims.rows() && fakePSF.cols() > ims.cols() )
    {
        eigenImageT cfake( ims.rows(), ims.cols() );
        cutPaddedImage( cfake, fakePSF, 0.5 * ( fakePSF.rows() - ims.rows() ) );
        fakePSF = cfake;
    }
 
    if( fakePSF.rows() != ims.rows() || fakePSF.cols() != ims.cols() )
    {
        mxThrowException(
            err::sizeerr,
            "ADIobservation::injectFake",
            "fake PSF has different dimensions and can't be sized properly (is it even in rows and cols?)" );
    }
 
    /*** Now shift to the separation and PA, scale, apply contrast, and inject ***/
    // allocate shifted fake psf
    eigenImageT shiftFake( fakePSF.rows(), fakePSF.cols() );
 
    realT ang, dx, dy;
 
    ang = math::dtor( -1 * PA ) + derotAngle;
 
    dx = sep * RDISepScale * sin( ang );
    dy = sep * RDISepScale * cos( ang );
 
    imageShift( shiftFake, fakePSF, dx, dy, cubicConvolTransform<realT>() );
 
    ims.image( image_i ) = ims.image( image_i ) + shiftFake * scale * RDIFluxScale * contrast;
 
    return 0;
}
 
template <typename realT, class derotFunctObj>
void ADIobservation<realT, derotFunctObj>::makeMaskCube()
{
    if( this->m_mask.rows() != this->m_Nrows || this->m_mask.cols() != this->m_Ncols )
    {
        // clang-format off
        std::string message = "Mask is not the same size as images.\n";
                   message += "    Mask:   rows=" + std::to_string(this->m_mask.rows()) + "\n";
                   message += "            cols=" + std::to_string(this->m_mask.cols()) + "\n";
                   message += "    Images: rows=" + std::to_string(this->m_Nrows) + "\n";
                   message += "            cols=" + std::to_string(this->m_Ncols) + "\n";
        // clang-format on
 
        mxThrowException( err::invalidconfig, "ADIobservation<realT, derotFunctObj>::makeMaskCube", message );
    }
 
    this->m_maskCube.resize( this->m_Nrows, this->m_Ncols, this->m_Nims );
 
#pragma omp parallel
    {
        eigenImageT rm;
 
#pragma omp for
        for( int i = 0; i < this->m_Nims; ++i )
        {
            rotateMask( rm, this->m_mask, m_derotF.derotAngle( i ) );
            this->m_maskCube.image( i ) = rm;
        }
    }
 
    ioutils::createDirectories( this->m_auxDataDir );
    std::ofstream fout( this->m_auxDataDir + "angles.dat" );
    for( int i = 0; i < this->m_Nims; ++i )
    {
        fout << m_derotF.derotAngle( i ) << "\n";
    }
    fout.close();
 
    fits::fitsFile<realT> ff;
    ff.write( this->m_auxDataDir + "maskCube.fits", this->m_maskCube );
}
 
template <typename _realT, class _derotFunctObj>
void ADIobservation<_realT, _derotFunctObj>::derotate()
{
    t_derotate_begin = sys::get_curr_time();
 
    for( size_t n = 0; n < this->m_psfsub.size(); ++n )
    {
#pragma omp parallel
        {
            eigenImageT rotim;
            realT derot;
 
#pragma omp for
            for( int i = 0; i < this->m_psfsub[n].planes(); ++i )
            {
                derot = m_derotF.derotAngle( i );
                if( derot != 0 )
                {
                    imageRotate( rotim, this->m_psfsub[n].image( i ), derot, cubicConvolTransform<realT>() );
                    this->m_psfsub[n].image( i ) = rotim;
                }
            }
        }
    }
 
    t_derotate_end = sys::get_curr_time();
}
 
// If fakeFileName == "" or skipPreProcess == true then use the structure of propagated values
 
template <typename _realT, class _derotFunctObj>
void ADIobservation<_realT, _derotFunctObj>::stdFitsHeader( fits::fitsHeader *head )
{
    if( head == 0 )
        return;
 
    head->append( "", fits::fitsCommentType(), "----------------------------------------" );
    head->append( "", fits::fitsCommentType(), "mx::ADIobservation parameters:" );
    head->append( "", fits::fitsCommentType(), "----------------------------------------" );
 
    head->append( "POSTMEDS", m_postMedSub, "median subtraction after processing" );
 
    if( m_fakeFileName != "" )
    {
        head->append( "FAKEFILE", m_fakeFileName, "name of fake planet PSF file" );
 
        if( m_fakeScaleFileName != "" )
        {
            head->append( "FAKESCFL", m_fakeScaleFileName, "name of fake planet scale file name" );
        }
 
        std::stringstream str;
 
        if( m_fakeSep.size() > 0 )
        {
            for( size_t nm = 0; nm < m_fakeSep.size() - 1; ++nm )
            {
                str << m_fakeSep[nm] << ",";
            }
            str << m_fakeSep[m_fakeSep.size() - 1];
            head->append<char *>( "FAKESEP", (char *)str.str().c_str(), "separation of fake planets" );
        }
 
        if( m_fakePA.size() > 0 )
        {
            str.str( "" );
            for( size_t nm = 0; nm < m_fakePA.size() - 1; ++nm )
            {
                str << m_fakePA[nm] << ",";
            }
            str << m_fakePA[m_fakePA.size() - 1];
            head->append<char *>( "FAKEPA", (char *)str.str().c_str(), "PA of fake planets" );
        }
 
        if( m_fakeContrast.size() > 0 )
        {
            str.str( "" );
            for( size_t nm = 0; nm < m_fakeContrast.size() - 1; ++nm )
            {
                str << m_fakeContrast[nm] << ",";
            }
            str << m_fakeContrast[m_fakeContrast.size() - 1];
            head->append<char *>( "FAKECONT", (char *)str.str().c_str(), "Contrast of fake planets" );
        }
    }
}
 
template <typename realT>
class ADIDerotator;
 
extern template class ADIobservation<float, ADIDerotator<float>>;
extern template class ADIobservation<double, ADIDerotator<double>>;
 
} // namespace improc
} // namespace mx
 
#endif //__ADIobservation_hpp__