Kernels.C

00001 /*@file Image/Kernels.C Functions to construct various kinds of
00002   filter kernels
00003 */
00004 
00005 // //////////////////////////////////////////////////////////////////// //
00006 // The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 by the //
00007 // University of Southern California (USC) and the iLab at USC.         //
00008 // See http://iLab.usc.edu for information about this project.          //
00009 // //////////////////////////////////////////////////////////////////// //
00010 // Major portions of the iLab Neuromorphic Vision Toolkit are protected //
00011 // under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
00012 // in Visual Environments, and Applications'' by Christof Koch and      //
00013 // Laurent Itti, California Institute of Technology, 2001 (patent       //
00014 // pending; application number 09/912,225 filed July 23, 2001; see      //
00015 // http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
00016 
00017 
00018 
00019 // //////////////////////////////////////////////////////////////////// //
00020 // This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
00021 //                                                                      //
00022 // The iLab Neuromorphic Vision C++ Toolkit is free software; you can   //
00023 // redistribute it and/or modify it under the terms of the GNU General  //
00024 // Public License as published by the Free Software Foundation; either  //
00025 // version 2 of the License, or (at your option) any later version.     //
00026 //                                                                      //
00027 // The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope  //
00028 // that it will be useful, but WITHOUT ANY WARRANTY; without even the   //
00029 // implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR      //
00030 // PURPOSE.  See the GNU General Public License for more details.       //
00031 //                                                                      //
00032 // You should have received a copy of the GNU General Public License    //
00033 // along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
00034 // to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
00035 // Boston, MA 02111-1307 USA.                                           //
00036 // //////////////////////////////////////////////////////////////////// //
00037 //
00038 // Primary maintainer for this file: Rob Peters <rjpeters@klab.caltech.edu>
00039 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/Image/Kernels.C $
00040 // $Id: Kernels.C 13757 2010-08-04 22:11:39Z siagian $
00041 //
00042 
00043 #include "Image/Kernels.H"
00044 
00045 #include "Image/Image.H"
00046 #include "Image/MathOps.H"  // for mean()
00047 #include "Image/CutPaste.H" // for crop()
00048 #include "Util/Assert.H"
00049 #include "Util/MathFunctions.H"
00050 #include "rutz/compat_cmath.h"
00051 #include <cmath>
00052 
00053 
00054 // ######################################################################
00055 template <class T>
00056 Image<T> dogFilterHmax(const float theta, const float gamma, const int size, const float div)
00057 {
00058   // resize the data buffer
00059 
00060   Image<T> dest(size, size, NO_INIT);
00061 
00062   // change the angles in degree to the those in radian : rotation degree
00063   float rtDeg = M_PI / 180.0F * theta;
00064 
00065   // calculate constants
00066   float lambda = size*2.0F/div;
00067   float sigma = lambda*0.8F;
00068   float sigq = sigma*sigma;
00069   int center    = (int)ceil(size/2.0F);
00070   int filtSizeL = center-1;
00071   int filtSizeR = size-filtSizeL-1;
00072   typename Image<T>::iterator aptr = dest.beginw();
00073 
00074   // for DOG operation : to give orientation, it uses omit y-directional
00075   // component
00076   for (int x = -filtSizeL; x <= filtSizeR; x ++)
00077     for (int y = -filtSizeL; y <= filtSizeR; y ++) {
00078       if(sqrt(x*x +y*y)>size/2.0F){
00079         *aptr++ = 0;
00080       }
00081       else{
00082         float rtX =  y * cos(rtDeg) - x * sin(rtDeg);
00083         float rtY = y * sin(rtDeg) + x * cos(rtDeg);
00084         *aptr++ = T(exp(-(rtX*rtX + gamma*gamma*rtY*rtY)/(2.0F*sigq)) *
00085                     cos(2*M_PI*rtX/lambda));
00086       }
00087     }
00088   return dest;
00089 }
00090 
00091 // ######################################################################
00092 template <class T>
00093 Image<T> dogFilter(const float stddev, const float theta, const int halfsize)
00094 {
00095   // resize the data buffer
00096   int size = halfsize;
00097   if (size <= 0) size = int(ceil(stddev * sqrt(7.4F)));
00098   Image<T> dest(2 * size + 1, 2 * size + 1, NO_INIT);
00099 
00100   // change the angles in degree to the those in radian : rotation degree
00101   float rtDeg = M_PI / 180.0F * theta;
00102 
00103   // calculate constants
00104   float sigq = stddev * stddev;
00105   typename Image<T>::iterator aptr = dest.beginw();
00106 
00107   // for DOG operation : to give orientation, it uses omit y-directional
00108   // component
00109   for (int x = -size; x <= size; x ++)
00110     for (int y = -size; y <= size; y ++) {
00111       float rtX = x * cos(rtDeg) + y * sin(rtDeg);
00112       float rtY = -x * sin(rtDeg) + y * cos(rtDeg);
00113       *aptr++ = T((((rtX * rtX)) / sigq - 1.0F)/sigq *
00114                   exp(-(rtX * rtX + rtY * rtY) / (2.0F * sigq)));
00115     }
00116   return dest;
00117 }
00118 
00119 
00120 
00121 // ######################################################################
00122 template <class T>
00123 Image<T> dogFilter(const float stddev, const int halfsize)
00124 {
00125   float sFact = 4.0F;
00126 
00127   // resize the data buffer
00128   int size = halfsize;
00129   if (size <= 0) size = int(ceil(2*sFact*stddev));
00130   Image<T> dest(2 * size + 1, 2 * size + 1, NO_INIT);
00131 
00132   // calculate constants
00133   float NC = 1.0F/(2.0F*M_PI*stddev * stddev);
00134   float NS = 1.0F/(sFact*sFact) * NC;
00135   typename Image<T>::iterator aptr = dest.beginw();
00136 
00137   // for DOG operation : to give orientation, it uses omit y-directional
00138   // component
00139   for (int x = -size; x <= size; x ++)
00140     for (int y = -size; y <= size; y ++) {
00141       float xy = -1.0F * (x*x + y*y)/(2.0F*stddev*stddev);
00142       *aptr++ = T( NS*exp((1.0F/(sFact*sFact))*xy) - NC*exp(xy));
00143     }
00144   return dest;
00145 }
00146 
00147 // ######################################################################
00148 template <class T>
00149 Image<T> dogFilterHmax(const float stddev, const float theta,
00150                        const int cBegin, const int cEnd)
00151 {
00152   const Image<T> srcFilt = dogFilter<T>(stddev, theta);
00153   const Rectangle cropRegion =
00154     Rectangle::tlbrI(cBegin, cBegin, cEnd-1, cEnd-1);
00155   const Image<T> cropped = crop(srcFilt, cropRegion, false);
00156   const Image<T> diff = cropped - T(mean(cropped));
00157   return diff / T(sum(squared(diff)));
00158 }
00159 
00160 // ######################################################################
00161 template <class T>
00162 Image<T> gaborFilter(const float stddev, const float period,
00163                      const float phase, const float theta,
00164                      const float bg, const float ampl)
00165 {
00166   // figure the proper size for the result
00167   int size = int(ceil(stddev * sqrt(-2.0F * log(exp(-5.0F)))));
00168 
00169   Image<T> result(2 * size + 1, 2 * size + 1, NO_INIT);
00170 
00171   // change the angles in degree to the those in radians:
00172   float psi = M_PI / 180.0F * phase;
00173 
00174   float rtDeg = M_PI / 180.0F * theta;
00175 
00176   // calculate constants:
00177   float omega = (2.0F * M_PI) / period;
00178   float co = cos(rtDeg), si = sin(rtDeg);
00179   float sigq = 2.0F * stddev * stddev;
00180   typename Image<T>::iterator aptr = result.beginw();
00181   float a = float(ampl), b = float(bg);
00182 
00183   // compute gabor:
00184   for (int y = -size; y <= size; y ++)
00185     for (int x = -size; x <= size; x ++)
00186       *aptr++ = T(a * cos(omega * (x * co + y * si) + psi) *
00187                   exp(-(x * x + y * y) / sigq) + b);
00188   return result;
00189 }
00190 
00191 // ######################################################################
00192 template <class T>
00193 Image<T> gaborFilter(const float scale, const float theta)
00194 {
00195   // figure the proper size for the result
00196   int size = int(scale * 12+0.5f);
00197 
00198   Image<T> result(2 * size + 1, 2 * size + 1, NO_INIT);
00199 
00200   // calculate constants:
00201   float cosTheta = cos(theta + M_PI/2), sinTheta = sin(theta + M_PI/2);
00202   float normConst = 50/M_PI/scale/scale;
00203   typename Image<T>::iterator aptr = result.beginw();
00204 
00205   // compute gabor:
00206   for (int y = -size; y <= size; ++y)
00207     for (int x = -size; x <= size; ++x)
00208     {
00209       if (x*x+y*y > size*size)
00210         *aptr++ = T(); //zero out side the circle
00211       else
00212       {
00213         float xx = (x*cosTheta + y * sinTheta)/scale;
00214         float yy = (y*cosTheta - x * sinTheta)/scale;
00215         *aptr++ = T(exp(-(4.0F*xx*xx+yy*yy)/100.0F)/normConst);
00216       }
00217     }
00218 
00219 
00220   return result;
00221 }
00222 
00223 
00224 // ######################################################################
00225 //pixRGB version for colored gabor patches
00226 
00227 Image<PixRGB<byte> > gaborFilterRGB(const float stddev, const float freq,
00228                                     const float theta,const float hueShift)
00229 {
00230 
00231 
00232   // figure the proper size for the result
00233    int size = int(ceil(stddev * sqrt(-2.0F * log(exp(-5.0F)))));
00234 
00235   typedef std::complex<float> complexf;
00236 
00237   complexf gauss, sinu, gResult, ctemp;
00238   complexf i(0.0, 1.0);
00239 
00240   Image<PixHSV<float> > resultHSV(2 * size + 1, 2 * size + 1, NO_INIT);
00241 
00242 
00243   // change the angles in degree to the those in radians:
00244   float rtDeg = M_PI / 180.0F * theta;
00245   Image<PixHSV<float> >::iterator aptr = resultHSV.beginw();
00246 
00247 
00248   float tempf = 0.0;
00249   PixHSV<float> tempPix;
00250 
00251   int totalcnt = (2*size + 1)*(2*size + 1);
00252   std::vector<float> hVals(totalcnt), vVals(totalcnt);
00253 
00254   int cnt=0;
00255   for (int y = -size; y <= size; y ++)
00256    for (int x = -size; x <= size; x ++)
00257         {
00258           gauss   = std::complex<float>(exp(-(x * x + y * y) /(stddev*stddev)),0.0);
00259           ctemp   = std::complex<float>(0.0,(2.0*M_PI*freq)*(x*cos(rtDeg) + y*sin(rtDeg)));
00260           sinu    = std::exp(ctemp);
00261           gResult = gauss*sinu;
00262 
00263           tempf = arg(gResult);
00264 
00265           tempf = tempf * 180/M_PI; //convert back to degs
00266           hVals[cnt]   = tempf;
00267           vVals[cnt++] = real(gResult)*real(gResult);
00268 
00269            *aptr++ = tempPix;
00270 
00271           }
00272 
00273 
00274     //find min max of hues and values and rescale h from 0-360 and v from 0-255
00275   cnt = 0;
00276   float minH = hVals[0], maxH = hVals[0], minV = vVals[0], maxV = vVals[0];
00277   while(cnt!=totalcnt)
00278     {
00279       if(hVals[cnt] < minH)
00280         minH = hVals[cnt];
00281       else if(hVals[cnt] > maxH)
00282         maxH = hVals[cnt];
00283 
00284 
00285       if(vVals[cnt] < minV)
00286         minV = vVals[cnt];
00287       else if(vVals[cnt] > maxV)
00288         maxV = vVals[cnt];
00289 
00290 
00291       cnt++;
00292     }
00293 
00294   aptr = resultHSV.beginw();
00295 
00296   for(int i = 0; i <= totalcnt ; i++)
00297     {
00298       hVals[i] = ((hVals[i] - minH)/(maxH - minH)) * 360;
00299       hVals[i] = fabs(hVals[i] + hueShift);
00300 
00301       if(hVals[i] > 360)
00302         hVals[i] = hVals[i] - 360;
00303 
00304       tempPix.setH(hVals[i]);
00305       tempPix.setS(100.0) ;
00306       tempPix.setV(vVals[i]);
00307       *aptr++ = tempPix;
00308     }
00309 
00310   LINFO("stddev:%f freq:%f theta:%f hueShift:%f",stddev,freq,theta,hueShift);
00311 
00312   Image<PixRGB<float> > resultfloat = resultHSV;
00313   Image<PixRGB<byte> > result = resultfloat*256;
00314 
00315   return result;
00316 
00317 }
00318 
00319 // ######################################################################
00320 template <class T>
00321 Image<T> gaborFilter2(const float stddev, const float period,
00322                       const float phase, const float theta,
00323                       const float sigMod = 1.0F,
00324                       const float amplitude = 128.0F)
00325 {
00326   // figure the proper size for the result
00327   int size = int(ceil(stddev * sqrt(-2.0F * log(exp(-5.0F)))));
00328 
00329   Image<T> result(2 * size + 1, 2 * size + 1, NO_INIT);
00330 
00331   // change the angles in degree to the those in radians:
00332   float psi = M_PI / 180.0F * phase;
00333   float rtDeg = M_PI / 180.0F * theta;
00334 
00335   // calculate constants:
00336   float omega = (2.0F * M_PI) / period;
00337   float co = cos(rtDeg), si = sin(rtDeg);
00338   float sigq = sigMod * stddev * stddev;
00339   typename Image<T>::iterator aptr = result.beginw();
00340 
00341   // compute gabor:
00342   for (int y = -size; y <= size; y ++)
00343     for (int x = -size; x <= size; x ++)
00344       *aptr++ = T(amplitude*(cos(omega * (x * co + y * si) + psi) *
00345                              exp(-(x * x + y * y) / sigq)));
00346   // who's your daddy?
00347   return result;
00348 }
00349 
00350 // ######################################################################
00351 // Produces a Gabor kernel with optionally unequal major+minor axis lengths.
00352 Image<float> gaborFilter3(const float major_stddev, const float minor_stddev,
00353                           const float period, const float phase,
00354                           const float theta, int size)
00355 {
00356   const float max_stddev =
00357     major_stddev > minor_stddev ? major_stddev : minor_stddev;
00358 
00359   // figure the proper size for the result
00360   if (size == -1) size = int(ceil(max_stddev * sqrt(-2.0F * log(exp(-5.0F)))));
00361   else size = size/2;
00362   //const int size = int(ceil(max_stddev * 2));
00363 
00364   Image<float> result(2 * size + 1, 2 * size + 1, NO_INIT);
00365 
00366   // change the angles in degree to the those in radians:
00367   const float psi = M_PI / 180.0F * phase;
00368   const float rtDeg = M_PI / 180.0F * theta;
00369 
00370   // calculate constants:
00371   const float omega = (2.0F * M_PI) / period;
00372   const float co = cos(rtDeg), si = sin(rtDeg);
00373   const float major_sigq = 2.0F * major_stddev * major_stddev;
00374   const float minor_sigq = 2.0F * minor_stddev * minor_stddev;
00375   Image<float>::iterator aptr = result.beginw();
00376 
00377   // compute gabor:
00378   for (int y = -size; y <= size; ++y)
00379     for (int x = -size; x <= size; ++x)
00380       {
00381         const float major = x*co + y*si;
00382         const float minor = x*si - y*co;
00383         *aptr++ = float(cos(omega * major + psi)
00384                         * exp(-(major*major) / major_sigq)
00385                         * exp(-(minor*minor) / minor_sigq));
00386       }
00387 
00388   // make the result have (mean == 0.0)
00389   return (result - mean(result));
00390 }
00391 
00392 // ######################################################################
00393 
00394 template <class T>
00395 Image<T> gaussian2D(const float stddev, const float sigMod = 1.0F,
00396                     const float amplitude = 255.0F)
00397 {
00398   // figure the proper size for the result
00399   int size = int(ceil(stddev * sqrt(-2.0F * log(exp(-5.0F)))));
00400   //float norm = 40.010587F;
00401   Image<T> result(2 * size + 1, 2 * size + 1, NO_INIT);
00402   T resultX,resultY;
00403 
00404   typename Image<T>::iterator aptr = result.beginw();
00405   for (int y = -size; y <= size; y++)
00406     {
00407       resultY = T(/*((1.0F/(stddev*sqrt(2.0F*M_PI)))*/1.0F*
00408                   exp(-((y*y)/(stddev*stddev*sigMod))));
00409       for (int x = -size; x <= size; x++)
00410         {
00411           resultX = T(/*((1.0F/(stddev*sqrt(2.0F*M_PI)))*/1.0F*
00412                       exp(-((x*x)/(stddev*stddev*sigMod))));
00413           *aptr++ = T((resultX*resultY)*amplitude);
00414         }
00415     }
00416   return result;
00417 }
00418 
00419 // ######################################################################
00420 template <class T>
00421 Image<T> gaussianBlob(const Dims& dims, const Point2D<int>& center,
00422                       const float sigmaX, const float sigmaY)
00423 {
00424   Image<T> ret(dims, NO_INIT);
00425   const int w = dims.w(), h = dims.h();
00426 
00427   const float fac = 0.5f / (M_PI * sigmaX * sigmaY);
00428   const float vx = -0.5f / (sigmaX * sigmaX);
00429   const float vy = -0.5f / (sigmaY * sigmaY);
00430   for (int jj = 0; jj < h; jj ++)
00431     {
00432       float vydy2 = float(jj - center.j); vydy2 *= vydy2 * vy;
00433       for (int ii = 0; ii < w; ii ++)
00434         {
00435           float dx2 = float(ii - center.i); dx2 *= dx2;
00436 
00437           ret.setVal(ii, jj, clamped_convert<T>(fac * expf(vx * dx2 + vydy2)));
00438         }
00439     }
00440   return ret;
00441 }
00442 
00443 // ######################################################################
00444 template <class T>
00445 Image<T> gaussianBlobUnnormalized(const Dims& dims, const Point2D<int>& center,
00446                                   const float sigmaX, const float sigmaY)
00447 {
00448   Image<T> ret(dims, NO_INIT);
00449   const int w = dims.w(), h = dims.h();
00450 
00451   const float vx = -0.5f / (sigmaX * sigmaX);
00452   const float vy = -0.5f / (sigmaY * sigmaY);
00453   for (int jj = 0; jj < h; jj ++)
00454     {
00455       float vydy2 = float(jj - center.j); vydy2 *= vydy2 * vy;
00456       for (int ii = 0; ii < w; ii ++)
00457         {
00458           float dx2 = float(ii - center.i); dx2 *= dx2;
00459 
00460           ret.setVal(ii, jj, clamped_convert<T>(expf(vx * dx2 + vydy2)));
00461         }
00462     }
00463   return ret;
00464 }
00465 
00466 // ######################################################################
00467 namespace
00468 {
00469   // Helper function for binomialKernel():
00470   double binomial(int n, int k)
00471   {
00472     /*
00473       binomial coefficient refresher:
00474 
00475       k:    0   1   2   3   4   5   6   7   8
00476       N:  +------------------------------------
00477       0   |   1
00478       1   |   1   1
00479       2   |   1   2   1
00480       3   |   1   3   3   1
00481       4   |   1   4   6   4   1
00482       5   |   1   5  10  10   5   1
00483       6   |   1   6  15  20  15   6   1
00484       7   |   1   7  21  35  35  21   7   1
00485       8   |   1   8  28  56  70  56  28   8   1
00486 
00487       n-choose-k := n! / k! (n-k)!
00488 
00489       == n*(n-1)*...(n-(k-1)) / k*(k-1)*...*1
00490 
00491       gamma(n+1) == n! for integral n
00492 
00493       lgamma(x) := ln(gamma(x))
00494 
00495       By doing the factorial multiplications as additions in the
00496       log-domain, we avoid numerical overflow in intermediate results.
00497     */
00498     return floor(0.5 + exp(lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1)));
00499   }
00500 }
00501 
00502 Image<float> binomialKernel(const int sz)
00503 {
00504   ASSERT(sz > 0);
00505 
00506   Image<float> result(sz, 1, NO_INIT);
00507 
00508   const int N = sz-1;
00509 
00510   const double div = pow(2.0, double(N));
00511 
00512   // compute a series of N-choose-k, where N == (sz-1), with k running
00513   // from 0 to N, inclusive
00514   for (int k = 0; k <= N; ++k)
00515     {
00516       result.setVal(k, 0, float(binomial(N, k) / div));
00517     }
00518 
00519   return result;
00520 }
00521 
00522 
00523 // ######################################################################
00524 template <class T>
00525 Image<T> grating(const int width, const int height,
00526                  const float period, const float phase, const float theta)
00527 {
00528   Image<T> result(width, height, NO_INIT);
00529 
00530   // change the angles in degree to the those in radian
00531   float psi = M_PI / 180.0F * phase;
00532   float rtDeg = M_PI / 180.0F * theta;
00533 
00534   // calculate constants
00535   float omega = (2.0F * M_PI) / period;
00536   float co = cos(rtDeg), si = sin(rtDeg);
00537   typename Image<T>::iterator aptr = result.beginw();
00538 
00539   // for gabor operation
00540   for (int y = 0; y < height; y ++)
00541     for (int x = 0; x < width; x ++)
00542       *aptr++ = T(cos(omega * (x * co + y * si) + psi));
00543   return result;
00544 }
00545 
00546 // ######################################################################
00547 template <class T>
00548 Image<T> gaussian(const float coeff, const float sigma,
00549                   const int maxhw, const float threshperc)
00550 {
00551   // determine size: keep only values larger that threshperc*max (here max=1)
00552   int hw = (int)(sigma * sqrt(-2.0F * log(threshperc / 100.0F)));
00553 
00554   // if kernel turns out to be too large, cut it off:
00555   if (maxhw > 0 && hw > maxhw) hw = maxhw;
00556 
00557   // allocate image for result:
00558   Image<T> result(2 * hw + 1, 1, NO_INIT);
00559 
00560   // if coeff is given as 0, compute it from sigma:
00561   float c = coeff;
00562   if (coeff == 0.0f) c = 1.0f / (sigma * sqrtf(2.0f * float(M_PI)));
00563 
00564   // build both halves simultaneously
00565   result.setVal(hw, T(c));
00566   const float sig22 = - 0.5F / (sigma * sigma);
00567   for (int i = 1; i <= hw; ++i)
00568     {
00569       T val = T(c * exp(float(i * i) * sig22));
00570       result.setVal(hw + i, val);
00571       result.setVal(hw - i, val);
00572     }
00573   return result;
00574 }
00575 
00576 // ######################################################################
00577 template <class T>
00578 Image<T> longRangeExcFilter(const float factor, const float orient)
00579 {
00580   int siz = 9;         // kernel size
00581   float peak = 3.0F;   // [3] excitation is max 'peak' pixels away from center
00582   float sigmaR = 1.0F; // [1] radial sigma, in pixels
00583   float sigmaT = 5.0F; // [15] angular sigma, in degrees
00584 
00585   Image<T> result(siz, siz, NO_INIT); int mid = (siz - 1) / 2;
00586   float o = orient * M_PI / 180.0F, st = sigmaT * M_PI / 180.0F;
00587   st = 2.0F * st * st; float sr = 2.0F * sigmaR * sigmaR;
00588 
00589   typename Image<T>::iterator aptr = result.beginw();
00590   for (int j = 0; j < siz; j ++)
00591     for (int i = 0; i < siz; i ++)
00592       {
00593         float x = (float)(i - mid), y = (float)(j - mid);
00594         float r = sqrt(x * x + y * y) - peak, t = atan2(y, x);
00595         float odiff = t - o;
00596         while (odiff < -M_PI_2) odiff += M_PI;
00597         while (odiff > M_PI_2) odiff -= M_PI;
00598         *aptr++ = T(factor * exp(- r*r/sr - odiff*odiff/st));
00599       }
00600   return result;
00601 }
00602 
00603 // ######################################################################
00604 template <class T>
00605 Image<T> fixationMask(const Dims& dims, const Point2D<int>& fixation,
00606                       const float pixperdeg, const T maxval, const float sigma)
00607 {
00608   Image<T> img(dims, ZEROS);
00609   return fixationMask(img, fixation, pixperdeg, maxval, sigma);
00610 }
00611 
00612 // ######################################################################
00613 template <class T>
00614 Image<T> fixationMask(const Image<T>& mask, const Point2D<int>& fixation,
00615                       const float pixperdeg, const T maxval, const float sigma)
00616 {
00617   Image<T> img = mask;
00618   typename Image<T>::iterator src = img.beginw();
00619   float mvf = float(maxval); int w = img.getWidth(), h = img.getHeight();
00620   float sigfac = 0.5f / (sigma * sigma);
00621   int fi = fixation.i, fj = fixation.j;
00622   for (int j = 0; j < h; j ++)
00623     for (int i = 0; i < w; i ++)
00624       {
00625         float ang = sqrt((i-fi) * (i-fi) + (j-fj) * (j-fj)) / pixperdeg;
00626         float val = exp(- ang * ang * sigfac);
00627         *src++ += T(mvf * val);  // range 0..maxval
00628       }
00629   return img;
00630 }
00631 
00632 // ######################################################################
00633 Image<byte> twofiftyfives(Dims d)
00634 {
00635   Image<byte> result(d,NO_INIT);
00636   result.clear(255);
00637   return result;
00638 }
00639 
00640 // ######################################################################
00641 Image<byte> twofiftyfives(int w, int h)
00642 {
00643   return twofiftyfives(Dims(w,h));
00644 }
00645 
00646 // ######################################################################
00647 Image<byte> twofiftyfives(int w)
00648 {
00649   return twofiftyfives(Dims(w,w));
00650 }
00651 
00652 // Include the explicit instantiations
00653 #include "inst/Image/Kernels.I"
Generated on Sun May 8 08:40:55 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3