test-canny.C

Go to the documentation of this file.
00001 /*! @file BeoSub/test-canny.C [put description here] */
00002 
00003 /*******************************************************************************
00004 * --------------------------------------------
00005 *(c) 2001 University of South Florida, Tampa
00006 * Use, or copying without permission prohibited.
00007 * PERMISSION TO USE
00008 * In transmitting this software, permission to use for research and
00009 * educational purposes is hereby granted.  This software may be copied for
00010 * archival and backup purposes only.  This software may not be transmitted
00011 * to a third party without prior permission of the copyright holder. This
00012 * permission may be granted only by Mike Heath or Prof. Sudeep Sarkar of
00013 * University of South Florida (sarkar@csee.usf.edu). Acknowledgment as
00014 * appropriate is respectfully requested.
00015 *
00016 *  Heath, M., Sarkar, S., Sanocki, T., and Bowyer, K. Comparison of edge
00017 *    detectors: a methodology and initial study, Computer Vision and Image
00018 *    Understanding 69 (1), 38-54, January 1998.
00019 *  Heath, M., Sarkar, S., Sanocki, T. and Bowyer, K.W. A Robust Visual
00020 *    Method for Assessing the Relative Performance of Edge Detection
00021 *    Algorithms, IEEE Transactions on Pattern Analysis and Machine
00022 *    Intelligence 19 (12),  1338-1359, December 1997.
00023 *  ------------------------------------------------------
00024 *
00025 * PROGRAM: canny_edge
00026 * PURPOSE: This program implements a "Canny" edge detector. The processing
00027 * steps are as follows:
00028 *
00029 *   1) Convolve the image with a separable gaussian filter.
00030 *   2) Take the dx and dy the first derivatives using [-1,0,1] and [1,0,-1]'.
00031 *   3) Compute the magnitude: sqrt(dx*dx+dy*dy).
00032 *   4) Perform non-maximal suppression.
00033 *   5) Perform hysteresis.
00034 *
00035 * The user must input three parameters. These are as follows:
00036 *
00037 *   sigma = The standard deviation of the gaussian smoothing filter.
00038 *   tlow  = Specifies the low value to use in hysteresis. This is a
00039 *           fraction (0-1) of the computed high threshold edge strength value.
00040 *   thigh = Specifies the high value to use in hysteresis. This fraction (0-1)
00041 *           specifies the percentage point in a histogram of the gradient of
00042 *           the magnitude. Magnitude values of zero are not counted in the
00043 *           histogram.
00044 *
00045 * NAME: Mike Heath
00046 *       Computer Vision Laboratory
00047 *       University of South Floeida
00048 *       heath@csee.usf.edu
00049 *
00050 * DATE: 2/15/96
00051 *
00052 * Modified: 5/17/96 - To write out a floating point RAW headerless file of
00053 *                     the edge gradient "up the edge" where the angle is
00054 *                     defined in radians counterclockwise from the x direction.
00055 *                     (Mike Heath)
00056 *******************************************************************************/
00057 
00058 //
00059 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/BeoSub/test-canny.C $
00060 // $Id: test-canny.C 9412 2008-03-10 23:10:15Z farhan $
00061 // Primary maintainer for this file: Zack Gossman <gossman@usc.edu>
00062 //
00063 
00064 
00065 /**************************TODO*************************
00066     -NOTE: switching to std vecotrs BREAKS powell, even if changed throughout. Note sure why
00067 
00068     -MAKE a tester execuatable to spool camera values (similar to those displayed at the start of test-multigrab) continually, checking white balance, exposure and gain in order to manually clibrate the cameras at wet test time. Auto set is not working well. NOTE that every color setting is currently rcognizing white! this is bad! something wrong with the V range? <-seems fixable using a mix of changing V and hardware stuff
00069 
00070     -SUGGESTION: instead of returning final recognized position data using a struct, inseat store data in shape class and return using reference. (allows shapes with different dimensions to return useful data, unlike current standard 5 dims)
00071 *******************************************************/
00072 
00073 #include "BeoSub/CannyModel.H"
00074 #include "BeoSub/hysteresis.H"
00075 #include "Component/ModelManager.H"
00076 #include "Devices/FrameGrabberFactory.H"
00077 #include "GUI/XWindow.H"
00078 #include "Image/ColorOps.H"
00079 #include "Image/DrawOps.H"
00080 #include "Image/FilterOps.H"
00081 #include "Image/Image.H"
00082 #include "Image/MathOps.H"
00083 #include "Image/MorphOps.H"     // for erodeImg(), closeImg()
00084 #include "Image/Transforms.H"
00085 #include "Raster/Raster.H"
00086 #include "Util/MathFunctions.H" // for squareOf()
00087 #include "Util/Timer.H"
00088 #include "Util/Types.H"
00089 #include "Util/log.H"
00090 #include "VFAT/segmentImageTrackMC.H"
00091 #include "rutz/shared_ptr.h"
00092 
00093 #include <algorithm> // for std::max()
00094 #include <cmath>
00095 #include <cstdio>
00096 #include <cstdlib>
00097 #include <cstring>
00098 
00099 //canny
00100 #define BOOSTBLURFACTOR 90.0
00101 //powell
00102 #define TOL 2.0e-4
00103 #define ITMAX 200
00104 //brent...
00105 #define ITMAXB 100
00106 #define CGOLD 0.3819660
00107 #define ZEPS 1.0e-10
00108 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
00109 #define SIGN(a,b) ((b)>=0.0?fabs(a):-fabs(a))
00110 //mnbrak
00111 #define GOLD 1.618034
00112 #define GLIMIT 100.0
00113 #define TINY 1.0e-20
00114 
00115 //CAMERA STUFF
00116 
00117 
00118 #define FREE_ARG char*
00119 //END CAMERA STUFF
00120 
00121 
00122 int canny(unsigned char *image, int rows, int cols, float sigma,
00123          float tlow, float thigh, unsigned char **edge, const char* fname);
00124 void gaussian_smooth(unsigned char *image, int rows, int cols, float sigma,
00125         short int **smoothedim);
00126 void make_gaussian_kernel(float sigma, float **kernel, int *windowsize);
00127 void derrivative_x_y(short int *smoothedim, int rows, int cols,
00128         short int **delta_x, short int **delta_y);
00129 void magnitude_x_y(short int *delta_x, short int *delta_y, int rows, int cols,
00130         short int **magnitude);
00131 void radian_direction(short int *delta_x, short int *delta_y, int rows,
00132     int cols, float **dir_radians, int xdirtag, int ydirtag);
00133 void radian_direction(short int *delta_x, short int *delta_y, int rows,
00134     int cols, float **dir_radians, int xdirtag, int ydirtag);
00135 double angle_radians(double x, double y);
00136 void grabImage(Image<PixRGB<byte> >* image);
00137 
00138 Image<float> distMap;
00139 rutz::shared_ptr<XWindow> xwin, win;
00140 const bool debugmode = true;
00141 
00142 //Moved to global
00143 
00144 bool fromFile = true;//boolean to check if input is from file or from camera
00145 
00146 const char *infilename = NULL;  /* Name of the input image */
00147 const char *dirfilename = NULL; /* Name of the output gradient direction image */
00148 const char *shapeArg = NULL;       /* Shape for recognition*/
00149 const char *colorArg = NULL;       /* Color for tracking*/
00150 char outfilename[128];    /* Name of the output "edge" image */
00151 char composedfname[128];  /* Name of the output "direction" image */
00152 unsigned char *edge;      /* The output edge image */
00153 float sigma,              /* Standard deviation of the gaussian kernel. */
00154     tlow,               /* Fraction of the high threshold in hysteresis. */
00155     thigh;              /* High hysteresis threshold control. The actual
00156                            threshold is the (100 * thigh) percentage point
00157                            in the histogram of the magnitude of the
00158                            gradient image that passes non-maximal
00159                            suppression. */
00160 
00161 int imgW, imgH, ww, hh;
00162 Image< PixRGB<byte> > colorImg;
00163 Image<byte> grayImg;
00164 
00165 Image< PixRGB<byte> > display;
00166 Image<byte> temp;
00167 Image<PixRGB<byte> > Aux;
00168 Image<PixH2SV2<float> > H2SVimage;
00169 
00170 // instantiate a model manager (for camera input):
00171 ModelManager manager("Canny Tester");
00172 // Instantiate our various ModelComponents:
00173 nub::soft_ref<FrameIstream>
00174 gb(makeIEEE1394grabber(manager, "cannycam", "cc"));
00175 /// all this should not be in global scope, but in main()!
00176 
00177 segmentImageTrackMC<float,unsigned int,4> *segmenter;
00178 
00179 //! Mean color to track
00180 std::vector<float> color(4,0.0F);
00181 
00182 bool runCanny(bool fromFile, bool useColor, const char* shapeArg); //FIX: change later so this returns struct with found info (position, etc)
00183 
00184 /************************************************************
00185 Optimizer function Prototypes (pt2)
00186 ************************************************************/
00187 
00188 //back to being ripped from a book
00189 
00190 //Note: powell and linmin now use a ShapeModel optimizee class. brent and mnbrak use func ptr to f1dim, which has copy of the optimizee
00191 bool powell(double p[], double **xi, int n, double ftol,
00192             int *iter, double *fret, rutz::shared_ptr<ShapeModel>& optimizee);
00193 int ncom;
00194 double *pcom,*xicom;
00195 
00196 double brent(double ax, double bx, double cx,
00197              double (*f)(double, rutz::shared_ptr<ShapeModel>&), double tol,
00198              double *xmin, rutz::shared_ptr<ShapeModel>& shape);
00199 void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
00200             double *fc, double (*func)(double, rutz::shared_ptr<ShapeModel>&),
00201             rutz::shared_ptr<ShapeModel>& shape);
00202 double f1dim(double x, rutz::shared_ptr<ShapeModel>& shape);
00203 void linmin(double p[], double xi[], int n, double *fret,
00204             rutz::shared_ptr<ShapeModel>& optimizee);
00205 
00206 
00207 ///////////////////////FROM NRUTIL/////////////////////////
00208 
00209 
00210 double *nrVector(long nl, long nh)
00211 /* allocate a double vector with subscript range v[nl..nh] */
00212 {
00213         double *v;
00214 
00215         v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
00216         if (!v) printf("allocation failure in nrVector()");
00217         return v-nl+1;
00218 }
00219 
00220 void free_nrVector(double *v, long nl, long nh)
00221 /* free a double vector allocated with nrVector() */
00222 {
00223         free((FREE_ARG) (v+nl-1));
00224 }
00225 
00226 //////////////////////////MAIN/////////////////////////////
00227 
00228 int main(int argc, char *argv[])
00229 {
00230   //  MYLOGVERB = 0; //sets extern value to suppress unwanted output from other modules
00231 
00232   //set math details. For now, harcoded
00233   sigma = .5;
00234   tlow = .5;
00235   thigh = .5;
00236 
00237   bool useColor = true;
00238 
00239   manager.addSubComponent(gb);
00240 
00241 
00242   // set the camera number (in IEEE1394 lingo, this is the
00243   // "subchannel" number):
00244   gb->setModelParamVal("FrameGrabberSubChan", 0);
00245 
00246 
00247 
00248   /****************************************************************************
00249    * Get the command line arguments.
00250    ****************************************************************************/
00251   if(argc < 3){
00252     fprintf(stderr,"\n<USAGE> %s image shape color \n",argv[0]);
00253     fprintf(stderr,"\n      image:      An image to process. Must be in PGM format.\n");
00254     fprintf(stderr,"                  Type 'none' for camera input.\n");
00255     fprintf(stderr,"      shape:      Shape on which to run recognition\n");
00256     fprintf(stderr,"                  Candidates: Rectangle, Square, Circle, Octagon.\n");
00257     fprintf(stderr,"      color:       Color to track\n");
00258     fprintf(stderr,"                  Candidates: Blue, Yellow, none (for no color tracking).\n");
00259     exit(1);
00260   }
00261   infilename = argv[1];
00262   shapeArg = argv[2];
00263   colorArg = argv[3];
00264 
00265   printf("READ: 1: %s 2: %s 3: %s\n", infilename, shapeArg, colorArg);
00266 
00267   if(!strcmp(infilename, "none")){
00268     fromFile = false;
00269     infilename = "../../bin/Camera.pgm";//hacky~ should let user say where they want file placed
00270   }
00271 
00272   //what's this?
00273   if(argc == 6) dirfilename = infilename;
00274   else dirfilename = NULL;
00275 
00276 
00277   /***************************************************************************
00278    * Perform color exclusion. Using grayImg for now, should likely use separate image later to preserve step output.
00279    *NOTE: a lot of code is repeated in file and camera input, because up in-if/else declaration of segmenter
00280    **************************************************************************/
00281 
00282 
00283   //COLORTRACKING DECS: note that much of these may not need to be in this file
00284   //Colors MUST use H2SV2 pixel values! use test-sampleH2SV2 to sample needed values! --Z
00285 
00286   //0 = H1 (0-1), 1=H2 (0-1), 2=S (0-1), 3=V (0-1)
00287   if(!strcmp(colorArg, "Red")){
00288     //RED
00289     color[0] = 0.70000F; color[1] = 0.94400F; color[2] = 0.67200F; color[3] = 0.96900F;
00290   }
00291 
00292   else if(!strcmp(colorArg, "Green")){
00293     //Green
00294     color[0] = 0.53000F; color[1] = 0.36000F; color[2] = 0.44600F; color[3] = 0.92000F;
00295   }
00296 
00297   else if(!strcmp(colorArg, "Orange")){
00298     //"ORANGE"
00299     color[0] = 0.80000F; color[1] = 0.80000F; color[2] = 0.57600F; color[3] = 0.96800F;
00300   }
00301   else if(!strcmp(colorArg, "Blue")){
00302     //BLUE
00303     color[0] = 0.27500F; color[1] = 0.52500F; color[2] = 0.44500F; color[3] = 0.97200F;
00304   }
00305   else if(!strcmp(colorArg, "Yellow")){
00306     //YELLOW
00307     color[0] = 0.90000F; color[1] = 0.62000F; color[2] = 0.54000F; color[3] = 0.97000F;
00308   }
00309   else{
00310     useColor = false;
00311   }
00312 
00313    //! +/- tollerance value on mean for track
00314    std::vector<float> std(4,0.0F);
00315    //NOTE that the saturation tolerance is important (if it gos any higher thn this, it will nearly always recognize white!)
00316    std[0] = 0.30000; std[1] = 0.30000; std[2] = 0.44500; std[3] = 0.45000;
00317 
00318    //! normalizer over color values (highest value possible)
00319    std::vector<float> norm(4,0.0F);
00320    norm[0] = 1.0F; norm[1] = 1.0F; norm[2] = 1.0F; norm[3] = 1.0F;
00321 
00322    //! how many standard deviations out to adapt, higher means less bias
00323    //The lower these are, the more strict recognition will be in subsequent frames
00324    //TESTED AND PROVED do NOT change unless you're SURE
00325    std::vector<float> adapt(4,0.0F);
00326    //adapt[0] = 3.5F; adapt[1] = 3.5F; adapt[2] = 3.5F; adapt[3] = 3.5F;
00327    adapt[0] = 5.0F; adapt[1] = 5.0F; adapt[2] = 5.0F; adapt[3] = 5.0F;
00328 
00329    //! highest value for color adaptation possible (hard boundry)
00330    std::vector<float> upperBound(4,0.0F);
00331    upperBound[0] = color[0] + 0.45F; upperBound[1] = color[1] + 0.45F; upperBound[2] = color[2] + 0.55F; upperBound[3] = color[3] + 0.55F;
00332 
00333    //! lowest value for color adaptation possible (hard boundry)
00334    std::vector<float> lowerBound(4,0.0F);
00335    lowerBound[0] = color[0] - 0.45F; lowerBound[1] = color[1] - 0.45F; lowerBound[2] = color[2] - 0.55F; lowerBound[3] = color[3] - 0.55F;
00336    //END DECS
00337 
00338    //set up image details
00339    if(fromFile){
00340      colorImg = Raster::ReadRGB(infilename);
00341      imgW = colorImg.getWidth();
00342      imgH = colorImg.getHeight();
00343    }
00344    else{
00345      /**************************************************
00346        CAMERA INPUT
00347         -Code ripped from test-grab. Forces to IEEE1394. has short stream for buffering and reads in pic from camera, then starts powell minimization, which displays in new window.
00348      **************************************************/
00349 
00350      imgH = gb->getHeight();
00351      imgW = gb->getWidth();
00352 
00353  }
00354 
00355    segmenter = new segmentImageTrackMC<float,unsigned int,4> (imgW*imgH);
00356    //color tracking stuff
00357    segmenter->SITsetTrackColor(&color,&std,&norm,&adapt,&upperBound,&lowerBound);
00358 
00359    ww = imgW/4;
00360    hh = imgH/4;
00361 
00362    segmenter->SITsetFrame(&ww,&hh);
00363 
00364    // Set display colors for output of tracking. Strictly aesthetic //TEST
00365    //And definitely not needed, although for some reason it's required --Z
00366    segmenter->SITsetCircleColor(0,255,0);
00367    segmenter->SITsetBoxColor(255,255,0,0,255,255);
00368    segmenter->SITsetUseSmoothing(true,10);
00369 
00370    //method provided by Nathan to turn off adaptation (useful for dynamic environments)
00371    segmenter->SITtoggleColorAdaptation(false);
00372    //enable single-frame tracking
00373    segmenter->SITtoggleCandidateBandPass(false);
00374 
00375    Image<PixRGB<byte> > Aux;
00376    Aux.resize(100,450,true);
00377    Image<byte> temp;
00378    Image< PixRGB<byte> > display = colorImg;
00379 
00380    if(fromFile){   //for input from file
00381 
00382      if(useColor){
00383        segmenter->SITtrackImage(colorImg,&display,&Aux,true);
00384        temp = segmenter->SITreturnCandidateImage();//get image we'll be using for edge detection
00385 
00386        grayImg = luminance(temp);
00387      }
00388      else{
00389        grayImg = luminance(colorImg);
00390      }
00391    }
00392    else{
00393      // get ready for main loop:
00394      win.reset( new XWindow(gb->peekDims(), -1, -1, "test-canny window") );
00395      //color tracking stuffages
00396      segmenter->SITsetTrackColor(&color,&std,&norm,&adapt,&upperBound,&lowerBound);
00397      ww = imgW/4;
00398      hh = imgH/4;
00399      segmenter->SITsetFrame(&ww,&hh);
00400      // Set display colors for output of tracking. Strictly asthetic //TEST
00401      //And definitely not needed, although for soe reason it's required --Z
00402      segmenter->SITsetCircleColor(0,255,0);
00403      segmenter->SITsetBoxColor(255,255,0,0,255,255);
00404      segmenter->SITsetUseSmoothing(true,10);
00405 
00406      Aux.resize(100,450,true);
00407 
00408    }
00409 
00410    if(fromFile){
00411      runCanny(fromFile, useColor, shapeArg);
00412    }
00413    else{
00414 
00415      // let's get all our ModelComponent instances started:
00416      manager.start();
00417 
00418      // get the frame grabber to start streaming:
00419      gb->startStream();
00420 
00421      colorImg = gb->readRGB();
00422 
00423      xwin.reset(new XWindow(colorImg.getDims()));
00424      xwin->setPosition(colorImg.getWidth(), 0);
00425      while(1){
00426        runCanny(fromFile, useColor, shapeArg);
00427      }
00428 
00429     manager.stop();
00430    }
00431 }
00432 
00433 bool runCanny(bool fromFile, bool useColor, const char* shapeArg){
00434 
00435   if(!fromFile){
00436     int i;
00437     for(i = 0; i < 30; i++) {//give camera holder time to aim
00438 
00439       colorImg = gb->readRGB();
00440       win->drawImage(colorImg);
00441     }
00442 
00443     //Get candidate image representing blobs of desired color
00444     H2SVimage = colorImg;//changed since transfer to H2SV
00445     display = colorImg;
00446     segmenter->SITtrackImageAny(H2SVimage,&display,&Aux,true);
00447     temp = segmenter->SITreturnCandidateImage();//get image we'll be using for edge detection
00448 
00449     //set up structuring elements (disk)
00450     Image<byte> se1(5, 5, ZEROS);
00451     Point2D<int> t1(2,2);
00452     drawDisk(se1, t1, 2, static_cast<byte>(255));
00453 
00454     Image<byte> se2(9, 9, ZEROS);
00455     Point2D<int> t2(4,4);
00456     drawDisk(se2, t2, 3, static_cast<byte>(255));
00457 
00458     //Erode image to prevent noise NOTE: isn't working well, erodes too much, probably because temp is so lowrez
00459     //temp = erodeImg(temp, se1);
00460 
00461     //Close candidate image to prevent holes
00462     temp = closeImg(temp, se2);
00463 
00464     if(useColor){
00465       grayImg = (Image<byte>)(quickInterpolate(temp,4));//note that this is very low rez
00466     }
00467     else{
00468       grayImg = luminance(colorImg);
00469     }
00470     //END CAMERA
00471   }
00472    /****************************************************************************
00473     * Perform the edge detection.
00474     ****************************************************************************/
00475   /*
00476    if(dirfilename != NULL){
00477      if(!fromFile){
00478        sprintf(composedfname, "%s_s_%3.2f_l_%3.2f_h_%3.2f.fim", infilename,
00479                sigma, tlow, thigh);
00480      }
00481      else{
00482        sprintf(composedfname, "%s_s_%3.2f_l_%3.2f_h_%3.2f.fim", "Camera",
00483                sigma, tlow, thigh);
00484      }
00485      dirfilename = composedfname;
00486    }
00487   */
00488    canny(grayImg.getArrayPtr(), grayImg.getHeight(), grayImg.getWidth(), sigma, tlow, thigh, &edge, dirfilename);
00489 
00490 
00491    /****************************************************************************
00492     * Write out the edge image to a file.
00493     *****************************************************************************/
00494    /*
00495    if(!fromFile){
00496      sprintf(outfilename, "%s_s_%3.2f_l_%3.2f_h_%3.2f.pgm", infilename, sigma, tlow, thigh);
00497    }
00498    else{
00499      sprintf(outfilename, "%s_s_%3.2f_l_%3.2f_h_%3.2f.pgm", "Camera", sigma, tlow, thigh);
00500    }
00501    */
00502    Image<unsigned char> cannyCharImg(edge, grayImg.getWidth(), grayImg.getHeight());
00503 
00504    Image<float> cannyImg = cannyCharImg;
00505 
00506    Raster::WriteGray(cannyImg, outfilename, RASFMT_PNM);//edge image
00507 
00508    xwin->drawImage((Image< PixRGB<byte> >)grayImg);
00509 
00510    Point2D<int> cent = centroid(cannyCharImg);
00511    distMap = chamfer34(cannyImg, 5000.0f);
00512    //END EDGE DETECTION
00513 
00514    //PERFORM SHAPE DETECTION
00515    rutz::shared_ptr<ShapeModel> shape;
00516    double* p = NULL;
00517 
00518    int iter = 0;
00519    double fret = 0.0;
00520    bool shapeFound = false;
00521 
00522    //NOTE: following bracket set is to define scope outside which the shape model will be destroyed
00523    {
00524      /***********************************************
00525       *Set values for the model and declare it here
00526         Common Values:
00527            x_center: initial central x position of the model
00528            y_center: initial central y position of the model
00529               -in the future, try multiple arrangements of these,
00530                     perhaps in some meta-omptimizer, since noisy
00531                     images are sensitive to these values
00532            height/width/etc: initial dimensions of the model.
00533                Be careful to note if ratio is to be fixed or free!
00534            alpha: initial angular orientation of the model
00535      **********************************************/
00536 
00537      if(!strcmp(shapeArg, "Rectangle")){
00538        //rectangle
00539        //NOTE: the order of dimensions is changed. This is a test. May need to change thm back
00540 
00541        p = (double*)calloc(6, sizeof(double));
00542        p[1] = (double)(cent.i); //Xcenter
00543        p[2] = (double)(cent.j); //Ycenter
00544        p[4] = 170.0f; // Width
00545        p[5] = 140.0f; // Height
00546        p[3] = 0.0; //alpha
00547 
00548        shape.reset(new RectangleShape(120.0, p, true));
00549      }
00550      else if(!strcmp(shapeArg, "Square")){
00551        //try it with a square
00552 
00553        p = (double*)calloc(5, sizeof(double));
00554        p[1] = (double)(cent.i); //Xcenter
00555        p[2] = (double)(cent.j); //Ycenter
00556        p[3] = 100.0f; // Height
00557        p[4] = 0.0; // alpha
00558 
00559        shape.reset(new SquareShape(100.0, p, true));
00560      }
00561      else if(!strcmp(shapeArg, "Octagon")){
00562        //try it with an octagon
00563 
00564        p = (double*)calloc(5, sizeof(double));
00565        p[1] = (double)(cent.i); //Xcenter
00566        p[2] = (double)(cent.j); //Ycenter
00567        p[3] = 80.0f; // Height
00568        p[4] = 0.0; // alpha
00569 
00570        shape.reset(new OctagonShape(80.0, p, true));
00571      }
00572      //    else if(!strcmp(shapeArg, "Parallel")){
00573      //        //parallel lines (SUCKS for now)
00574      //        p[1] = (double)(cent.i); //Xcenter
00575      //        p[2] = (double)(cent.j); //Ycenter
00576      //        p[3] = 50.0f; // Width
00577      //        p[4] = 160.0f; // Height
00578      //        p[5] = 0.0; //alpha
00579 
00580      //        shape = new ParallelShape(40.0, true);
00581      //      }
00582      else if(!strcmp(shapeArg, "Circle")){
00583        //circle
00584        p = (double*)calloc(4, sizeof(double));
00585        p[1] = (double)(cent.i); //Xcenter
00586        p[2] = (double)(cent.j); //Ycenter
00587        p[3] = 50.0f; // Radius
00588 
00589 
00590        shape.reset(new CircleShape(40.0, p, true));
00591      }
00592      else{
00593        printf("Cannot run shape recognition without a shape to recognize! Returning...\n");
00594        shape.reset(new CircleShape(9999.0, p, false));
00595        return(false);
00596      }
00597 
00598      double** xi = (double**)calloc(shape->getNumDims(), sizeof(double));
00599 
00600      //NOTE: HERE is where the initial direction matrix is built. Change these values to change axis assignment
00601      for(int i=0; i < shape->getNumDims(); i++) {
00602        xi[i] = (double*)calloc(shape->getNumDims(), sizeof(double));
00603        xi[i][i] = 1.0;
00604      }
00605      //TEST axes
00606      xi[0][0] = 0.5;
00607      xi[0][1] = 0.5;
00608 
00609      //here, we will be using powell classes, and, instead of passing in calcDist, will be passing in a shape.
00610      shapeFound = powell(shape->getDimensions(), xi, shape->getNumDims(), 0.5, &iter, &fret, shape);
00611      if(shapeFound){
00612        p = shape->getDimensions();//TEST
00613        printf("Shape found!\n");
00614        printf("Final dimensions are -- ");
00615        for(int w = 1; w <=shape->getNumDims(); w++){
00616          printf("%d: %f ", w, p[w]);//outputs final values
00617        }
00618        printf("Certainty (lower is better): %f\n", fret);
00619      }else{
00620        printf("Shape not found -- certainty was %f\n", fret);
00621      }
00622    }
00623    //commenting out annoying outputs -- Z
00624 
00625    //printf("%f %f %f %f %f %f\n", p[0], p[1], p[2], p[3], p[4], fret);
00626    //sprintf(outfilename, "%s_d_%3.2f_l_%3.2f_h_%3.2f.pgm", infilename, sigma, tlow, thigh);
00627    //Raster::WriteFloat(distMap, FLOAT_NORM_0_255, outfilename);//distance transfer
00628 
00629    return(shapeFound);
00630 }
00631 
00632 void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
00633             double *fc, double (*func)(double, rutz::shared_ptr<ShapeModel>& ),
00634             rutz::shared_ptr<ShapeModel>& shape)
00635 {
00636         double ulim,u,r,q,fu,dum;
00637 
00638         *fa=(*func)(*ax, shape);
00639         *fb=(*func)(*bx, shape);
00640         if (*fb > *fa) {
00641                 SHFT(dum,*ax,*bx,dum)
00642                 SHFT(dum,*fb,*fa,dum)
00643         }
00644         *cx=(*bx)+GOLD*(*bx-*ax);
00645         *fc=(*func)(*cx, shape);
00646         while (*fb > *fc) {
00647                 r=(*bx-*ax)*(*fb-*fc);
00648                 q=(*bx-*cx)*(*fb-*fa);
00649                 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
00650                         (2.0*SIGN(std::max(fabs(q-r),TINY),q-r));
00651                 ulim=(*bx)+GLIMIT*(*cx-*bx);
00652                 if ((*bx-u)*(u-*cx) > 0.0) {
00653                         fu=(*func)(u, shape);
00654                         if (fu < *fc) {
00655                                 *ax=(*bx);
00656                                 *bx=u;
00657                                 *fa=(*fb);
00658                                 *fb=fu;
00659                                 return;
00660                         } else if (fu > *fb) {
00661                                 *cx=u;
00662                                 *fc=fu;
00663                                 return;
00664                         }
00665                         u=(*cx)+GOLD*(*cx-*bx);
00666                         fu=(*func)(u, shape);
00667                 } else if ((*cx-u)*(u-ulim) > 0.0) {
00668                         fu=(*func)(u, shape);
00669                         if (fu < *fc) {
00670                                 SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
00671                                 SHFT(*fb,*fc,fu,(*func)(u, shape))
00672                         }
00673                 } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
00674                         u=ulim;
00675                         fu=(*func)(u, shape);
00676                 } else {
00677                         u=(*cx)+GOLD*(*cx-*bx);
00678                         fu=(*func)(u, shape);
00679                 }
00680                 SHFT(*ax,*bx,*cx,u)
00681                 SHFT(*fa,*fb,*fc,fu)
00682         }
00683 }
00684 
00685 double brent(double ax, double bx, double cx, double (*f)(double, rutz::shared_ptr<ShapeModel>&), double tol,
00686         double *xmin, rutz::shared_ptr<ShapeModel>& shape)
00687 {
00688         int iter;
00689         double a,b,d=0.0,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
00690         double e=0.0;
00691 
00692         a=(ax < cx ? ax : cx);
00693         b=(ax > cx ? ax : cx);
00694         x=w=v=bx;
00695         fw=fv=fx=(*f)(x, shape);
00696         for (iter=1;iter<=ITMAXB;iter++) {
00697                 xm=0.5*(a+b);
00698                 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
00699                 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
00700                         *xmin=x;
00701                         return fx;
00702                 }
00703                 if (fabs(e) > tol1) {
00704                         r=(x-w)*(fx-fv);
00705                         q=(x-v)*(fx-fw);
00706                         p=(x-v)*q-(x-w)*r;
00707                         q=2.0*(q-r);
00708                         if (q > 0.0) p = -p;
00709                         q=fabs(q);
00710                         etemp=e;
00711                         e=d;
00712                         if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
00713                                 d=CGOLD*(e=(x >= xm ? a-x : b-x));
00714                         else {
00715                                 d=p/q;
00716                                 u=x+d;
00717                                 if (u-a < tol2 || b-u < tol2)
00718                                         d=SIGN(tol1,xm-x);
00719                         }
00720                 } else {
00721                         d=CGOLD*(e=(x >= xm ? a-x : b-x));
00722                 }
00723                 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
00724                 fu=(*f)(u, shape);
00725                 if (fu <= fx) {
00726                         if (u >= x) a=x; else b=x;
00727                         SHFT(v,w,x,u)
00728                         SHFT(fv,fw,fx,fu)
00729                 } else {
00730                         if (u < x) a=u; else b=u;
00731                         if (fu <= fw || w == x) {
00732                                 v=w;
00733                                 w=u;
00734                                 fv=fw;
00735                                 fw=fu;
00736                         } else if (fu <= fv || v == x || v == w) {
00737                                 v=u;
00738                                 fv=fu;
00739                         }
00740                 }
00741         }
00742         printf("Too many iterations in brent");
00743         *xmin=x;
00744         return fx;
00745 }
00746 
00747 //Changed to run dim on ShapeModel
00748 double f1dim(double x, rutz::shared_ptr<ShapeModel>& shape)
00749 {
00750         int j;
00751         double f, *xt;
00752         xt=nrVector(1,ncom);
00753 
00754         for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
00755         f=shape->calcDist(xt, distMap);
00756         free_nrVector(xt,1,ncom);
00757         return f;
00758 }
00759 
00760 //Changed to accept ShapeModel optimizee class
00761 void linmin(double p[], double xi[], int n, double *fret,
00762             rutz::shared_ptr<ShapeModel>& optimizee)
00763 {
00764         int j;
00765         double xx,xmin,fx,fb,fa,bx,ax;
00766 
00767         ncom=n;
00768         pcom=nrVector(1,n);
00769         xicom=nrVector(1,n);
00770         for (j=1;j<=n;j++) {
00771                 pcom[j]=p[j];
00772                 xicom[j]=xi[j];
00773         }
00774 
00775         ax=0.0;
00776         xx=1.0;
00777 
00778         mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim,optimizee);
00779         *fret=brent(ax,xx,bx,f1dim,TOL,&xmin, optimizee);
00780         for (j=1;j<=n;j++) {
00781                 xi[j] *= xmin;
00782                 p[j] += xi[j];
00783         }
00784         free_nrVector(xicom,1,n);
00785         free_nrVector(pcom,1,n);
00786 }
00787 
00788 //Changed to accept ShapeModel optimizee class
00789 bool powell(double p[], double **xi, int n, double ftol,
00790             int *iter, double *fret, rutz::shared_ptr<ShapeModel>& optimizee)
00791 {
00792         int i,ibig,j;
00793         double del,fp,fptt,t, *pt, *ptt, *xit;
00794 
00795         pt=nrVector(1,n);
00796         ptt=nrVector(1,n);
00797         xit=nrVector(1,n);
00798 
00799         fptt = optimizee->getThreshold(); //hacky, to prevent warning
00800 
00801         *fret=optimizee->calcDist(p, distMap);
00802         for (j=1;j<=n;j++) pt[j]=p[j];
00803         for (*iter=1;;++(*iter)) {
00804                 fp=(*fret);
00805                 ibig=0;
00806                 del=0.0;
00807                 for (i=1;i<=n;i++) {
00808                         for (j = 1; j <= n; j++){
00809                           //NOTE: This stupid change using -1 for the xi vector is necessary to handle weirdness of the nrVectors
00810                           xit[j] = xi[j-1][i-1];
00811                         }
00812                         fptt=(*fret);
00813                         linmin(p,xit,n,fret,optimizee);
00814                         if (fabs(fptt-(*fret)) > del) {
00815                                 del=fabs(fptt-(*fret));
00816                                 ibig=i;
00817                         }
00818                 }
00819                 if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
00820 
00821                         //If cost is below threshold, then return that the shape was found
00822                         //perhaps this is not in the most ideal spot. --Z
00823                         if(fptt < optimizee->getThreshold()){/*if cost estimate is below threshold, shape is present in
00824                                            image (very rudimentary, need to take dimensions into account)*/
00825                           return(true);//test
00826                         }
00827                         //otherwise, return that it was not
00828                         else{return(false);}//test
00829                 }
00830                 if (*iter == ITMAX){
00831                   return(false);//test
00832                   printf("powell exceeding maximum iterations.");
00833                 }
00834                 for (j=1;j<=n;j++) {
00835                         ptt[j]=2.0*p[j]-pt[j];
00836                         xit[j]=p[j]-pt[j];
00837                         pt[j]=p[j];
00838                 }
00839                 fptt=optimizee->calcDist(ptt, distMap);
00840                 if (fptt < fp) {
00841                   t=2.0*(fp-2.0*(*fret)+fptt)*squareOf(fp-(*fret)-del)-del*squareOf(fp-fptt);
00842                   if (t < 0.0) {
00843                     linmin(p,xit,n,fret,optimizee);
00844                     for (j=0;j<n;j++) {
00845                       xi[j][ibig]=xi[j][n];
00846                       xi[j][n]=xit[j];
00847                     }
00848                   }
00849                 }
00850 
00851         }
00852 }
00853 
00854 /*******************************************************************************
00855 * PROCEDURE: canny
00856 * PURPOSE: To perform canny edge detection.
00857 * NAME: Mike Heath
00858 * DATE: 2/15/96
00859 //Pradeep: returns the centroid of the "white" pixels
00860 *******************************************************************************/
00861 int canny(unsigned char *image, int rows, int cols, float sigma,
00862          float tlow, float thigh, unsigned char **edge, const char* fname)
00863 {
00864    FILE *fpdir=NULL;          /* File to write the gradient image to.     */
00865    unsigned char *nms;        /* Points that are local maximal magnitude. */
00866    short int *smoothedim,     /* The image after gaussian smoothing.      */
00867              *delta_x,        /* The first devivative image, x-direction. */
00868              *delta_y,        /* The first derivative image, y-direction. */
00869              *magnitude;      /* The magnitude of the gadient image.      */
00870    float *dir_radians=NULL;   /* Gradient direction image.                */
00871 
00872    /****************************************************************************
00873    * Perform gaussian smoothing on the image using the input standard
00874    * deviation.
00875    ****************************************************************************/
00876    gaussian_smooth(image, rows, cols, sigma, &smoothedim);
00877 
00878    /****************************************************************************
00879    * Compute the first derivative in the x and y directions.
00880    ****************************************************************************/
00881    derrivative_x_y(smoothedim, rows, cols, &delta_x, &delta_y);
00882 
00883    /****************************************************************************
00884    * This option to write out the direction of the edge gradient was added
00885    * to make the information available for computing an edge quality figure
00886    * of merit.
00887    ****************************************************************************/
00888    if(fname != NULL){
00889       /*************************************************************************
00890       * Compute the direction up the gradient, in radians that are
00891       * specified counteclockwise from the positive x-axis.
00892       *************************************************************************/
00893       radian_direction(delta_x, delta_y, rows, cols, &dir_radians, -1, -1);
00894 
00895       /*************************************************************************
00896       * Write the gradient direction image out to a file.
00897       *************************************************************************/
00898       if((fpdir = fopen(fname, "wb")) == NULL){
00899          fprintf(stderr, "Error opening the file %s for writing.\n", fname);
00900          exit(1);
00901       }
00902       fwrite(dir_radians, sizeof(float), rows*cols, fpdir);
00903       fclose(fpdir);
00904       free(dir_radians);
00905    }
00906 
00907    /****************************************************************************
00908    * Compute the magnitude of the gradient.
00909    ****************************************************************************/
00910    magnitude_x_y(delta_x, delta_y, rows, cols, &magnitude);
00911 
00912    /****************************************************************************
00913    * Perform non-maximal suppression.
00914    ****************************************************************************/
00915    if((nms = (unsigned char *) calloc(rows*cols,sizeof(unsigned char)))==NULL){
00916       fprintf(stderr, "Error allocating the nms image.\n");
00917       exit(1);
00918    }
00919 
00920    non_max_supp(magnitude, delta_x, delta_y, rows, cols, nms);
00921 
00922    /****************************************************************************
00923    * Use hysteresis to mark the edge pixels.
00924    ****************************************************************************/
00925    if((*edge=(unsigned char *)calloc(rows*cols,sizeof(unsigned char))) ==NULL){
00926       fprintf(stderr, "Error allocating the edge image.\n");
00927       exit(1);
00928    }
00929    //printf("\n%d %d %d %d %d %d\n", magnitude, nms, rows, cols, tlow, thigh);
00930      int centroid = apply_hysteresis(magnitude, nms, rows, cols, tlow, thigh, *edge);
00931    /****************************************************************************
00932    * Free all of the memory that we allocated except for the edge image that
00933    * is still being used to store out result.
00934    ****************************************************************************/
00935    free(smoothedim);
00936    free(delta_x);
00937    free(delta_y);
00938    free(magnitude);
00939    free(nms);
00940    return centroid;
00941 }
00942 
00943 /*******************************************************************************
00944 * Procedure: radian_direction
00945 * Purpose: To compute a direction of the gradient image from component dx and
00946 * dy images. Because not all derriviatives are computed in the same way, this
00947 * code allows for dx or dy to have been calculated in different ways.
00948 *
00949 * FOR X:  xdirtag = -1  for  [-1 0  1]
00950 *         xdirtag =  1  for  [ 1 0 -1]
00951 *
00952 * FOR Y:  ydirtag = -1  for  [-1 0  1]'
00953 *         ydirtag =  1  for  [ 1 0 -1]'
00954 *
00955 * The resulting angle is in radians measured counterclockwise from the
00956 * xdirection. The angle points "up the gradient".
00957 *******************************************************************************/
00958 void radian_direction(short int *delta_x, short int *delta_y, int rows,
00959     int cols, float **dir_radians, int xdirtag, int ydirtag)
00960 {
00961    int r, c, pos;
00962    float *dirim=NULL;
00963    double dx, dy;
00964 
00965    /****************************************************************************
00966    * Allocate an image to store the direction of the gradient.
00967    ****************************************************************************/
00968    if((dirim = (float *) calloc(rows*cols, sizeof(float))) == NULL){
00969       fprintf(stderr, "Error allocating the gradient direction image.\n");
00970       exit(1);
00971    }
00972    *dir_radians = dirim;
00973 
00974    for(r=0,pos=0;r<rows;r++){
00975       for(c=0;c<cols;c++,pos++){
00976          dx = (double)delta_x[pos];
00977          dy = (double)delta_y[pos];
00978 
00979          if(xdirtag == 1) dx = -dx;
00980          if(ydirtag == -1) dy = -dy;
00981 
00982          dirim[pos] = (float)angle_radians(dx, dy);
00983       }
00984    }
00985 }
00986 
00987 /*******************************************************************************
00988 * FUNCTION: angle_radians
00989 * PURPOSE: This procedure computes the angle of a vector with components x and
00990 * y. It returns this angle in radians with the answer being in the range
00991 * 0 <= angle <2*PI.
00992 *******************************************************************************/
00993 double angle_radians(double x, double y)
00994 {
00995    double xu, yu, ang;
00996 
00997    xu = fabs(x);
00998    yu = fabs(y);
00999 
01000    if((xu == 0) && (yu == 0)) return(0);
01001 
01002    ang = atan(yu/xu);
01003 
01004    if(x >= 0){
01005       if(y >= 0) return(ang);
01006       else return(2*M_PI - ang);
01007    }
01008    else{
01009       if(y >= 0) return(M_PI - ang);
01010       else return(M_PI + ang);
01011    }
01012 }
01013 
01014 /*******************************************************************************
01015 * PROCEDURE: magnitude_x_y
01016 * PURPOSE: Compute the magnitude of the gradient. This is the square root of
01017 * the sum of the squared derivative values.
01018 * NAME: Mike Heath
01019 * DATE: 2/15/96
01020 *******************************************************************************/
01021 void magnitude_x_y(short int *delta_x, short int *delta_y, int rows, int cols,
01022         short int **magnitude)
01023 {
01024    int r, c, pos, sq1, sq2;
01025 
01026    /****************************************************************************
01027    * Allocate an image to store the magnitude of the gradient.
01028    ****************************************************************************/
01029    if((*magnitude = (short *) calloc(rows*cols, sizeof(short))) == NULL){
01030       fprintf(stderr, "Error allocating the magnitude image.\n");
01031       exit(1);
01032    }
01033 
01034    for(r=0,pos=0;r<rows;r++){
01035       for(c=0;c<cols;c++,pos++){
01036          sq1 = (int)delta_x[pos] * (int)delta_x[pos];
01037          sq2 = (int)delta_y[pos] * (int)delta_y[pos];
01038          (*magnitude)[pos] = (short)(0.5 + sqrt((float)sq1 + (float)sq2));
01039       }
01040    }
01041 
01042 }
01043 
01044 /*******************************************************************************
01045 * PROCEDURE: derrivative_x_y
01046 * PURPOSE: Compute the first derivative of the image in both the x any y
01047 * directions. The differential filters that are used are:
01048 *
01049 *                                          -1
01050 *         dx =  -1 0 +1     and       dy =  0
01051 *                                          +1
01052 *
01053 * NAME: Mike Heath
01054 * DATE: 2/15/96
01055 *******************************************************************************/
01056 void derrivative_x_y(short int *smoothedim, int rows, int cols,
01057         short int **delta_x, short int **delta_y)
01058 {
01059    int r, c, pos;
01060 
01061    /****************************************************************************
01062    * Allocate images to store the derivatives.
01063    ****************************************************************************/
01064    if(((*delta_x) = (short *) calloc(rows*cols, sizeof(short))) == NULL){
01065       fprintf(stderr, "Error allocating the delta_x image.\n");
01066       exit(1);
01067    }
01068    if(((*delta_y) = (short *) calloc(rows*cols, sizeof(short))) == NULL){
01069       fprintf(stderr, "Error allocating the delta_x image.\n");
01070       exit(1);
01071    }
01072 
01073    /****************************************************************************
01074    * Compute the x-derivative. Adjust the derivative at the borders to avoid
01075    * losing pixels.
01076    ****************************************************************************/
01077    for(r=0;r<rows;r++){
01078       pos = r * cols;
01079       (*delta_x)[pos] = smoothedim[pos+1] - smoothedim[pos];
01080       pos++;
01081       for(c=1;c<(cols-1);c++,pos++){
01082          (*delta_x)[pos] = smoothedim[pos+1] - smoothedim[pos-1];
01083       }
01084       (*delta_x)[pos] = smoothedim[pos] - smoothedim[pos-1];
01085    }
01086 
01087    /****************************************************************************
01088    * Compute the y-derivative. Adjust the derivative at the borders to avoid
01089    * losing pixels.
01090    ****************************************************************************/
01091    for(c=0;c<cols;c++){
01092       pos = c;
01093       (*delta_y)[pos] = smoothedim[pos+cols] - smoothedim[pos];
01094       pos += cols;
01095       for(r=1;r<(rows-1);r++,pos+=cols){
01096          (*delta_y)[pos] = smoothedim[pos+cols] - smoothedim[pos-cols];
01097       }
01098       (*delta_y)[pos] = smoothedim[pos] - smoothedim[pos-cols];
01099    }
01100 }
01101 
01102 /*******************************************************************************
01103 * PROCEDURE: gaussian_smooth
01104 * PURPOSE: Blur an image with a gaussian filter.
01105 * NAME: Mike Heath
01106 * DATE: 2/15/96
01107 *******************************************************************************/
01108 void gaussian_smooth(unsigned char *image, int rows, int cols, float sigma,
01109         short int **smoothedim)
01110 {
01111    int r, c, rr, cc,     /* Counter variables. */
01112       windowsize,        /* Dimension of the gaussian kernel. */
01113       center;            /* Half of the windowsize. */
01114    float *tempim,        /* Buffer for separable filter gaussian smoothing. */
01115          *kernel,        /* A one dimensional gaussian kernel. */
01116          dot,            /* Dot product summing variable. */
01117          sum;            /* Sum of the kernel weights variable. */
01118 
01119    /****************************************************************************
01120    * Create a 1-dimensional gaussian smoothing kernel.
01121    ****************************************************************************/
01122    make_gaussian_kernel(sigma, &kernel, &windowsize);
01123    center = windowsize / 2;
01124 
01125    /****************************************************************************
01126    * Allocate a temporary buffer image and the smoothed image.
01127    ****************************************************************************/
01128    if((tempim = (float *) calloc(rows*cols, sizeof(float))) == NULL){
01129       fprintf(stderr, "Error allocating the buffer image.\n");
01130       exit(1);
01131    }
01132    if(((*smoothedim) = (short int *) calloc(rows*cols,
01133          sizeof(short int))) == NULL){
01134       fprintf(stderr, "Error allocating the smoothed image.\n");
01135       exit(1);
01136    }
01137 
01138    /****************************************************************************
01139    * Blur in the x - direction.
01140    ****************************************************************************/
01141    for(r=0;r<rows;r++){
01142       for(c=0;c<cols;c++){
01143          dot = 0.0;
01144          sum = 0.0;
01145          for(cc=(-center);cc<=center;cc++){
01146             if(((c+cc) >= 0) && ((c+cc) < cols)){
01147                dot += (float)image[r*cols+(c+cc)] * kernel[center+cc];
01148                sum += kernel[center+cc];
01149             }
01150          }
01151          tempim[r*cols+c] = dot/sum;
01152       }
01153    }
01154 
01155    /****************************************************************************
01156    * Blur in the y - direction.
01157    ****************************************************************************/
01158    for(c=0;c<cols;c++){
01159       for(r=0;r<rows;r++){
01160          sum = 0.0;
01161          dot = 0.0;
01162          for(rr=(-center);rr<=center;rr++){
01163             if(((r+rr) >= 0) && ((r+rr) < rows)){
01164                dot += tempim[(r+rr)*cols+c] * kernel[center+rr];
01165                sum += kernel[center+rr];
01166             }
01167          }
01168          (*smoothedim)[r*cols+c] = (short int)(dot*BOOSTBLURFACTOR/sum + 0.5);
01169       }
01170    }
01171 
01172    free(tempim);
01173    free(kernel);
01174 }
01175 
01176 /*******************************************************************************
01177 * PROCEDURE: make_gaussian_kernel
01178 * PURPOSE: Create a one dimensional gaussian kernel.
01179 * NAME: Mike Heath
01180 * DATE: 2/15/96
01181 *******************************************************************************/
01182 void make_gaussian_kernel(float sigma, float **kernel, int *windowsize)
01183 {
01184    int i, center;
01185    float x, fx, sum=0.0;
01186 
01187    *windowsize = int(1 + 2 * ceil(2.5 * sigma));
01188    center = (*windowsize) / 2;
01189 
01190    if((*kernel = (float *) calloc((*windowsize), sizeof(float))) == NULL){
01191       fprintf(stderr, "Error callocing the gaussian kernel array.\n");
01192       exit(1);
01193    }
01194 
01195    for(i=0;i<(*windowsize);i++){
01196       x = (float)(i - center);
01197       fx = pow(2.71828, -0.5*x*x/(sigma*sigma)) / (sigma * sqrt(6.2831853));
01198       (*kernel)[i] = fx;
01199       sum += fx;
01200    }
01201 
01202    for(i=0;i<(*windowsize);i++) (*kernel)[i] /= sum;
01203 
01204 }
Generated on Sun May 8 08:40:20 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3