[Debian-astro-commits] [gyoto] 84/221: Changed the anglekind API and "Rectilinear" option, better suited for photo-realistic ray-tracing.

Thibaut Jean-Claude Paumard thibaut at moszumanska.debian.org
Fri May 22 20:52:35 UTC 2015


This is an automated email from the git hooks/post-receive script.

thibaut pushed a commit to branch master
in repository gyoto.

commit 409df1d00db643a9b9c0bd9aa7bad3ef33c313f3
Author: Thibaut Paumard <paumard at users.sourceforge.net>
Date:   Thu Nov 6 17:17:01 2014 +0100

    Changed the anglekind API and "Rectilinear" option, better suited for photo-realistic ray-tracing.
---
 include/GyotoDefs.h     |  1 +
 include/GyotoScreen.h   |  9 +++++--
 lib/Screen.C            | 62 ++++++++++++++++++++++++++++++++++++++-----------
 yorick/gyoto.i          | 13 ++++++++---
 yorick/gyoto_Screen.C   |  5 ++--
 yorick/matte-painting.i | 32 +++++++++++++++++--------
 6 files changed, 91 insertions(+), 31 deletions(-)

diff --git a/include/GyotoDefs.h b/include/GyotoDefs.h
index 8271162..d2bf69b 100644
--- a/include/GyotoDefs.h
+++ b/include/GyotoDefs.h
@@ -550,6 +550,7 @@ namespace Gyoto {
 # define getDelta0          delta0
 # define setAlpha0          alpha0
 # define getAlpha0          alpha0
+# define setAnglekind       anglekind
 # define getTmin            tMin
 # define setTmin            tMin
 # define getTime            time
diff --git a/include/GyotoScreen.h b/include/GyotoScreen.h
index 4828812..d1a926c 100644
--- a/include/GyotoScreen.h
+++ b/include/GyotoScreen.h
@@ -183,7 +183,10 @@ class Gyoto::Screen : protected Gyoto::SmartPointee {
   double distance_; ///< Distance to the observer in m
   double dmax_; ///< Maximum distance from which the photons are launched (geometrical units) 
 
-  int anglekind_; ///< Screen angles kind (0: equatorial, 1: spherical)
+  enum anglekind_e { equatorial_angles=0, rectilinear=1, spherical_angles=2};
+  typedef int anglekind_t;
+
+  anglekind_t anglekind_; ///< Screen angles kind (0: equatorial, 1: spherical)
   
   /**
    * The angles are position angle of the line of nodes (North of
@@ -430,7 +433,9 @@ class Gyoto::Screen : protected Gyoto::SmartPointee {
   double delta0(std::string unit);
 
   /// Set Screen::anglekind_
-  void setAnglekind(int);
+  void anglekind(int);
+  void anglekind(std::string);
+  std::string anglekind() const;
 
   /// Get Screen::npix_
   size_t resolution();
diff --git a/lib/Screen.C b/lib/Screen.C
index 5b129f7..1bef1fc 100644
--- a/lib/Screen.C
+++ b/lib/Screen.C
@@ -47,7 +47,7 @@ using namespace Gyoto;
 Screen::Screen() : 
   //tobs_(0.), fov_(M_PI*0.1), tmin_(0.), npix_(1001),
   tobs_(0.), fov_(M_PI*0.1), npix_(1001), mask_(NULL), mask_filename_(""),
-  distance_(1.), dmax_(GYOTO_SCREEN_DMAX), anglekind_(0),
+  distance_(1.), dmax_(GYOTO_SCREEN_DMAX), anglekind_(equatorial_angles),
   alpha0_(0.), delta0_(0.),
   gg_(NULL), spectro_(NULL),
   freq_obs_(1.),
@@ -362,34 +362,42 @@ void Screen::spectrometer(SmartPointer<Spectrometer::Generic> spr) { spectro_=sp
 SmartPointer<Spectrometer::Generic> Screen::spectrometer() const { return spectro_; }
 
 void Screen::getRayCoord(const size_t i, const size_t j, double coord[]) const {
-  const double delta= fov_/double(npix_);
   double xscr, yscr;
 # if GYOTO_DEBUG_ENABLED
   GYOTO_DEBUG << "(i=" << i << ", j=" << j << ", coord)" << endl;
 # endif
-  if (anglekind_){
+  if (anglekind_ == Screen::spherical_angles){
     /*
       GYOTO screen labelled by spherical
       angles a and b (see Fig. in user guide)
      */
     xscr = double(i-1)*fov_/(2.*double(npix_-1));
-    yscr = double(j-1)*2.*M_PI/double(npix_-1);
-    getRayCoord(xscr,M_PI-yscr, coord);
+    yscr = M_PI-(double(j-1)*2.*M_PI/double(npix_-1));
     // NB: here xscr and yscr are the spherical angles
     // a and b ; the b->pi-b transformation boils down
     // to performing X->-X, just as below for equat angles.
-  }else{
+  }else if (anglekind_ == Screen::equatorial_angles) {
     /*
       GYOTO screen labelled by equatorial
       angles alpha and delta (see Fig. in user guide)
      */
+    const double delta= fov_/double(npix_);
     yscr=delta*(double(j)-double(npix_+1)/2.);
-    xscr=delta*(double(i)-double(npix_+1)/2.);
-    getRayCoord(-xscr, yscr, coord); 
+    xscr=-delta*(double(i)-double(npix_+1)/2.);
     // transforming X->-X (X being coord along e_1 observer vector)
     // this is due to the orientation convention of the screen 
     // (cf InitialConditions.pdf)
+  } else if (anglekind_ == Screen::rectilinear) {
+    if (fov_ >= M_PI)
+      throwError("Rectilinear projection requires fov_ < M_PI");
+    const double xfov=2.*tan(fov_*0.5);
+    const double delta= xfov/double(npix_);
+    yscr=delta*(double(j)-double(npix_+1)/2.);
+    xscr=-delta*(double(i)-double(npix_+1)/2.);
+  } else {
+    throwError("Unrecognized anglekind_");
   }
+  getRayCoord(xscr, yscr, coord); 
 }
 
 void Screen::getRayCoord(double alpha, double delta,
@@ -419,7 +427,7 @@ void Screen::getRayCoord(double alpha, double delta,
   double spherical_angle_a,
     spherical_angle_b;
   
-  if (anglekind_){
+  if (anglekind_ == spherical_angles){
     /*
       GYOTO screen labelled by spherical
       angles a and b (see Fig. in user guide)
@@ -427,7 +435,10 @@ void Screen::getRayCoord(double alpha, double delta,
 
     spherical_angle_a = alpha;
     spherical_angle_b = delta;
-  }else{
+  } else if (anglekind_ == rectilinear) {
+    spherical_angle_a = atan(sqrt(alpha*alpha+delta*delta));
+    spherical_angle_b = atan2(delta, alpha);
+  } else if (anglekind_ == equatorial_angles) {
     /*
       GYOTO screen labelled by equatorial
       angles alpha and delta (see Fig. in user guide)
@@ -1002,7 +1013,29 @@ void Screen::delta0(double fov, const string &unit) {
   delta0(fov);
 }
 
-void Screen::setAnglekind(int kind) { anglekind_ = kind; }
+void Screen::anglekind(int kind) { anglekind_ = kind; }
+void Screen::anglekind(std::string skind) {
+  if      (skind=="EquatorialAngles") anglekind_=equatorial_angles;
+  else if (skind=="SphericalAngles")  anglekind_=spherical_angles;
+  else if (skind=="Rectilinear")      anglekind_=rectilinear;
+  else throwError("Invalid string value for anglekind_");
+}
+
+std::string Screen::anglekind() const {
+  switch (anglekind_) {
+  case equatorial_angles:
+    return "EquatorialAngles";
+    break;
+  case spherical_angles:
+    return "SphericalAngles";
+    break;
+  case rectilinear:
+    return "Rectilinear";
+    break;
+  default:
+    throwError("Invalid integer value for Screen::anglekind_");
+  }
+}
 
 size_t Screen::resolution() { return npix_; }
 void Screen::resolution(size_t n) {
@@ -1097,7 +1130,7 @@ void Screen::fillElement(FactoryMessenger *fmp) {
     delete child; child = NULL;
   }
   fmp -> setParameter ("FreqObs", freq_obs_);
-  fmp -> setParameter ( anglekind_? "SphericalAngles" : "EquatorialAngles");
+  fmp -> setParameter (anglekind());
   if (mask_filename_!="") fmp->setParameter("Mask",
 					    (mask_filename_.compare(0,1,"!") ?
 					     mask_filename_ :
@@ -1179,8 +1212,9 @@ SmartPointer<Screen> Screen::Subcontractor(FactoryMessenger* fmp) {
       delta0 = atof(tc); delta0_found=1;
     }
     else if (name=="FreqObs") scr -> freqObs(atof(tc), unit);
-    else if (name=="SphericalAngles")  scr -> setAnglekind(1);
-    else if (name=="EquatorialAngles") scr -> setAnglekind(0);
+    else if (name=="SphericalAngles" ||
+	     name=="EquatorialAngles" ||
+	     name=="Rectilinear")  scr -> anglekind(name);
 #   ifdef GYOTO_USE_CFITSIO
     else if (name=="Mask")        scr -> fitsReadMask(fmp->fullPath(content));
 #   endif
diff --git a/yorick/gyoto.i b/yorick/gyoto.i
index 1dca40b..ccec27e 100644
--- a/yorick/gyoto.i
+++ b/yorick/gyoto.i
@@ -878,7 +878,8 @@ func gyoto_matte_paint(set, paint, kind=, yaw=, pitch=, roll=)
 
   bg=paint(theta, phi);
 
-  bg *= mask0(-,,);
+  if (structof(bg)==char) {bg*=mask0(-,,);}
+  else bg *= mask0;
   
   return bg;
 }
@@ -976,9 +977,15 @@ func gyoto_painters_panorama_eval(theta, phi, mask=) {
   x=long(i0 + phi*phi_scale) ;
   y=long(j0 + (theta-0.5*pi)*theta_scale) ;
 
-  if (anyof(x<=0 | y<=0 | x>nphi | y>ntheta))
+  // Rounding algo + errors can lead to an overflow by one pixel, but
+  // more would be fishy
+  if (anyof(x<-1 | y<-1 | x>(nphi+1) | y>(ntheta+1)))
     error, "something fishy in pixel computation";
-
+  if (numberof((ind=where(x==0)))) x(ind)=1;
+  if (numberof((ind=where(y==0)))) y(ind)=1;
+  if (numberof((ind=where(x==nphi+1)))) x(ind)=nphi;
+  if (numberof((ind=where(x==ntheta+1)))) x(ind)=ntheta;
+  
   dd=dimsof(phi);
   
   bg=array(char, 3, dd(2), dd(3));
diff --git a/yorick/gyoto_Screen.C b/yorick/gyoto_Screen.C
index c7bea2a..89479ff 100644
--- a/yorick/gyoto_Screen.C
+++ b/yorick/gyoto_Screen.C
@@ -39,7 +39,7 @@ extern "C" {
       "unit",
       "metric",
       "time","fov","resolution", "mask", "maskwrite",
-      "alpha0", "delta0",
+      "alpha0", "delta0", "anglekind",
       "distance", "dmax", "inclination", "paln", "argument",
       "freqobs",
       "projection", "observerpos",
@@ -49,7 +49,7 @@ extern "C" {
       "xmlwrite", "clone",
       0
     };
-    YGYOTO_WORKER_INIT1(Screen, Screen, knames, 26)
+    YGYOTO_WORKER_INIT1(Screen, Screen, knames, 27)
 
     YGYOTO_WORKER_SET_UNIT;
     YGYOTO_WORKER_GETSET_OBJECT2(metric,Metric);
@@ -100,6 +100,7 @@ extern "C" {
 
     YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(alpha0);
     YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(delta0);
+    YGYOTO_WORKER_GETSET_STRING2(anglekind);
     YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(distance);
     YGYOTO_WORKER_GETSET_DOUBLE2(dMax);
     YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(inclination);
diff --git a/yorick/matte-painting.i b/yorick/matte-painting.i
index 3a0b2cc..fdc4438 100644
--- a/yorick/matte-painting.i
+++ b/yorick/matte-painting.i
@@ -29,18 +29,18 @@ restore, gyoto;
 
 // Choose a metric, set its mass:
 earth_mass=5.9722e24; // mass of the Earth, in kg
-kerrbl=KerrBL(spin=0.999999, mass=1*earth_mass, unit="kg");
+kerrbl=KerrBL(spin=0., mass=1*earth_mass, unit="kg");
 kerrks=KerrKS(mass=kerrbl.mass(), spin=kerrbl.spin());
 minkowski=Metric("Minkowski");
-minkowski, mass=kerrbl.mass(), setparameter="Cartesian";
-metric=kerrks;
+minkowski, mass=kerrbl.mass(), setparameter="Spherical";
+metric=kerrbl;
 
 // Define the screen (=camera):
-res=512; // resolution of the output image
+res=128; // resolution of the output image
 scr=Screen(inclination=pi/2, paln=pi, argument=pi/2,
-           metric=metric, resolution=res);
+           metric=metric, resolution=res, anglekind="Rectilinear");
 scr,distance=2, unit="m"; // camera is at 1m from the BH
-scr, fov=1.0;
+scr, fov=pi/2.;
 
 // Use an empty Astrobj, just set rmax to 10 meters:
 ao=Astrobj("Complex", rmax=10./metric.unitlength());
@@ -50,17 +50,29 @@ sc=Scenery(metric=metric, astrobj=ao, screen=scr, nthreads=8);
 
 ///// Second, choose a painter to paint the sky:
 // To use a jpeg file containing a full-sky, 360°x180° panorama:
-painter=gyoto.painters.mk_panorama(img=jpeg_read("456185645_e56abcc2cd_o.jpg")(,,::-1));
+// download file:
+//   http://farm1.staticflickr.com/192/456185667_adde9d2f8e_o_d.jpg
+//painter=gyoto.painters.mk_panorama(img=jpeg_read("456185645_e56abcc2cd_o.jpg")(,,::-1));
 // To use a p-mode-like pattern:
-//painter=gyoto.painters.mk_p_mode(ntheta=80, nphi=80);
+painter=gyoto.painters.mk_p_mode(ntheta=80, nphi=80);
+
 
 ///// Third, simply call the appropriate function:
-ipct=sc(,,impactcoords=);
-bg=matte_paint(ipct, painter, kind=metric.coordkind(), yaw=1);
+// matte_paint() can be called directly on sc, but to specify only a
+// region or to call matte_paint repeatedly (e.g. on distinct
+// painters, of to fiddle with yaw, pitch & roll), we can precompute
+// ipct:
+ny=res*10/16; // use a 16/10 aspect ratio
+y1=(res-ny)/2+1;
+y2=res+1-y1;
+ipct=sc(,y1:y2,impactcoords=);
+//ipct=fits_read("ipct.fits");
+bg=matte_paint(ipct, painter, kind=metric.coordkind(), yaw=-1.62, roll=-0.32);
 
 ///// Fourth, display the image:
 fma;
 pli, bg;
+limits, square=1;
 
 ///// Fifth, possibly, save it as a JPEG:
 //jpeg_write, "file_out.jpg", bg(,,::-1);

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-astro/packages/gyoto.git



More information about the Debian-astro-commits mailing list