panorama-projection.c 15.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
/* This file is an image processing operation for GEGL
 *
 * GEGL is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 3 of the License, or (at your option) any later version.
 *
 * GEGL 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
14
 * License along with GEGL; if not, see <https://www.gnu.org/licenses/>.
15
 *
16
 * Copyright 2014, 2018 Øyvind Kolås <pippin@gimp.org>
17 18 19 20
 */

#include <math.h>

21

22 23
#ifdef GEGL_PROPERTIES

24 25 26 27
property_double (pan, _("Pan"), 0.0)
  description   (_("Horizontal camera panning"))
  value_range   (-360.0, 360.0)
  ui_meta       ("unit", "degree")
28
  ui_meta       ("direction", "cw")
29 30 31 32 33 34

property_double (tilt, _("Tilt"), 0.0)
  description   (_("Vertical camera panning"))
  value_range   (-180.0, 180.0)
  ui_range      (-180.0, 180.0)
  ui_meta       ("unit", "degree")
35
  ui_meta       ("direction", "cw")
36 37 38 39

property_double (spin, _("Spin"), 0.0)
  description   (_("Spin angle around camera axis"))
  value_range   (-360.0, 360.0)
40
  ui_meta       ("direction", "cw")
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

property_double (zoom, _("Zoom"), 100.0)
  description   (("Zoom level"))
  value_range   (0.01, 1000.0)

property_int    (width, _("Width"), -1)
  description   (_("output/rendering width in pixels, -1 for input width"))
  value_range   (-1, 10000)
  ui_meta       ("role", "output-extent")
  ui_meta       ("axis", "x")

property_int    (height, _("Height"), -1)
  description   (_("output/rendering height in pixels, -1 for input height"))
  value_range   (-1, 10000)
  ui_meta       ("role", "output-extent")
  ui_meta       ("axis", "y")

58
property_boolean(inverse, _("Inverse transform"), FALSE)
59
  description   (_("Do the inverse mapping, useful for touching up zenith, nadir or other parts of panorama."))
60

61
property_enum   (sampler_type, _("Resampling method"),
62 63
                  GeglSamplerType, gegl_sampler_type, GEGL_SAMPLER_NEAREST)
  description   (_("Image resampling method to use, for good results with double resampling when retouching panoramas, use nearest to generate the view and cubic or better for the inverse transform back to panorama."))
64

65 66
#else

67
#define GEGL_OP_FILTER
68
#define GEGL_OP_NAME     panorama_projection
69
#define GEGL_OP_C_SOURCE panorama-projection.c
70 71 72

#include "config.h"
#include <glib/gi18n-lib.h>
73
#include "gegl-op.h"
74

75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
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;
91 92
  float in_width;
  float in_height;
93
  void (*mapfun) (Transform *transform, float x, float  y, float *lon, float *lat);
94
  int   reverse;
95 96
  int   do_spin;
  int   do_zoom;
97 98 99 100 101 102 103 104 105 106 107 108 109
};

/* 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;

110
  y -= 0.5f;
111 112
  x -= transform->xoffset;

113
  if (transform->do_spin)
114 115 116 117 118 119
  {
    float tx = x, ty = y;
    x = tx * transform->cos_spin - ty * transform->sin_spin;
    y = ty * transform->cos_spin + tx * transform->sin_spin;
  }

120 121 122 123 124 125
  if (transform->do_zoom)
  {
    x /= transform->zoom;
    y /= transform->zoom;
  }

126 127 128 129 130 131 132 133 134
  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);

135
  if (longtitude < 0.0f)
136 137 138 139 140 141 142 143 144 145 146
    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)
{
147
  float cos_c, sin_lat, cos_lat, cos_lon_minus_pan;
148

149 150 151 152 153 154
  lat = lat * M_PI - M_PI/2;
  lon = lon * (M_PI * 2);

  sin_lat = sinf (lat);
  cos_lat = cosf (lat);
  cos_lon_minus_pan = cosf (lon - transform->pan);
155

156
  cos_c = (transform->sin_tilt * sin_lat +
157
           transform->cos_tilt * cos_lat *
158 159
            cos_lon_minus_pan);

160 161 162 163 164 165 166
  if (cos_c <= 0.01f)
  {
    *x = -.1f;
    *y = -.1f;
    return;
  }

167 168 169
  *x = ((cos_lat * sin (lon - transform->pan)) / cos_c);
  *y = ((transform->cos_tilt * sin_lat -
         transform->sin_tilt * cos_lat * cos_lon_minus_pan) / cos_c);
170

171 172 173 174 175 176
  if (transform->do_zoom)
  {
    *x *= transform->zoom;
    *y *= transform->zoom;
  }
  if (transform->do_spin)
177 178 179 180 181
  {
    float tx = *x, ty = *y;
    *x = tx * transform->cos_negspin - ty * transform->sin_negspin;
    *y = ty * transform->cos_negspin + tx * transform->sin_negspin;
  }
182 183 184 185

  *x += transform->xoffset;
  *y += 0.5f;

186 187 188 189 190
}

static void prepare_transform (Transform *transform,
                               float pan, float spin, float zoom, float tilt,
                               float width, float height,
191 192
                               float input_width, float input_height,
                               int inverse)
193
{
194
  float xoffset = 0.5f;
195
  transform->reverse = inverse;
196 197
  if (inverse)
    transform->mapfun = gnomonic_ll2xy;
198
  else
199
    transform->mapfun = gnomonic_xy2ll;
200 201 202

  pan  = pan / 360 * M_PI * 2;
  spin = spin / 360 * M_PI * 2;
203
  zoom = zoom / 100.0f;
204 205 206 207 208 209 210
  tilt = tilt / 360 * M_PI * 2;

  while (pan > M_PI)
    pan -= 2 * M_PI;

  if (width <= 0 || height <= 0)
  {
211
    width = input_height;
212
    height = width;
213
    xoffset = ((input_width - height)/height) / 2 + 0.5f;
214 215 216 217 218
  }
  else
  {
    float orig_width = width;
    width = height;
219
    xoffset = ((orig_width - height)/height) / 2 + 0.5f;
220 221
  }

222 223


224 225 226
  transform->do_spin = fabs (spin) > 0.000001 ? 1 : 0;
  transform->do_zoom = fabs (zoom-1.0) > 0.000001 ? 1 : 0;

227 228 229
  transform->pan         = pan;
  transform->tilt        = tilt;
  transform->spin        = spin;
230
  transform->zoom        = zoom;
231 232 233
  transform->xoffset     = xoffset;
  transform->sin_tilt    = sinf (tilt);
  transform->cos_tilt    = cosf (tilt);
234 235 236 237
  transform->sin_spin    = sinf (spin);
  transform->cos_spin    = cosf (spin);
  transform->sin_negspin = sinf (-spin);
  transform->cos_negspin = cosf (-spin);
238 239
  transform->width       = width;
  transform->height      = height;
240 241 242 243 244 245 246 247 248 249 250
  transform->in_width    = input_width;
  transform->in_height   = input_height;

  if (inverse)
  {
    float tmp;
#define swap(a,b) tmp = a; a = b; b= tmp;
    swap(transform->width, transform->in_width);
    swap(transform->height, transform->in_height);
#undef swap
  }
251 252
}

253 254 255
static void
prepare (GeglOperation *operation)
{
256
  const Babl *space = gegl_operation_get_source_space (operation, "input");
257 258 259
  GeglProperties *o = GEGL_PROPERTIES (operation);
  const Babl *format;
  if (o->sampler_type == GEGL_SAMPLER_NEAREST)
260
    format = babl_format_with_space ("RGBA float", space);
261
  else
262
    format = babl_format_with_space ("RaGaBaA float", space);
263 264 265 266 267

  gegl_operation_set_format (operation, "input", format);
  gegl_operation_set_format (operation, "output", format);
}

268 269 270
static GeglRectangle
get_bounding_box (GeglOperation *operation)
{
271
  GeglProperties *o = GEGL_PROPERTIES (operation);
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
  GeglRectangle result = {0,0,0,0};

  if (o->width <= 0 || o->height <= 0)
  {
     GeglRectangle *in_rect;
     in_rect = gegl_operation_source_get_bounding_box (operation, "input");
     if (in_rect)
       {
          result = *in_rect;
       }
     else
     {
       result.width = 320;
       result.height = 200;
     }
  }
  else
  {
    result.width = o->width;
    result.height = o->height;
  }
  return result;
}

296
static void prepare_transform2 (Transform *transform,
297 298
                                GeglOperation *operation,
                                gint level)
299
{
300
  gint factor = 1 << level;
301
  GeglProperties *o = GEGL_PROPERTIES (operation);
302 303
  GeglRectangle in_rect = *gegl_operation_source_get_bounding_box (operation, "input");

304
  prepare_transform (transform,
305
                     o->pan, o->spin, o->zoom, o->tilt,
306
                     o->width / factor, o->height / factor,
307
                     in_rect.width, in_rect.height,
308
                     o->inverse);
309

310 311
}

312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
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");
  return result;
}

static gboolean
process (GeglOperation       *operation,
         GeglBuffer          *input,
         GeglBuffer          *output,
         const GeglRectangle *result,
         gint                 level)
{
328
  GeglProperties *o = GEGL_PROPERTIES (operation);
329
  Transform           transform;
330
  GeglSampler        *sampler;
331
  gint                factor       = 1 << level;
332
  GeglBufferIterator *it;
333 334
  GeglBufferMatrix2   scale_matrix;
  GeglBufferMatrix2  *scale        = NULL;
335
  gint                sampler_type = o->sampler_type;
336
  const Babl         *format_io    = gegl_operation_get_format (operation, "output");
337
  GeglSamplerGetFun   getfun;
338

339 340
  level = 0;
  factor = 1;
341
  prepare_transform2 (&transform, operation, level);
342

343 344
  if (level)
    sampler_type = GEGL_SAMPLER_NEAREST;
345

346
  if (transform.reverse)
347
  {
348
    /* artifacts have been observed with these samplers */
349 350 351
    if (sampler_type == GEGL_SAMPLER_NOHALO ||
        sampler_type == GEGL_SAMPLER_LOHALO)
      sampler_type = GEGL_SAMPLER_CUBIC;
352
  }
353

354
  if (sampler_type != GEGL_SAMPLER_NEAREST &&
355
      !(o->inverse == FALSE && abs(o->tilt < 33)))
356 357
    /* skip the computation of sampler neighborhood scale matrix in cases where
     * we are unlikely to be scaling down */
358
    scale = &scale_matrix;
359 360

  sampler = gegl_buffer_sampler_new_at_level (input, format_io, sampler_type, 0);
361
  getfun = gegl_sampler_get_fun (sampler);
362

363
  {
364 365
    float   ud = ((1.0f/transform.width)*factor);
    float   vd = ((1.0f/transform.height)*factor);
366
    int abyss_mode = transform.reverse ? GEGL_ABYSS_NONE : GEGL_ABYSS_LOOP;
367 368

    it = gegl_buffer_iterator_new (output, result, level, format_io,
369
                                   GEGL_ACCESS_WRITE, GEGL_ABYSS_NONE, 1);
370 371 372

    while (gegl_buffer_iterator_next (it))
      {
373
        GeglRectangle *roi = &it->items[0].roi;
374 375
        gint i;
        gint n_pixels = it->length;
376
        gint x = roi->width; /* initial x                   */
377

378
        float   u0 = (((roi->x*factor * 1.0f)/transform.width));
379 380
        float   u, v;

381
        float *out = it->items[0].data;
382 383

        u = u0;
384
        v = ((roi->y*factor * 1.0/transform.height));
385 386 387 388 389 390

        if (scale)
          {
            for (i=0; i<n_pixels; i++)
              {
                float cx, cy;
391 392 393 394 395 396 397 398
/* we need our own jacobian matrix approximator,
 * since we do not operate on pixel values
 */
#define gegl_sampler_compute_scale2(matrix, x, y)        \
{                                                        \
  float ax, ay, bx, by;                                  \
  gegl_unmap(x + 0.5 * ud, y, ax, ay);                   \
  gegl_unmap(x - 0.5 * ud, y, bx, by);                   \
399 400
  matrix.coeff[0][0] = (ax - bx);   \
  matrix.coeff[1][0] = (ay - by);  \
401 402
  gegl_unmap(x, y + 0.5 * ud, ax, ay);                   \
  gegl_unmap(x, y - 0.5 * ud, bx, by);                   \
403 404
  matrix.coeff[0][1] = (ax - bx);   \
  matrix.coeff[1][1] = (ay - by);  \
405 406 407 408 409
}

#define gegl_unmap(xx,yy,ud,vd) {                                  \
                  float rx, ry;                                    \
                  transform.mapfun (&transform, xx, yy, &rx, &ry); \
410
                  ud = rx;vd = ry;}
411
                gegl_sampler_compute_scale2 (scale_matrix, u, v);
412
                gegl_unmap(u,v, cx, cy);
413
#undef gegl_unmap
414

415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
#if 1
                if (scale_matrix.coeff[0][0] > 0.5f)
                  scale_matrix.coeff[0][0] = (scale_matrix.coeff[0][0]-1.0) * transform.in_width;
                else if (scale_matrix.coeff[0][0] < -0.5f) 
                  scale_matrix.coeff[0][0] = (scale_matrix.coeff[0][0]+1.0) * transform.in_width;
                else
                  scale_matrix.coeff[0][0] *= transform.in_width;

                if (scale_matrix.coeff[0][1] > 0.5f)
                  scale_matrix.coeff[0][1] = (scale_matrix.coeff[0][1]-1.0) * transform.in_width;
                else if (scale_matrix.coeff[0][1] < -0.5f)
                  scale_matrix.coeff[0][1] = (scale_matrix.coeff[0][1]+1.0) * transform.in_width;
                else
                  scale_matrix.coeff[0][1] *= transform.in_width;

                scale_matrix.coeff[1][1] *= transform.in_height;
                scale_matrix.coeff[1][0] *= transform.in_height;
432 433 434 435 436
#endif
                getfun (sampler,
                        cx * transform.in_width + 0.5f,
                        cy * transform.in_height + 0.5f,
                        scale, out, abyss_mode);
437

438 439 440
                out += 4;

                /* update x, y and u,v coordinates */
441
                x--;
442
                u+=ud;
443
                if (x == 0)
444
                  {
445
                    x = roi->width;
446 447 448
                    u = u0;
                    v += vd;
                  }
449
              }
450 451 452 453
          }
        else
          {
            for (i=0; i<n_pixels; i++)
454
              {
455
                float cx, cy;
456
                transform.mapfun (&transform, u, v, &cx, &cy);
457 458 459 460
                getfun (sampler,
                        cx * transform.in_width + 0.5f,
                        cy * transform.in_height + 0.5f,
                        scale, out, abyss_mode);
461 462 463
                out += 4;

                /* update x, y and u,v coordinates */
464
                x--;
465
                u+=ud;
466
                if (x <= 0)
467
                  {
468
                    x = roi->width;
469 470
                    u = u0;
                    v += vd;
471
                  }
472
              }
473 474 475 476
          }
      }
  }

477 478 479 480
  g_object_unref (sampler);
  return TRUE;
}

481 482 483 484 485 486 487 488 489 490
static gchar *composition = "<?xml version='1.0' encoding='UTF-8'?>"
    "<gegl>"
    "<node operation='gegl:panorama-projection' width='200' height='200'/>"
    "<node operation='gegl:load'>"
    "  <params>"
    "    <param name='path'>standard-panorama.png</param>"
    "  </params>"
    "</node>"
    "</gegl>";

491
static void
492
gegl_op_class_init (GeglOpClass *klass)
493 494 495 496 497 498 499 500
{
  GeglOperationClass         *operation_class;
  GeglOperationFilterClass *filter_class;
  operation_class = GEGL_OPERATION_CLASS (klass);
  filter_class  = GEGL_OPERATION_FILTER_CLASS (klass);

  filter_class->process = process;
  operation_class->prepare = prepare;
501
  operation_class->threaded = TRUE;
502
  operation_class->get_bounding_box = get_bounding_box;
503 504 505
  operation_class->get_required_for_output = get_required_for_output;

  gegl_operation_class_set_keys (operation_class,
506
    "name",                  "gegl:panorama-projection",
507
    "title",                 _("Panorama Projection"),
508
    "reference-composition", composition,
509
    "reference-hash",        "3ab9831053ff0a9e32623ecc8a148e67",
510
    "position-dependent",    "true",
511
    "categories" ,           "map",
512
    "description", _("Do panorama viewer rendering mapping or its inverse for an equirectangular input image. (2:1 ratio containing 360x180 degree panorama)."),
513 514 515
    NULL);
}
#endif