GRALE
imageplane.h
1 /*
2 
3  This file is a part of GRALE, a library to facilitate the simulation
4  and inversion of gravitational lenses.
5 
6  Copyright (C) 2008-2012 Jori Liesenborgs
7 
8  Contact: jori.liesenborgs@gmail.com
9 
10  This program is free software; you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation; either version 2 of the License, or
13  (at your option) any later version.
14 
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23 
24 */
25 
26 #ifndef GRALE_IMAGEPLANE_H
27 
28 #define GRALE_IMAGEPLANE_H
29 
30 #include "graleconfig.h"
31 #include "vector2d.h"
32 #include "gridfunction.h"
33 //#include "randomnumbergenerator.h"
34 //#include "simpleuniformdistribution.h"
35 #include "lensplane.h"
36 #include <enut/networklayeraddress.h>
37 #include <serut/serializationinterface.h>
38 #include <stdint.h>
39 #include <stdio.h>
40 #include <vector>
41 #include <list>
42 #include <string>
43 
44 namespace grale
45 {
46 
47 class GravitationalLens;
48 class LensPlane;
49 class SourcePlane;
50 class OpenCLKernel;
51 
52 class GRALE_IMPORTEXPORT ImagePlane : public errut::ErrorBase
53 {
54 public:
55  ImagePlane(const GravitationalLens *pLens, SourcePlane *sourceplane);
56  virtual ~ImagePlane();
57  bool init(Vector2D<double> bottomleft,Vector2D<double> topright,int xpoints,int ypoints, OpenCLKernel *pCLKernel = 0);
58  bool init(Vector2D<double> bottomleft,Vector2D<double> topright,int xpoints,int ypoints,
59  const nut::NetworkLayerAddress &serverAddress, uint16_t serverPort);
60  bool init(LensPlane *lensplane);
61  bool isInit() const { if (m_betas.size() == 0) return false; return true; }
62  const GravitationalLens *getLens() const { return m_pLens; }
63  SourcePlane *getSourcePlane() { return sourceplane; }
64 
65  int getNumXPoints() const { return numx; }
66  int getNumYPoints() const { return numy; }
67  int getNumXPixels() const { return numx-1; }
68  int getNumYPixels() const { return numy-1; }
69  Vector2D<double> getBottomLeft() const { return bottomleft; }
70  Vector2D<double> getTopRight() const { return topright; }
71  double getXStep() const { return xstep; }
72  double getYStep() const { return ystep; }
73 
74  double getImageIntensityAccurate(int xPixel, int yPixel, int subSamples = 10) const;
75  double getSourceIntensityAccurate(int xPixel, int yPixel, int subSamples = 10) const;
76 
77  double getInverseMagnification(int xpos,int ypos) const { double axx = m_alphaxx[xpos+ypos*numx]; double ayy = m_alphayy[xpos+ypos*numx]; double axy = m_alphaxy[xpos+ypos*numx]; return (1.0-axx)*(1.0-ayy)-axy*axy; }
78  void getInverseMagnificationComponents(int xpos,int ypos, double &kappa, double &gamma1, double &gamma2) const { double axx = m_alphaxx[xpos+ypos*numx]; double ayy = m_alphayy[xpos+ypos*numx]; double axy = m_alphaxy[xpos+ypos*numx]; kappa = 0.5*(axx+ayy); gamma1 = 0.5*(axx-ayy); gamma2 = axy; }
79 
80  const std::list<Vector2D<double> > &getCriticalLinePoints() const { return criticallines; }
81  const std::list<Vector2D<double> > &getCausticPoints() const { return caustics; }
82  const std::list<std::list<Vector2D<double> > *> &getCriticalLineSegments() const { return criticallines2; }
83  const std::list<std::list<Vector2D<double> > *> &getCausticSegments() const { return caustics2; }
84 
85  Vector2D<double> getIndexCoordinate(int xpos, int ypos) const;
86  Vector2D<double> getPixelCoordinate(int xPixel, int yPixel) const;
87 
88  Vector2D<double> traceIndexCoordinate(int xpos, int ypos) const;
89  bool getIndices(Vector2D<double> v,int *xpos,int *ypos) const;
90 
91  bool traceBeta(Vector2D<double> beta, std::list<Vector2D<double> > &thetaPoints, int numIterations = 4, bool check = true) const;
92  bool traceTheta(Vector2D<double> theta, Vector2D<double> *pBeta) const;
93  bool traceThetaApprox(Vector2D<double> theta, Vector2D<double> *pBeta) const;
94  double getInverseMagnification(Vector2D<double> theta) const;
95  double getInverseMagnificationApprox(Vector2D<double> theta) const;
96 
97  bool write(serut::SerializationInterface &si,bool writesources = false) const;
98  static bool read(serut::SerializationInterface &si,ImagePlane **ip,std::string &errstr);
99  static bool load(const std::string &fname,ImagePlane **ip,std::string &errstr);
100  bool save(const std::string &fname,bool writesources = false) const;
101 protected:
102  virtual void setFeedbackStatus(const std::string &msg) { }
103  virtual void setFeedbackPercentage(int pct) { }
104 private:
105  class InternalLensPlane : public LensPlane
106  {
107  public:
108  InternalLensPlane(const GravitationalLens *pLens, ImagePlane *pImgPlane) : LensPlane(pLens)
109  {
110  m_pImgPlane = pImgPlane;
111  }
112 
113  ~InternalLensPlane()
114  {
115  }
116 
117  void setFeedbackStatus(const std::string &msg)
118  {
119  m_pImgPlane->setFeedbackStatus(msg);
120  }
121 
122  virtual void setFeedbackPercentage(int pct)
123  {
124  m_pImgPlane->setFeedbackPercentage(pct);
125  }
126  private:
127  ImagePlane *m_pImgPlane;
128  };
129 
130  ImagePlane() /* : m_uniDist(&m_rndGen) */ { }
131  static double getTime();
132 
133  void searchCritCaust(double *inversemagnifications);
134  bool refinePosition(Vector2D<double> beta, Vector2D<double> startTheta, Vector2D<double> &theta, int numIterations) const;
135  bool calculateBeta(Vector2D<double> theta, Vector2D<double> *pBeta) const;
136  Vector2D<double> approximateBeta(int xIdx, int yIdx, Vector2D<double> theta) const;
137  double getTriangleIntensity(Vector2D<double> p1, Vector2D<double> p2, Vector2D<double> p3,
138  double axx1, double axx2, double axx3,
139  double ayy1, double ayy2, double ayy3,
140  double axy1, double axy2, double axy3,
141  int numPoints, bool useMagnification) const;
142 
143  GravitationalLens *m_pLens;
144  SourcePlane *sourceplane;
145 
146  std::vector<Vector2D<double> > m_betas;
147  std::vector<double> m_alphaxx;
148  std::vector<double> m_alphayy;
149  std::vector<double> m_alphaxy;
150  GridFunction *m_pAlphaxxFunction;
151  GridFunction *m_pAlphayyFunction;
152  GridFunction *m_pAlphaxyFunction;
153  std::list<Vector2D<double> > criticallines;
154  std::list<Vector2D<double> > caustics;
155  std::list<std::list<Vector2D<double> > *> criticallines2;
156  std::list<std::list<Vector2D<double> > *> caustics2;
157 
158  Vector2D<double> bottomleft;
159  Vector2D<double> topright;
160  double xstep,ystep;
161  int numx,numy;
162 
163  //RandomNumberGenerator m_rndGen;
164  //SimpleUniformDistribution m_uniDist;
165 };
166 
167 inline Vector2D<double> ImagePlane::getIndexCoordinate(int xpos,int ypos) const
168 {
169  return Vector2D<double>((((double)xpos))*xstep+bottomleft.getX(),(((double)ypos))*ystep+bottomleft.getY());
170 }
171 
172 inline Vector2D<double> ImagePlane::getPixelCoordinate(int xPixel, int yPixel) const
173 {
174  return Vector2D<double>((((double)xPixel))*xstep+bottomleft.getX(),(((double)yPixel))*ystep+bottomleft.getY())+Vector2D<double>(xstep/2.0, ystep/2.0);
175 }
176 
177 inline Vector2D<double> ImagePlane::traceIndexCoordinate(int xpos,int ypos) const
178 {
179  return m_betas[xpos+ypos*numx];
180 }
181 
182 inline bool ImagePlane::getIndices(Vector2D<double> v,int *xpos,int *ypos) const
183 {
184  int x,y;
185 
186  x = (int)(((v.getX()-bottomleft.getX())/xstep)+0.5);
187  y = (int)(((v.getY()-bottomleft.getY())/ystep)+0.5);
188 
189  if (x < 0 || y < 0 || x >= numx || y >= numy)
190  return false;
191  *xpos = x;
192  *ypos = y;
193  return true;
194 }
195 
196 inline double ImagePlane::getInverseMagnificationApprox(Vector2D<double> theta) const
197 {
198  double axx = (*m_pAlphaxxFunction)(theta);
199  double ayy = (*m_pAlphayyFunction)(theta);
200  double axy = (*m_pAlphaxyFunction)(theta);
201  double inverseMagnification = ABS((1.0-axx)*(1.0-ayy)-axy*axy);
202 
203  return inverseMagnification;
204 }
205 
206 } // end namespace
207 
208 #endif // GRALE_IMAGEPLANE_H
209 
Generated on Wed Dec 19 2012 09:09:23 for GRALE by   doxygen 1.8.1.1