src/panolib.c

   1 /* Panorama_Tools       -       Generate, Edit and Convert Panoramic Images
   2    Copyright (C) 1998,1999 - Helmut Dersch  der@fh-furtwangen.de
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2, or (at your option)
   7    any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  */
  17
  18 /*------------------------------------------------------------*/
  19 #include <stdlib.h>
  20 #include <stdint.h>
  21 #include <assert.h>
  22
  23 #include "filter.h"
  24 #include "panolib.h"
  25 #include "utils_math.h"
  26 #include "utils_image.h"
  27
  28 #define DEG_TO_RAD(x) ((x) * 2.0 * M_PI / 360.0 )
  29
  30 // Lookup Tables for Trig-functions and interpolator
  31
  32 #define NATAN 2048
  33 #define NSQRT 2048
  34
  35 static void matrix_matrix_mult( double m1[3][3],double m2[3][3],double result[3][3])
  36 {
  37         int i,k;
  38
  39         for(i=0;i<3;i++)
  40                 for(k=0; k<3; k++)
  41                         result[i][k] = m1[i][0] * m2[0][k] + m1[i][1] * m2[1][k] + m1[i][2] * m2[2][k];
  42 } /* void matrix_matrix_mult */
  43
  44 /* Set matrix elements based on Euler angles a, b, c */
  45 static void set_transformation_matrix (double m[3][3],
  46                 double a, double b, double c)
  47 {
  48         double mx[3][3], my[3][3], mz[3][3], dummy[3][3];
  49
  50
  51         // Calculate Matrices;
  52
  53         mx[0][0] = 1.0 ;                                mx[0][1] = 0.0 ;                                mx[0][2] = 0.0;
  54         mx[1][0] = 0.0 ;                                mx[1][1] = cos(a) ;                     mx[1][2] = sin(a);
  55         mx[2][0] = 0.0 ;                                mx[2][1] =-mx[1][2] ;                   mx[2][2] = mx[1][1];
  56
  57         my[0][0] = cos(b);                              my[0][1] = 0.0 ;                                my[0][2] =-sin(b);
  58         my[1][0] = 0.0 ;                                my[1][1] = 1.0 ;                                my[1][2] = 0.0;
  59         my[2][0] = -my[0][2];                   my[2][1] = 0.0 ;                                my[2][2] = my[0][0];
  60
  61         mz[0][0] = cos(c) ;                     mz[0][1] = sin(c) ;                     mz[0][2] = 0.0;
  62         mz[1][0] =-mz[0][1] ;                   mz[1][1] = mz[0][0] ;                   mz[1][2] = 0.0;
  63         mz[2][0] = 0.0 ;                                mz[2][1] = 0.0 ;                                mz[2][2] = 1.0;
  64
  65         /* Calculate `m = mz * mx * my' */
  66
  67         matrix_matrix_mult( mz, mx,     dummy);
  68         matrix_matrix_mult( dummy, my, m);
  69 } /* void SetMatrix */
  70
  71 #if 0
  72 static double lanczos (double x)
  73 {
  74         static double table[257];
  75
  76         assert ((x >= 0.0) && (x <= 1.0));
  77
  78         if (table[0] != 1.0)
  79         {
  80                 int i;
  81
  82                 table[0] = 1.0;
  83                 for (i = 1; i <= 256; i++)
  84                 {
  85                         double t = ((double) i) / 256.0;
  86                         table[i] = 2 * sin (M_PI * t) * sin (M_PI_2 * t) / (M_PI * M_PI * t * t);
  87                 }
  88         }
  89
  90         return (table[(int) (0.5 + 256.0 * x)]);
  91 } /* double lanczos */
  92 #endif
  93
  94 static int copy_pixel (ui_image_t *dest, const ui_image_t *src,
  95                 int x_dest, int y_dest,
  96                 double x_src_fp, double y_src_fp,
  97                 interpolator_t interp)
  98 {
  99         uint32_t pixel_dest = (y_dest * dest->width) + x_dest;
 100
 101         assert (x_src_fp >= 0.0);
 102
 103         if (y_src_fp < 0.0)
 104                 y_src_fp = 0.0;
 105
 106         if (interp == BILINEAR)
 107         {
 108                 int x_src_left;
 109                 int x_src_right;
 110                 int y_src_top;
 111                 int y_src_bottom;
 112
 113                 double x_right_frac;
 114                 double y_bottom_frac;
 115
 116                 int i;
 117
 118                 x_src_left = (int) x_src_fp;
 119                 x_src_right = (x_src_left + 1) % src->width;
 120
 121                 y_src_top  = (int) y_src_fp;
 122                 y_src_bottom = y_src_top + 1;
 123                 if (y_src_bottom >= src->height)
 124                         y_src_bottom = src->height - 1;
 125
 126                 x_right_frac = x_src_fp - x_src_left;
 127                 y_bottom_frac = y_src_fp - y_src_top;
 128
 129                 /*
 130                  * The polynomial function
 131                  *   2x^3 - 3x^2 + 1
 132                  * is similar to the Lanczos function, but much easier to compute. It's
 133                  * applied here to improve sharpness.
 134                  */
 135                 {
 136                         double t = x_src_fp - x_src_left;
 137                         x_right_frac = -2.0 * t * t * t + 3.0 * t * t;
 138
 139                         t = y_src_fp - y_src_top;
 140                         y_bottom_frac = -2.0 * t * t * t + 3.0 * t * t;
 141                 }
 142
 143                 assert ((x_right_frac >= 0.0) && (x_right_frac <= 1.0));
 144                 assert ((y_bottom_frac >= 0.0) && (y_bottom_frac <= 1.0));
 145
 146                 for (i = 0; i < 3; i++)
 147                 {
 148                         uint8_t values[2][2];
 149                         double value_left;
 150                         double value_right;
 151                         double value_final;
 152
 153                         values[0][0] = src->data[i][(y_src_top * src->width) + x_src_left];
 154                         values[0][1] = src->data[i][(y_src_top * src->width) + x_src_right];
 155                         values[1][0] = src->data[i][(y_src_bottom * src->width) + x_src_left];
 156                         values[1][1] = src->data[i][(y_src_bottom * src->width) + x_src_right];
 157
 158                         value_left = (1.0 - y_bottom_frac) * values[0][0]
 159                                 + y_bottom_frac * values[1][0];
 160                         value_right = (1.0 - y_bottom_frac) * values[0][1]
 161                                 + y_bottom_frac * values[1][1];
 162
 163                         value_final = (1.0 - x_right_frac) * value_left + x_right_frac * value_right;
 164
 165                         assert ((value_final >= 0.0) && (value_final <= 255.0));
 166
 167                         dest->data[i][pixel_dest] = (uint8_t) (value_final + 0.5);
 168                 }
 169         }
 170         else /* if (interp == NNEIGHBOUR) */
 171         {
 172                 int x_src = (int) (x_src_fp + 0.5) % src->width;
 173                 int y_src = (int) (y_src_fp + 0.5) % src->height;
 174
 175                 if ((x_src < 0) || (x_src >= src->width)
 176                                 || (y_src < 0) || (y_src >= src->height))
 177                 {
 178                         dest->data[0][pixel_dest] = 0;
 179                         dest->data[1][pixel_dest] = 0;
 180                         dest->data[2][pixel_dest] = 0;
 181                 }
 182                 else
 183                 {
 184                         uint32_t pixel_src = (y_src * src->width) + x_src;
 185
 186                         dest->data[0][pixel_dest] = src->data[0][pixel_src];
 187                         dest->data[1][pixel_dest] = src->data[1][pixel_src];
 188                         dest->data[2][pixel_dest] = src->data[2][pixel_src];
 189                 }
 190         }
 191
 192         return (0);
 193 } /* int copy_pixel */
 194
 195 //    Main transformation function. Destination image is calculated using transformation
 196 //    Function "func". Either all colors (color = 0) or one of rgb (color =1,2,3) are
 197 //    determined. If successful, TrPtr->success = 1. Memory for destination image
 198 //    must have been allocated and locked!
 199
 200 int pl_extract_view (ui_image_t *view, const ui_image_t *pano,
 201                 double pitch, double yaw, double fov, interpolator_t interp)
 202 {
 203         int x_dest, y_dest; // Loop through destination image
 204
 205         int dest_width_left = view->width  / 2 ;
 206         int dest_height_top = view->height / 2 ;
 207         int src_width_left =  pano->width   / 2 ;
 208         int src_height_top =  pano->height  / 2 ;
 209
 210         double v[3];
 211         int     x_min, x_max, y_min, y_max;
 212
 213         /* The transformation matrix */
 214         double tm[3][3];
 215
 216         double dist_r;
 217         double dist_e;
 218
 219         { /* What the fuck does this? -octo */
 220                 double a = DEG_TO_RAD (fov);
 221                 double b = 2.0 * M_PI; /* DEG_TO_RAD (360.0) */
 222
 223                 dist_r = ((double) view->width) / (2.0 * tan (a / 2.0));
 224                 dist_e = ((double) pano->width) / b;
 225         }
 226
 227         set_transformation_matrix (tm, DEG_TO_RAD (pitch), DEG_TO_RAD (yaw), 0.0);
 228
 229         x_min = -dest_width_left; x_max = view->width - dest_width_left;
 230         y_min = -dest_height_top; y_max = view->height - dest_height_top;
 231
 232         for(y_dest = y_min; y_dest < y_max; y_dest++)
 233         {
 234                 for(x_dest = x_min; x_dest < x_max; x_dest++)
 235                 {
 236                         double x_src_fp;
 237                         double y_src_fp;
 238
 239                         v[0] = tm[0][0] * x_dest + tm[1][0] * y_dest + tm[2][0] * dist_r;
 240                         v[1] = tm[0][1] * x_dest + tm[1][1] * y_dest + tm[2][1] * dist_r;
 241                         v[2] = tm[0][2] * x_dest + tm[1][2] * y_dest + tm[2][2] * dist_r;
 242
 243                         x_src_fp = dist_e * atan2 (v[0], v[2]);
 244                         y_src_fp = dist_e * atan2 (v[1], sqrt (v[2] * v[2] + v[0] * v[0]));
 245
 246                         copy_pixel (view, pano,
 247                                         dest_width_left + x_dest, dest_height_top + y_dest,
 248                                         src_width_left + x_src_fp, src_height_top + y_src_fp,
 249                                         interp); /* FIXME */
 250                 }
 251         }
 252
 253         return (0);
 254 }
 255
 256 /*
 257  * vim: set tabstop=4 softtabstop=4 shiftwidth=4 :
 258  */