[gegl] panorama-projection: add inverse transforms to core

[Øyvind Kolås <ok src gnome org>]

commit 408ba8efc941edcdff2dd8524303bdef0db59a4f
Author: Øyvind Kolås <pippin gimp org>
Date:   Sat May 3 23:58:15 2014 +0200

    panorama-projection: add inverse transforms to core
    
    The core transformation machinery is now abstracted to be work as an isolated
    unit of code.

 operations/common/panorama-projection.c |  399 +++++++++++++++++++++---------
 1 files changed, 279 insertions(+), 120 deletions(-)
---
diff --git a/operations/common/panorama-projection.c b/operations/common/panorama-projection.c
index 8045724..37ad374 100644
--- a/operations/common/panorama-projection.c
+++ b/operations/common/panorama-projection.c
@@ -14,7 +14,6 @@
  * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
  *
  * Copyright 2014 Øyvind Kolås <pippin gimp org>
- *
  */
 
 #include <math.h>
@@ -35,7 +34,7 @@ gegl_chant_int (height, _("height"), -1, 10000, -1, _("input buffer height, -1 f
 gegl_chant_boolean (little_planet, _("little planet"), FALSE, _("use the pan/tilt location as center for a 
stereographic/little planet projection."))
 
 gegl_chant_enum (sampler_type, _("Sampler"), GeglSamplerType, gegl_sampler_type,
-                 GEGL_SAMPLER_CUBIC, _("Image resampling method to use"))
+                 GEGL_SAMPLER_NEAREST, _("Image resampling method to use"))
 
 #else
 
@@ -48,6 +47,224 @@ gegl_chant_enum (sampler_type, _("Sampler"), GeglSamplerType, gegl_sampler_type,
 
 //#include "speedmath.inc"
 
+typedef struct _Transform Transform;
+struct _Transform
+{
+  float pan;
+  float tilt;
+  float sin_tilt;
+  float cos_tilt;
+  float sin_spin;
+  float cos_spin;
+  float sin_negspin;
+  float cos_negspin;
+  float zoom;
+  float spin;
+  float xoffset;
+  float width;
+  float height;
+  void (*xy2ll) (Transform *transform, float x, float  y, float *lon, float *lat);
+  void (*ll2xy) (Transform *transform, float lon, float  lat, float *x, float *y);
+  int   do_spin;
+  int   do_zoom;
+};
+
+/* formulas from:
+ * http://mathworld.wolfram.com/GnomonicProjection.html
+ */
+static void inline
+gnomonic_xy2ll (Transform *transform, float x, float y,
+                float *lon, float *lat)
+{
+  float p, c;
+  float longtitude, latitude;
+  float sin_c, cos_c;
+
+  if (transform->do_spin)
+  {
+    float tx = x, ty = y;
+    x = tx * transform->cos_spin - ty * transform->sin_spin;
+    y = ty * transform->cos_spin + tx * transform->sin_spin;
+  }
+
+  if (transform->do_zoom)
+  {
+    x /= transform->zoom;
+    y /= transform->zoom;
+  }
+
+  p = sqrtf (x*x+y*y);
+  c = atan2f (p, 1);
+
+  sin_c = sinf(c);
+  cos_c = cosf(c);
+
+  latitude = asinf (cos_c * transform->sin_tilt + ( y * sin_c * transform->cos_tilt) / p);
+  longtitude = transform->pan + atan2f (x * sin_c, p * transform->cos_tilt * cos_c - y * transform->sin_tilt 
* sin_c);
+
+  if (longtitude < 0)
+    longtitude += M_PI * 2;
+
+  *lon = (longtitude / (M_PI * 2));
+  *lat = ((latitude + M_PI/2) / M_PI);
+}
+
+static void inline
+gnomonic_ll2xy (Transform *transform,
+                float lon, float lat,
+                float *x, float *y)
+{
+  float cos_c;
+  float sin_lon = sinf (lon);
+  float cos_lon = cosf (lon);
+  float cos_lat_minus_pan = cosf (lat - transform->pan);
+
+  cos_c = (transform->sin_tilt * sin_lon +
+           transform->cos_tilt * cos (lon) *
+            cos_lat_minus_pan);
+
+  *x = ((cos_lon * sin (lat - transform->pan)) / cos_c);
+  *y = ((transform->cos_tilt * sin_lon -
+         transform->sin_tilt * cos_lon * cos_lat_minus_pan) / cos_c);
+
+  if (transform->do_zoom)
+  {
+    *x *= transform->zoom;
+    *y *= transform->zoom;
+  }
+  if (transform->do_spin)
+  {
+    float tx = *x, ty = *y;
+    *x = tx * transform->cos_negspin - ty * transform->sin_negspin;
+    *y = ty * transform->cos_negspin + tx * transform->sin_negspin;
+  }
+}
+
+static void inline
+stereographic_ll2xy (Transform *transform,
+                     float lon, float lat,
+                     float *x,  float *y)
+{
+  float k;
+  float sin_lon = sinf (lon);
+  float cos_lon = cosf (lon);
+  float cos_lat_minus_pan = cosf (lat - transform->pan);
+
+  k = 2/(1 + transform->sin_tilt * sin_lon +
+            transform->cos_tilt * cos (lon) *
+            cos_lat_minus_pan);
+
+  *x = k * ((cos_lon * sin (lat - transform->pan)));
+  *y = k * ((transform->cos_tilt * sin_lon -
+        transform->sin_tilt * cos_lon * cos_lat_minus_pan));
+
+  if (transform->do_zoom)
+  {
+    *x *= transform->zoom;
+    *y *= transform->zoom;
+  }
+
+  if (transform->do_spin)
+  {
+    float tx = *x, ty = *y;
+    *x = tx * transform->cos_negspin - ty * transform->sin_negspin;
+    *y = ty * transform->cos_negspin + tx * transform->sin_negspin;
+  }
+}
+
+static void inline
+stereographic_xy2ll (Transform *transform,
+                     float x, float y,
+                     float *lon, float *lat)
+{
+  float p, c;
+  float longtitude, latitude;
+  float sin_c, cos_c;
+
+  if (transform->do_spin)
+  {
+    float tx = x, ty = y;
+    x = tx * transform->cos_spin - ty * transform->sin_spin;
+    y = ty * transform->cos_spin + tx * transform->sin_spin;
+  }
+
+  if (transform->do_zoom)
+  {
+    x /= transform->zoom;
+    y /= transform->zoom;
+  }
+
+  p = sqrtf (x*x+y*y);
+  c = 2 * atan2f (p / 2, 1);
+
+  sin_c = sinf (c);
+  cos_c = cosf (c);
+
+  latitude = asinf (cos_c * transform->sin_tilt + ( y * sin_c * transform->cos_tilt) / p);
+  longtitude = transform->pan + atan2f ( x * sin_c, p * transform->cos_tilt * cos_c - y * 
transform->sin_tilt * sin_c);
+
+  if (longtitude < 0)
+    longtitude += M_PI * 2;
+
+  *lon = (longtitude / (M_PI * 2));
+  *lat = ((latitude + M_PI/2) / M_PI);
+}
+
+static void prepare_transform (Transform *transform,
+                               float pan, float spin, float zoom, float tilt,
+                               int little_planet,
+                               float width, float height,
+                               float input_width, float input_height)
+{
+  float xoffset = 0.5;
+  transform->xy2ll = gnomonic_xy2ll;
+  transform->ll2xy = gnomonic_ll2xy;
+
+  pan  = pan / 360 * M_PI * 2;
+  spin = spin / 360 * M_PI * 2;
+  zoom = little_planet?zoom / 1000.0:zoom / 100.0;
+  tilt = tilt / 360 * M_PI * 2;
+
+  while (pan > M_PI)
+    pan -= 2 * M_PI;
+
+  if (width <= 0 || height <= 0)
+  {
+    width = input_height; 
+    height = width;
+    xoffset = ((input_width - height)/height) / 2 + 0.5;
+  }
+  else
+  {
+    float orig_width = width;
+    width = height;
+    xoffset = ((orig_width - height)/height) / 2 + 0.5;
+  }
+
+  if (little_planet)
+  {
+    transform->xy2ll = stereographic_xy2ll;
+    transform->ll2xy = stereographic_ll2xy;
+  }
+
+  transform->do_spin = fabs (spin) > 0.000001 ? 1 : 0;
+  transform->do_zoom = fabs (zoom-1.0) > 0.000001 ? 1 : 0;
+
+  transform->pan         = pan;
+  transform->tilt        = tilt;
+  transform->spin        = spin;
+  transform->zoom        = zoom;
+  transform->xoffset     = xoffset;
+  transform->sin_tilt    = sinf (tilt);
+  transform->cos_tilt    = cosf (tilt);
+  transform->sin_spin    = sinf (spin);
+  transform->cos_spin    = cosf (spin);
+  transform->sin_negspin = sinf (-spin);
+  transform->cos_negspin = cosf (-spin);
+  transform->width       = width;
+  transform->height      = height;
+}
+
 static void
 prepare (GeglOperation *operation)
 {
@@ -82,76 +299,31 @@ get_bounding_box (GeglOperation *operation)
     result.width = o->width;
     result.height = o->height;
   }
-
   return result;
 }
 
+
+static void prepare_transform2 (Transform *transform,
+                                GeglOperation *operation)
+{
+  GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+  GeglRectangle in_rect = *gegl_operation_source_get_bounding_box (operation, "input");
+
+  prepare_transform (transform, 
+                     o->pan, o->spin, o->zoom, o->tilt,
+                     o->little_planet, o->width, o->height,
+                     in_rect.width, in_rect.height);
+}
+
 static GeglRectangle
 get_required_for_output (GeglOperation       *operation,
                          const gchar         *input_pad,
                          const GeglRectangle *region)
 {
   GeglRectangle result = *gegl_operation_source_get_bounding_box (operation, "input");
-  /* XXX: do inverse computations; and use their bounding box */
-
   return result;
 }
 
-/* formulas from:
- * http://mathworld.wolfram.com/GnomonicProjection.html
- */
-
-static void inline
-gnomonic_calc_long_lat (float x,   float  y,
-                        float pan, float spin,
-                        float sin_tilt, float cos_tilt,
-                        float *in_long, float *in_lat)
-{
-  float p, c;
-  float longtitude, latitude;
-  float sin_c, cos_c;
-
-  p = sqrtf (x*x+y*y);
-  c = atan2f (p, 1);
-
-  sin_c = sinf(c);
-  cos_c = cosf(c);
-
-  latitude = asinf (cos_c * sin_tilt + ( y * sin_c * cos_tilt) / p);
-  longtitude = pan + atan2f (x * sin_c, p * cos_tilt * cos_c - y * sin_tilt * sin_c);
-
-  if (longtitude < 0)
-    longtitude += M_PI * 2;
-
-  *in_long = (longtitude / (M_PI * 2));
-  *in_lat = ((latitude + M_PI/2) / M_PI);
-}
-
-static void inline
-stereographic_calc_long_lat (float x,   float y,
-                             float pan, float spin,
-                             float sin_tilt, float  cos_tilt,
-                             float *in_long, float *in_lat)
-{
-  float p, c;
-  float longtitude, latitude;
-  float sin_c, cos_c;
-
-  p = sqrtf (x*x+y*y);
-  c = 2 * atan2f (p / 2, 1);
-
-  sin_c = sinf (c);
-  cos_c = cosf (c);
-
-  latitude = asinf (cos_c * sin_tilt + ( y * sin_c * cos_tilt) / p);
-  longtitude = pan + atan2f ( x * sin_c, p * cos_tilt * cos_c - y * sin_tilt * sin_c);
-
-  if (longtitude < 0)
-    longtitude += M_PI * 2;
-
-  *in_long = (longtitude / (M_PI * 2));
-  *in_lat = ((latitude + M_PI/2) / M_PI);
-}
 
 static gboolean
 process (GeglOperation       *operation,
@@ -161,67 +333,27 @@ process (GeglOperation       *operation,
          gint                 level)
 {
   GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
-  float pan  = o->pan / 360 * M_PI * 2;
-  float spin = o->spin / 360 * M_PI * 2;
-  float zoom = o->zoom / 100.0;
-  float tilt = o->tilt / 360 * M_PI * 2;
-  float width;
-  float height;
+  Transform           transform;
   const Babl         *format_io;
   GeglSampler        *sampler;
   GeglBufferIterator *it;
-  float sin_tilt, cos_tilt;
-  float cos_spin, sin_spin;
   gint  index_in, index_out;
-  float xoffset = 0.5;
   GeglRectangle in_rect = *gegl_operation_source_get_bounding_box (operation, "input");
   GeglMatrix2  scale_matrix;
   GeglMatrix2 *scale = NULL;
-  void (*calc_long_lat) (float x, float  y,
-                         float pan, float spin,
-                         float sin_tilt, float cos_tilt,
-                         float *in_long, float *in_lat) =
-    gnomonic_calc_long_lat;
+
+  prepare_transform2 (&transform, operation);
 
   format_io = babl_format ("RaGaBaA float");
   sampler = gegl_buffer_sampler_new (input, format_io, o->sampler_type);
 
-  while (pan > M_PI)
-    pan -= M_PI;
-
-  sin_tilt = sinf (tilt);
-  cos_tilt = cosf (tilt);
-  sin_spin = sinf (spin);
-  cos_spin = cosf (spin);
-
-  width = o->width;
-  height = o->height;
-
-  if (width <= 0 || height <= 0)
-  {
-    width = in_rect.height;
-    height = width;
-    xoffset = ((in_rect.width - height)/height) / 2 + 0.5;
-  }
-  else
-  {
-    width = height;
-    xoffset = ((o->width - height)/height) / 2 + 0.5;
-  }
-
-  if (o->little_planet)
-  {
-    calc_long_lat = stereographic_calc_long_lat;
-    zoom /= 10;
-  }
-
   if (o->sampler_type == GEGL_SAMPLER_NOHALO ||
       o->sampler_type == GEGL_SAMPLER_LOHALO)
     scale = &scale_matrix;
 
     {
-      float   ud = ((1.0/width))/zoom;
-      float   vd = ((1.0/height))/zoom;
+      float   ud = ((1.0/transform.width));
+      float   vd = ((1.0/transform.height));
       it = gegl_buffer_iterator_new (output, result, level, format_io, GEGL_BUFFER_WRITE, GEGL_ABYSS_NONE);
       index_out = 0;
 
@@ -232,14 +364,14 @@ process (GeglOperation       *operation,
           gint    x = it->roi->x; /* initial x                   */
           gint    y = it->roi->y; /*           and y coordinates */
 
-          float   u0 = ((x/width) - xoffset)/zoom;
+          float   u0 = ((x/transform.width) - transform.xoffset);
           float   u, v;
 
           float *in = it->data[index_in];
           float *out = it->data[index_out];
 
           u = u0;
-          v = ((y/height) - 0.5) / zoom; 
+          v = ((y/transform.height) - 0.5); 
 
           if (scale)
             {
@@ -248,12 +380,7 @@ process (GeglOperation       *operation,
                   float cx, cy;
 #define gegl_unmap(xx,yy,ud,vd) { \
                   float rx, ry;\
-                  calc_long_lat (\
-                     xx * cos_spin - yy * sin_spin,\
-                     yy * cos_spin + xx * sin_spin,\
-                     pan, spin,\
-                     sin_tilt, cos_tilt,\
-                     &rx, &ry);\
+                  transform.xy2ll (&transform, xx, yy, &rx, &ry);\
                   ud = rx;vd = ry;}
                   gegl_sampler_compute_scale (scale_matrix, u, v);
                   gegl_unmap(u,v, cx, cy);
@@ -283,14 +410,10 @@ process (GeglOperation       *operation,
                   {
                     float cx, cy;
 
-                    calc_long_lat (
-                        u * cos_spin - v * sin_spin,
-                        v * cos_spin + u * sin_spin,
-                        pan, spin,
-                        sin_tilt, cos_tilt,
-                        &cx, &cy);
+                    transform.xy2ll (&transform, u, v, &cx, &cy);
 
-                    gegl_sampler_get (sampler, cx * in_rect.width, cy * in_rect.height,
+                    gegl_sampler_get (sampler,
+                                      cx * in_rect.width, cy * in_rect.height,
                                       scale, out, GEGL_ABYSS_LOOP);
                     in  += 4;
                     out += 4;
@@ -311,6 +434,42 @@ process (GeglOperation       *operation,
     }
   g_object_unref (sampler);
 
+#if 0
+  {
+    float t;
+    float lat0  = 0;
+    float lon0 = 0;
+    float lat1  = 0.5;
+    float lon1 = 0.5;
+    int i = 0;
+    guchar pixel[4] = {255,0,0,255};
+
+    for (t = 0; t < 1.0; t+=0.01, i++)
+    {
+      float lat = lat0 * (1.0 - t) + lat1 * t;
+      float lon = lon0 * (1.0 - t) + lon1 * t;
+      float x, y;
+      float xx, yy;
+      GeglRectangle prect = {0,0,1,1};
+
+      ll2xy (&transform, lon, lat, &x, &y);
+
+      x += xoffset;
+      y += 0.5;
+
+      x *= width;
+      y *= height;
+
+      prect.x = floor (x);
+      prect.y = floor (y);
+      prect.width = 1;
+      prect.height = 1;
+
+      gegl_buffer_set (output, &prect, 0, babl_format ("R'G'B' u8"), pixel, 8);
+    }
+  }
+#endif
+
   return TRUE;
 }