wiki:Dc3bPsfApi
Last modified 10 years ago Last modified on 09/17/2009 10:27:54 PM

Proposed Changes to Psf/Kernel? API for DC3b

None of the stuff below is complete - it's just pseudocode highlighting the differences from existing code and specifying how UC Davis would like to use some of these classes. It will probably be useful to have source:/DMS/afw/trunk/include/lsst/afw/math/Kernel.h open when looking at the code below.

The main use case for kernel/psf redesign is for model convolution. This is mathematically a symmetric ooperation, so it does not make any more sense for kernels to know about models than it does for models to know about kernels. Model convolution requires internal knowledge about both models and kernels, but this can be simplified by rephrasing it as two questions: what type of convolution operation to perform, and what parts of the model need to be convolved. As we can only conceive of a very limited number of types of convolutions, the first question is easier to answer. The proposal below utilizes the GoF visitor design pattern, with each convolution type represented as different visitor subclass. The seconds question is addressed in the Model::convolve method of each model subclass, in which it is then the model's responsibility to apply the convolution to specific parts of itself, without exposing the details of what sort of things actually get convolved.

The new interface involves an as-yet-undefined FourierCutout object we expect will behave more or less like the Image class. It will encapsulate the packing system FFTW uses for real-to-complex transforms. We'll make a more concrete proposal for what it should have later.

ConvolutionVisitor

/** Just a typedef for use below. */
typedef Eigen::MatrixXd CovarianceMatrix;


/** 
 *  All model convolution goes through a visitor that defines one of
 *  the (we hope very few) different ways to convolve a model.  We can
 *  imagine no more than five of these at present (only one of which
 *  will see action in DC3b).
 *
 *  Both a model and PSF will be asked which convolution methods they
 *  support; a PSF may support many, but a particular model class will
 *  only ever support one.  The free function setModelPsf, defined
 *  below, acts as a kind of visitor factory, checking if any of the
 *  convolution methods supported by an instance of a psf are
 *  compatible with a model, and asking the Psf to create a
 *  convolution visitor if possible.  This visitor is then passed to
 *  the Model, providing it all the information it needs to convolve
 *  itself.
 *
 *  The "local linear combination of functions" discussed at the last
 *  meeting is a combination of special cases of two of these visitors
 *  (IMAGE and FOURIER).  A linear combination is less fundamental
 *  than we thought (we actually just need to compute the derivative
 *  of a kernel with respect to its local parameters).  Splitting
 *  regular image convolution from Fourier-space convolution also
 *  seemed like a good idea.
 */
class ConvolutionVisitor {
    boost:::shared_ptr<CovarianceMatrix> _covariance;    
public:

    /**
     * The types of convolution are one-hot encoded so that they can be used
     * as a bit mask by kernels to list all the convolution types it supports.
     * \note the number of types of convolution is expected to remain small.
     * We may find the types defined below are already more than we will 
     * ever need.   
     */
    enum TypeFlag { IMAGE=0x1, FOURIER=0x2, MULTIGAUSSIAN=0x4, SHAPELET=0x8 };

    virtual void visit(Model & model) = 0;


    /**
     *  This allows Psf to set a covariance matrix for a visitor after
     *  it has been constructed by a Kernel (which has no knowledge of
     *  uncertainty, and hence initializes the covariance to zero).
     */
    void setCovarianceMatrix(boost::shared_ptr<CovarianceMatrix> covariance) {
        _covariance = covariance;
    }
    boost::shared_ptr<CovarianceMatrix> getCovariance() const {
        return _covariance;
    }
};

/**
 *  ConvolutionVisitor corresponding to models that construct their
 *  unconvolved selves as a regular images and convolve in real space.
 */
class ImageConvolutionVisitor : public ConvolutionVisitor {
public:

    /**
     *  Create an image of the kernel.
     */
    void getImage(Image::Ptr image);

    /**
     *  Create images of the derivative of the kernel w.r.t. its local
     *  parameters (evaluated at the values of those parameters).
     *
     *  For a LinearCombinationKernel, these are just the images of
     *  the KernelList.
     */
    void getDerivative(std::vector<Image::Ptr> & derivative);

    /**
     *  Get the local parameters of the kernel.
     */
    std::vector<double> getParameters() const;

    virtual void visit(Model & model) { model.convolve(*this); }
};

/**
 *  ConvolutionVisitor corresponding to models that construct
 *  their unconvolved selves in Fourier space.
 */
class FourierConvolutionVisitor : public ConvolutionVisitor {
public:

    /**
     *  Create a Fourier-space image of the kernel.
     */
    void getFourierImage(FourierCutout & fourierCutout);


    /**
     *  Create Fourier-space images of the derivative of the kernel
     *  w.r.t. its local parameters (evaluated at the values of those
     *  parameters).
     *
     *  For a LinearCombinationKernel, these are just the Fourier
     *  images of the KernelList.
     */
    void getDerivative(std::vector<FourierCutout::Ptr> derivative);

    /**
     *  Get the local parameters of the kernel.
     */
    std::vector<double> getParameters() const;

    virtual void visit(Model & model) { model.convolve(*this); }
};

/**
 *  ConvolutionVisitor corresponding to models that are sums of
 *  Gaussians and want to convolve themselves with Kernels that are
 *  sums of Gaussians.
 *
 *  We probably won't implement this until it's actually needed; it's
 *  only necessary for Gaussian-based models, not Gaussian-based PSFs.
 */
class MultiGaussianConvolutionVisitor : public ConvolutionVisitor {
public:
    
    /**
     *  Get the local parameters of the kernel.
     *
     *  In this case that would be a flattened sequence of
     *  (amplitude,xx,yy,xy) tuples or something.
     */
    std::vector<double> getParameters() const;

    virtual void visit(Model & model) { model.convolve(*this); }
};


/**
 *  ConvolutionVisitor corresponding to models that are represented
 *  in shapelet space.
 *
 *  We probably won't implement this until it's actually needed; it's
 *  only necessary for shapelet based models, not shapelet-based PSFs.
 */
class ShapeletConvolutionVisitor : public ConvolutionVisitor {
public:

    /**
     *  Get the scale factor of the shapelet expansion.
     */
    double getScale() const;

    /**
     *  Get the local parameters of the kernel.
     *
     *  In this case that would be the shapelet coefficients.
     */
    std::vector<double> getParameters() const;

    virtual void visit(Model & model) { model.convolve(*this); }
};

Kernel

/**
 *  This is the base Kernel class, with only three new member functions.
 *  All the new functions have default implementations.
 */
class Kernel {
public:

    /**
     *  Compute a Fourier-space representation of the Kernel.
     *
     *  A default implementation in the base class will just FFT the
     *  output of computeImage.
     */
    virtual void computeFourierImage(
        FourierCutout::Ptr fourierImage,
        bool const doNormalize,
        double const x = 0;
        double const y = 0
    ) const;
    
    /**
     *  Return a bitwise OR of all supported convolution methods.
     *
     *  All Kernels will support at least IMAGE and FOURIER convolution.
     */
    virtual ConvolutionVisitor::TypeFlag getConvolutionTypeFlags() const { 
        return ConvolutionVisitor::IMAGE | ConvolutionVisitor::FOURIER;
    }

    /**
     *  Return a ConvolutionVisitor that matches the type requested,
     *  at the given location.
     *
     *  The default implementation would support the creation of IMAGE
     *  and FOURIER visitors with zeroed derivatives.
     *  LinearCombinationKernel (and possibly AnalyticKernel) can
     *  override to provide versions with derivatives.  Gaussian-based
     *  analytic kernels would of course support  MultiGaussianVisitor,
     *  and someday a ShapeletKernel would be capable of providing a 
     *  ShapeletVisitor.  
     */
    virtual ConvolutionVisitor::Ptr computeConvolutionVisitor(
        ConvolutionVisitor::TypeFlag visitorType,
        Point const & location,
    ) const;

};

PSF

/**
 *  Psf can (for now) be a concrete class that simply owns a kernel.
 *  Once we figure out how to compute uncertainty information for a
 *  Kernel, we can store that here too.
 *
 *  In fact, since we won't be implementing either of the two things
 *  Psf will do that Kernel doesn't, we don't need Psf yet for
 *  DC3b, but here's how it might look.
 */
class Psf {
public:
    
    typedef double PsfColor; ///< not used for now

    /**
     *  Return a bitwise OR of all supported convolution methods.
     *
     *  All Psfs will support at least IMAGE and FOURIER convolution.
     */
    ConvolutionVisitor::TypeFlag getConvolutionTypeFlags() const { 
        return _kernel->getConvolutionFlags();
    }

    /**
     *  Get a Kernel for a given PsfColor.
     */
    Kernel::Ptr getKernel(PsfColor color=PsfColor()) const;

    /**
     *  Return a ConvolutionVisitor that matches the type requested at
     *  the given location and color.
     */
    virtual ConvolutionVisitor::Ptr computeConvolutionVisitor(
        ConvolutionVisitor::TypeFlag visitorType,
        Point const & location,
        PsfColor color=PsfColor()
    ) const {
        ConvolutionVisitor::Ptr visitor = getKernel(color)->computeConvolutionVisitor(visitorType,location);
        // ... compute covariance matrix ... 
        visitor.setCovariance(/* covariance matrix */);
        return visitor;
    }
};

Model

/**
 *  The Model base class below is a convolved model that is
 *  "assigned" a Psf/Kernel by the free function setModelPsf().  Some
 *  of the model's behavior (which is defined by a separate interface
 *  relating to the fitting algorithm) may be undefined unless the
 *  Model's Psf has been set.  Once the Psf is set, the model can be
 *  convolved by calling its "convolve()" member function.
 *
 *  There are two major reasons for this design, as opposed to a
 *  "Kernel::convolve(Model)" or free function
 *  "convolve(Kernel,Model)" design:
 * 
 *  - In the course of fitting, a single Model object will be asked to
 *  update its parameters and re-evaluate itself (including
 *  convolution) many times without ever changing the Psf.  This
 *  design does not require potentially expensive convolution
 *  intermediates to be computed at every iteration.
 *
 *  - We do not presently have a clear idea of all the different
 *  model-fitting algorithms our API will support.  They may involve
 *  convolving not only the models, but their derivatives.  They may
 *  involve fitting to Fourier- or shapelet-transformed data.  And
 *  they may or may not want to propogate PSF uncertainty into the
 *  model uncertainty.  These complexities make it difficult for
 *  models to expose a simple interface to any external convolver.
 *  Moreover, the particular way in which we convolve affects not just
 *  the speed of the operation, but also the accuracy of the
 *  convolution. This makes it crucial for models to specify the mode
 *  of their convolution.
 */
class Model {
public:
    /**
     *  Whether the Model has been provided with a Psf.
     */
    virtual bool hasPsf() const = 0;

    /**
     *  For our purposes, this is the location the Psf will be evaluated at.
     *
     *  This might vary in the course of fitting the Model, but probably not
     *  enough to merit updating the Psf (checking for such an occurrence is
     *  something the fitting algorithm would have to be responsible for.
     */
    virtual Point getCenter() const = 0;

    /**
     *  This just defines the color of the object for color-dependent
     *  Psf purposes.  Same notes as getCenter() apply.
     */
    virtual Psf::PsfColor getPsfColor() const = 0;
private:

    /**
     *  Return the convolution method this model supports (only one allowed).
     */
    virtual ConvolutionVisitor::TypeFlag getConvolutionType() const = 0;

    /**
     *  Model subclasses will override one of the following (the one that
     *  matches the return of getConvolutionType()).  
     *  The others will throw exceptions.
     */ 
    virtual void convolve(ImageConvolutionVisitor & visitor);
    virtual void convolve(FourierConvolutionVisitor & visitor);
    virtual void convolve(MultiGaussianConvolutionVisitor & visitor);
    virtual void convolve(ShapeletConvolutionVisitor & visitor);

};

/**
 *  Free function to set the Psf or kernel of a Model. 
 *  Checks for compatibility between the convolution type(s) supported
 *  by the psf/kernel, and the type of convolution supported by the model.
 *  If a compatible convolution type is found, the appropriate
 *  ConvolutionVisitor is created and set in the model.
 *
 *  There's actually nothing at present to keep this from being a member function
 *  of Model, but since it is where the double-dispatch starts we put it outside 
 *  the class for now.
 */
void setModelPsf(Psf::ConstPtr psf, Model & model) {
    int visitorType = psf->getConvolutionTypeFlags() & model.getConvolutionType();
    if (visitorType == 0) throw ConvolutionImpossible();
    ConvolutionVisitor::Ptr visitor = psf->computeConvolutionVisitor(
        (ConvolutionVisitor::TypeFlag)visitorType
        Model.getCenter(),
        Model.getPsfColor()
    );
    visitor->visit(model);
}

void setModelKernel(Kernel::ConstPtr kernel, Model & model) {
    int visitorType = kernel->getConvolutionTypeFlags() & model->getConvolutionType();
    if (visitorType == 0) throw ConvolutionImpossible();
    ConvolutionVisitor::Ptr visitor = kernel->computeConvolutionVisitor(
        (ConvolutionVisitor::TypeFlag)visitorType
        Model.getCenter(),
        Model.getPsfColor()
    );
    visitor->visit(model);
    if(visitorType != 0)
}