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 * See \ref im_color.h
00038 * \ingroup util */
00039
00040
00041 /** Does Zero Order Decimation (Mean).
00042 * \ingroup math */
00043 template <class T, class TU>
00044 inline T imZeroOrderDecimation(int width, int height, T *map, float xl, float yl, float box_width, float box_height, TU Dummy)
00045 {
00046 int x0,x1,y0,y1;
00047 (void)Dummy;
00048
00049 x0 = (int)floor(xl - box_width/2.0 - 0.5) + 1;
00050 y0 = (int)floor(yl - box_height/2.0 - 0.5) + 1;
00051 x1 = (int)floor(xl + box_width/2.0 - 0.5);
00052 y1 = (int)floor(yl + box_height/2.0 - 0.5);
00053
00054 if (x0 == x1) x1++;
00055 if (y0 == y1) y1++;
00056
00057 x0 = x0<0? 0: x0>width-1? width-1: x0;
00058 y0 = y0<0? 0: y0>height-1? height-1: y0;
00059 x1 = x1<0? 0: x1>width-1? width-1: x1;
00060 y1 = y1<0? 0: y1>height-1? height-1: y1;
00061
00062 TU Value;
00063 int Count = 0;
00064
00065 Value = 0;
00066
00067 for (int y = y0; y <= y1; y++)
00068 {
00069 for (int x = x0; x <= x1; x++)
00070 {
00071 Value += map[y*width+x];
00072 Count++;
00073 }
00074 }
00075
00076 if (Count == 0)
00077 {
00078 Value = 0;
00079 return (T)Value;
00080 }
00081
00082 return (T)(Value/(float)Count);
00083 }
00084
00085 /** Does Bilinear Decimation.
00086 * \ingroup math */
00087 template <class T, class TU>
00088 inline T imBilinearDecimation(int width, int height, T *map, float xl, float yl, float box_width, float box_height, TU Dummy)
00089 {
00090 int x0,x1,y0,y1;
00091 (void)Dummy;
00092
00093 x0 = (int)floor(xl - box_width/2.0 - 0.5) + 1;
00094 y0 = (int)floor(yl - box_height/2.0 - 0.5) + 1;
00095 x1 = (int)floor(xl + box_width/2.0 - 0.5);
00096 y1 = (int)floor(yl + box_height/2.0 - 0.5);
00097
00098 if (x0 == x1) x1++;
00099 if (y0 == y1) y1++;
00100
00101 x0 = x0<0? 0: x0>width-1? width-1: x0;
00102 y0 = y0<0? 0: y0>height-1? height-1: y0;
00103 x1 = x1<0? 0: x1>width-1? width-1: x1;
00104 y1 = y1<0? 0: y1>height-1? height-1: y1;
00105
00106 TU Value, LineValue;
00107 float LineNorm, Norm, dxr, dyr;
00108
00109 Value = 0;
00110 Norm = 0;
00111
00112 for (int y = y0; y <= y1; y++)
00113 {
00114 dyr = yl - (y+0.5f);
00115 if (dyr < 0) dyr *= -1;
00116
00117 LineValue = 0;
00118 LineNorm = 0;
00119
00120 for (int x = x0; x <= x1; x++)
00121 {
00122 dxr = xl - (x+0.5f);
00123 if (dxr < 0) dxr *= -1;
00124
00125 LineValue += map[y*width+x] * dxr;
00126 LineNorm += dxr;
00127 }
00128
00129 Value += LineValue * dyr;
00130 Norm += dyr * LineNorm;
00131 }
00132
00133 if (Norm == 0)
00134 {
00135 Value = 0;
00136 return (T)Value;
00137 }
00138
00139 return (T)(Value/Norm);
00140 }
00141
00142 /** Does Zero Order Interpolation (Nearest Neighborhood).
00143 * \ingroup math */
00144 template <class T>
00145 inline T imZeroOrderInterpolation(int width, int height, T *map, float xl, float yl)
00146 {
00147 int x0 = (int)(xl-0.5f);
00148 int y0 = (int)(yl-0.5f);
00149 x0 = x0<0? 0: x0>width-1? width-1: x0;
00150 y0 = y0<0? 0: y0>height-1? height-1: y0;
00151 return map[y0*width + x0];
00152 }
00153
00154 /** Does Bilinear Interpolation.
00155 * \ingroup math */
00156 template <class T>
00157 inline T imBilinearInterpolation(int width, int height, T *map, float xl, float yl)
00158 {
00159 int x0, y0, x1, y1;
00160 float t, u;
00161
00162 if (xl < 0.5)
00163 {
00164 x1 = x0 = 0;
00165 t = 0;
00166 }
00167 else if (xl > width-0.5)
00168 {
00169 x1 = x0 = width-1;
00170 t = 0;
00171 }
00172 else
00173 {
00174 x0 = (int)(xl-0.5f);
00175 x1 = x0+1;
00176 t = xl - (x0+0.5f);
00177 }
00178
00179 if (yl < 0.5)
00180 {
00181 y1 = y0 = 0;
00182 u = 0;
00183 }
00184 else if (yl > height-0.5)
00185 {
00186 y1 = y0 = height-1;
00187 u = 0;
00188 }
00189 else
00190 {
00191 y0 = (int)(yl-0.5f);
00192 y1 = y0+1;
00193 u = yl - (y0+0.5f);
00194 }
00195
00196 T fll = map[y0*width + x0];
00197 T fhl = map[y0*width + x1];
00198 T flh = map[y1*width + x0];
00199 T fhh = map[y1*width + x1];
00200
00201 return (T)((fhh - flh - fhl + fll) * u * t +
00202 (fhl - fll) * t +
00203 (flh - fll) * u +
00204 fll);
00205 }
00206
00207 /** Does Bicubic Interpolation.
00208 * \ingroup math */
00209 template <class T, class TU>
00210 inline T imBicubicInterpolation(int width, int height, T *map, float xl, float yl, TU Dummy)
00211 {
00212 int x0, y0, x1, y1, x2, y2;
00213 float t, u;
00214 (void)Dummy;
00215
00216 if (xl < 0.5)
00217 {
00218 x2 = x1 = x0 = 0;
00219 t = 0;
00220 }
00221 else if (xl > width-0.5)
00222 {
00223 x2 = x1 = x0 = width-1;
00224 t = 0;
00225 }
00226 else
00227 {
00228 x0 = (int)(xl-0.5f);
00229 x1 = x0-1;
00230 if (x1 < 0) x1 = 0;
00231 x2 = x0+2; // 4 pixels (-1 0 +1 +2)
00232 if (x2 > width-1) x2 = width-1;
00233 t = xl - (x0+0.5f);
00234 }
00235
00236 if (yl < 0.5)
00237 {
00238 y2 = y1 = y0 = 0;
00239 u = 0;
00240 }
00241 else if (yl > height-0.5)
00242 {
00243 y2 = y1 = y0 = height-1;
00244 u = 0;
00245 }
00246 else
00247 {
00248 y0 = (int)(yl-0.5f);
00249 y1 = y0-1;
00250 if (y1 < 0) y1 = 0;
00251 y2 = y0+2; // 4 pixels (-1 0 +1 +2)
00252 if (y2 > height-1) y2 = height-1;
00253 u = yl - (y0+0.5f);
00254 }
00255
00256 float CX[4], CY[4];
00257
00258 // Optimize calculations
00259 {
00260 float x, x2, x3;
00261
00262 #define C0 (-x3 + 2.0f*x2 - x)
00263 #define C1 ( x3 - 2.0f*x2 + 1.0f)
00264 #define C2 (-x3 + x2 + x)
00265 #define C3 ( x3 - x2)
00266
00267 x = t;
00268 x2 = x*x; x3 = x2*x;
00269 CX[0] = C0; CX[1] = C1; CX[2] = C2; CX[3] = C3;
00270
00271 x = u;
00272 x2 = x*x; x3 = x2*x;
00273 CY[0] = C0; CY[1] = C1; CY[2] = C2; CY[3] = C3;
00274
00275 #undef C0
00276 #undef C1
00277 #undef C2
00278 #undef C3
00279 }
00280
00281 TU LineValue, Value;
00282 float LineNorm, Norm;
00283
00284 Value = 0;
00285 Norm = 0;
00286
00287 for (int y = y1; y <= y2; y++)
00288 {
00289 LineValue = 0;
00290 LineNorm = 0;
00291
00292 for (int x = x1; x <= x2; x++)
00293 {
00294 LineValue += map[y*width+x] * CX[x-x1];
00295 LineNorm += CX[x-x1];
00296 }
00297
00298 Value += LineValue * CY[y-y1];
00299 Norm += CY[y-y1] * LineNorm;
00300 }
00301
00302 if (Norm == 0)
00303 {
00304 Value = 0;
00305 return (T)Value;
00306 }
00307
00308 Value = (Value/Norm);
00309
00310 int size = sizeof(T);
00311 if (size == 1)
00312 return (T)(Value<=(TU)0? (TU)0: Value<=(TU)255? Value: (TU)255);
00313 else
00314 return (T)(Value);
00315 }
00316
00317 /** Calculates minimum and maximum values.
00318 * \ingroup math */
00319 template <class T>
00320 inline void imMinMax(const T *map, int count, T& min, T& max)
00321 {
00322 min = *map++;
00323 max = min;
00324 for (int i = 1; i < count; i++)
00325 {
00326 T value = *map++;
00327
00328 if (value > max)
00329 max = value;
00330 else if (value < min)
00331 min = value;
00332 }
00333 }
00334
00335 #endif