Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   Related Pages   Examples  

FloatKernel.cpp

00001 /* Copyright (c) 2001 C. Grigorescu */
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    /* initialization */
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    /* initialization */
00269   for (j=0; j<height; j++)
00270     for (i=0; i<width; i++)
00271       pixels[0][j*width+i] = 0.0;
00272   
00273   /* transform the phase in radians */
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 }