Filters.cpp

Go to the documentation of this file.
00001 
00004 // This file is part of the SaliencyToolbox - Copyright (C) 2006-2007
00005 // by Dirk B. Walther and the California Institute of Technology.
00006 // See the enclosed LICENSE.TXT document for the license agreement. 
00007 // More information about this project is available at: 
00008 // http://www.saliencytoolbox.net
00009 
00010 #include "mexLog.h"
00011 #include "Image.h"
00012 #include "Filters.h"
00013 
00014 #include <algorithm>
00015 #include <cmath>
00016 #include <limits>
00017 
00018 // ######################################################################
00019 // kernel: 1 5 10 10 5 1
00020 Image lowPass6yDecY(const Image& src)
00021 {
00022   ASSERT(ecxStr);
00023   const int w = src.getWidth(), hs = src.getHeight();
00024   const float ecw = ecx / w;
00025   int hr = hs / 2;
00026   if (hr == 0) hr = 1;
00027   
00028   Image result(w,hr);
00029   Image::iterator rptr = result.beginw();
00030   Image::const_iterator sptr = src.begin();
00031   
00032   if (hs <= 1)
00033     result = src;
00034   else if (hs == 2)
00035     for (int x = 0; x < w; ++x)
00036       {
00037         // use kernel [1 1]^T / 2
00038         *rptr++ = (sptr[0] + sptr[1]) / 2.0;
00039         sptr += 2;
00040       }
00041   else if (hs == 3)
00042     for (int x = 0; x < w; ++x)
00043       {
00044         // use kernel [1 2 1]^T / 4
00045         *rptr++ = (sptr[0] + sptr[1] * 2.0 + sptr[2]) / 4.0;
00046         sptr += 3;
00047       }
00048   else // general case with hs >= 4
00049     for (int x = 0; x < w; ++x)
00050       {
00051         // top most point - use kernel [10 10 5 1]^T / 26
00052         *rptr++ = ((sptr[0] + sptr[1]) * 10.0 + 
00053                     sptr[2] * 5.0 + sptr[3]) / 26.0;
00054         //++sptr;
00055         
00056         // general case
00057         int y;
00058         for (y = 0; y < (hs - 5); y += 2)
00059           {
00060             // use kernel [1 5 10 10 5 1]^T / 32
00061             *rptr++ = ((sptr[1] + sptr[4])  *  5.0 +
00062                        (sptr[2] + sptr[3])  * 10.0 +
00063                        (sptr[0] + sptr[5])) / 32.0;
00064             sptr += 2;
00065           }
00066         
00067         // find out how to treat the bottom most point
00068         if (y == (hs - 5))
00069           {
00070             // use kernel [1 5 10 10 5]^T / 31
00071              *rptr++ = ((sptr[1] + sptr[4])  *  5.0 +
00072                         (sptr[2] + sptr[3])  * 10.0 +
00073                          sptr[0])            / 31.0;
00074             sptr += 5;
00075           }
00076         else
00077           {
00078             // use kernel [1 5 10 10]^T / 26
00079             *rptr++ = ( sptr[0] + sptr[1]  *  5.0 +
00080                        (sptr[2] + sptr[3]) * 10.0) / 26.0;
00081             sptr += 4;
00082           }
00083         }
00084             
00085    return result;      
00086 }
00087 
00088 // ######################################################################
00089 // kernel: 1 5 10 10 5 1
00090 Image lowPass6xDecX(const Image& src)
00091 {
00092   ASSERT(ecxStr);
00093   const int ws = src.getWidth(), h = src.getHeight();
00094   const float ecw = ecx / ws;
00095   const int h2 = h * 2, h3 = h * 3, h4 = h * 4, h5 = h * 5;
00096   int wr = ws / 2;
00097   if (wr == 0) wr = 1;
00098   
00099   Image result(wr,h);
00100   Image::iterator rptr = result.beginw();
00101   Image::const_iterator sptr = src.begin();
00102 
00103   if (ws <= 1)
00104     result = src;
00105   else if (ws == 2)
00106     for (int y = 0; y < h; ++y)
00107       {
00108         // use kernel [1 1] / 2
00109         *rptr++ = (sptr[0] + sptr[h]) / 2.0;
00110         ++sptr;
00111       }
00112   else if (ws == 3)
00113     for (int y = 0; y < h; ++y)
00114       {
00115         // use kernel [1 2 1] / 4
00116         *rptr++ = (sptr[0] + sptr[h] * 2.0 + sptr[h2]) / 4.0;
00117         ++sptr;
00118       }
00119   else // general case for ws >= 4
00120     {
00121       // left most point - use kernel [10 10 5 1] / 26
00122       for (int y = 0; y < h; ++y)
00123         {
00124           *rptr++ = ((sptr[0] + sptr[h]) * 10.0 + 
00125                       sptr[h2] * 5.0 + sptr[h3]) / 26.0;
00126           ++sptr;
00127         }
00128       sptr -= h;
00129       
00130       // general case
00131       int x;
00132       for (x = 0; x < (ws - 5); x += 2)
00133         {
00134           for (int y = 0; y < h; ++y)
00135             {
00136               // use kernel [1 5 10 10 5 1] / 32
00137               *rptr++ = ((sptr[h]  + sptr[h4])  *  5.0 +
00138                          (sptr[h2] + sptr[h3])  * 10.0 +
00139                          (sptr[0]  + sptr[h5])) / 32.0;
00140               ++sptr;
00141             }
00142           sptr += h;
00143         }
00144         
00145       // find out how to treat the right most point
00146       if (x == (ws - 5))
00147         for (int y = 0; y < h; ++y)
00148           {
00149             // use kernel [1 5 10 10 5] / 31
00150             *rptr++ = ((sptr[h]  + sptr[h4])  *  5.0 +
00151                        (sptr[h2] + sptr[h3])  * 10.0 +
00152                         sptr[0]) / 31.0;
00153             ++sptr;
00154           }
00155       else
00156         for (int y = 0; y < h; ++y)
00157           {
00158             // use kernel [1 5 10 10] / 26
00159             *rptr++ = ( sptr[0]  + sptr[h]   * 5.0 + 
00160                        (sptr[h2] + sptr[h3]) * 10.0) / 26.0;
00161             ++sptr;
00162           }
00163     }
00164   return result;
00165 }
00166 
00167 // ######################################################################
00168 Image conv2PreserveEnergy(const Image& src, const Image& f)
00169 {
00170   ASSERT(src.isInitialized());
00171   ASSERT(ecxStr);
00172 
00173   const int sw = src.getWidth();
00174   const int sh = src.getHeight();
00175   const float ecw = ecx / sw;
00176 
00177   Image::const_iterator filter = f.begin();
00178   const int fw = f.getWidth();
00179   const int fh = f.getHeight();
00180 
00181   ASSERT((fw & 1) && (fh & 1));
00182 
00183   Image result(sw, sh);
00184   Image::const_iterator sptr = src.begin();
00185   Image::iterator rptr = result.beginw();
00186 
00187   const int fend = fw * fh - 1;
00188   const int fw2 = (fw - 1) / 2;
00189   const int fh2 = (fh - 1) / 2;
00190 
00191   const int sSkip = sh - fh;
00192 
00193   for (int rx = 0; rx < sw; ++rx)
00194     {
00195       // Determine if we're safely inside the image in the x-direction:
00196       const bool isXclean = ((rx >= fw2) && (rx <  (sw - fw2)));
00197       for (int ry = 0; ry < sh; ++ry, ++rptr)
00198         {
00199           // Determine if we're safely inside the image in the y-direction:
00200           const bool isYclean = ((ry >= fh2) && (ry <  (sh - fh2)));
00201           
00202           if (isXclean && isYclean)
00203             {
00204               // well inside the image: use straight-forward convolution
00205               float rval = 0.0f;
00206               Image::const_iterator fptr = filter+fend;
00207               Image::const_iterator sptr2 = sptr + sh*(rx-fw2) + ry - fh2;
00208 
00209               for (int fx = 0; fx < fw; ++fx)
00210                 {
00211                   for (int fy = 0; fy < fh; ++fy)
00212                     rval += (*sptr2++) * (*fptr--);
00213 
00214                   sptr2 += sSkip;
00215                 }
00216               *rptr = rval;
00217               continue;
00218             }
00219           else
00220             {
00221               // near the border: assume that the value missing pixels
00222               // is equal to the average over the present pixels
00223               // to minimize artifacts at the edge.
00224 
00225               float rval = 0.0f;
00226               float sSum = 0.0f;
00227               int sCount = 0;
00228               float fsum_skipped = 0.0f;
00229 
00230               for (int fx = 0; fx < fw; ++fx)
00231                 {
00232                   const int sx = rx + fx - fw2;
00233                   if (sx >= 0 && sx < sw)
00234                     {
00235                       for (int fy = 0; fy < fh; ++fy)
00236                         {
00237                           const float fil = filter[fend - fx * fh - fy];
00238                           const int sy = ry + fy - fh2;
00239                           if (sy >= 0 && sy < sh)
00240                             {
00241                               // here we have the pixel present: 
00242                               // store it and count it to compute the mean
00243                               const float sVal = sptr[sx * sh + sy];
00244                               rval += sVal * fil;
00245                               sSum += sVal;
00246                               ++sCount;
00247                             }
00248                           else
00249                             {
00250                               // this is accumulataing the filter values 
00251                               // over missing pixels
00252                               fsum_skipped += fil;
00253                             }
00254                         }
00255                     }
00256                   else
00257                     {
00258                       // this is accumulating the filter values
00259                       // over an entire column of mixxing pixels
00260                       for (int fy = 0; fy < fh; ++fy)
00261                         fsum_skipped += filter[fend - fx * fh - fy];
00262                     }
00263                 }
00264               // compute the average over the present pixels
00265               const float sAvg = sSum / sCount;
00266               
00267               // add the that average x accumulated filter values to the result
00268               *rptr = rval + (fsum_skipped * sAvg);
00269             }
00270         }
00271     }
00272   return result;
00273 }

Generated on Fri Sep 7 13:09:49 2007 for SaliencyToolbox by  doxygen 1.5.2