[Debian-astro-commits] [gyoto] 82/221: put most of the matte_painting machinery in gyoto.i
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 185edcb99782a5ed0138a3ca4a53966245403a8f
Author: Thibaut Paumard <paumard at users.sourceforge.net>
Date: Tue Nov 4 18:35:04 2014 +0100
put most of the matte_painting machinery in gyoto.i
---
yorick/gyoto.i | 277 +++++++++++++++++++++++++++++++++++++++++++++++
yorick/gyoto_namespace.i | 6 +-
yorick/matte-painting.i | 129 +++-------------------
3 files changed, 300 insertions(+), 112 deletions(-)
diff --git a/yorick/gyoto.i b/yorick/gyoto.i
index 1ba747a..d187152 100644
--- a/yorick/gyoto.i
+++ b/yorick/gyoto.i
@@ -763,6 +763,283 @@ BROKEN
return data;
}
+func gyoto_matte_paint(set, paint, kind=, yaw=, pitch=, roll=)
+/* DOCUMENT bg = gyoto.matte_paint(set, paint[, coordkind=coordkind])
+
+ Vizualize lensing effects on a painted background. Only geometrical
+ distorsion effects are rendered.
+
+ gyoto.matte_paint() will ray-trace the Scenery SET, compute the
+ origin direction of each Photon that do originate from infinite,
+ and pass these directions to PAINT.
+
+ INPUTS:
+ set: a Gyoto Scenery or impact coordinates as computed with
+ sc(,,impactcoords=);
+ paint: a function-like object that will perform the actual
+ painting. PAINT will be called as:
+ bg = paint(theta, phi, mask=mask);
+ where THETA and PHI denote the direction from which each Photon
+ originates, as seen from the center of the coordinate system,
+ and mask is 0 for Photons that do not exist (e.g. photons that
+ would fall into the black-hole.
+
+ OUTPUT:
+ The rendered image.
+
+ */
+{
+ restore, gyoto;
+ coordkind=kind
+ if (is_Scenery(set)) {
+ data=sc(,,impactcoords=)(9:,,);
+ coordkind=sc.metric().coordkind();
+ } else {
+ if (is_void(coordkind))
+ error, "KIND must be provided if SET is not a Scenery";
+ if (dimsof(set)(2)==8) data=set;
+ else data=set(9:,,);
+ }
+
+
+ // Prepare rotation matrix to transform from compact object
+ // coordinates to painter coordinates. Yaw, pitch and roll are the
+ // angles defining the orientation of the painter relative to the
+ // metric coordinate system.
+ // First, R is identity.
+ R=array(0., 3, 3);
+ for (i=1; i<=3; ++i) R(i, i)=1.;
+ // Then, apply each rotation in order.
+ // yaw, about Z axis
+ if (!is_void(yaw)) R = gyoto.rotation(3, yaw);
+ // pitch, about rotated Z axis
+ if (!is_void(pitch)) R = gyoto.rotation(1, pitch)(,+)*R(+,);
+ // roll, about rotated X axis
+ if (!is_void(roll)) R = gyoto.rotation(2, roll)(,+)*R(+,);
+
+ // Transform the photon's 4-velocity into Cartesian coordinates
+ if (coordkind==gyoto.coordkind.spherical) {
+
+ t=data(1,,);
+ r=data(2,,);
+ theta=data(3,,);
+ phi=data(4,,);
+ tdot=data(5,,);
+ rp=rdot=data(6,,);
+ thp=thetadot=data(7,,);
+ php=phidot=data(8,,);
+
+ ind=where(tdot);
+ taup=tdot; taup(ind)=1./tdot(ind);
+ rp =rdot *taup;
+ php=phidot *taup;
+ thp=thetadot*taup;
+
+ sth=sin(theta);
+ cth=cos(theta);
+ sph=sin(phi);
+ cph=cos(phi);
+
+ ur= [sth*cph, sth*sph, cth];
+ uth=[cth*cph, cth*sph, -sth];
+ uph=[ -sph, cph, 0];
+
+ durdph=sth*uph;
+ durdth=uth;
+
+ durdt=php*durdph+thp*durdth;
+
+ v=rp*ur + r*durdt;
+
+ mask0=char(r<1e200);
+
+ } else {
+ tdot=data(5,,);
+ xdot=data(6,,);
+ ydot=data(7,,);
+ zdot=data(8,,);
+ v=[xdot, ydot, zdot]/tdot(,,-);
+ mask0=char(tdot<1e200);
+ }
+
+ // Express this in the painter base:
+ V=R(,+)*v(,,+);
+
+ // transform to spherical coordinates
+ vproj2=V(1,,)^2+V(2,,)^2;
+
+ vr=sqrt(vproj2+V(3,,)^2);
+ vph=atan(V(2,,), V(1,,));
+ vth=atan(sqrt(vproj2), V(3,,));
+
+ theta=pi-vth;
+ phi=vph-pi;
+ if (numberof((ind=where(phi<-pi)))) phi(ind)+=2.*pi;
+
+ bg=paint(theta, phi);
+
+ bg *= mask0(-,,);
+
+ return bg;
+}
+
+if (is_func(use))
+ painters_junk=save(painters_junk,
+ gyoto_painters_p_mode_eval,
+ gyoto_painters_mk_p_mode,
+ gyoto_painters_panorama_eval,
+ gyoto_painters_mk_panorama);
+
+func gyoto_painters_p_mode_eval(theta, phi, mask=) {
+ use, ntheta, nphi;
+ bg=sin(ntheta*theta)*sin(nphi*phi)+1.;
+ if (!is_void(mask)) bg*=mask;
+ return bg;
+}
+
+func gyoto_painters_mk_p_mode(ntheta=, nphi=)
+/* DOCUMENT painter = gyoto.painters.mk_p_mode([ntheta=ntheta, [nphi=nphi])
+
+ Make a painter for gyoto.matte_paint(). In this painter, the sky
+ is mainted in a p-mode spherical harmonic fashion, with m=nphi and
+ l=nphi+ntheta.
+
+ SEE ALSO: gyoto.matte_paint
+*/
+ {
+ if (is_void(ntheta)) ntheta=1;
+ if (is_void(nphi)) nphi=1;
+ return closure(save(ntheta=ntheta, nphi=nphi,
+ painter_eval=gyoto.painters.p_mode_eval),
+ painter_eval);
+}
+
+func gyoto_rotation(axis, angle)
+/* DOCUMENT R=gyoto.rotation(axis, angle)
+
+ Yield rotation matrix about axis AXIS (Ox=1, Oy=2, Oz=3) by
+ amount ANGLE.
+
+ If vector V holds the coordinates of a point in original the
+ coordinate base,
+ R(,+)*v(+)
+ holds coordinates of same point in rotated base while
+ R(+,)*v(+)
+ holds coordinates of rotated vector in original base.
+
+ In both cases, if S and T are two rotations, R=T(,+)*S(+,) if the
+ rotation obtained by first applying S, then T.
+
+ SEE ALSO:
+ */
+{
+
+ ca=cos(angle);
+ sa=sin(angle);
+
+ R2=[[ca, -sa], [sa, ca]];
+
+ R3=array(0., 3, 3);
+ for (i=1; i<=3; ++i) R3(i, i)=1.;
+
+ plane=(axis+[1, 2])%3;
+
+ if (numberof((ind=where(plane==0)))) plane(ind)=3;
+
+ R3(plane, plane)=R2;
+
+ return R3;
+}
+
+func gyoto_painters_panorama_eval(theta, phi, mask=) {
+ use, img, phi_fov, theta_fov;
+ dmatte=dimsof(img);
+ nphi=dmatte(3);
+ ntheta=dmatte(4);
+
+ if (is_void(phi_fov)) phi_fov=2.*pi;
+ if (is_void(theta_fov)) theta_fov=pi;
+
+
+ // Thats in pixel per radian. Minus signs because theta increases
+ // from top to bottom and phi from right to left, as seen from
+ // inside the sphere.
+ phi_scale=-nphi/phi_fov;
+ theta_scale=-ntheta/theta_fov;
+
+ // Below, +0.5 so that for a 4x4 grid, the center is at 2.5, and
+ // another +0.5 so that the long() conversion yields the closest
+ // approximation rahter than the bottom truncation.
+ i0=nphi*0.5 +1.;
+ j0=ntheta*0.5 +1.;
+
+ 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))
+ error, "something fishy in pixel computation";
+
+ dd=dimsof(phi);
+
+ bg=array(char, 3, dd(2), dd(3));
+ for(k=1; k<=3; ++k) {
+ for (i=1; i<=dd(2); ++i) {
+ for (j=1; j<=dd(3); ++j) {
+ bg(k,i,j)=img(k,x(i, j),y(i, j));
+ }
+ }
+ }
+
+ return bg;
+}
+
+
+func gyoto_painters_mk_panorama(img=, phi_fov=, theta_fov=)
+/* DOCUMENT painter = gyoto.painters.mk_picture(img=img)
+
+ Make a painter for gyoto.matte_paint(). In this painter, the sky
+ is painted using an image. The image should be a 360°x180°
+ panorama rendered in equirectagular projection. If the image does
+ not cover the full celestial sphere, use phi_fov and theta_fov to
+ specify the field of view in each direction. img should be an
+ array(char, 3, nphi, ntheta).
+
+ SEE ALSO: gyoto.matte_paint, .painters.mk_p_mode
+*/
+{
+ d=dimsof(img);
+ nx=d(3);
+ ny=d(4);
+ if (is_void(phi_fov)) phi_fov=2.*pi;
+ if (is_void(theta_fov)) theta_fov=pi;
+ return closure(save(img, phi_fov, theta_fov,
+ painter_eval=gyoto.painters.panorama_eval),
+ painter_eval);
+}
+
+local gyoto_painters;
+/* DOCUMENT namespace gyoto.painters
+
+ Contains painters to be used with gyoto.matte_paint.
+
+ EXAMPLE:
+ sc=gyoto.Scenery("file.xml");
+ painter=gyoto.painters.mk_p_mode(ntheta=80, nphi=80);
+ bg=gyoto.matte_paint(sc, painter);
+
+ SEE ALSO: gyoto.matte_paint, .painters.mk_p_mode
+ */
+
+if (is_func(use)) {
+ gyoto_painters=save(p_mode_eval=gyoto_painters_p_mode_eval,
+ mk_p_mode=gyoto_painters_mk_p_mode,
+ picture_eval=gyoto_painters_picture_eval,
+ mk_picture=gyoto_painters_mk_picture,
+ panorama_eval=gyoto_painters_panorama_eval,
+ mk_panorama=gyoto_painters_mk_panorama);
+ restore, painters_junk;
+ }
+
// PHOTON CLASS
extern gyoto_Photon;
/* DOCUMENT photon = gyoto.Photon([filename], [members=values])
diff --git a/yorick/gyoto_namespace.i b/yorick/gyoto_namespace.i
index 4e1dbe0..6c016a6 100644
--- a/yorick/gyoto_namespace.i
+++ b/yorick/gyoto_namespace.i
@@ -13,6 +13,7 @@ gyoto=save(
Scenery=gyoto_Scenery,
Scenery_rayTrace=gyoto_Scenery_rayTrace,
+ matte_paint=gyoto_matte_paint,
Photon=gyoto_Photon,
Metric=gyoto_Metric,
Astrobj=gyoto_Astrobj,
@@ -46,7 +47,10 @@ gyoto=save(
SUN_RADIUS=GYOTO_SUN_RADIUS,
KPC=GYOTO_KPC,
- coordkind=save(unspecified=0, cartesian=1, spherical=2)
+ coordkind=save(unspecified=0, cartesian=1, spherical=2),
+
+ painters=gyoto_painters,
+ rotation=gyoto_rotation
);
diff --git a/yorick/matte-painting.i b/yorick/matte-painting.i
index 88c6538..3a0b2cc 100644
--- a/yorick/matte-painting.i
+++ b/yorick/matte-painting.i
@@ -25,33 +25,22 @@
#include "jpeg.i"
restore, gyoto;
-matte=[];
-// To use a jpeg file as matte:
-// - uncomment the following line;
-// - set the path to the file you want;
-// - set x0, y0 and fov0 below.
-//matte=jpeg_read("file_in.jpg")(,,::-1);
-
-// Pixel coordinates of the center of the zone of interest in the
-// JPEG image:
-x0=2000;
-y0=1500;
-// Size in pixel of the zone of the JPEG file that should map the
-// screen field-of-view:
-fov0=3000;
+///// First, define an appropriate Scenery:
// Choose a metric, set its mass:
earth_mass=5.9722e24; // mass of the Earth, in kg
-kerrbl=KerrBL(spin=0., mass=0.1*earth_mass, unit="kg");
+kerrbl=KerrBL(spin=0.999999, mass=1*earth_mass, unit="kg");
kerrks=KerrKS(mass=kerrbl.mass(), spin=kerrbl.spin());
minkowski=Metric("Minkowski");
minkowski, mass=kerrbl.mass(), setparameter="Cartesian";
-metric=kerrbl;
+metric=kerrks;
// Define the screen (=camera):
-res=128; // resolution of the output image
-scr=Screen(inclination=pi/2, paln=pi, metric=metric, resolution=res);
-scr,distance=1, unit="m"; // camera is at 1m from the BH
+res=512; // resolution of the output image
+scr=Screen(inclination=pi/2, paln=pi, argument=pi/2,
+ metric=metric, resolution=res);
+scr,distance=2, unit="m"; // camera is at 1m from the BH
+scr, fov=1.0;
// Use an empty Astrobj, just set rmax to 10 meters:
ao=Astrobj("Complex", rmax=10./metric.unitlength());
@@ -59,101 +48,19 @@ ao=Astrobj("Complex", rmax=10./metric.unitlength());
// Build a Scenery with this Metric, Screen and Astrobj:
sc=Scenery(metric=metric, astrobj=ao, screen=scr, nthreads=8);
-// Compute the last coordinates of the "escaping" photons
-// (Actually coming from infinity)
-data=sc(,,impactcoords=)(9:,,);
-
-// Transform the photon's 4-velocity into equatorial coordinates
-if (sc.metric().coordkind()==gyoto.coordkind.spherical) {
-
- t=data(1,,);
- r=data(2,,);
- theta=data(3,,);
- phi=data(4,,);
- tdot=data(5,,);
- rp=rdot=data(6,,);
- thp=thetadot=data(7,,);
- php=phidot=data(8,,);
-
- ind=where(tdot);
- taup=tdot; taup(ind)=1./tdot(ind);
- rp =rdot *taup;
- php=phidot *taup;
- thp=thetadot*taup;
-
- sth=sin(theta);
- cth=cos(theta);
- sph=sin(phi);
- cph=cos(phi);
-
- ur= [sth*cph, sth*sph, cth];
- uth=[cth*cph, cth*sph, -sth];
- uph=[ -sph, cph, 0];
-
- durdph=sth*uph;
- durdth=uth;
-
- durdt=php*durdph+thp*durdth;
-
- v=rp*ur + r*durdt;
-
- mask0=r<1e200;
-
- } else {
- tdot=data(5,,);
- xdot=data(6,,);
- ydot=data(7,,);
- zdot=data(8,,);
- v=[xdot, ydot, zdot]/tdot(,,-);
- mask0=tdot<1e200;
- }
-
-vproj2=v(,,1)^2+v(,,2)^2
+///// 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));
+// To use a p-mode-like pattern:
+//painter=gyoto.painters.mk_p_mode(ntheta=80, nphi=80);
-vr=sqrt(vproj2+v(,,3)^2);
-vph=atan(v(,,2), v(,,1));
-vth=atan(sqrt(vproj2), v(,,3));
+///// Third, simply call the appropriate function:
+ipct=sc(,,impactcoords=);
+bg=matte_paint(ipct, painter, kind=metric.coordkind(), yaw=1);
-ind=where(r>1e308);
-
-// A simple analytical matte:
-bg=sin(80*vph)*sin(80*vth);
-if (numberof(ind)) bg(ind)=bg(min, min);
-
-// If a JPEG was provided:
-if (!is_void(matte)) {
-
- dmatte=dimsof(matte);
-
- phi0=(vph(0,0)+vph(1,1))*0.5;
- theta0=(vth(0,0)+vth(1,1))*0.5;
- Dphi=(vph(0,0)-vph(1,1));
- Dtheta=(vth(0,0)-vth(1,1));
-
- x=long((vph-phi0)/Dphi*fov0)+x0;
- y=long((vth-theta0)/Dtheta*fov0)+y0;
-
- if (anyof((mask=(x<1)))) x(where(mask))=1;
- if (anyof((mask=(y<1)))) y(where(mask))=1;
- if (anyof((mask=(x>dmatte(3))))) x(where(mask))=dmatte(3);
- if (anyof((mask=(y>dmatte(4))))) y(where(mask))=dmatte(4);
-
- bg=array(char, 3, res, res);
- for(k=1; k<=3; ++k) {
- for (i=1; i<=res; ++i) {
- for (j=1; j<=res; ++j) {
- bg(k,i,j)=matte(k,x(i, j),y(i, j));
- }
- }
- bg(k,,) *= mask0;
- }
-
- }
-
-// Display the image:
+///// Fourth, display the image:
fma;
pli, bg;
-
-// Possibly, save it as a JPEG:
+///// 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