im_math.h

Go to the documentation of this file.
00001 /** \file
00002  * \brief Math Utilities
00003  *
00004  * See Copyright Notice in im_lib.h
00005  */
00006 
00007 #ifndef __IM_MATH_H
00008 #define __IM_MATH_H
00009 
00010 #include <math.h>
00011 #include "im_util.h"
00012 
00013 #ifdef IM_DEFMATHFLOAT
00014 inline float acosf(float _X) {return ((float)acos((double)_X)); }
00015 inline float asinf(float _X) {return ((float)asin((double)_X)); }
00016 inline float atanf(float _X) {return ((float)atan((double)_X)); }
00017 inline float atan2f(float _X, float _Y) {return ((float)atan2((double)_X, (double)_Y)); }
00018 inline float ceilf(float _X) {return ((float)ceil((double)_X)); }
00019 inline float cosf(float _X)  {return ((float)cos((double)_X)); }
00020 inline float coshf(float _X) {return ((float)cosh((double)_X)); }
00021 inline float expf(float _X) {return ((float)exp((double)_X)); }
00022 inline float fabsf(float _X) {return ((float)fabs((double)_X)); }
00023 inline float floorf(float _X) {return ((float)floor((double)_X)); }
00024 inline float fmodf(float _X, float _Y) {return ((float)fmod((double)_X, (double)_Y)); }
00025 inline float logf(float _X) {return ((float)log((double)_X)); }
00026 inline float log10f(float _X) {return ((float)log10((double)_X)); }
00027 inline float powf(float _X, float _Y) {return ((float)pow((double)_X, (double)_Y)); }
00028 inline float sinf(float _X) {return ((float)sin((double)_X)); }
00029 inline float sinhf(float _X) {return ((float)sinh((double)_X)); }
00030 inline float sqrtf(float _X) {return ((float)sqrt((double)_X)); }
00031 inline float tanf(float _X) {return ((float)tan((double)_X)); }
00032 inline float tanhf(float _X) {return ((float)tanh((double)_X)); }
00033 #endif
00034 
00035 /** \defgroup math Math Utilities
00036  * \par
00037  * When converting between continuous and discrete use: \n
00038  * Continuous = Discrete + 0.5         [Reconstruction/Interpolation] \n
00039  * Discrete = Round(Continuous - 0.5)  [Sampling/Quantization] \n
00040  * \par
00041  * Notice that must check 0-max limits when converting from Continuous to Discrete.
00042  * \par
00043  * When converting between discrete and discrete use: \n
00044  * integer src_size, dst_len, src_i, dst_i            \n
00045  * real factor = (real)(dst_size)/(real)(src_size)    \n
00046  * dst_i = Round(factor*(src_i + 0.5) - 0.5)
00047  * \par
00048  * See \ref im_math.h
00049  * \ingroup util */
00050 
00051 
00052 /** Round a real to the nearest integer.
00053  * \ingroup math */
00054 inline int imRound(float x) 
00055 {
00056   return (int)(x < 0? x-0.5f: x+0.5f);
00057 }
00058 inline int imRound(double x) 
00059 {
00060   return (int)(x < 0? x-0.5: x+0.5);
00061 }
00062 
00063 /** Converts between two discrete grids.
00064  * factor is "dst_size/src_size".
00065  * \ingroup math */
00066 inline int imResample(int x, float factor)
00067 {
00068   float xr = factor*(x + 0.5f) - 0.5f;
00069   return (int)(xr < 0? xr-0.5f: xr+0.5f);  /* Round */
00070 }
00071 
00072 /** Does Zero Order Decimation (Mean).
00073  * \ingroup math */
00074 template <class T, class TU>
00075 inline T imZeroOrderDecimation(int width, int height, T *map, float xl, float yl, float box_width, float box_height, TU Dummy)
00076 {
00077   int x0,x1,y0,y1;
00078   (void)Dummy;
00079 
00080   x0 = (int)floor(xl - box_width/2.0 - 0.5) + 1;
00081   y0 = (int)floor(yl - box_height/2.0 - 0.5) + 1;
00082   x1 = (int)floor(xl + box_width/2.0 - 0.5);
00083   y1 = (int)floor(yl + box_height/2.0 - 0.5);
00084 
00085   if (x0 == x1) x1++;
00086   if (y0 == y1) y1++;
00087 
00088   x0 = x0<0? 0: x0>width-1? width-1: x0;
00089   y0 = y0<0? 0: y0>height-1? height-1: y0;
00090   x1 = x1<0? 0: x1>width-1? width-1: x1;
00091   y1 = y1<0? 0: y1>height-1? height-1: y1;
00092 
00093   TU Value;
00094   int Count = 0;
00095 
00096   Value = 0;
00097 
00098   for (int y = y0; y <= y1; y++)
00099   {
00100     for (int x = x0; x <= x1; x++)
00101     {
00102       Value += map[y*width+x];
00103       Count++;
00104     }
00105   }
00106 
00107   if (Count == 0)
00108   {
00109     Value = 0;
00110     return (T)Value;
00111   }
00112 
00113   return (T)(Value/(float)Count);
00114 }
00115 
00116 /** Does Bilinear Decimation.
00117  * \ingroup math */
00118 template <class T, class TU>
00119 inline T imBilinearDecimation(int width, int height, T *map, float xl, float yl, float box_width, float box_height, TU Dummy)
00120 {
00121   int x0,x1,y0,y1;
00122   (void)Dummy;
00123 
00124   x0 = (int)floor(xl - box_width/2.0 - 0.5) + 1;
00125   y0 = (int)floor(yl - box_height/2.0 - 0.5) + 1;
00126   x1 = (int)floor(xl + box_width/2.0 - 0.5);
00127   y1 = (int)floor(yl + box_height/2.0 - 0.5);
00128 
00129   if (x0 == x1) x1++;
00130   if (y0 == y1) y1++;
00131 
00132   x0 = x0<0? 0: x0>width-1? width-1: x0;
00133   y0 = y0<0? 0: y0>height-1? height-1: y0;
00134   x1 = x1<0? 0: x1>width-1? width-1: x1;
00135   y1 = y1<0? 0: y1>height-1? height-1: y1;
00136 
00137   TU Value, LineValue;
00138   float LineNorm, Norm, dxr, dyr;
00139 
00140   Value = 0;
00141   Norm = 0;
00142 
00143   for (int y = y0; y <= y1; y++)
00144   {
00145     dyr = yl - (y+0.5f);
00146     if (dyr < 0) dyr *= -1;
00147 
00148     LineValue = 0;
00149     LineNorm = 0;
00150 
00151     for (int x = x0; x <= x1; x++)
00152     {
00153       dxr = xl - (x+0.5f);
00154       if (dxr < 0) dxr *= -1;
00155 
00156       LineValue += map[y*width+x] * dxr;
00157       LineNorm += dxr;
00158     }
00159 
00160     Value += LineValue * dyr;
00161     Norm += dyr * LineNorm;
00162   }
00163 
00164   if (Norm == 0)
00165   {
00166     Value = 0;
00167     return (T)Value;
00168   }
00169 
00170   return (T)(Value/Norm);
00171 }
00172 
00173 /** Does Zero Order Interpolation (Nearest Neighborhood).
00174  * \ingroup math */
00175 template <class T>
00176 inline T imZeroOrderInterpolation(int width, int height, T *map, float xl, float yl)
00177 {
00178   int x0 = imRound(xl-0.5f);
00179   int y0 = imRound(yl-0.5f);
00180   x0 = x0<0? 0: x0>width-1? width-1: x0;
00181   y0 = y0<0? 0: y0>height-1? height-1: y0;
00182   return map[y0*width + x0];
00183 }
00184 
00185 /** Does Bilinear Interpolation.
00186  * \ingroup math */
00187 template <class T>
00188 inline T imBilinearInterpolation(int width, int height, T *map, float xl, float yl)
00189 {
00190   int x0, y0, x1, y1;
00191   float t, u;
00192 
00193   if (xl < 0.5)
00194   {
00195     x1 = x0 = 0; 
00196     t = 0;
00197   }
00198   else if (xl > width-0.5)
00199   {
00200     x1 = x0 = width-1;
00201     t = 0;
00202   }
00203   else
00204   {
00205     x0 = (int)(xl-0.5f);
00206     x1 = x0+1;
00207     t = xl - (x0+0.5f);
00208   }
00209 
00210   if (yl < 0.5)
00211   {
00212     y1 = y0 = 0; 
00213     u = 0;
00214   }
00215   else if (yl > height-0.5)
00216   {
00217     y1 = y0 = height-1;
00218     u = 0;
00219   }
00220   else
00221   {
00222     y0 = (int)(yl-0.5f);
00223     y1 = y0+1;
00224     u = yl - (y0+0.5f);
00225   }
00226 
00227   T fll = map[y0*width + x0];
00228   T fhl = map[y0*width + x1];
00229   T flh = map[y1*width + x0];
00230   T fhh = map[y1*width + x1];
00231 
00232   return (T)((fhh - flh - fhl + fll) * u * t +
00233                          (fhl - fll) * t +
00234                          (flh - fll) * u +
00235                                 fll);
00236 }
00237 
00238 /** Does Bicubic Interpolation.
00239  * \ingroup math */
00240 template <class T, class TU>
00241 inline T imBicubicInterpolation(int width, int height, T *map, float xl, float yl, TU Dummy)
00242 {
00243   int X[4], Y[4];
00244   float t, u;
00245   (void)Dummy;
00246 
00247   if (xl > width-0.5)
00248   {
00249     X[3] = X[2] = X[1] = width-1;
00250     X[0] = X[1]-1;
00251     t = 0;
00252   }
00253   else
00254   {
00255     X[1] = (int)(xl-0.5f);
00256     if (X[1] < 0) X[1] = 0;
00257 
00258     X[0] = X[1]-1;
00259     X[2] = X[1]+1;
00260     X[3] = X[1]+2;
00261 
00262     if (X[0] < 0) X[0] = 0;
00263     if (X[3] > width-1) X[3] = width-1;
00264 
00265     t = xl - (X[1]+0.5f);
00266   }
00267 
00268   if (yl > height-0.5)
00269   {
00270     Y[3] = Y[2] = Y[1] = height-1;
00271     Y[0] = Y[1]-1;
00272     u = 0;
00273   }
00274   else
00275   {
00276     Y[1] = (int)(yl-0.5f);
00277     if (Y[1] < 0) Y[1] = 0;
00278 
00279     Y[0] = Y[1]-1;
00280     Y[2] = Y[1]+1;
00281     Y[3] = Y[1]+2;
00282 
00283     if (Y[0] < 0) Y[0] = 0;
00284     if (Y[3] > height-1) Y[3] = height-1;
00285 
00286     u = yl - (Y[1]+0.5f);
00287   }
00288 
00289   float CX[4], CY[4];
00290 
00291   // Optimize calculations
00292   {
00293     float c, c2, c3;
00294 
00295 #define C0 (-c3 + 2.0f*c2 - c)
00296 #define C1 ( c3 - 2.0f*c2 + 1.0f)
00297 #define C2 (-c3 + c2 + c)
00298 #define C3 ( c3 - c2)
00299 
00300     c = t;
00301     c2 = c*c; c3 = c2*c;
00302     CX[0] = C0; CX[1] = C1; CX[2] = C2; CX[3] = C3;
00303 
00304     c = u;
00305     c2 = c*c; c3 = c2*c;
00306     CY[0] = C0; CY[1] = C1; CY[2] = C2; CY[3] = C3;
00307 
00308 #undef C0
00309 #undef C1
00310 #undef C2
00311 #undef C3
00312   }
00313 
00314   TU LineValue, Value;
00315   float LineNorm, Norm;
00316 
00317   Value = 0;
00318   Norm = 0;
00319 
00320   for (int y = 0; y < 4; y++)
00321   {
00322     LineValue = 0;
00323     LineNorm = 0;
00324 
00325     for (int x = 0; x < 4; x++)
00326     {
00327       LineValue += map[Y[y]*width+X[x]] * CX[x];
00328       LineNorm += CX[x];
00329     }
00330 
00331     Value += LineValue * CY[y];
00332     Norm += CY[y] * LineNorm;
00333   }
00334 
00335   if (Norm == 0)
00336   {
00337     Value = 0;
00338     return (T)Value;
00339   }
00340 
00341   Value = (Value/Norm);
00342 
00343   int size = sizeof(T); 
00344   if (size == 1)
00345     return (T)(Value<=(TU)0? (TU)0: Value<=(TU)255? Value: (TU)255);
00346   else
00347     return (T)(Value);
00348 }
00349 
00350 /** Calculates minimum and maximum values.
00351  * \ingroup math */
00352 template <class T> 
00353 inline void imMinMax(const T *map, int count, T& min, T& max)
00354 {
00355   min = *map++;
00356   max = min;
00357   for (int i = 1; i < count; i++)
00358   {
00359     T value = *map++;
00360 
00361     if (value > max)
00362       max = value;
00363     else if (value < min)
00364       min = value;
00365   }
00366 }
00367 
00368 #endif

Generated on Thu Oct 1 11:40:01 2009 for IM by  doxygen 1.6.1