From: Florian Forster Date: Thu, 16 Aug 2007 13:36:01 +0000 (+0200) Subject: Moved the square-root-function to `src/utils_math.c'. X-Git-Url: https://git.verplant.org/?a=commitdiff_plain;h=0d448cc6928c5e448bc6b42dfb39fc4c193d4733;p=libopano.git Moved the square-root-function to `src/utils_math.c'. --- diff --git a/src/panolib.c b/src/panolib.c index 58d190b..9a6a12d 100644 --- a/src/panolib.c +++ b/src/panolib.c @@ -17,6 +17,7 @@ /*------------------------------------------------------------*/ #include "filter.h" +#include "utils_math.h" // Lookup Tables for Trig-functions and interpolator @@ -79,31 +80,33 @@ static void bil( unsigned char *dst, unsigned char **rgb, -// Extract image from pano in TrPtr->src -// using parameters in prefs (ignore image parameters -// in TrPtr) - +/* + * Extract image from pano in TrPtr->src using parameters in prefs (ignore + * image parameters in TrPtr) + */ void PV_ExtractStill( TrformStr *TrPtr ) { - double a,b; // field of view in rad + /* field of view in rad */ + double a; + double b; + double p[2]; double mt[3][3]; int mi[3][3],i,k; - TrPtr->success = 1; - a = DEG_TO_RAD( TrPtr->dest->hfov ); // field of view in rad b = DEG_TO_RAD( TrPtr->src->hfov ); - SetMatrix( DEG_TO_RAD( TrPtr->dest->pitch ), - DEG_TO_RAD( TrPtr->dest->yaw ), - 0.0 , - mt, - 1 ); + /* 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); p[0] = (double)TrPtr->dest->width/ (2.0 * tan( a / 2.0 ) ); - p[1] = (double)TrPtr->src->width/ b; + p[1] = (double) TrPtr->src->width / b; for(i=0; i<3; i++){ for(k=0; k<3; k++){ @@ -127,13 +130,12 @@ void PV_ExtractStill( TrformStr *TrPtr ) void PV_transForm( TrformStr *TrPtr, int dist_r, int dist_e, int mt[3][3]) { - register int x, y; // Loop through destination image + int x_dest, y_dest; // Loop through destination image unsigned char *dest,*src,*sry, *dst;// Source and destination image data long cy; // rownum in destimage int dx,dy; - int xc; // Cartesian Coordinates of point ("target") in Destination image - int xs, ys; + int x_src, y_src; unsigned char *rgb[2] , cdata[16]; // Image data handed to sampler @@ -155,15 +157,11 @@ void PV_transForm( TrformStr *TrPtr, int dist_r, int dist_e, int mt[3][3]) int x_min, x_max, y_min, y_max; int dr1, dr2, dr3; - + 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 @@ -174,74 +172,72 @@ void PV_transForm( TrformStr *TrPtr, int dist_r, int dist_e, int mt[3][3]) if( TrPtr->interpolator == _bilinear ) { - for(y = y_min; y < y_max; y++, cy+=TrPtr->dest->bytesPerLine) + for(y_dest = y_min; y_dest < y_max; y_dest++, cy+=TrPtr->dest->bytesPerLine) { dst = dest + cy + 1; - for(x = x_min; x < x_max; x++, dst+=4) + for(x_dest = x_min; x_dest < x_max; x_dest++, dst+=4) { - v[0] = mt[0][0] * x + mt[1][0] * y + dr1; - v[1] = mt[0][1] * x + mt[1][1] * y + dr2; - v[2] = mt[0][2] * x + mt[1][2] * y + dr3; + 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; v[0] = v[0] >> 8; v[2] = v[2] >> 8; - xs = dist_e * PV_atan2( v[0], v[2] ) / NATAN ; - - ys = dist_e * PV_atan2( v[1], PV_sqrt( abs(v[2]), abs(v[0]) ) ) / NATAN ; + x_src = dist_e * PV_atan2( v[0], v[2] ) / NATAN ; + y_src = dist_e * PV_atan2( v[1], utils_sqrt2 ( abs(v[2]), abs(v[0]) ) ) / NATAN ; - dx = xs & 255; dy = ys & 255; // fraction - + dx = x_src & 255; dy = y_src & 255; // fraction - xs = (xs >> 8) + sw2; - ys = (ys >> 8) + sh2; + x_src = (x_src >> 8) + sw2; + y_src = (y_src >> 8) + sh2; - if( ys >= 0 && ys < miy && xs >= 0 && xs < mix ) // all interpolation pixels inside image + if( y_src >= 0 && y_src < miy && x_src >= 0 && x_src < mix ) // all interpolation pixels inside image // (most pixels) { - sry = src + ys * BytesPerLine + xs * 4; + sry = src + y_src * BytesPerLine + x_src * 4; rgb[0] = sry; rgb[1] = sry + BytesPerLine; } else // edge pixels { - xc = xs; + int x_copy = x_src; rgb[0] = cdata; - if( ys < 0 ) + if( y_src < 0 ) sry = src; - else if( ys > miy ) + else if( y_src > miy ) sry = src + miy * BytesPerLine; else - sry = src + ys * BytesPerLine; + sry = src + y_src * BytesPerLine; - if( xs < 0 ) xs = mix; - if( xs > mix) xs = 0; - *(long*)rgb[0] = *(long*)(sry + xs*4); + if( x_src < 0 ) x_src = mix; + if( x_src > mix) x_src = 0; + *(long*)rgb[0] = *(long*)(sry + x_src*4); - xs = xc+1; - if( xs < 0 ) xs = mix; - if( xs > mix) xs = 0; - *(long*)(rgb[0]+4) = *(long*)(sry + xs*4); + x_src = x_copy+1; + if( x_src < 0 ) x_src = mix; + if( x_src > mix) x_src = 0; + *(long*)(rgb[0]+4) = *(long*)(sry + x_src*4); rgb[1] = cdata+8; - ys+=1; - if( ys < 0 ) + y_src+=1; + if( y_src < 0 ) sry = src; - else if( ys > miy ) + else if( y_src > miy ) sry = src + miy * BytesPerLine; else - sry = src + ys * BytesPerLine; - xs = xc; - if( xs < 0 ) xs = mix; - if( xs > mix) xs = 0; - *(long*)rgb[1] = *(long*)(sry + xs*4); - xs = xc+1; - if( xs < 0 ) xs = mix; - if( xs > mix) xs = 0; - *(long*)(rgb[1]+4) = *(long*)(sry + xs*4); + sry = src + y_src * BytesPerLine; + x_src = x_copy; + if( x_src < 0 ) x_src = mix; + if( x_src > mix) x_src = 0; + *(long*)rgb[1] = *(long*)(sry + x_src*4); + x_src = x_copy+1; + if( x_src < 0 ) x_src = mix; + if( x_src > mix) x_src = 0; + *(long*)(rgb[1]+4) = *(long*)(sry + x_src*4); } @@ -251,33 +247,33 @@ void PV_transForm( TrformStr *TrPtr, int dist_r, int dist_e, int mt[3][3]) } else if( TrPtr->interpolator == _nn ) { - for(y = y_min; y < y_max; y++, cy+=TrPtr->dest->bytesPerLine) + for(y_dest = y_min; y_dest < y_max; y_dest++, cy+=TrPtr->dest->bytesPerLine) { dst = dest + cy; - for(x = x_min; x < x_max; x++, dst+=4) + for(x_dest = x_min; x_dest < x_max; x_dest++, dst+=4) { - v[0] = mt[0][0] * x + mt[1][0] * y + dr1; - v[1] = mt[0][1] * x + mt[1][1] * y + dr2; - v[2] = mt[0][2] * x + mt[1][2] * y + dr3; + 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; v[0] = v[0] >> 8; v[2] = v[2] >> 8; - xs = dist_e * PV_atan2( v[0], v[2] ) / NATAN ; + x_src = dist_e * PV_atan2( v[0], v[2] ) / NATAN ; - ys = dist_e * PV_atan2( v[1], PV_sqrt( abs(v[2]), abs(v[0]) ) ) / NATAN ; + y_src = dist_e * PV_atan2( v[1], utils_sqrt2( abs(v[2]), abs(v[0]) ) ) / NATAN ; - dx = xs & 255; dy = ys & 255; // fraction + dx = x_src & 255; dy = y_src & 255; // fraction - xs = (xs >> 8) + sw2; - ys = (ys >> 8) + sh2; + x_src = (x_src >> 8) + sw2; + y_src = (y_src >> 8) + sh2; - if( xs < 0 ) xs = 0; - if( xs > mix ) xs = mix; - if( ys < 0) ys = 0; - if( ys > miy ) ys = miy; + if( x_src < 0 ) x_src = 0; + if( x_src > mix ) x_src = mix; + if( y_src < 0) y_src = 0; + if( y_src > miy ) y_src = miy; - *(long*)dst = *(long*)(src + ys * BytesPerLine + xs * 4); + *(long*)dst = *(long*)(src + y_src * BytesPerLine + x_src * 4); } } } @@ -286,7 +282,7 @@ void PV_transForm( TrformStr *TrPtr, int dist_r, int dist_e, int mt[3][3]) return; } - +#if 0 void matrix_inv_mult( double m[3][3], double vector[3] ) { register int i; @@ -299,6 +295,7 @@ void matrix_inv_mult( double m[3][3], double vector[3] ) 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 @@ -321,25 +318,23 @@ void SetMatrix( double a, double b, double c , double m[3][3], int cl ) 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 */ @@ -547,3 +542,6 @@ int extractJPEG( fullPath *image, fullPath *jpeg ) } +/* + * vim: set tabstop=4 softtabstop=4 shiftwidth=4 : + */ diff --git a/src/utils_math.c b/src/utils_math.c new file mode 100644 index 0000000..09b424f --- /dev/null +++ b/src/utils_math.c @@ -0,0 +1,141 @@ +/* + * Copyright (C) 1998,1999 Helmut Dersch + * Copyright (C) 2007 Florian octo Forster + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 2 of the License, or (at your option) + * any later version. + * + * This program 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 General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., 51 + * Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + **/ + +#include +#include +#include + +#include "utils_math.h" + +// Lookup Tables for Trig-functions and interpolator + +#define ATAN_TABLE_SIZE 4096 +#define SQRT_TABLE_SIZE 4096 + +static int *atan_table = NULL; +static int *sqrt_table = NULL; + +static int setup_sqrt_table (void) +{ + double z; + int i; + + sqrt_table = (int *) malloc (SQRT_TABLE_SIZE * sizeof (int)); + if (sqrt_table == NULL) + return (-1); + + /* Calculate the table entries. */ + for (i = 0; i < SQRT_TABLE_SIZE; i++) + { + z = ((double) i) / ((double) (SQRT_TABLE_SIZE - 1)); + sqrt_table[i] = (int) (sqrt( 1.0 + z*z ) * 4096.0 + .5); + } + + return (0); +} /* int setup_sqrt_table */ + +/* + * Calculate + * 256 * sqrt(x1^2 + x2^2) + * efficiently + */ +int utils_sqrt2 (int x1, int x2) +{ + int ret; + + if (sqrt_table == NULL) + if (setup_sqrt_table () != 0) + return (-1); + + if ((x1 == 0) && (x2 == 0)) + { + ret = 0; + } + else if (x1 > x2) + { + ret = x1 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x2 / x1] / 16; + } + else /* (x2 < x1) */ + { + ret = x2 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x1 / x2] / 16; + } + + return (ret); +} /* int utils_sqrt2 */ + +static int setup_atan_table (void) +{ + int i; + double z; + + atan_table = (int *) malloc (ATAN_TABLE_SIZE * sizeof (int)); + if (atan_table == NULL) + return (-1); + + for (i = 0; i < (ATAN_TABLE_SIZE - 1); i++) + { + z = ((double) i) / ((double) (ATAN_TABLE_SIZE - 1)); + atan_table[i] = atan ( z / (1.0 - z ) ) * ATAN_TABLE_SIZE * 256; + } + atan_table[ATAN_TABLE_SIZE - 1] = M_PI_4 * ATAN_TABLE_SIZE * 256; + + return 0; +} /* int setup_atan_table */ + +int utils_atan2 (int y, int x) +{ + if (atan_table == NULL) + if (setup_atan_table () != 0) + return (-1); + + // return atan2(y,x) * 256*ATAN_TABLE_SIZE; + if( x > 0 ) + { + if( y > 0 ) + { + return atan_table[(int)( ATAN_TABLE_SIZE * y / ( x + y ))]; + } + else + { + return -atan_table[ (int)(ATAN_TABLE_SIZE * (-y) / ( x - y ))]; + } + } + + if( x == 0 ) + { + if( y > 0 ) + return (int)(256*ATAN_TABLE_SIZE * M_PI_2); + else + return -(int)(256*ATAN_TABLE_SIZE * M_PI_2); + } + + if( y < 0 ) + { + return atan_table[(int)( ATAN_TABLE_SIZE * y / ( x + y ))] - (int)(M_PI*256*ATAN_TABLE_SIZE); + } + else + { + return -atan_table[ (int)(ATAN_TABLE_SIZE * (-y) / ( x - y ))] + (int)(M_PI*256*ATAN_TABLE_SIZE); + } + +} + +/* + * vim: set tabstop=8 softtabstop=2 shiftwidth=2 : + */ diff --git a/src/utils_math.h b/src/utils_math.h new file mode 100644 index 0000000..08e2831 --- /dev/null +++ b/src/utils_math.h @@ -0,0 +1,20 @@ +/* + * Copyright (C) 2007 Florian octo Forster + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 2 of the License, or (at your option) + * any later version. + * + * This program 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 General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., 51 + * Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + **/ + +int utils_sqrt2 (int x1, int x2); +int utils_atan2 (int y, int x);