00001
00002
00003 #include <math.h>
00004 #include <fstream.h>
00005 #include "imlib_io.h"
00006 #include "FloatKernel.h"
00007
00008 void FloatKernel::makeRectangularPulse(float max_amplitude, int support)
00009 {
00010 int i, j, pos_x, pos_y, x_coord, y_coord, index;
00011 int env_x, env_y;
00012
00013 kernel_type = RECTANGULAR_PULSE;
00014 for (j=0; j<height; j++)
00015 for (i=0; i<width; i++)
00016 pixels[0][j*width+i] = 0.0;
00017
00018 pos_x = (int)(width/2+width%2);
00019 pos_y = (int)(height/2+height%2);
00020
00021 if (support > width/2) support = width/2;
00022 if (support > height/2) support = height/2;
00023
00024 env_x = (int)(support/2 + support%2);
00025 env_y = (int)(support/2 + support%2);
00026 for (j=-env_y;j<env_y;j++)
00027 for (i=-env_x;i<env_x;i++) {
00028 x_coord = pos_x + i;
00029 y_coord = pos_y + j;
00030 if (y_coord >= height)
00031 y_coord = pos_y + j - height;
00032 if (y_coord < 0)
00033 y_coord = pos_y + j + height;
00034 index = (int)(y_coord*width + x_coord);
00035 if ((index >=0) && (index < width*height))
00036 if ((i == -env_x) || (i == env_x-1) || (j == -env_y) || (j == env_y-1))
00037 pixels[0][index] = max_amplitude/2;
00038 else
00039 pixels[0][index] = max_amplitude;
00040 }
00041 }
00042
00043 void FloatKernel::makeTriangularPulse(float max_amplitude, int support)
00044 {
00045 int i, j, pos_x, pos_y, x_coord, y_coord, index;
00046 int env_x, env_y;
00047 float val_x, val_y;
00048
00049 kernel_type = TRIANGULAR_PULSE;
00050 for (j=0; j<height; j++)
00051 for (i=0; i<width; i++)
00052 pixels[0][j*width+i] = 0.0;
00053
00054 pos_x = (int)(width/2+width%2);
00055 pos_y = (int)(height/2+height%2);
00056
00057 if (support > width/2) support = width/2;
00058 if (support > height/2) support = height/2;
00059
00060 env_x = (int)(support/2 + support%2);
00061 env_y = (int)(support/2 + support%2);
00062 for (j=-env_y;j<=0;j++)
00063 for (i=-env_x;i<=0;i++) {
00064 x_coord = pos_x + i - 1;
00065 y_coord = pos_y + j - 1;
00066 val_x = (float)(env_x+i)/(float)(env_x);
00067 val_y = (float)(env_y+j)/(float)(env_y);
00068 if (y_coord >= height)
00069 y_coord = pos_y + j - height;
00070 if (y_coord < 0)
00071 y_coord = pos_y + j + height;
00072 index = (int)(y_coord*width + x_coord);
00073 if ((index >=0) && (index < width*height)) {
00074 pixels[0][index] = max_amplitude * MIN(val_x,val_y);
00075 }
00076 }
00077 for (j=0;j<env_y;j++)
00078 for (i=-env_x;i<=0;i++) {
00079 x_coord = pos_x + i - 1;
00080 y_coord = pos_y + j;
00081 val_x = (float)(env_x+i)/(float)(env_x);
00082 val_y = (float)(-env_y+j)/(float)(-env_y);
00083 if (y_coord >= height)
00084 y_coord = pos_y + j - height;
00085 if (y_coord < 0)
00086 y_coord = pos_y + j + height;
00087 index = (int)(y_coord*width + x_coord);
00088 if ((index >=0) && (index < width*height)) {
00089 pixels[0][index] = max_amplitude * MIN(val_x,val_y);
00090 }
00091 }
00092 for (j=-env_y;j<=0;j++)
00093 for (i=0;i<env_x;i++) {
00094 x_coord = pos_x + i;
00095 y_coord = pos_y + j - 1;
00096 val_x = (float)(-env_x+i)/(float)(-env_x);
00097 val_y = (float)(env_y+j)/(float)(env_y);
00098 if (y_coord >= height)
00099 y_coord = pos_y + j - height;
00100 if (y_coord < 0)
00101 y_coord = pos_y + j + height;
00102 index = (int)(y_coord*width + x_coord);
00103 if ((index >=0) && (index < width*height)) {
00104 pixels[0][index] = max_amplitude * MIN(val_x,val_y);
00105 }
00106 }
00107 for (j=0;j<env_y;j++)
00108 for (i=0;i<env_x;i++) {
00109 x_coord = pos_x + i;
00110 y_coord = pos_y + j;
00111 val_x = (float)(-env_x+i)/(float)(-env_x);
00112 val_y = (float)(-env_y+j)/(float)(-env_y);
00113 if (y_coord >= height)
00114 y_coord = pos_y + j - height;
00115 if (y_coord < 0)
00116 y_coord = pos_y + j + height;
00117 index = (int)(y_coord*width + x_coord);
00118 if ((index >=0) && (index < width*height)) {
00119 pixels[0][index] = max_amplitude * MIN(val_x,val_y);
00120 }
00121 }
00122 }
00123
00124 void FloatKernel::makeGauss(float sigma)
00125 {
00126 kernel_type = GAUSS;
00127 makeGauss(sigma,1.0,0.0);
00128 }
00129
00130 void FloatKernel::makeGauss(float sigma, float aspect, float orientation)
00131 {
00132 kernel_type = GAUSS;
00133 int i, j, index, x_coord, y_coord;
00134 int pos_x, pos_y, env_x, env_y;
00135 float co;
00136 float cth, sth, x1, y1;
00137 float gauss;
00138
00139 find_limits(sigma, aspect, orientation);
00140
00141
00142 for (j=0; j<height; j++)
00143 for (i=0; i<width; i++)
00144 pixels[0][j*width+i] = 0.0;
00145
00146 env_x = (int)(h_m_x);
00147 env_y = (int)(h_m_y);
00148 pos_x = (int)(width/2+width%2);
00149 pos_y = (int)(height/2);
00150
00151 cth = cos(PI*orientation/180.0);
00152 sth = sin(PI*orientation/180.0);
00153 co = 1/(sqrt(2*PI)*sigma);
00154
00155 for (j=0; j<height; j++)
00156 for (i=0; i<width; i++)
00157 pixels[0][j*width+i] = 0.0;
00158
00159 for (j=-env_y;j<env_y;j++)
00160 for (i=-env_x;i<env_x;i++) {
00161 x1 = i*cth - j*sth;
00162 y1 = j*cth + i*sth;
00163 gauss = co*exp((-(x1*x1)-(aspect*y1)*(aspect*y1))/(2*sigma*sigma));
00164 x_coord = pos_x + i;
00165 y_coord = pos_y + j;
00166 if (y_coord >= height)
00167 y_coord = pos_y + j - height;
00168 if (y_coord < 0)
00169 y_coord = pos_y + j + height;
00170 index = (int)(y_coord*width + x_coord);
00171 if ((index >=0) && (index < width*height))
00172 pixels[0][index] = gauss;
00173 }
00174 }
00175
00176 void FloatKernel::makeDoG(float sigma, float ratio)
00177 {
00178 kernel_type = DOG;
00179 makeDoG(sigma, 1.0, 0.0, ratio);
00180 }
00181
00182 void FloatKernel::makeDoG(float sigma, float aspect, float orientation, float ratio)
00183
00184 {
00185 int i, j, env_x, env_y, x_coord, y_coord, index;
00186 int pos_x, pos_y;
00187 float l1_norm, l2_norm, gauss1, gauss2, gauss;
00188 float co1, co2, cth, sth, x1, y1;
00189
00190 kernel_type = DOG;
00191 find_limits(sigma, aspect, orientation);
00192
00193 env_x = (int)(h_m_x);
00194 env_y = (int)(h_m_y);
00195 pos_x = (int)(width/2+width%2);
00196 pos_y = (int)(height/2);
00197
00198 for (j=0; j<height; j++)
00199 for (i=0; i<width; i++)
00200 pixels[0][j*width+i] = 0.0;
00201
00202 l1_norm = 0;
00203 l2_norm = 0;
00204
00205 cth = cos(PI*orientation/180.0);
00206 sth = sin(PI*orientation/180.0);
00207 co1 = 1/(sqrt(2*PI)*sigma);
00208 co2 = 1/(sqrt(2*PI)*ratio*sigma);
00209
00210 for (j=-env_y;j<env_y;j++)
00211 for (i=-env_x;i<env_x;i++) {
00212 x1 = i*cth - j*sth;
00213 y1 = j*cth + i*sth;
00214 gauss1 = co1 * exp((-(x1*x1)-(aspect*y1)*(aspect*y1))/(2*sigma*sigma));
00215 gauss2 = co2 * exp((-(x1*x1)-(aspect*y1)*(aspect*y1))/
00216 (2*ratio*ratio*sigma*sigma));
00217 gauss = gauss1 - gauss2;
00218 x_coord = pos_x + i;
00219 y_coord = pos_y + j;
00220 if (y_coord >= height)
00221 y_coord = pos_y + j - height;
00222 if (y_coord < 0)
00223 y_coord = pos_y + j + height;
00224 index = (int)(y_coord*width + x_coord);
00225 if ((index >=0) && (index < width*height)) {
00226 pixels[0][index] = gauss;
00227 l1_norm += gauss;
00228 l2_norm += gauss*gauss;
00229 }
00230 }
00231 l1_norm = l1_norm/(width*height);
00232 l2_norm = sqrt(l2_norm);
00233
00234 for (j=-env_y;j<env_y;j++)
00235 for (i=-env_x;i<env_x;i++) {
00236 index = (int)((pos_y+j)*width + pos_x+i);
00237 if ((index >=0) && (index < width*height)) {
00238 pixels[0][index] = pixels[0][index]-l1_norm;
00239 pixels[0][index] = pixels[0][index]/l2_norm;
00240 }
00241 }
00242 }
00243
00244 void FloatKernel::makeGabor(float sigma, float aspect, float orientation,
00245 float lambda, float phase)
00246 {
00247 int i, j, env_x, env_y, index, x_coord, y_coord;
00248 float cth, sth, ph, x1, y1;
00249 int pos_x, pos_y;
00250 float gauss;
00251 float l1_norm, l2_norm, co;
00252
00253 kernel_type = GABOR;
00254 find_limits(sigma, aspect, orientation);
00255
00256 env_x = (int)(h_m_x);
00257 env_y = (int)(h_m_y);
00258 pos_x = (int)(width/2+width%2);
00259 pos_y = (int)(height/2);
00260
00261 l1_norm = 0.0;
00262 l2_norm = 0.0;
00263
00264 cth = cos(PI*orientation/180.0);
00265 sth = sin(PI*orientation/180.0);
00266 co = 1/(sqrt(2*PI)*sigma);
00267
00268
00269 for (j=0; j<height; j++)
00270 for (i=0; i<width; i++)
00271 pixels[0][j*width+i] = 0.0;
00272
00273
00274 ph = PI*phase/180.0;
00275
00276 for (j=-env_y;j<env_y;j++)
00277 for (i=-env_x;i<env_x;i++) {
00278 x1 = i*cth - j*sth;
00279 y1 = j*cth + i*sth;
00280 gauss = co*exp((-(x1*x1)-(aspect*y1)*(aspect*y1))/(2*sigma*sigma));
00281 x_coord = pos_x + i;
00282 y_coord = pos_y + j;
00283 if (y_coord >= height)
00284 y_coord = pos_y + j - height;
00285 if (y_coord < 0)
00286 y_coord = pos_y + j + height;
00287 index = (int)(y_coord*width + x_coord);
00288 if ((index >=0) && (index < width*height)) {
00289 pixels[0][index] = gauss * cos(2*PI*x1/lambda + ph);
00290 l1_norm = l1_norm + pixels[0][index];
00291 l2_norm = l2_norm + pixels[0][index]*pixels[0][index];
00292 }
00293 }
00294 l1_norm = l1_norm/(width*height);
00295 l2_norm = sqrt(l2_norm);
00296
00297 for (j=-env_y;j<env_y;j++)
00298 for (i=-env_x;i<env_x;i++) {
00299 index = (int)((pos_y+j)*width + pos_x+i);
00300 if ((index >=0) && (index < width*height)) {
00301 pixels[0][index] = pixels[0][index]-l1_norm;
00302 pixels[0][index] = pixels[0][index]/l2_norm;
00303 }
00304 }
00305 }
00306
00307 void FloatKernel::makeLoG(float sigma)
00308 {
00309 kernel_type = LOG;
00310 makeLoG(sigma,1.0,0.0);
00311 }
00312
00313 void FloatKernel::makeLoG(float sigma, float aspect, float orientation)
00314 {
00315 int i, j, env_x, env_y, index, x_coord, y_coord;
00316 float cth, sth, x1, y1;
00317 int pos_x, pos_y;
00318 float gauss;
00319 float l1_norm, l2_norm, co;
00320
00321 kernel_type = LOG;
00322 find_limits(sigma, aspect, orientation);
00323
00324 env_x = (int)(h_m_x);
00325 env_y = (int)(h_m_y);
00326 pos_x = (int)(width/2+width%2);
00327 pos_y = (int)(height/2);
00328 cth = cos(PI*orientation/180.0);
00329 sth = sin(PI*orientation/180.0);
00330
00331 l1_norm = 0.0;
00332 l2_norm = 0.0;
00333
00334 for (j=-env_y;j<env_y;j++)
00335 for (i=-env_x;i<env_x;i++) {
00336 x1 = i*cth - j*sth;
00337 y1 = j*cth + i*sth;
00338 co = 1-(x1*x1+y1*y1)/(2*sigma*sigma);
00339 gauss = co*exp((-(x1*x1)-(aspect*y1)*(aspect*y1))/(2*sigma*sigma));
00340 x_coord = pos_x + i;
00341 y_coord = pos_y + j;
00342 if (y_coord >= height)
00343 y_coord = pos_y + j - height;
00344 if (y_coord < 0)
00345 y_coord = pos_y + j + height;
00346 index = (int)(y_coord*width + x_coord);
00347 if ((index >=0) && (index < width*height)) {
00348 pixels[0][index] = gauss;
00349 l1_norm = l1_norm + pixels[0][index];
00350 l2_norm = l2_norm + pixels[0][index]*pixels[0][index];
00351 }
00352 }
00353 l1_norm = l1_norm/(width*height);
00354 l2_norm = sqrt(l2_norm);
00355
00356 for (j=-env_y;j<env_y;j++)
00357 for (i=-env_x;i<env_x;i++) {
00358 index = (int)((pos_y+j)*width + pos_x+i);
00359 if ((index >=0) && (index < width*height)) {
00360 pixels[0][index] = pixels[0][index]-l1_norm;
00361 pixels[0][index] = pixels[0][index]/l2_norm;
00362 }
00363 }
00364 }
00365
00366 void FloatKernel::find_limits(float sigma, float aspect, float orientation)
00367 {
00368 float proj_1_x, proj_2_x, proj_1_y, proj_2_y;
00369
00370 h_m_x = 7 * sigma;
00371 h_m_y = 7 * (h_m_x)/aspect;
00372
00373 if (orientation <= 90)
00374 {
00375 proj_1_x = h_m_x * cos(PI*orientation/180.0);
00376 proj_2_y = h_m_y * cos(PI*orientation/180.0);
00377 }
00378 else
00379 {
00380 proj_1_x = h_m_x * -cos(PI*orientation/180.0);
00381 proj_2_y = h_m_y * -cos(PI*orientation/180.0);
00382 }
00383
00384 proj_2_x = h_m_x * sin(PI*orientation/180.0);
00385 proj_1_y = h_m_y * sin(PI*orientation/180.0);
00386
00387 if (proj_1_x > proj_1_y)
00388 h_m_x = proj_1_x;
00389 else
00390 h_m_x = proj_1_y;
00391 if (proj_2_x > proj_2_y)
00392 h_m_y = proj_2_x;
00393 else
00394 h_m_y = proj_2_y;
00395
00396 if (h_m_x >= width/2)
00397 h_m_x = width/2+width%2;
00398 if (h_m_y >= height/2)
00399 h_m_y = height/2+height%2;
00400 }