The Machine Perception Toolbox

[Introduction]- [News]- [Download]- [Screenshots]- [Manual (pdf)]- [Forums]- [API Reference]- [Repository ]

 

LGN.cpp

Go to the documentation of this file.

#include <iostream.h>
#include <mpisearch.h>
#include <Magick++.h>

using namespace Magick;
using namespace std;



class CenterSurround : public MPISearchObjectDetector< float > { 
public:
  CenterSurround();
  ~CenterSurround() {};
  void search(RImage< float> &pixels); 
};

CenterSurround::CenterSurround(){
  // Because we are working with discrete pixels, not real-valued
  // points, we run into off-by-one issues.  To compute the integral
  // between points a and b in 1-D continuous space, given the
  // integrated function f(x), we would compute f(b) - f(a).  In
  // words, you want to look up the integral at b, then subtract off
  // everything to the left of a.  In discrete pixel space, since
  // pixels take up space, in order to subtract off everything to the
  // left of a, you have to go to the pixel at a-1 so you don't
  // subtract off the integral *at* pixel a.  In our 2-D case, this
  // means that we always have to widen our features by one, up and to
  // the left.  The result is that the total size of a patch is one
  // plus the width of the window of interest in each direction.  Lets
  // take an example: To get the sum of image pixels within the
  // rectangle with upper left hand corner at (1,1) and lower right
  // hand corner at (2,2), the corners used in the *integral-image*
  // feature are (0,0), (0,2), (2,0), and (2,2).  So the total extent
  // of the feature, in terms of the distance (inclusive) between the
  // pixels we will do lookups on in the *integral* image, is 3, not
  // 2.  It may be confusing at first to see that the patch width is
  // searchWidth +1, so until you get used to it, just keep in mind
  // that this off-by-one business can be a little confusing, so be
  // careful!

  // That said, we want our basic center-surround to be of size 5 x 5.
  // Thus, the smallest valid integral image patch would be 6 x 6 (if
  // we had logic to check if we went off the edge of the image at
  // every pixel access, we could make it 5 x 5, but that would really
  // slow things down).
  data.patch_width = 5 + 1;
  data.patch_height = 5 + 1;  
  int searchWidth = data.patch_width - 1, searchHeight = data.patch_height - 1;
  
  // Now, allocate space for the features
  data.numfeatures = 2;
  data.features = new Feature[data.numfeatures]; 

  // The first feature is a rectangle over the entire patch -- this is
  // used in DC subtraction Note that we're specifying indexes into
  // pixels from 0,0 to 5,5 -- which requires a 6x6 patch of valid
  // pixels.
  data.features[0].numcorners = 4;  
  data.features[0].corners = new Corner[data.features[0].numcorners]; 
  data.features[0].corners[0].x = 0;           data.features[0].corners[0].y = 0;            data.features[0].corners[0].value =  1;
  data.features[0].corners[1].x = searchWidth; data.features[0].corners[1].y = 0;            data.features[0].corners[1].value = -1;
  data.features[0].corners[2].x = 0;           data.features[0].corners[2].y = searchHeight; data.features[0].corners[2].value = -1;
  data.features[0].corners[3].x = searchWidth; data.features[0].corners[3].y = searchHeight; data.features[0].corners[3].value =  1;
  
  int innerWidth = (int)(searchWidth*.20), outerWidth = searchWidth-innerWidth;
  int innerHeight = (int)(searchWidth*.20), outerHeight = searchHeight-innerHeight;
  data.features[1].numcorners = 8;
  // Outer rect
  data.features[1].corners = new Corner[data.features[1].numcorners]; 
  data.features[1].corners[0].x = 0;           data.features[1].corners[0].y = 0;            data.features[1].corners[0].value =  1;
  data.features[1].corners[1].x = searchWidth; data.features[1].corners[1].y = 0;            data.features[1].corners[1].value = -1;
  data.features[1].corners[2].x = 0;           data.features[1].corners[2].y = searchHeight; data.features[1].corners[2].value = -1;
  data.features[1].corners[3].x = searchWidth; data.features[1].corners[3].y = searchHeight; data.features[1].corners[3].value =  1;
  // Inner rect. In LGN X cells, the gain of the inner rectangler is
  // about -3.52 times the gain of the outer rectangle. Since "value"
  // is an integer, we use -4
  data.features[1].corners[4].x = innerWidth; data.features[1].corners[4].y = innerHeight; data.features[1].corners[4].value = -4;
  data.features[1].corners[5].x = outerWidth; data.features[1].corners[5].y = innerHeight; data.features[1].corners[5].value =  4;
  data.features[1].corners[6].x = innerWidth; data.features[1].corners[6].y = outerHeight; data.features[1].corners[6].value =  4;
  data.features[1].corners[7].x = outerWidth; data.features[1].corners[7].y = outerHeight; data.features[1].corners[7].value = -4;
  
}

void CenterSurround::search(RImage< float> &pixels){

  // stream.images is an array of pointers to the RImage<T> base
  // class.  We cast stream.images[0] as a pointer to its derived
  // class RIntegral, and call the integrate method.  Javier: Why not
  // let RImage have an integrate method itself?  Ian: Because RImage
  // is the *base* class from which RIntegral and RDerivative are
  // derived.  We do this because we want to be able to use lists of
  // images of different types.  Thus we use a list of pointers to the
  // base class, even though in this case we are always looking at
  // RIntegrals.  In the future, I think a good idea might be to use
  // the mixed list type from the Boost Library, which allows mixed
  // types in the list.  Then we wouldn't have to typescast the
  // pointer.  If we do that, then perhaps we should also think about
  // using Boost shared_pointers, which are reference counted, so we
  // don't have to worry about memory management as much.  Currently
  // memory is managed pretty well by hand, but as libmpisearch
  // evolves, use of shared_pointers would make development of less
  // error-prone code a lot easier.  The only possible downside would
  // be a potential speed hit from shared_pointer, but I think this
  // could be minimized, even down to *no* effective penalty.
  static_cast<RIntegral< float >* >(stream.images[0])->integrate(pixels);
  
  float a,b,c,d,s,mean;
  int scale_index, x, y;
  float scale_factor, area;
  unsigned int rect_ind = 0, cent_surr_ind = 1; 
  int numrectcorners = data.features[rect_ind].numcorners;
  int numcent_surrcorners = data.features[cent_surr_ind].numcorners;
  int numWindows = 0;


    // An image pyramid consists of scales and scales consist of
    // windows.  The outer for loop below iterates through scales. The
    // inner loop iterates over windows within each scale.




  // Get begin and end iterators for the outer loop.  mpi is a
  // datamember of the class MPISearchStream and it is of type
  // MPIImagePyramid.  MPIImagePyramid is a container class
  // representing all the patches at all scales in an image.

   
  MPIImagePyramid< float >::const_iterator scale = stream.mpi->begin();
  MPIImagePyramid< float >::const_iterator last_scale = stream.mpi->end();
  
  for( ; scale != last_scale; ++scale){

    // Retreive the cached info about this scale

    scale_index = scale.getScale(scale_factor);

    //The class CornerCache has 4 members: scaledCornerX,
    //scaledCornerY, scaledIndex, value 


    // Javier: What are CornerCache
    //used for?  stream is of type MPISearchStream. corners is a
    //private datamember of MPISearchStream. Corners is of type
    //vector<CornerCache<T>**>

    CornerCache< float > **corners = stream.corners[scale_index];
    CornerCache< float > *rect_corners = corners[rect_ind];
    CornerCache< float > *cs_corners = corners[cent_surr_ind];
  
    

    // stream.fns[i] is a vector which counts, at each scale i, the
    // total number of additions or subtraction of a pixel by each
    // feature. This can be used to remove the DC component of the
    // filter, i.e, we want for the filter to have zero output to
    // patches that have pixels of constant value.  For instance,
    // suppose a feature operates on a 5x5 patch by adding the 25
    // pixels of the patch and then substracting the center pixel. If
    // presented with a 5x5 patch in which all the pixels have a
    // constant intensity x then the output of the filter would be
    // (25-1) * x . If we want for the filter to have zero output to
    // constant valued patches we would then have to substract (25 -1)
    // * mean, where mean is the average value of the 25 pixels in the
    // patch.


    // Below this comment, we store the inverse of the number of
    // pixels in the 0th feature ahead of time, which helps us compute
    // the mean.  Why the inverse?  So we don't have to issue a divide
    // at each step, rather we can issue a mult, which is much faster
    // on most CPUs.
    
    // Javier: Should not we take into consideration the weights of
    // the feature componets. In our case, for example, the

    float one_over_rect_area = 1.0 / stream.fns[scale_index][rect_ind];



    // Finally, we will have to scale the image to between 0 and 1 at
    // the end due to ImageMagick's assumptions about how float images
    // are represented.
    float minval = 99999, maxval= -99999;
    



    // A pyramid consists of scales and scales consist of windows. We
    // now iterate over the windows within each scale. 



    MPIScaledImage< float >::const_iterator window = (*scale).begin(), last_window = (*scale).end();
    for( ; window != last_window; ++window, ++numWindows){      
      // First, compute sum of entire patch, so we can remove DC component.
      // since we know this is just a square, we don't need any for loops
      a = window.getPixel0( rect_corners[0].scaledIndex ) * rect_corners[0].value;
      b = window.getPixel0( rect_corners[1].scaledIndex ) * rect_corners[1].value;
      c = window.getPixel0( rect_corners[2].scaledIndex ) * rect_corners[2].value;
      d = window.getPixel0( rect_corners[3].scaledIndex ) * rect_corners[3].value;
      mean = a + b + c + d;
      mean *= one_over_rect_area;
      
      // Now, compute the center surround feature. In this case, we
      // could compute all 8 corners explicitly as above, but insetad
      // lets do it in a more generic way, so that we don't have to
      // know anything about the feature a priori. Done this way, we
      // can compute an arbitrary feature stored as feature 1 of
      // "data".  We could generalize this even more by looping over
      // data.numfeatures instead of specifying the center_surround
      // feature as we do here.  This way, you could compute N
      // features in an identical way.  To see this in action, look at
      // MPISearchObjectDetector.classifyWindow
      //
      // Instead of looping over every corner, we unroll this loop a
      // little, which clues the compiler in to the fact that the
      // value of each corner (or at least, each group of 4 corners)
      // is independent.  This allows it to schedule instructions in
      // such a way that achieved a 40% speed improvement on a
      // PowerMac G4 using gcc 3.1 due to efficient use of the
      // pipeline. This sets the stage for vectorization using AltiVec
      // or SSD, if someone ever wants to try this.
      s = 0;
      for(int corner = 0; corner < numcent_surrcorners; corner+=4) {
        a = window.getPixel0( cs_corners[corner].scaledIndex   ) * cs_corners[corner].value;
        b = window.getPixel0( cs_corners[corner+1].scaledIndex ) * cs_corners[corner+1].value;
        c = window.getPixel0( cs_corners[corner+2].scaledIndex ) * cs_corners[corner+2].value;
        d = window.getPixel0( cs_corners[corner+3].scaledIndex ) * cs_corners[corner+3].value;
        s += a + b + c + d;
      }
      // Now, remove the DC component. This is done by counting the
      // number of times a pixel was added or subtracted ( stored in
      // stream.fns ), multiplying by the mean, then subtracting it
      // from the total.
      s -= mean*stream.fns[scale_index][cent_surr_ind];

      // Now set every pixel of the corresponding block in our output image to s
      for(int i = 0; i < scale_factor; ++i)
        for(int j = 0; j < scale_factor; ++j)
          window.setPixel(1, j, i, s);

      // Finally, because of the way we are displaying the image, we
      // need to get min and max vals.
      if(s < minval)
        minval = s;
      if(s > maxval)
        maxval = s;
      window.getCoords(x,y);
    }
    // We're finished, so now display the result.  
    // First, make sure all values are between 0 and 1;
    float range = maxval - minval; 
    float one_over_range = 1.0/range;
    unsigned int numpixels = stream.images[1]->width*stream.images[1]->height;
    float * p = stream.images[1]->array;
    for(unsigned int i = 0; i < numpixels ; ++i)
      *p = (*p++ - minval)*one_over_range;
    // Now, copy the pixels into an ImageMagick image.
    Image img((const unsigned int)stream.images[1]->width,(const unsigned int)stream.images[1]->height,
              "I",FloatPixel,stream.images[1]->array);

    // Crop and scale the image for visibility.
    const Geometry cropregion(x+scale_factor,y+scale_factor);
    img.crop(cropregion);
    const Geometry newsize(stream.images[1]->width-1,stream.images[1]->height-1);
    img.scale(newsize);
    img.type(GrayscaleType);
    img.fontPointsize(40);
    //Geometry(width, height, xoffset, yoffset);
    char tmpstring[20];
    sprintf(tmpstring,"Scale %d", (int)scale_factor);
    img.annotate(tmpstring, Geometry(50,50, 40, 40) );
    img.display();
    // img.write("Output2.jpg");
  }
}

int main()
{
  // Load the image
   Image img("../../Media/ExampleImage.jpg");
  //  Image img("../../Media/spock.gif");
  // Write the word "Hello" on the image
  //   img.fontPointsize(100);
   //  Geometry(width, height, xoffset, yoffset)    
   //  img.annotate("Hello2", Geometry(400,400, 0, 220) );

  img.display();

// Construct an RImage, which is expected by MPISearch
// This constructor sets the member variables for width, height, and
// number of pixels and allocates memory for the member variable
// "array", which is a vector of size npixels

  RImage< float > pixels(img.columns(), img.rows()); 

  // Copy the pixels into an RImage, which is expected by MPISearch
  img.write(0,0, img.columns(), img.rows(), "I", FloatPixel, pixels.array); 
// The 'I' tells ImageMagick to take the intensity channel.  I think
// there are better ways to do this, but I haven't looked into them.


  CenterSurround CS;

  // stream is a data member of CenterSurround, of type
  // MPISearchStream MPISearchStream holds all the memory and caches
  // for the class MPISearchObjectDetector, which Provides an
  // architecture for performing operations on every sub-window of an
  // image at multiple scales.

  // CS.inistream calls stream.init(width,height,data,WINSHIFT)
  // were data is a data member of CenterSurround  of type FeatureData
  // stream.init 



  CS.initStream(pixels.width, pixels.height, 1.0f);

// Now run over the image at multiple scales and compute the mean
// pixel at that size


  CS.search(pixels);
  
  return 0;
}

---