panolib.c: Use something similar to Lanczos resampling for sharper images.
[libopano.git] / src / panolib.c
index e9c04ef..58948d3 100644 (file)
 #include <assert.h>
 
 #include "filter.h"
+#include "panolib.h"
 #include "utils_math.h"
+#include "utils_image.h"
 
+#define DEG_TO_RAD(x) ((x) * 2.0 * M_PI / 360.0 )
 
 // Lookup Tables for Trig-functions and interpolator
 
 #define NATAN 2048
 #define NSQRT 2048
 
-int *atan_LU;
-int *sqrt_LU;
-int *mweights[256];
+static void matrix_matrix_mult( double m1[3][3],double m2[3][3],double result[3][3])
+{
+       int i,k;
+       
+       for(i=0;i<3;i++)
+               for(k=0; k<3; k++)
+                       result[i][k] = m1[i][0] * m2[0][k] + m1[i][1] * m2[1][k] + m1[i][2] * m2[2][k];
+} /* void matrix_matrix_mult */
 
+/* Set matrix elements based on Euler angles a, b, c */
+static void set_transformation_matrix (double m[3][3],
+               double a, double b, double c)
+{
+       double mx[3][3], my[3][3], mz[3][3], dummy[3][3];
+       
 
-void   matrix_matrix_mult      ( double m1[3][3],double m2[3][3],double result[3][3]);
-void   PV_transForm( TrformStr *TrPtr, double dist_r, double dist_e, int mt[3][3]);
-int    PV_atan2(int y, int x);
-int    PV_sqrt( int x1, int x2 );
+       // Calculate Matrices;
 
+       mx[0][0] = 1.0 ;                                mx[0][1] = 0.0 ;                                mx[0][2] = 0.0;
+       mx[1][0] = 0.0 ;                                mx[1][1] = cos(a) ;                     mx[1][2] = sin(a);
+       mx[2][0] = 0.0 ;                                mx[2][1] =-mx[1][2] ;                   mx[2][2] = mx[1][1];
+       
+       my[0][0] = cos(b);                              my[0][1] = 0.0 ;                                my[0][2] =-sin(b);
+       my[1][0] = 0.0 ;                                my[1][1] = 1.0 ;                                my[1][2] = 0.0;
+       my[2][0] = -my[0][2];                   my[2][1] = 0.0 ;                                my[2][2] = my[0][0];
+       
+       mz[0][0] = cos(c) ;                     mz[0][1] = sin(c) ;                     mz[0][2] = 0.0;
+       mz[1][0] =-mz[0][1] ;                   mz[1][1] = mz[0][0] ;                   mz[1][2] = 0.0;
+       mz[2][0] = 0.0 ;                                mz[2][1] = 0.0 ;                                mz[2][2] = 1.0;
 
-/*
- * Extract image from pano in TrPtr->src using parameters in prefs (ignore
- * image parameters in TrPtr)
- */
-void PV_ExtractStill( TrformStr *TrPtr )
-{
-       /* field of view in rad */
-       double a;
-       double b;
+       /* Calculate `m = mz * mx * my' */
 
-       double      p[2];
-       double          mt[3][3];
-       int             mi[3][3],i,k;
+       matrix_matrix_mult( mz, mx,     dummy);
+       matrix_matrix_mult( dummy, my, m);
+} /* void SetMatrix */
 
-       a =      DEG_TO_RAD( TrPtr->dest->hfov );       // field of view in rad         
-       b =      DEG_TO_RAD( TrPtr->src->hfov );
+#if 0
+static double lanczos (double x)
+{
+       static double table[257];
 
-       /* Set up the transformation matrix `mt' using Euler angles (somehow..) */
-       SetMatrix (DEG_TO_RAD (TrPtr->dest->pitch), /* alpha */
-                       DEG_TO_RAD (TrPtr->dest->yaw), /* beta */
-                       0.0, /* gamma */
-                       mt, /* output */
-                       1);
+       assert ((x >= 0.0) && (x <= 1.0));
 
+       if (table[0] != 1.0)
+       {
+               int i;
 
-       p[0] = (double)TrPtr->dest->width/ (2.0 * tan( a / 2.0 ) );
-       p[1] = (double) TrPtr->src->width / b;
-       
-       for(i=0; i<3; i++){
-               for(k=0; k<3; k++){
-                       mi[i][k] = 256 * mt[i][k];
+               table[0] = 1.0;
+               for (i = 1; i <= 256; i++)
+               {
+                       double t = ((double) i) / 256.0;
+                       table[i] = 2 * sin (M_PI * t) * sin (M_PI_2 * t) / (M_PI * M_PI * t * t);
                }
        }
 
+       return (table[(int) (0.5 + 256.0 * x)]);
+} /* double lanczos */
+#endif
 
-       PV_transForm( TrPtr, p[0], p[1], mi);
-       return;
-}
-
-typedef enum interpolator_e
-{
-       NNEIGHBOUR,
-       BILINEAR
-} interpolator_t;
-
-static int copy_pixel (Image *dest, Image *src,
+static int copy_pixel (ui_image_t *dest, const ui_image_t *src,
                int x_dest, int y_dest,
                double x_src_fp, double y_src_fp,
                interpolator_t interp)
 {
-       uint32_t pixel_value;
-       uint32_t *src_data  = (uint32_t *) (*src->data);
-       uint32_t *dest_data = (uint32_t *) (*dest->data);
+       uint32_t pixel_dest = (y_dest * dest->width) + x_dest;
+
+       assert (x_src_fp >= 0.0);
 
-       interp = NNEIGHBOUR;
+       if (y_src_fp < 0.0)
+               y_src_fp = 0.0;
 
-       if (interp == NNEIGHBOUR)
+       if (interp == BILINEAR)
        {
-               int x_src = (int) (x_src_fp + 0.5) % src->width;
-               int y_src = (int) (y_src_fp + 0.5) % src->height;
+               int x_src_left;
+               int x_src_right;
+               int y_src_top;
+               int y_src_bottom;
 
-               if ((x_src < 0) || (x_src >= src->width)
-                               || (y_src < 0) || (y_src >= src->height))
-                       pixel_value = 0x000000FF;
-               else
-                       pixel_value = src_data[(y_src * src->width) + x_src];
-       }
+               double x_right_frac;
+               double y_bottom_frac;
 
-       dest_data[(y_dest * dest->width) + x_dest] = pixel_value;
+               int i;
 
-       return (0);
-} /* int copy_pixel */
+               x_src_left = (int) x_src_fp;
+               x_src_right = (x_src_left + 1) % src->width;
 
-//    Main transformation function. Destination image is calculated using transformation
-//    Function "func". Either all colors (color = 0) or one of rgb (color =1,2,3) are
-//    determined. If successful, TrPtr->success = 1. Memory for destination image
-//    must have been allocated and locked!
-
-void PV_transForm( TrformStr *TrPtr, double dist_r, double dist_e, int mt[3][3])
-{
-       int x_dest, y_dest;                     // Loop through destination image
-       unsigned char           *dest,*src;// Source and destination image data
-
-       // Variables used to convert screen coordinates to cartesian coordinates
+               y_src_top  = (int) y_src_fp;
+               y_src_bottom = y_src_top + 1;
+               if (y_src_bottom >= src->height)
+                       y_src_bottom = src->height - 1;
 
-               
-       int                             dest_width_left = TrPtr->dest->width  / 2 ;  
-       int                             dest_height_top = TrPtr->dest->height / 2 ;
-       int                             src_width_left =  TrPtr->src->width   / 2 ;
-       int                             src_height_top =  TrPtr->src->height  / 2 ;
-       
-       int                             v[3];
-       int                                     x_min, x_max, y_min, y_max;
+               x_right_frac = x_src_fp - x_src_left;
+               y_bottom_frac = y_src_fp - y_src_top;
 
-       int                                     dr1, dr2, dr3;
+               /*
+                * The polynomial function
+                *   2x^3 - 3x^2 + 1
+                * is similar to the Lanczos function, but much easier to compute. It's
+                * applied here to improve sharpness. 
+                */
+               {
+                       double t = x_src_fp - x_src_left;
+                       x_right_frac = -2.0 * t * t * t + 3.0 * t * t;
 
-       dr1 = mt[2][0] * dist_r;
-       dr2 = mt[2][1] * dist_r;
-       dr3 = mt[2][2] * dist_r;
-       
-       dest = *TrPtr->dest->data;
-       src  = *TrPtr->src->data; // is locked
+                       t = y_src_fp - y_src_top;
+                       y_bottom_frac = -2.0 * t * t * t + 3.0 * t * t;
+               }
 
-       x_min = -dest_width_left; x_max = TrPtr->dest->width - dest_width_left;
-       y_min = -dest_height_top; y_max = TrPtr->dest->height - dest_height_top;
+               assert ((x_right_frac >= 0.0) && (x_right_frac <= 1.0));
+               assert ((y_bottom_frac >= 0.0) && (y_bottom_frac <= 1.0));
 
-       for(y_dest = y_min; y_dest < y_max; y_dest++)
-       {
-               for(x_dest = x_min; x_dest < x_max; x_dest++)
+               for (i = 0; i < 3; i++)
                {
-                       double x_src_fp;
-                       double y_src_fp;
+                       uint8_t values[2][2];
+                       double value_left;
+                       double value_right;
+                       double value_final;
 
-                       v[0] = mt[0][0] * x_dest + mt[1][0] * y_dest + dr1;
-                       v[1] = mt[0][1] * x_dest + mt[1][1] * y_dest + dr2;
-                       v[2] = mt[0][2] * x_dest + mt[1][2] * y_dest + dr3;
+                       values[0][0] = src->data[i][(y_src_top * src->width) + x_src_left];
+                       values[0][1] = src->data[i][(y_src_top * src->width) + x_src_right];
+                       values[1][0] = src->data[i][(y_src_bottom * src->width) + x_src_left];
+                       values[1][1] = src->data[i][(y_src_bottom * src->width) + x_src_right];
 
-                       v[0] = v[0] >> 8; v[2] = v[2] >> 8;
-                       v[1] = v[1] >> 8;
+                       value_left = (1.0 - y_bottom_frac) * values[0][0]
+                               + y_bottom_frac * values[1][0];
+                       value_right = (1.0 - y_bottom_frac) * values[0][1]
+                               + y_bottom_frac * values[1][1];
 
-                       x_src_fp = dist_e * atan2 (v[0], v[2]);
-                       y_src_fp = dist_e * atan2 (v[1], sqrt (v[2] * v[2] + v[0] * v[0]));
+                       value_final = (1.0 - x_right_frac) * value_left + x_right_frac * value_right;
 
-                       copy_pixel (TrPtr->dest, TrPtr->src,
-                                       dest_width_left + x_dest, dest_height_top + y_dest,
-                                       src_width_left + x_src_fp, src_height_top + y_src_fp,
-                                       NNEIGHBOUR);
+                       assert ((value_final >= 0.0) && (value_final <= 255.0));
+
+                       dest->data[i][pixel_dest] = (uint8_t) (value_final + 0.5);
                }
        }
-       
-       return;
-}
-
-#if 0
-void matrix_inv_mult( double m[3][3], double vector[3] )
-{
-       register int i;
-       register double v0 = vector[0];
-       register double v1 = vector[1];
-       register double v2 = vector[2];
-       
-       for(i=0; i<3; i++)
+       else /* if (interp == NNEIGHBOUR) */
        {
-               vector[i] = m[0][i] * v0 + m[1][i] * v1 + m[2][i] * v2;
-       }
-}
-#endif
-
-// Set matrix elements based on Euler angles a, b, c
-
-void SetMatrix( double a, double b, double c , double m[3][3], int cl )
-{
-       double mx[3][3], my[3][3], mz[3][3], dummy[3][3];
-       
-
-       // Calculate Matrices;
-
-       mx[0][0] = 1.0 ;                                mx[0][1] = 0.0 ;                                mx[0][2] = 0.0;
-       mx[1][0] = 0.0 ;                                mx[1][1] = cos(a) ;                     mx[1][2] = sin(a);
-       mx[2][0] = 0.0 ;                                mx[2][1] =-mx[1][2] ;                   mx[2][2] = mx[1][1];
-       
-       my[0][0] = cos(b);                              my[0][1] = 0.0 ;                                my[0][2] =-sin(b);
-       my[1][0] = 0.0 ;                                my[1][1] = 1.0 ;                                my[1][2] = 0.0;
-       my[2][0] = -my[0][2];                   my[2][1] = 0.0 ;                                my[2][2] = my[0][0];
-       
-       mz[0][0] = cos(c) ;                     mz[0][1] = sin(c) ;                     mz[0][2] = 0.0;
-       mz[1][0] =-mz[0][1] ;                   mz[1][1] = mz[0][0] ;                   mz[1][2] = 0.0;
-       mz[2][0] = 0.0 ;                                mz[2][1] = 0.0 ;                                mz[2][2] = 1.0;
-
-       /* Calculate `m = mz * mx * my' */
-
-       if( cl )
-               matrix_matrix_mult( mz, mx,     dummy);
-       else
-               matrix_matrix_mult( mx, mz,     dummy);
-       matrix_matrix_mult( dummy, my, m);
-} /* void SetMatrix */
-
-void matrix_matrix_mult( double m1[3][3],double m2[3][3],double result[3][3])
-{
-       register int i,k;
-       
-       for(i=0;i<3;i++)
-               for(k=0; k<3; k++)
-                       result[i][k] = m1[i][0] * m2[0][k] + m1[i][1] * m2[1][k] + m1[i][2] * m2[2][k];
-} /* void matrix_matrix_mult */
+               int x_src = (int) (x_src_fp + 0.5) % src->width;
+               int y_src = (int) (y_src_fp + 0.5) % src->height;
 
-#if 0
-static int orig_PV_atan2(int y, int x)
-{
-       // return atan2(y,x) * 256*NATAN;
-       if( x > 0 )
-       {
-               if( y > 0 )
+               if ((x_src < 0) || (x_src >= src->width)
+                               || (y_src < 0) || (y_src >= src->height))
                {
-                       return  atan_LU[(int)( NATAN * y / ( x + y ))];
+                       dest->data[0][pixel_dest] = 0;
+                       dest->data[1][pixel_dest] = 0;
+                       dest->data[2][pixel_dest] = 0;
                }
                else
                {
-                       return -atan_LU[ (int)(NATAN * (-y) / ( x - y ))];
+                       uint32_t pixel_src = (y_src * src->width) + x_src;
+
+                       dest->data[0][pixel_dest] = src->data[0][pixel_src];
+                       dest->data[1][pixel_dest] = src->data[1][pixel_src];
+                       dest->data[2][pixel_dest] = src->data[2][pixel_src];
                }
        }
 
-       if( x == 0 )
-       {
-               if( y > 0 )
-                       return  (int)(256*NATAN*PI / 2.0);
-               else
-                       return  -(int)(256*NATAN*PI / 2.0);
-       }
-       
-       if( y < 0 )
-       {
-               return  atan_LU[(int)( NATAN * y / ( x + y ))] - (int)(PI*256*NATAN);
-       }
-       else
-       {
-               return -atan_LU[ (int)(NATAN * (-y) / ( x - y ))] + (int)(PI*256*NATAN);
-       }
-       
-} /* int orig_PV_atan2 */
-#endif
+       return (0);
+} /* int copy_pixel */
 
-int PV_atan2(int y, int x)
-{
-       double ret = atan2 (y, x) * NATAN * 256.0;
-       return ((int) ret);
-} /* int PV_atan2 */
+//    Main transformation function. Destination image is calculated using transformation
+//    Function "func". Either all colors (color = 0) or one of rgb (color =1,2,3) are
+//    determined. If successful, TrPtr->success = 1. Memory for destination image
+//    must have been allocated and locked!
 
-int SetUpAtan()
+int pl_extract_view (ui_image_t *view, const ui_image_t *pano,
+               double pitch, double yaw, double fov, interpolator_t interp)
 {
-       int i;
-       double dz = 1.0 / (double)NATAN;
-       double z = 0.0;
-       
-       atan_LU = (int*) malloc( (NATAN+1) * sizeof( int ));
-       
-       if( atan_LU == NULL )
-               return -1;
-               
-       for( i=0; i< NATAN; i++, z+=dz )
-               atan_LU[i] = atan( z / (1.0 - z ) ) * NATAN * 256;
-               
-       atan_LU[NATAN] = PI/4.0 * NATAN * 256;
-       
-       // Print a test
-#if 0  
-       for(i = -10; i< 10; i++)
-       {
-               int k;
-               for(k=-10; k<10; k++)
-               {
-                       printf("i =  %d  k = %d   atan2(i,k) = %g    LUatan(i,k) = %g diff = %g\n", i,k,atan2(i,k), 
-                               (double)PV_atan2(i,k) / (256*NATAN) , atan2(i,k) - (double)PV_atan2(i,k) / (256*NATAN));
-               }
-       }
-       exit(0);
-#endif 
-       return 0;
-}
+       int x_dest, y_dest; // Loop through destination image
 
-int SetUpSqrt()
-{
-       int i;
-       double dz = 1.0 / (double)NSQRT;
-       double z = 0.0;
-       
-       sqrt_LU = (int*) malloc( (NSQRT+1) * sizeof( int ));
+       int dest_width_left = view->width  / 2 ;  
+       int dest_height_top = view->height / 2 ;
+       int src_width_left =  pano->width   / 2 ;
+       int src_height_top =  pano->height  / 2 ;
        
-       if( sqrt_LU == NULL )
-               return -1;
-               
-       for( i=0; i< NSQRT; i++, z+=dz )
-               sqrt_LU[i] = sqrt( 1.0 + z*z ) * 256 * NSQRT;
-               
-       sqrt_LU[NSQRT] = sqrt(2.0) * 256 * NSQRT;
-       
-       return 0;
-}
+       double v[3];
+       int     x_min, x_max, y_min, y_max;
 
-int SetUpMweights()
-{
-       int i,k;
-       
-       for(i=0; i<256; i++)
-       {
-               mweights[i] = (int*)malloc( 256 * sizeof(int) );
-               if( mweights[i] == NULL ) return -1;
-       }
-       for(i=0; i<256; i++)
-       {
-               for(k=0; k<256; k++)
-               {
-                       mweights[i][k] = i*k;
-               }
-       }
-       
-       return 0;
-}
+       /* The transformation matrix */
+       double tm[3][3];
 
+       double dist_r;
+       double dist_e;
 
-int PV_sqrt( int x1, int x2 )
-{
-       if( x1 > x2 )
-       {
-               return  x1 * sqrt_LU[ NSQRT * x2 /  x1 ] / NSQRT;
-       }
-       else
-       {
-               if( x2 == 0 ) return 0;
-               return x2 * sqrt_LU[ NSQRT * x1 /  x2 ] / NSQRT;
-       }
-}
+       { /* What the fuck does this? -octo */
+               double a = DEG_TO_RAD (fov);
+               double b = 2.0 * M_PI; /* DEG_TO_RAD (360.0) */
 
+               dist_r = ((double) view->width) / (2.0 * tan (a / 2.0));
+               dist_e = ((double) pano->width) / b;
+       }
 
-#define ID_0 0xff
-#define ID_1 0xd8
-#define ID_2 0xff
-#define ID_3 0xe0
-#define ID_4 0x00
-#define ID_5 0x10
-#define ID_6 0x4a
-#define ID_7 0x46
-#define ID_8 0x49
-#define ID_9 0x46
-#define ID_10 0x00
+       set_transformation_matrix (tm, DEG_TO_RAD (pitch), DEG_TO_RAD (yaw), 0.0);
 
-#define ID_LENGTH 11
+       x_min = -dest_width_left; x_max = view->width - dest_width_left;
+       y_min = -dest_height_top; y_max = view->height - dest_height_top;
 
-// Find last jpeg inside image; create and copy to file jpeg
-int extractJPEG( fullPath *image, fullPath *jpeg )
-{
-       file_spec                               fnum;
-       long                                    count;
-       unsigned char*                  im;
-       int                                     i, idx = -1;
-       unsigned char                   ch;
-               
-       if( myopen( image, read_bin, fnum ) )
-               return -1;
-       
-       
-       count = 1; i=0; // Get file length
-       
-       while( count == 1 )
+       for(y_dest = y_min; y_dest < y_max; y_dest++)
        {
-               myread(  fnum, count, &ch ); 
-               if(count==1) i++;
-       }
-       myclose(fnum);
+               for(x_dest = x_min; x_dest < x_max; x_dest++)
+               {
+                       double x_src_fp;
+                       double y_src_fp;
 
-       count = i;
-       
-       im = (UCHAR*)malloc( count );
-       if( im == NULL )
-       {
-               PrintError("Not enough memory");
-               return -1;
-       }
-       
-       if( myopen( image, read_bin, fnum ) )
-               return -1;
-               
-       myread(fnum,count,im);
-       myclose(fnum);
-       
-       if( i != count )
-               return -1;
-       
-       count -= ID_LENGTH;
-               
-       for(i=0; i<count; i++)
-       {
-               if( im[i] == ID_0 && im[i+1] == ID_1 && im[i+2] == ID_2 && im[i+3] == ID_3 
-                                       && im[i+4] == ID_4 && im[i+5] == ID_5 && im[i+6] == ID_6 && im[i+7] == ID_7
-                                       && im[i+8] == ID_8 && im[i+9] == ID_9 && im[i+10] == ID_10)
-                               idx = i;
-       }
+                       v[0] = tm[0][0] * x_dest + tm[1][0] * y_dest + tm[2][0] * dist_r;
+                       v[1] = tm[0][1] * x_dest + tm[1][1] * y_dest + tm[2][1] * dist_r;
+                       v[2] = tm[0][2] * x_dest + tm[1][2] * y_dest + tm[2][2] * dist_r;
 
-       if( idx == -1 ) // No jpeg found
-       {
-               free(im);
-               return -1;
+                       x_src_fp = dist_e * atan2 (v[0], v[2]);
+                       y_src_fp = dist_e * atan2 (v[1], sqrt (v[2] * v[2] + v[0] * v[0]));
+
+                       copy_pixel (view, pano,
+                                       dest_width_left + x_dest, dest_height_top + y_dest,
+                                       src_width_left + x_src_fp, src_height_top + y_src_fp,
+                                       interp); /* FIXME */
+               }
        }
        
-       count = count + ID_LENGTH - idx;
-       
-       mycreate( jpeg, 'GKON','JPEG');
-       if( myopen( jpeg, write_bin, fnum ) )
-               return -1;
-       mywrite( fnum, count, im+idx );
-       free( im );
-       myclose( fnum );
-       return 0;
+       return (0);
 }
-               
 
 /*
  * vim: set tabstop=4 softtabstop=4 shiftwidth=4 :