backprojecttreebase.h

Go to the documentation of this file.

/*

  This file is a part of GRALE, a library to facilitate the simulation
  and inversion of gravitational lenses.

  Copyright (C) 2008-2012 Jori Liesenborgs

  Contact: jori.liesenborgs@gmail.com
  
  This program 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 2 of the License, or
  (at your option) any later version.
  
  This program 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 this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  
*/

#ifndef GRALE_BACKPROJECTTREEBASE_H

#define GRALE_BACKPROJECTTREEBASE_H

#include "graleconfig.h"
#include "projectedimagesinterface.h"
#include <vector>
#include <list>

namespace grale
{

class ImagesDataExtended;
class GridNodeBase;

class GRALE_IMPORTEXPORT BackProjectTreeBase : public ProjectedImagesInterface
{
public:
        ~BackProjectTreeBase();

        void calculateDeflection(const GridNodeBase *pRootNode, bool intensities);
        void calculateBetas(float massFactor);
        void calculateMagnifications(float massFactor);

        double getAngularScale() const                                                                  { return m_angularScale; }
        const Vector2D<float> *getBetas(int sourceNumber) const                                         { return m_betaPointers[sourceNumber][0]; }
        const Vector2D<float> *getBetas(int sourceNumber, int imageNumber) const                        { return m_betaPointers[sourceNumber][imageNumber]; }
        const Vector2D<float> *getThetas(int sourceNumber) const                                        { return m_thetaPointers[sourceNumber][0]; }
        const Vector2D<float> *getThetas(int sourceNumber, int imageNumber) const                       { return m_thetaPointers[sourceNumber][imageNumber]; }
        const float *getInverseMagnifications(int sourceNumber) const                                   { return m_invMagPointers[sourceNumber][0]; }
        const float *getInverseMagnifications(int sourceNumber, int imageNumber) const                  { return m_invMagPointers[sourceNumber][imageNumber]; }
        const float *getDerivativesXX(int sourceNumber) const                                           { return 0; }
        const float *getDerivativesXX(int sourceNumber, int imageNumber) const                          { return 0; }
        const float *getDerivativesYY(int sourceNumber) const                                           { return 0; }
        const float *getDerivativesYY(int sourceNumber, int imageNumber) const                          { return 0; }
        const float *getDerivativesXY(int sourceNumber) const                                           { return 0; }
        const float *getDerivativesXY(int sourceNumber, int imageNumber) const                          { return 0; }
        const float *getShearComponents1(int sourceNumber) const                                        { return 0; }
        const float *getShearComponents1(int sourceNumber, int imageNumber) const                       { return 0; }
        const float *getShearComponents2(int sourceNumber) const                                        { return 0; }
        const float *getShearComponents2(int sourceNumber, int imageNumber) const                       { return 0; }
        const float *getConvergence(int sourceNumber) const                                             { return 0; }
        const float *getConvergence(int sourceNumber, int imageNumber) const                            { return 0; }

        float getTimeDelay(int sourceNumber, int imageNumber, int pointNumber, 
                           Vector2D<float> beta) const; // in days!
protected:
        BackProjectTreeBase();

        void init(const std::list<ImagesDataExtended *> &images,
                  Vector2D<double> center, double areaSize, int firstLevelSubdiv, int maxDepth, bool useWeights);
        virtual void getDeflectionInformation(double squareSize, Vector2D<double> theta, 
                                              Vector2D<double> &alpha, 
                                              double &alphaxx, double &alphayy, double &alphaxy) = 0;
private:
        void processNode(int numSubdiv, int level, int offset, const GridNodeBase *pNode, bool intensities);
        void fillInRows(int numSubdiv, int level, int offset, int maxDepth, Vector2D<double> center, double areaSize, const std::vector<Vector2D<double> > &unscaledThetas, const std::vector<double> &distanceFractions);

        std::vector<std::vector<std::vector<Vector2D<float> > > > m_deflectionRows;
        std::vector<std::vector<std::vector<float> > > m_magnificationRows[3];

        std::vector<Vector2D<float> > m_betas, m_thetas, m_subThetas;
        std::vector<std::vector<Vector2D<float> * > > m_betaPointers, m_thetaPointers;

        std::vector<float> m_invMag, m_subMag0, m_subMag1, m_subMag2;
        std::vector<std::vector<float *> > m_invMagPointers;

        std::vector<float> m_tmpBuf, m_tmpBuf2;
        std::vector<float> m_oneVector;
        std::vector<Vector2D<float> > m_zeroVector;

        double m_angularScale;
        int m_firstLevelSubdiv;
        int m_totalPoints;
        bool m_useWeights;
};

} // end namespace

#endif // GRALE_BACKPROJECTTREEBASE_H