digimon.cpp

Package [view]: cutter.zip
Upload User: happy94
Upload Date: 2009-01-08
Package Size: 28k
Code Size: 112k
Category: Special Effects
Development Platform: C++
  1. /***************************************************************************
  2.  * digimon.c
  3.  * Note: Images are oriented with pixel (0,0) in the top left corner.
  4.  ***************************************************************************/
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <string.h>
  8. #include <math.h>
  9. #include <iostream>
  10. #include <GL/glut.h>
  11. #include <limits.h>    /* This is where SHRT_MAX is defined */
  12. #include <time.h>
  13. #define TRUE  1
  14. #define FALSE 0
  15. #ifndef DBL_MAX
  16. #define DBL_MAX         1.7976931348623157E+308
  17. #endif
  18. //#define PIE 3.14159265358979323846
  19. #define PIE 3.141592653589
  20. #define TWO_PIE (2*PIE)
  21. //6.2831853072
  22. #define MAX_NUM_CLASSES 10
  23. #define HIST_DIM 256
  24. #define FEATURE_SIZE 128
  25. #define MAX_SIGNAL_SIZE 65536
  26. #define X_AXIS 1
  27. #define Y_AXIS 2
  28. #define SINE 1
  29. #define COSINE 2
  30. #define GRY 0
  31. #define RED 1
  32. #define GRN 2
  33. #define BLU 3
  34. #define ABS(a) ((a)>(-a)?(a):(-a))
  35. #define DBL_TOL .0000001
  36. #define MAX_NUM_POINTS 100000
  37. #define MAX_NUM_SEEDS 30
  38. #define WAVE_SIZE 256
  39. #define WAVE_MOD 0xff
  40. #define BOUNDARY_PIXEL 255
  41. #define SHAPE_REP_SIZE 64
  42. #define NUM_GROUPS 4
  43. class Image;
  44. class Complex;
  45. class Cutter;
  46. //GLOBAL VARAIBLES
  47. int cut_xres = 1, cut_yres = 1;
  48. Image *cut_image = NULL;
  49. Cutter *cutter = NULL;
  50. int g_delay;
  51. int button_down = FALSE;
  52. //Function prototypes
  53. void display();
  54. void getNeighbor4(int x, int y, int *nx, int *ny, int id);
  55. void getNeighbor8(int x, int y, int *nx, int *ny, int id);
  56. void getColorFromID(int id, unsigned char *r,  unsigned char *g,  unsigned char *b );
  57. void convolve1D(double *s1, int size1,
  58.                 double *s2, int size2,
  59.                 double *s3, int *size3);
  60. Image* convolve2D(Image *img1, Image *img2);
  61. Image* sobelFilter(Image *img2);
  62. void medianFilter1D(int range, double *s1, int size1, double *s2, int *size2);
  63. void enhanceEdge1D (double *s1, int size1,
  64.                     double *s2, int size2,
  65.                     double *s3, int *size3);
  66. void slowSort(double *a, int size);
  67. void threshold(double *a, int size, double val, double bot, double top);
  68. void threshBlock2D(double *a, int w, int h, int range, double val, double bot, double top, double *b);
  69. void FourierTransform1D(double *f, int size, Complex *F);
  70. void invFourierTransform1D(Complex *F, int size, double *f);
  71. void FFT1D(double *f, int size, Complex *F);
  72. void FFT1Drec(Complex *F, int size, int mid, int dir);
  73. void FFT2D(double *f, int size, Complex *F);
  74. void invFFT2D(Complex *F, int size, double *f);
  75. void invFFT1D(Complex *F, int size, double *f);
  76. int  findPowerOf2(int size);
  77. void createSinusoid(int size, double w, double *s, double mag, int flag);
  78. void createSinusoid2D(int size, double freq, double *s, double mag, int axis, int flag);
  79. void createSobel(int size, double *sobel);
  80. void plot1D(double *a, int size, char *plot_name, char *filename);
  81. void plotComplex1D(Complex *a, int size, char *plot_name, char *filename);
  82. //Complex Number stuff
  83. Complex add(Complex c1, Complex c2);
  84. Complex sub(Complex c1, Complex c2);
  85. Complex mult(Complex c1, Complex c2);
  86. Complex div(Complex c1, Complex c2);
  87. void add(double *s1, int size1, double *s2, int size2, double *s3, int *size3);
  88. void mult(double *s1, int size1, double *s2, int size2, double *s3, int *size3);
  89. void mult(Complex *s1, int size1, Complex *s2, int size2, Complex *s3, int *size3);
  90. void div(Complex *s1, int size1, Complex *s2, int size2, Complex *s3, int *size3);
  91. //More functions
  92. int  boxThreshold(Image *im, Image *out_im);
  93. void copyDouble(double *dest, double *src, int size);
  94. void dilate(double *in, int w, int h, double *out, double *filter, int size, double tol);
  95. void erode(double *in, int w, int h, double *out, double *filter, int size, double tol);
  96. double getAverage(double *s, int size, double zero_value, double tol);
  97. double getPercentValue(double *s, int size, double val, double tol);
  98. double getDiff(double *s1, double *s2, int size);
  99. double getAveDiff(double *s1, double *s2, int size);
  100. void invert(Complex *s1, int size1);
  101. void rotate2D(double *s, int n, double *d, double theta);
  102. void swapMagPhase(Complex *c1, Complex *c2, int size);
  103. void removeInterference(Complex *D_Data, int size);
  104. void enhanceEdge2D (Complex *c1, int size1,
  105.                     Complex *c2, int size2,
  106.                     Complex *c3, int *size3);
  107. void medianFilter2D(int range, double *s1, int w, int h, double *s2);
  108. void normalize(double *d, int size);
  109. void writeUB(double *s, int size, char *filename);
  110. void butterFilter(Complex *c, int size, double cut_off, double n);
  111. void filterFourier2D(Complex *c, int size, int begin, int end);
  112. //Other Functions
  113. void fwrite_long(unsigned long number, FILE *stream);
  114. void fwrite_short(unsigned short number, FILE *stream);
  115. void write_wav(char *filename, double *sample, unsigned long sample_len, unsigned long freq);
  116. double *create_wav(int sample_len, double freq_scale, double noise_frac);
  117. //Classes
  118. //Classes
  119. //Classes
  120. //Classes
  121. class Complex{
  122.   public:
  123.     double r, i;
  124.     double getMag(){ return sqrt(r*r + i*i); };
  125.     double getPhase(){
  126.       if (i < DBL_TOL || r < DBL_TOL){
  127.         return atan2(0.0, 0.0);
  128.       }
  129.       return atan2(i, r);
  130.     };
  131. };
  132. class Image{
  133.   public:
  134.     unsigned char *data;
  135.     int height, width;
  136.     int type, max;
  137.     Image();
  138.     Image(Image *from_img);
  139.     Image(int w, int h, int type);
  140.     Image(double *signal, int size);
  141.     ~Image();
  142.     unsigned char colorShapes(unsigned char tol, int min_area);
  143.     void copy(Image *img);
  144.     void createInfoImage(int num_shapes);
  145.     void clear(unsigned char val);
  146.     void drawLine(int x0, int y0, int x1, int y1, unsigned char r, unsigned char g, unsigned char b);
  147.     void diff(Image *img);
  148.     void dilate(int size);
  149.     void erode(int size);
  150.     void predCode1(int flag);
  151.     void predCode2(int flag);
  152.     void equalize();
  153.     void flip();
  154.     void getCenterOfGravity(int *x, int *y, unsigned char tol);
  155.     int markCenterOfGravityColor(int *cx, int *cy, unsigned char i);
  156.     double markOrientation(int i, int cy, int cy);
  157.     Image* getHistImage(int flag);
  158.     Image* getEqualizeTableImage();
  159.     void getPDF(double *pdf, int n, double *max_val);
  160.     void getShapeColor(unsigned char *r, unsigned char *g, unsigned char *b, unsigned char *br, unsigned char *bg, unsigned char *bb, unsigned char tol);
  161.     double *getDouble(double offset);
  162.     void getItemBounds(unsigned char id, int *bot_x, int *bot_y, int *top_x, int *top_y);
  163.     unsigned char getBoundAverage(int bot_x, int bot_y, int top_x, int top_y);
  164.     void unionThreshRegion(Image *img, int bot_x, int bot_y, int top_x, int top_y, int range, unsigned char t);
  165.     void invert();
  166.     void makePowerOf2(int new_size);
  167.     void makeColor();
  168.     void makeGray();
  169.     void pad(int new_width, int new_height);
  170.     int read(char *filename);
  171.     void rotate(double deg);
  172.     void makeSinusoid(int size, double freq, double mag, int axis);
  173.     void makeMag(Complex *D_Data, int size);
  174.     void makePhase(Complex *D_Data, int size);
  175.     void mask(Image *m);
  176.     void signalToImage2D(double *d_data, int w, int h);
  177.     void swapQuads();
  178.     void threshold(unsigned char val);
  179.     int write(char *filename);
  180.     void seedFill(int x, int y,
  181.                   unsigned char r, unsigned char g, unsigned char b,
  182.                   unsigned char f_r, unsigned char f_g, unsigned char f_b, 
  183.                   unsigned char tol, int *area);
  184.   private:
  185. };
  186. //Cutter stuff
  187. class GraphNode{
  188.   /*
  189.    * 0,0   ->   width
  190.    *
  191.    * |
  192.    * v
  193.    * height
  194.    *
  195.    * Parent specification
  196.    *  3 2 1
  197.    *  4 - 0
  198.    *  5 6 7
  199.    */
  200.   public:
  201.     unsigned int cost;
  202.     unsigned char parent;
  203.     int total_cost;
  204.     char expanded;
  205.     int x, y;
  206.     GraphNode *next;
  207.     GraphNode(){
  208.       x = y = 0;
  209.       total_cost = 0;
  210.       cost = 0;
  211.       parent = 8;
  212.       expanded = FALSE;
  213.       next = NULL;
  214.     }
  215. };
  216. class PointList{
  217.   public:
  218. //  int curr_point;
  219.   int id;
  220.   int num_points;
  221.   int num_seeds;
  222.   int cx, cy;
  223.   int x[MAX_NUM_POINTS];
  224.   int y[MAX_NUM_POINTS];
  225.   double r[SHAPE_REP_SIZE];//for shape representation.
  226.   int seed[MAX_NUM_SEEDS];//a list of all the points in the x, y arrays that are
  227.                           //seeds.
  228.   PointList *next;
  229.   PointList(){
  230.     id = 0;
  231.     num_points = 0;
  232.     num_seeds = 0;
  233.     cx = cy = 0;
  234.     next = NULL;
  235.   }
  236.   ~PointList(){
  237.     if (next != NULL) delete next;
  238.   }
  239.   void createShapeRep();
  240.   double getDist(PointList *plist);
  241.   double getDist(double *r2);
  242.   void writeShapeRep(char *prefix);
  243. };
  244. class ContourMaker{
  245.   public:
  246.   PointList *plist_head;
  247.   int num_plist;
  248.   ContourMaker(Image *img){
  249.     Image im2;
  250.     plist_head = NULL;
  251.     createBoundaries(img);
  252.     //drawBoundaries(img);
  253.     matchBoundaries();
  254.     im2.copy(img);
  255.     im2.clear(0);
  256.     drawBoundaries(&im2);
  257.     im2.write("out_groups.ppm");
  258.     num_plist = 0;
  259.   }
  260.   ~ContourMaker(){
  261.     if (NULL != plist_head) delete plist_head;
  262.   }
  263.   //TOBECONT
  264.   void createBoundaries(Image *img);
  265.   void matchBoundaries();
  266.   void drawBoundaries(Image *img, PointList *plist);
  267.   void drawBoundaries(Image *img);//draws the boundaries onto this image
  268.                                   //the color of the boundary will be the
  269.                                   //group color.
  270. };
  271. class Cutter{
  272.   public:
  273.   GraphNode *graph;
  274.   GraphNode* wavefront[WAVE_SIZE];
  275.   PointList plist;
  276.   int best_node;
  277.   int wave_done;
  278.   int height, width;
  279.   int started;
  280.   Cutter(){ started = FALSE; height = width = 0; best_node = 0; wave_done = 0;}
  281.   Cutter(Image *img);
  282.   ~Cutter();
  283.   void addSeedPoint();
  284.   void addToWaveFront(int x, int y);
  285.   void clearGraph();
  286.   void deleteWaveFront();
  287.   void dumpCostGraph();
  288.   void dumpPaths();
  289.   void dumpWaveFront();
  290.   int  findInWaveFront(int x, int y);
  291.   int getCost(int x1, int y1, int x2, int y2);
  292.   void initCostGraph(Image *img);//initializes all the costs for the nodes in the graph;
  293.   void initWaveFront(int seed_x, int seed_y);
  294.   void removeSeedPoint();
  295.   void removeFromWaveFront(int x, int y);
  296.   void saveBoundary(int x, int y);
  297.   int  updateWaveFront();
  298.   void updatePointList(int x, int y);//add points to connect from last seed to
  299.                                      //this point in the graph
  300. };
  301. //POINTLIST
  302. //POINTLIST
  303. //POINTLIST
  304. //POINTLIST
  305. //Create the shape rep for each point list.
  306. void PointList::createShapeRep(){
  307.   double d_j, inc;
  308.   int i, j;
  309.   double max;
  310.   int dx, dy;
  311.   char fname[255];
  312.   //Use the R-S representation
  313.   //Radius as function of arc length.
  314.   //TOBECONT
  315.   if (num_points < SHAPE_REP_SIZE){
  316.     printf("Not enough points for my to be a valid shape for me.n");
  317.     return;
  318.   }
  319.   inc = (double)num_points/(double)SHAPE_REP_SIZE;
  320.   max = 0.0;
  321.   for (i=0; i<SHAPE_REP_SIZE; i++){
  322.     d_j = (double)i * inc;
  323.     j = (int) d_j;
  324. //    printf(" (i: %d, j: %d) ",i, j);
  325.     dx = x[j] - cx;
  326.     dy = y[j] - cy;
  327.     r[i] = sqrt(dx*dx + dy*dy);
  328.     if (max < r[i]) max = r[i];
  329.   }
  330. //  printf("num_points: %dn", num_points);
  331.   //Normalize the shape rep
  332.   max = 1.0/max;
  333.   for (i=0; i<SHAPE_REP_SIZE; i++){
  334.     r[i]*=max;
  335.   }
  336. //  writeShapeRep("out_shape_rep");
  337. }
  338. //Do corrolation to find the best match against the shape
  339. double PointList::getDist(PointList *plist){
  340.   int i, j;
  341.   double min, dist;
  342.   double *r2;
  343. //  if (id == plist->id) return -1; //same shape, but one is inside &&
  344.                                  // the other out
  345.   min = SHAPE_REP_SIZE*SHAPE_REP_SIZE;
  346.   r2 = plist->r;
  347.   for (i=0; i<SHAPE_REP_SIZE; i++){
  348.     dist = 0;
  349.     for (j=0; j<SHAPE_REP_SIZE; j++){
  350.       dist += fabs( r[j] - r2[(i+j)%SHAPE_REP_SIZE]);
  351.     }
  352.     if (dist < min) min = dist;
  353.   }
  354.   return min;
  355. }
  356. //Compare just with a double array
  357. double PointList::getDist(double *r2){
  358.   int i, j;
  359.   double min, dist;
  360.                                  
  361.   min = SHAPE_REP_SIZE*SHAPE_REP_SIZE;
  362.   for (i=0; i<SHAPE_REP_SIZE; i++){
  363.     dist = 0;
  364.     for (j=0; j<SHAPE_REP_SIZE; j++){
  365.       dist += fabs( r[j] - r2[(i+j)%SHAPE_REP_SIZE]);
  366.     }
  367.     if (dist < min) min = dist;
  368.   }
  369.   return min;
  370. }
  371. //Write out an image of each shape rep
  372. void PointList::writeShapeRep(char *prefix){
  373.   int i, j;
  374.   char fname[255];
  375.   //TOBECONT
  376.   Image *im;
  377.   im = new Image(SHAPE_REP_SIZE*2, SHAPE_REP_SIZE*2, 6);
  378.   im->clear(0);
  379.   for (i=0; i<SHAPE_REP_SIZE-1; i++){
  380.     j = i + 1;
  381.     im->drawLine(i<<1, SHAPE_REP_SIZE*2 - (int)(r[i]*SHAPE_REP_SIZE*2.0), j<<1, SHAPE_REP_SIZE*2 - (int)(r[j]*SHAPE_REP_SIZE*2.0), 255, 50, 64);
  382.   }
  383.   sprintf(fname, "%s%d.ppm", prefix, num_points);
  384.   im->write(fname);
  385.   delete im;
  386. }
  387. //CONTOURMAKER
  388. //CONTOURMAKER
  389. //CONTOURMAKER
  390. //CONTOURMAKER
  391. //Use the over the shoulder algorithm by Dr. Barrett.
  392. void ContourMaker::createBoundaries(Image *img){
  393.   int y, x, ind, ind2, i, j;
  394.   int num_shapes;
  395.   int height, width;
  396.   unsigned char bb, dir;
  397.   unsigned char *data;
  398.   int nx, ny, currx, curry, done;
  399.   PointList *new_plist;
  400.   int totalx, totaly;
  401.   Image *img_ptr;
  402.   Image im;
  403.   im.copy(img);
  404.   img_ptr = &im;
  405.   img_ptr->makeColor();
  406.   height = img_ptr->height;
  407.   width = img_ptr->width;
  408.   //First color the shapes a unique color
  409.   num_shapes = img_ptr->colorShapes(0, 10);
  410.   data = img_ptr->data;
  411.   //Get the background color
  412.   bb = 0;//data[2];
  413.   //TOBECONT
  414.   //Now find the contour
  415.   for(y=0; y<height; y++){
  416.     for(x=0; x<width-1; x++){
  417.       ind = (x+y*width)*3;
  418.       ind2 = (x+1+y*width)*3;
  419.       //We check the pixel after the current position
  420.       if ( data[ind+2] == bb &&
  421.            data[ind2+2] != bb &&
  422.            data[ind2+2] != BOUNDARY_PIXEL ){
  423.  //       printf("Boundary start: (%d, %d) curr:%d n:%dn", x, y, (int)data[ind+2], (int) data[ind2 +2]);
  424.         //We have found a shape, so make a contour
  425.         currx = x;
  426.         curry = y;
  427.         dir = 0;
  428.         done = FALSE;
  429.         totalx = totaly = 0;
  430.         new_plist = new PointList();
  431.         new_plist->next = plist_head;
  432.         plist_head = new_plist;
  433.         new_plist->id = data[ind2+2];
  434.         num_plist ++;
  435.         for (j=0; TRUE != done && j<100000; j++){
  436.           getNeighbor4(currx, curry, &nx, &ny, dir);
  437.           //While pixel in front is a region pixel, turn left
  438.           for(i=0;
  439.              (data[(nx + ny*width)*3+2] != bb &&
  440.               data[(nx + ny*width)*3+2] != BOUNDARY_PIXEL) && i<4;
  441.               i++){
  442.             dir++;
  443.             dir &= 3;
  444.             getNeighbor4(currx, curry, &nx, &ny, dir);
  445.           }
  446.           if (4==i){
  447.             printf("ERROR in contourMaker::createBoundaries(). Too many rotationsn");
  448.             return;
  449.           }
  450.           //Mark the current position as a boundary pixel
  451.           data[(currx + curry*width)*3+2] = BOUNDARY_PIXEL;
  452.           data[(currx + curry*width)*3] = 0;
  453.           data[(currx + curry*width)*3+1] = 0;
  454.           new_plist->x[new_plist->num_points] = currx;
  455.           new_plist->y[new_plist->num_points] = curry;
  456.           totalx += currx;
  457.           totaly += curry;
  458.           new_plist->num_points ++;
  459. //          printf("curr: (%d, %d)  nx: (%d, %d)n", currx, curry, nx, ny);
  460.           //Step
  461.           currx = nx;
  462.           curry = ny;
  463.           if ( currx == x && curry == y ){
  464.             done = TRUE;
  465.           }
  466.           else{
  467.             dir = (dir+3)&3;
  468.           }
  469.        //   if (j%10 == 0) display();  //For debug
  470.         } 
  471.         new_plist->cx = totalx/new_plist->num_points;
  472.         new_plist->cy = totaly/new_plist->num_points;
  473.         new_plist->createShapeRep();
  474.       }
  475.     }
  476.   }
  477.  
  478.   img_ptr->write("out_colored.pgm");
  479. }
  480. void ContourMaker::matchBoundaries(){
  481.   PointList *scan_i = NULL;
  482.   PointList *best_i = NULL;
  483. //  PointList **plist_ptr;
  484.   int i, j, k, best_g;
  485.   double group[NUM_GROUPS][SHAPE_REP_SIZE];
  486.   PointList pl[NUM_GROUPS];
  487.   double c, min, max_cost;
  488.   //TOBECONT
  489.   //Start by finding the most dissimilar shapes
  490.   min = SHAPE_REP_SIZE;
  491.   //Add in the first group to the list
  492.   copyDouble(group[0], plist_head->r, SHAPE_REP_SIZE);
  493.   for (k=1; k<NUM_GROUPS; k++){
  494.     best_i = plist_head;
  495.     max_cost = 0;
  496.     for (scan_i=plist_head->next; NULL != scan_i; scan_i = scan_i->next){
  497.       //compare against existing groups
  498.       c=0;
  499.       for (j=0; j<k; j++){
  500.         c += scan_i->getDist(group[j]);
  501.       }
  502.       if (c > max_cost){
  503.         best_i = scan_i;
  504.         max_cost = c;
  505.       }
  506.     }
  507.     copyDouble(group[k], best_i->r, SHAPE_REP_SIZE);
  508.   }
  509.   //Now assign plist to a group
  510.   for (scan_i=plist_head; NULL != scan_i; scan_i = scan_i->next){
  511.     best_g = 0;
  512.     min = SHAPE_REP_SIZE;
  513.     for (j=0; j<NUM_GROUPS; j++){
  514.       c = scan_i->getDist(group[j]);
  515.       if (c < min){
  516.         min = c;
  517.         best_g = j;
  518.       }
  519.     }
  520.     scan_i->id  = best_g;
  521.   }
  522.   for (k=0; k<NUM_GROUPS; k++){
  523.     copyDouble(pl[k].r, group[k], SHAPE_REP_SIZE);
  524.     pl[k].num_points = k;
  525. //    pl[k].writeShapeRep("out_group_");
  526.   }
  527.   //Now we have a cost graph we can use to match up the shapes with.
  528. //  delete [] cost;
  529. //  delete [] has_hole;
  530. //  delete [] plist_ptr;
  531. }
  532. void ContourMaker::drawBoundaries(Image *img){
  533.   unsigned char *data;
  534.   //int h = img->height;
  535.   int w = img->width;
  536.   PointList *scan_i;
  537.   int *x, *y;
  538.   int i;
  539.   unsigned char r, g, b;
  540.   img->makeColor();
  541.   data = img->data;
  542.   for (scan_i = plist_head; NULL != scan_i; scan_i = scan_i->next){
  543.     getColorFromID(scan_i->id, &r, &g, &b);
  544.     x = scan_i->x;
  545.     y = scan_i->y;
  546.     for (i=0; i<scan_i->num_points; i++){
  547.       data[(x[i] + y[i]*w)*3] = r;
  548.       data[(x[i] + y[i]*w)*3+1] = g;
  549.       data[(x[i] + y[i]*w)*3+2] = b;
  550.     } 
  551.   }
  552. }
  553. void ContourMaker::drawBoundaries(Image *img, PointList *plist){
  554.   unsigned char *data;
  555.   //int h = img->height;
  556.   int w = img->width;
  557.   int *x, *y;
  558.   int i;
  559.   unsigned char r, g, b;
  560.   r = 0; g = 0; b = 255;
  561.   img->makeColor();
  562.   data = img->data;
  563.   x = plist->x;
  564.   y = plist->y;
  565.   for (i=0; i<plist->num_points; i++){
  566.     data[(x[i] + y[i]*w)*3] = r;
  567.     data[(x[i] + y[i]*w)*3+1] = g;
  568.     data[(x[i] + y[i]*w)*3+2] = b;
  569.   } 
  570. }
  571. //CUTTER
  572. //CUTTER
  573. //CUTTER
  574. //CUTTER
  575. Cutter::Cutter(Image *img){
  576.   int i, j;
  577.   //int l, m;
  578.   img->makeGray();
  579.   height = img->height;
  580.   width = img->width;
  581.   graph = new GraphNode[height*width];
  582.   for (i=0; i<WAVE_SIZE; i++){
  583.     wavefront[i] = NULL;
  584.   }
  585.   for (j=0; j<height; j++){
  586.     for (i=0; i<width; i++){
  587.       graph[i + j*width].x = i;
  588.       graph[i + j*width].y = j;
  589.     }
  590.   }
  591.   best_node = 0;
  592.   wave_done = 0;
  593.   started = FALSE;
  594.   initCostGraph(img);
  595. }
  596. Cutter::~Cutter(){
  597.   delete [] graph;
  598. }
  599. void Cutter::addSeedPoint(){
  600.   plist.seed[plist.num_seeds] = plist.num_points-1;
  601.   plist.num_seeds ++;
  602. }
  603. void Cutter::removeSeedPoint(){
  604.   int seed_id;
  605.   
  606.   if (plist.num_seeds <= 0) return;
  607.   plist.num_seeds --;
  608.   if (plist.num_seeds == 0){
  609.     plist.num_points = 0;
  610.     started = FALSE;
  611.   }
  612.   else{
  613.     //Initialize the wave front for the previous seed
  614.     seed_id = plist.seed[plist.num_seeds-1];
  615.     initWaveFront(plist.x[seed_id], plist.y[seed_id]);
  616.     updateWaveFront();
  617.   }
  618. }
  619. void Cutter::removeFromWaveFront(int x, int y){
  620.   int done;
  621.   int loc;
  622.   GraphNode *prev, *curr, *next;
  623.   done = FALSE;
  624.   loc = (x + y*width) & WAVE_MOD;
  625.   for (prev = NULL, curr = wavefront[loc];
  626.        NULL != curr && !done;
  627.        prev = curr, curr = next){
  628.     next = curr->next;
  629.     if (curr->x == x && curr->y == y){
  630.       //Then remove it!
  631.       done = TRUE;
  632.       if (prev==NULL) wavefront[loc] = curr->next;
  633.       else prev->next = curr->next;
  634.       curr->next = NULL;
  635.     }
  636.   }
  637. }
  638. void Cutter::saveBoundary(int x, int y){
  639.   int    i;
  640.   FILE   *fpt;
  641.   if (!started) return;
  642.   if (plist.num_points == 0) return;
  643.   //Add a seed point here (x, y)
  644.   updatePointList(x, y);
  645.   addSeedPoint();
  646.   initWaveFront(x, y);
  647.   while( updateWaveFront() );
  648.   //Update the point list from the original seed point to this point.
  649.   updatePointList(plist.x[0], plist.y[0]);
  650.   //Save out the boundary minus the duplicated original point.
  651.   printf("saving "bound.txt"n");
  652.   fpt = fopen( "bound.txt", "wb" );
  653.   if ( !fpt ){
  654.     printf("Error: Could not create file: bound.txtn");
  655.     return;
  656.   }
  657.   plist.num_points--;//Drop the last duplicate point
  658.   fprintf( fpt, "%dn", plist.num_points);
  659.   for (i=0; i<plist.num_points; i++){
  660.     fprintf( fpt, "%d %dn", plist.x[i], plist.y[i]);
  661.   }
  662.   fclose( fpt );
  663.   plist.num_points ++;
  664.   removeSeedPoint();
  665. }
  666. void Cutter::addToWaveFront(int x, int y){
  667.   int total_cost;
  668.   int loc;
  669.   loc = x + y*width;
  670. //  printf("DBG002 addToWaveFront()n");
  671.   total_cost = graph[loc].total_cost;
  672.   //Add the new wave to the wave front
  673.   graph[loc].next = wavefront[ total_cost & WAVE_MOD ];
  674.   wavefront[total_cost & WAVE_MOD] = &graph[loc];
  675. }
  676. //Use the graph to add points from the last seed point to this point.
  677. void Cutter::updatePointList(int px, int py){
  678.   int parent, i, j;
  679.   int curr_x, curr_y;
  680.   int seed_x, seed_y, seed_id;
  681.   //This is a list that will be flipped later
  682.   int x[MAX_NUM_POINTS];
  683.   int y[MAX_NUM_POINTS];
  684.   //Need to make sure the wave front has been expanded enough to handle this
  685.   //point
  686.   if (!started) return;
  687.   if (plist.num_seeds < 1) return;
  688.   if (graph[px + py*width].expanded != TRUE) return;
  689.   //The last seed in the list is the seed whose wavefront is on the graph.
  690.   seed_id = plist.seed[plist.num_seeds-1];
  691.   seed_x = plist.x[seed_id];
  692.   seed_y = plist.y[seed_id];
  693.   i=0;
  694.   curr_x = px;
  695.   curr_y = py;
  696.   while((curr_x != seed_x || curr_y != seed_y) && i<(MAX_NUM_POINTS-seed_id-2)){
  697.     //Add point to the list
  698.     x[i] = curr_x;
  699.     y[i] = curr_y;
  700.     i++;
  701.     parent = graph[curr_x + curr_y*width].parent;
  702. //    printf("curr (%d, %d) seed (%d, %d) p:%d i:%dn", curr_x, curr_y, seed_x, seed_y, parent, i);
  703.     if (parent == 9){
  704.       printf("Error in Cutter::updatePointList(). seed found at (%d, %d) instead of (%d, %d)n", curr_x, curr_y, seed_x, seed_y);
  705.       curr_x = seed_x;
  706.       curr_y = seed_y;
  707.       //This is the seed point, end here.
  708.     }
  709.     else{
  710.       getNeighbor8(curr_x, curr_y, &curr_x, &curr_y, parent);
  711.     }
  712.   }
  713.   if (i>= MAX_NUM_POINTS-seed_id-2){
  714.     printf("Error in Cutter::updatePointList(). Not enough memory for pointlistn");
  715.   }
  716.   for (j=seed_id+1, i--; i>=0; i--, j++){
  717.     plist.x[j] = x[i];
  718.     plist.y[j] = y[i];
  719.   }
  720.   plist.num_points = j;
  721. //  printf("DBG006 MAX: %d  num_points: %dn", MAX_NUM_POINTS, j);
  722. }
  723. //Advance the wavefront forward
  724. int Cutter::updateWaveFront(){
  725.   int i, j;
  726.   GraphNode *wave;
  727.   int wave_x, wave_y, curr_x, curr_y;
  728.   int wave_loc, curr_loc;
  729.   int cost;
  730.   if (wave_done) return FALSE;
  731.   if (!started) return TRUE;
  732.   //First find the lowest cost node on the wavefront
  733.   //Advance the best_node pointer
  734.   for (j=0; j<WAVE_SIZE && wavefront[best_node] == NULL; j++){
  735.     best_node ++;
  736.     best_node &= WAVE_MOD;
  737.   }
  738.   if (j == WAVE_SIZE){
  739.     wave_done = TRUE;
  740.     return FALSE;
  741.   }
  742.   //Remove that node from the wavefront
  743.   wave = wavefront[best_node];
  744.   wavefront[best_node] = wave->next;
  745.   wave->next = NULL;
  746.   wave_x = wave->x;
  747.   wave_y = wave->y;
  748.   wave_loc = wave_x + wave_y*width;
  749.   //Now expand that node to the next level out
  750.   graph[wave_loc].expanded = TRUE;
  751.   //Look at all neighbors and Set and queue the best one not expanded
  752.   //Now put all the neighbors onto the wavefront.
  753.   for (i=0; i<8; i++){
  754.     getNeighbor8(wave_x, wave_y, &curr_x, &curr_y, i);
  755.     if (curr_x >= width||curr_x < 0 || curr_y >= height||curr_y < 0) continue;
  756.     curr_loc = curr_x  + curr_y*width;
  757.     if (graph[curr_loc].expanded == TRUE) continue;
  758.     cost = graph[wave_loc].total_cost + getCost(wave_x, wave_y, curr_x, curr_y);
  759.     if ( findInWaveFront(curr_x, curr_y) )
  760.      if ( cost < graph[curr_loc].total_cost )
  761.        removeFromWaveFront(curr_x, curr_y);
  762.     if ( !findInWaveFront(curr_x, curr_y) ){
  763.       graph[curr_loc].total_cost = cost;
  764.       graph[curr_loc].parent = (i+4)%8;
  765.       addToWaveFront(curr_x, curr_y);
  766.     }
  767.   }
  768.   return TRUE;
  769. }
  770. void Cutter::initCostGraph(Image *img){
  771.   int i=0, j;
  772.   int val1, val2, val3;
  773.   double max;
  774.   double cost;
  775.   unsigned char *data;
  776.   int p1, p2, p3;
  777.   int *lap;
  778.   Image lap2;
  779.   data = img->data;
  780.   if (height != img->height){
  781.     printf("ERROR: Cutter::initCostGraph()n");
  782.     return;
  783.   }
  784.   if (width != img->width){
  785.     printf("ERROR: Cutter::initCostGraph()n");
  786.     return;
  787.   }
  788.   max = 0;
  789.   for (j=0; j<height; j++){
  790.     if (j>0 && j < height-1){
  791.       for (i=0; i<width; i++){
  792.         if (i>0 && i<width-1){
  793.           //Calculate the gradient magnitude using the sobel filter
  794.           val1 =   - (int)data[i-1 + (j-1)*width] 
  795.                    - (int)data[i-1 + j*width]*2 
  796.                    - (int)data[i-1 + (j+1)*width] + 
  797.                    (int)data[i+1 + (j-1)*width] + 
  798.                    (int)data[i+1 + j*width]*2 + 
  799.                    (int)data[i+1 + (j+1)*width] ;
  800.           val2 =   - (int)data[i-1 + (j-1)*width]
  801.                    - (int)data[i +   (j-1)*width]*2
  802.                    - (int)data[i+1 + (j-1)*width] + 
  803.                    (int)data[i-1 + (j+1)*width] + 
  804.                    (int)data[i +   (j+1)*width]*2 + 
  805.                    (int)data[i+1 + (j+1)*width] ;
  806.           cost = sqrt(val1*val1 + val2*val2);
  807.           if (max < cost) max = cost;
  808.           graph[i + j*width].cost = (int)cost;
  809.         }
  810.         else{ graph[i + j*width].cost = 0;}
  811.       }
  812.     }
  813.     else{ graph[i + j*width].cost = 0;}
  814.   }
  815.   max = 128.0/max;
  816.   for (j=0; j<height; j++){
  817.     for (i=0; i<width; i++){
  818.       //Invert the cost graph
  819.       graph[i + j*width].cost = 128 - (int)((double)graph[i + j*width].cost * max);
  820.       if (graph[i + j*width].cost < 0){
  821.         printf("Error %d less than zeron", graph[i + j*width].cost);
  822.       }
  823.     }
  824.   }
  825.   //Now to calculate the laplacian stuff.
  826.           //Now calculate the laplacian
  827.   lap = new int[width*height];
  828.   for (i=0; i<width*height; i++){ lap[i] = 0; }
  829.   
  830.   for (j=1; j<height-1; j++){
  831.     for (i=1; i<width-1; i++){
  832.       lap[i + j*width] =
  833.                   (int)data[i + j*width]*4
  834.                  -(int)data[i +  (j-1)*width]
  835.                  -(int)data[i +  (j+1)*width]
  836.                  -(int)data[i-1 + j*width]
  837.                  -(int)data[i+1 + j*width];
  838.     }
  839.   }
  840.   lap2.copy(img);
  841.   lap2.clear(127);
  842.   //Now find the zero crossings.
  843.   for (j=0; j<height-1; j++){
  844.     for (i=0; i<width-1; i++){
  845.       p1 = lap[i +j*width];
  846.       p2 = lap[i+1 +j*width];
  847.       p3 = lap[i +(j+1)*width];
  848.       if (p1 == 0){
  849.         lap2.data[i+j*width] = 0;//zero additional cost, because of laplacian
  850.       }
  851.       else if (p1 < 0){
  852.         if (p2>0){
  853.           //zero crossing found
  854.           if (abs(p1) < abs(p2)){ lap2.data[i+j*width] = 0; }
  855.           else{ lap2.data[i+1+j*width] = 0; }
  856.         }
  857.         if (p3>0){
  858.           //zero crossing found
  859.           if (abs(p1) < abs(p3)){ lap2.data[i+j*width] = 0; }
  860.           else{ lap2.data[i+(j+1)*width] = 0; }
  861.         }
  862.       }
  863.       else {
  864.         if (p2<0){
  865.           //zero crossing found
  866.           if (abs(p1) < abs(p2)){ lap2.data[i+j*width] = 0; }
  867.           else{ lap2.data[i+1+j*width] = 0; }
  868.         }
  869.         if (p3<0){
  870.           //zero crossing found
  871.           if (abs(p1) < abs(p3)){ lap2.data[i+j*width] = 0; }
  872.           else{ lap2.data[i+(j+1)*width] = 0; }
  873.         }
  874.         //p1 > 0
  875.       }
  876.     }
  877.   }
  878.   lap2.write("out_lap2.pgm");
  879.   data = lap2.data;
  880.   for (j=0; j<height; j++){
  881.     for (i=0; i<width; i++){
  882.       //Invert the cost graph
  883.       graph[i + j*width].cost += (int)data[i +j*width];//128 - (int)((double)graph[i + j*width].cost * max);
  884.     }
  885.   }
  886.   delete [] lap;
  887.           //Add in the laplacian cost
  888.           //cost += abs(val3);
  889. }
  890. void Cutter::initWaveFront(int seed_x, int seed_y){
  891.   int loc;
  892.   loc = seed_x + seed_y*width;
  893.   deleteWaveFront();
  894.   clearGraph();
  895.   wave_done = FALSE;
  896.   started = TRUE;
  897.   //Add the first node the the wavefront
  898.   wavefront[0] = &graph[loc];
  899.   graph[loc].total_cost = 0;
  900.   best_node = 0;
  901.   graph[loc].parent = 9;
  902.   updateWaveFront();
  903. }
  904. void Cutter::clearGraph(){
  905.   int i;
  906.   for(i=0; i<width*height; i++){
  907.     graph[i].expanded = 0;
  908.     graph[i].total_cost = 0;
  909.     graph[i].parent = 8;
  910.     graph[i].next = NULL;
  911.   }
  912. }
  913. void Cutter::deleteWaveFront(){
  914.   int i;
  915.   for (i=0; i<WAVE_SIZE; i++){
  916.     wavefront[i] = NULL;
  917.   }
  918.   best_node = -1;
  919.   wave_done = FALSE;
  920. }
  921. void Cutter::dumpCostGraph(){
  922.   int i, j;
  923.   Image *im;
  924.   im = new Image(width, height, 5);
  925.   for (j=0; j<height; j++){
  926.     for (i=0; i<width; i++){
  927.       im->data[i + j*width] = graph[i + j*width].cost;
  928.     }
  929.   }
  930.   im->write("out_cost_graph.pgm");
  931.   delete im;
  932. }
  933. void Cutter::dumpWaveFront(){
  934.   int i, j;
  935.   Image *im;
  936.   double max_cost;
  937.   max_cost = 0;
  938.   for (j=0; j<height; j++){
  939.     for (i=0; i<width; i++){
  940.       if (graph[i + j*width].total_cost > max_cost)
  941.         max_cost = graph[i + j*width].total_cost;
  942.     }
  943.   }
  944.   max_cost = 255.0/max_cost;
  945.   im = new Image(width, height, 5);
  946.   for (j=0; j<height; j++){
  947.     for (i=0; i<width; i++){
  948.       im->data[i + j*width] = (unsigned char)(
  949.                                 (double)graph[i+j*width].total_cost*max_cost
  950.                               );
  951.     }
  952.   }
  953.   im->write("out_wavefront_image.pgm");
  954.   delete im;
  955. }
  956. void Cutter::dumpPaths(){
  957.   int j, i, parent;
  958.   int seed_x, seed_y;
  959.   int curr_x, curr_y, seed_id;
  960.   GraphNode *scan_wave;
  961.   seed_id = plist.seed[plist.num_seeds-1];
  962.   seed_x = plist.x[ seed_id ];
  963.   seed_y = plist.y[ seed_id ];
  964.   for (j=0; j<WAVE_SIZE; j++){
  965.     for (scan_wave = wavefront[j];
  966.          NULL != scan_wave;
  967.          scan_wave = scan_wave->next){
  968.       curr_x = scan_wave->x;
  969.       curr_y = scan_wave->y;
  970.       i=0;
  971.       while(curr_x != seed_x && curr_y != seed_y && i<(MAX_NUM_POINTS-seed_id-2)){
  972.         i++;
  973.         parent = graph[curr_x + curr_y*width].parent;
  974.         if (parent == 9){
  975.           printf("Error in dumpPaths()n");
  976.           curr_x = seed_x;
  977.           curr_y = seed_y;
  978.           //This is the seed point, end here.
  979.         }
  980.         else{
  981.           getNeighbor8(curr_x, curr_y, &curr_x, &curr_y, parent);
  982.         }
  983.       }
  984.       if ( i==(MAX_NUM_POINTS-seed_id-2) ){
  985.          printf("nBAD PATH STARTn");
  986.         curr_x = scan_wave->x;
  987.         curr_y = scan_wave->y;
  988.         i=0;
  989.         while(curr_x != seed_x && curr_y != seed_y && i<(MAX_NUM_POINTS-seed_id-2)){
  990.           printf("(%d, %d)  ", curr_x, curr_y );
  991.           i++;
  992.           parent = graph[curr_x + curr_y*width].parent;
  993.           if (parent == 9){
  994.             printf("Error in dumpPaths()n");
  995.             curr_x = seed_x;
  996.             curr_y = seed_y;
  997.             //This is the seed point, end here.
  998.           }
  999.           else{
  1000.             getNeighbor8(curr_x, curr_y, &curr_x, &curr_y, parent);
  1001.           }
  1002.         }
  1003.         printf("nBAD PATH ENDn");
  1004.       }
  1005.     }
  1006.   }
  1007. }
  1008. int Cutter::findInWaveFront(int x, int y){
  1009.   int loc;
  1010.   GraphNode *scan_wave;
  1011.   loc = graph[x + y*width].total_cost & WAVE_MOD;
  1012.   for (scan_wave = wavefront[loc];
  1013.        NULL != scan_wave;
  1014.        scan_wave = scan_wave->next){
  1015.     if (scan_wave->x == x && scan_wave->y == y){
  1016.       return TRUE;
  1017.     }
  1018.   }
  1019.   return FALSE;
  1020. }
  1021. //Cost from going from x1, y1, to x2, y2.
  1022. //Generally just the value at x2, y2.
  1023. int Cutter::getCost(int x1, int y1, int x2, int y2){
  1024.   return graph[x2 + y2*width].cost;
  1025. }
  1026. //IMAGE
  1027. //IMAGE
  1028. //IMAGE
  1029. //IMAGE
  1030. Image::Image(){
  1031.   data = NULL;
  1032.   height = width = max = type = 0;
  1033. }
  1034. Image::Image(Image *from_img){
  1035.   data = NULL;
  1036.   height = width = max = type = 0;
  1037.   copy (from_img);
  1038. }
  1039. Image::Image(int w, int h, int t){
  1040.   type = t;
  1041.   height = h;
  1042.   width = w;
  1043.   max = 255;
  1044.   if (2 == t || 5 == t)
  1045.     data = new unsigned char [w*h];
  1046.   if (3 == t || 6 == t)
  1047.     data = new unsigned char [w*h*3];
  1048. }
  1049. /*
  1050.  * Creates an image from this signal
  1051.  */
  1052. Image::Image(double *signal, int size){
  1053.   double max_val, min_val, range;
  1054.   int hist[MAX_SIGNAL_SIZE];
  1055.   int i, x, y, inc; if (size > 1024) size = 1024;
  1056.   type = 5;
  1057.   height = HIST_DIM;
  1058.   width = size;
  1059.   max = 255;
  1060.   data = new unsigned char [height*width];
  1061.   //Find the maximum value for the signal
  1062.   max_val = min_val = signal[0];
  1063.   for (i=0; i<size; i++){
  1064.     if (max_val < signal[i])
  1065.       max_val = signal[i];
  1066.     if (min_val > signal[i])
  1067.       min_val = signal[i];
  1068.   }
  1069.   if (max_val == min_val){
  1070.     max_val = 1.0;
  1071.     min_val = -1.0;
  1072.   }
  1073.     range = max_val - min_val;
  1074.   for (i=0; i<size; i++ ){
  1075.     hist[i] = (int)(( (signal[i]-min_val) /range)*255.0);
  1076.   }
  1077.     //color this white
  1078.   for (x=0; x<size; x++ ){
  1079.     for (y=0; y<HIST_DIM; y++ ){
  1080.       data[x + (HIST_DIM-y-1)*size] = 255;
  1081.     }
  1082.     }
  1083.     //Draw zero line
  1084.     y = (int) -((min_val / range)*255.0);
  1085.     if (y>=0){
  1086.         for (x=0; x<size; x++){
  1087.             data[x + (HIST_DIM-y-1)*size] = 127;
  1088.         }
  1089.     }
  1090.   for (x=1; x<size; x++){
  1091.         if (hist[x-1] > hist[x]) inc = -1;
  1092.         else inc = 1;
  1093.     for (y=hist[x-1]; y!=hist[x]; y+=inc ){
  1094.       //color this black
  1095.       data[x + (HIST_DIM-y-1)*size] = 0;
  1096.     }
  1097.     data[x + (HIST_DIM-hist[x]-1)*size] = 0;
  1098.   }
  1099.     
  1100. }
  1101. Image::~Image(){
  1102.   if (NULL != data) delete [] data;
  1103. }
  1104. /*
  1105.  * Purpose: Go through and color each connected (r,g,b) region in the image
  1106.  *          a different color.
  1107.  *
  1108.  *          returns the number of shapes colored.
  1109.  *
  1110.  *          The Id of the shape is equvalent to the Red component value.
  1111.  */
  1112. unsigned char Image::colorShapes(unsigned char tol, int min_area){
  1113.   int x, y, ind;
  1114.   unsigned char r, g, b, c_r, c_g, c_b;
  1115.   unsigned char br, bg, bb;
  1116.   unsigned char i=1;
  1117.   int area;
  1118.   makeColor();
  1119.   if (type != 6 && type != 3){
  1120.     printf("Error: Image::colorShapes(). must have color image to color shapesn");
  1121.     return 0;
  1122.   } 
  1123.   getShapeColor(&r, &g, &b, &br, &bg, &bb, 0);
  1124. //  printf("shape color %d, %d, %d n", (int)r, (int)g, (int)b);//DBG
  1125. //  printf("bg color %d, %d, %d min_area: %dn", (int)br, (int)bg, (int)bb, min_area);//DBG
  1126.   getColorFromID(i, &c_r, &c_g, &c_b);
  1127.   for(y=0; y<height; y++){
  1128.     for(x=0; x<width; x++){
  1129.       ind = (x+y*width)*3;
  1130.       if ( abs( (int)r-(int)data[ind] )   <= (int)tol &&
  1131.            abs( (int)g-(int)data[ind+1] ) <= (int)tol &&
  1132.            abs( (int)b-(int)data[ind+2] ) <= (int)tol){
  1133.         area = 0;
  1134.         seedFill(x, y, r, g, b, c_r, c_g, c_b, tol, &area);
  1135. //        printf("area: %dn", area);
  1136.         if (area < min_area){
  1137. //          printf("Area too small for shapen");
  1138.           //Area too small for shape
  1139.           area = 0;
  1140.           seedFill(x, y, c_r, c_g, c_b, br, bg, bb, tol, &area);
  1141.         }
  1142.         else{
  1143.           i++;
  1144.           getColorFromID(i, &c_r, &c_g, &c_b);
  1145.         }
  1146.       }
  1147.     }
  1148.   }
  1149.   return i-1;
  1150.   
  1151. }
  1152. /*
  1153.  * Purpose: print out and create image of info
  1154.  */
  1155. void Image::createInfoImage(int num_shapes){
  1156.   int i;
  1157.   int cx, cy;
  1158.   int area;
  1159.   double orient;
  1160.   for (i=1; i<=num_shapes; i++){
  1161.     area = markCenterOfGravityColor(&cx, &cy, (unsigned char) i);
  1162.     orient = markOrientation(i, cx, cy);
  1163.     printf("Shape %d| area: %d, CoG: (%d, %d), Oriet: %fn", i, area, cx, cy, orient);
  1164.   }
  1165.   //getArea
  1166.   //centerofgravity
  1167.   //orientation
  1168. }
  1169. void Image::clear(unsigned char val){
  1170.   int ind;
  1171.   if (2 == type || 5 == type){
  1172.     for (int y=0; y<height; y++ ){
  1173.       for (int x=0; x<width; x++){
  1174.         ind = y*width +x;
  1175.         data[ind] =  val;
  1176.       }
  1177.     }
  1178.   }
  1179.   else if (3 == type || 6 == type){
  1180.     for (int y=0; y<height; y++ ){
  1181.       for (int x=0; x<width; x++){
  1182.         ind = (y*width +x)*3;
  1183.         data[ind  ] =  val;
  1184.         data[ind+1] =  val;
  1185.         data[ind+2] =  val;
  1186.       }
  1187.     }
  1188.   }
  1189. }
  1190. void Image::copy(Image *from_img){
  1191.   int dim, loc;
  1192.   if (NULL != data) delete [] data;
  1193.   height = from_img->height;
  1194.   width  = from_img->width;
  1195.   type   = from_img->type;
  1196.   max    = from_img->max;
  1197.   if (3 == type || 6 == type) dim = width*height*3;
  1198.   else dim = width*height*1;
  1199.   //allocate mem, and copy
  1200.   data = new unsigned char [dim];
  1201.   for (loc=0; loc<dim; loc++){
  1202.     data[loc] = from_img->data[loc];
  1203.   }
  1204. }
  1205. /*
  1206.  * Purpose: Get the difference of the two images and assign it to this one.
  1207.  */
  1208. void Image::diff(Image *img){
  1209.   int x, y;
  1210.   if (height != img->height || width != img->width){
  1211.     printf("Error in Image::diff(Image *img). Images must have the same dimensionsn");
  1212.     return;
  1213.   }
  1214.   img->makeGray();
  1215.   makeGray();
  1216.   for(y=0; y<height; y++){
  1217.     for(x=0; x<width; x++){
  1218.       data[x + y*width] = (unsigned char) abs( (int)data[x + y*width] - (int)img->data[x + y*width]);
  1219.     }
  1220.   }
  1221. }
  1222. void Image::dilate(int size){
  1223.   int i, j, l, m;
  1224.   int offset;
  1225.   unsigned char *new_data;
  1226.   makeGray();
  1227.   offset = (size-1)/2;
  1228.   new_data = new unsigned char [height*width];
  1229.   for(i=0; i<height*width; i++){
  1230.     new_data[i] = 0;
  1231.   }
  1232.   for(i=offset; i<width-offset; i++){
  1233.     for(j=offset; j<height-offset; j++){
  1234.       //If the pixel is set.
  1235.       if(data[i + j*width] > 0){
  1236.         for(l=0; l<size; l++){
  1237.           for(m=0; m<size; m++){
  1238. //            if(filter[l + m*size] > tol)
  1239. //              new_data[i-offset+l+(j-offset+m)*width] = filter[l + m*size];
  1240.               new_data[i-offset+l+(j-offset+m)*width] = 255;//filter[l + m*size];
  1241.           }
  1242.         }
  1243.       }
  1244.     }
  1245.   }
  1246.   delete [] data;
  1247.   data = new_data;
  1248. }
  1249. void Image::erode(int size){
  1250.   int i, j, l, m;
  1251.   int offset;
  1252.   int set;
  1253.   unsigned char *new_data;
  1254.   makeGray();
  1255.   offset = (size-1)/2;
  1256.   new_data = new unsigned char [height*width];
  1257.   for(i=0; i<height*width; i++){
  1258.     new_data[i] = 0;
  1259.   }
  1260.   for(i=offset; i<width-offset; i++){
  1261.     for(j=offset; j<height-offset; j++){
  1262.       set = 1;
  1263.       for(l=0; l<size; l++){
  1264.         for(m=0; m<size; m++){
  1265.           if( data[i-offset+l+(j-offset+m)*width] == 0){
  1266.               set = 0;
  1267.               m = l = size;//break out
  1268.           }
  1269.         }
  1270.       }
  1271.       //Now set the pixel in the outgoing image if needs be.
  1272.       if (1 == set)
  1273.         new_data[i + j*width] = data[i + j*width];
  1274.       else
  1275.         new_data[i + j*width] = 0;
  1276.     }
  1277.   }
  1278.   delete data;
  1279.   data = new_data;
  1280. }
  1281. /*
  1282.  * Purpose: Makes the image a residual of a predicted encoding scheme.
  1283.  * -predictor
  1284.  *   a b c
  1285.  *   d e
  1286.  *       e = (a+b+c+d)/4;
  1287.  *  if flag == 1 => encode;
  1288.  *  if flag == 0 => decode;
  1289.  */
  1290. void Image::predCode1(int flag){
  1291.   int i, j;
  1292.   unsigned char *new_data;
  1293.   int pred, num_pred;
  1294.   if (flag == 1){
  1295.     new_data = new unsigned char[width*height];
  1296.   }
  1297.   else{
  1298.     new_data = data;
  1299.   }
  1300.   makeGray();
  1301.   for (j=0; j<height; j++){
  1302.     for (i=0; i<width; i++){
  1303.       pred = 0;
  1304.       num_pred = 0;
  1305.       if (i>0){
  1306.         pred += (int)data[i-1 + j*width];//d
  1307.         num_pred ++;
  1308.       }
  1309.       if (j>0){
  1310.         pred += (int)data[i + (j-1)*width];//b
  1311.         num_pred ++;
  1312.         if (i>0){
  1313.           pred += (int)data[i-1 + (j-1)*width];//a
  1314.           num_pred ++;
  1315.         }
  1316.         if (i<width-2){
  1317.           pred += (int)data[i+1 + (j-1)*width];//c
  1318.           num_pred ++;
  1319.         }
  1320.       }
  1321.       if (num_pred == 0){
  1322.         new_data[i + j*width] =data[i + j*width];
  1323.         continue;
  1324.       }
  1325.       pred = pred/num_pred;
  1326.       if (flag == 1){//encode
  1327.         if (data[i + j*width] < (unsigned char)pred)
  1328.           new_data[i + j*width] =(unsigned char)( 256 + (int)data[i + j*width] - pred);
  1329.         else
  1330.           new_data[i + j*width] = data[i + j*width] - (unsigned char)pred;
  1331.       }
  1332.       else{ //decode
  1333.         if( 256-(unsigned char)pred < data[i + j*width])
  1334.           new_data[i + j*width] = (unsigned char)((int)data[i + j*width] + pred - 256);
  1335.         else
  1336.           new_data[i + j*width] = data[i + j*width] + (unsigned char)pred;
  1337.       }
  1338.     }
  1339.   }
  1340.   if (flag == 1) delete data;
  1341.   data = new_data;
  1342. }
  1343. /*
  1344.  * Purpose: Makes the image a residual of a predicted encoding scheme.
  1345.  * -predictor
  1346.  *   a b c
  1347.  *   d e
  1348.  *       e = (a+b+c+d)/4;
  1349.  *  if flag == 1 => encode;
  1350.  *  if flag == 0 => decode;
  1351.  */
  1352. void Image::predCode2(int flag){
  1353.   int i, j;
  1354.   unsigned char *new_data;
  1355.   int pred, num_pred;
  1356.   if (flag == 1){
  1357.     new_data = new unsigned char[width*height];
  1358.   }
  1359.   else{
  1360.     new_data = data;
  1361.   }
  1362.   makeGray();
  1363.   for (j=0; j<height; j++){
  1364.     for (i=0; i<width; i++){
  1365.       pred = 0;
  1366.       num_pred = 0;
  1367.       if (i>1){
  1368.         //two step predictor
  1369.         pred += 2*(int)data[i-1 + j*width] - (int)data[i-2 + j*width];//d
  1370.         num_pred ++;
  1371.       }
  1372.       else if (i>0){
  1373.         pred += (int)data[i-1 + j*width];//d
  1374.         num_pred ++;
  1375.       }
  1376.       if (j>1){
  1377.         pred += 2*(int)data[i + (j-1)*width] - (int)data[i + (j-2)*width];//b
  1378.         num_pred ++;
  1379.         if (i>0){
  1380.           pred += 2*(int)data[i-1 + (j-1)*width] - (int)data[i-1 + (j-2)*width] ;//a
  1381.           num_pred ++;
  1382.         }
  1383.         if (i<width-2){
  1384.           pred += 2*(int)data[i+1 + (j-1)*width]- (int)data[i+1 + (j-2)*width];//c
  1385.           num_pred ++;
  1386.         }
  1387.       }
  1388.       else if (j>0){
  1389.         pred += (int)data[i + (j-1)*width];//b
  1390.         num_pred ++;
  1391.         if (i>0){
  1392.           pred += (int)data[i-1 + (j-1)*width];//a
  1393.           num_pred ++;
  1394.         }
  1395.         if (i<width-2){
  1396.           pred += (int)data[i+1 + (j-1)*width];//c
  1397.           num_pred ++;
  1398.         }
  1399.       }
  1400.       if (num_pred == 0){
  1401.         new_data[i + j*width] =data[i + j*width];
  1402.         continue;
  1403.       }
  1404.       if (pred < 0) pred = 0;
  1405.       pred = pred/num_pred;
  1406.       if (flag == 1){//encode
  1407.         if (data[i + j*width] < (unsigned char)pred)
  1408.           new_data[i + j*width] =(unsigned char)( 256 + (int)data[i + j*width] - pred);
  1409.         else
  1410.           new_data[i + j*width] = data[i + j*width] - (unsigned char)pred;
  1411.       }
  1412.       else{ //decode
  1413.         if( 256-(unsigned char)pred < data[i + j*width])
  1414.           new_data[i + j*width] = (unsigned char)((int)data[i + j*width] + pred - 256);
  1415.         else
  1416.           new_data[i + j*width] = data[i + j*width] + (unsigned char)pred;
  1417.       }
  1418.     }
  1419.   }
  1420.   if (flag == 1) delete data;
  1421.   data = new_data;
  1422. }
  1423. /*
  1424.  * Purpose: Equalizes a gray image
  1425.  */
  1426. void Image::equalize(){
  1427.   int i, j, x, y;
  1428.   double max;
  1429.   double lup[HIST_DIM];
  1430.   double *pdf = new double[HIST_DIM];
  1431.   //First check the image type
  1432.   if (type == 6 || type == 3){
  1433.     printf("Error: Image::equalize(). Image must be gray scale.n");
  1434.     return;
  1435.   } 
  1436.   //Second create the PDF (probability density function) for this image.
  1437.   getPDF(pdf, HIST_DIM, &max);
  1438.   //Third create a look up table for the color intensities transformation.
  1439.   for (i=0; i<HIST_DIM; i++){
  1440.     lup[i] = 0.0;
  1441.     for (j=0; j<=i; j++){
  1442.       lup[i] += pdf[j];
  1443.     }
  1444.     lup[i] *= (HIST_DIM-1);
  1445.     //printf("lup[%d]: %fn",i ,  lup[i]);
  1446.   }
  1447.   //Fourth now use the lookup table to get the new image.
  1448.   for (y=0; y<height; y++ ){
  1449.     for (x=0; x<width; x++){
  1450.       data[x + y*width] = (unsigned char) lup[data[x + y*width]];
  1451.     }
  1452.   }
  1453.   if (NULL != pdf) delete [] pdf;
  1454. }
  1455. /*
  1456.  * Purpose: flips the image vertically
  1457.  */
  1458. void Image::flip(){
  1459.   unsigned char r;
  1460.   int i, j, h_min_j, ind1, ind2;
  1461.   
  1462.   makeGray();
  1463.   for (j=0; j<height/2; j++){
  1464.     h_min_j = height - 1 - j;
  1465.     for (i=0; i<width; i++){
  1466.       ind1 = i + j*width;
  1467.       ind2 = i + h_min_j*width;
  1468.       r = data[ind1];
  1469.       data[ind1] = data[ind2];
  1470.       data[ind2] = r;
  1471.     }
  1472.   }
  1473.   /*
  1474.   unsigned char r, g, b;
  1475.   int i, j, h_min_j, ind1, ind2;
  1476.   makeColor();
  1477.   for (j=0; j<height/2; j++){
  1478.     h_min_j = height - 1 - j;
  1479.     for (i=0; i<width; i++){
  1480.       ind1 = (i + j*width)*3;
  1481.       ind2 = (i + h_min_j*width)*3;
  1482.       r = data[ind1];
  1483.       g = data[ind1+1];
  1484.       b = data[ind1+2];
  1485.       data[ind1] = data[ind2];
  1486.       data[ind1+1] = data[ind2+1];
  1487.       data[ind1+2] = data[ind2+2];
  1488.       data[ind2] = r;
  1489.       data[ind2+1] = g;
  1490.       data[ind2+2] = b;
  1491.     }
  1492.   }
  1493.   */
  1494. }
  1495. double Image::markOrientation(int i, int cx, int cy){
  1496.   double angle;
  1497.   double u11, u02, u20;
  1498.   int x, y;
  1499.   int total_x, total_y, count;
  1500.   makeColor();
  1501.   double vx, vy;
  1502.   total_x = total_y = count = 0;
  1503.   u11 = u20 = u02 = 0.0;
  1504.   for(y=0; y<height; y++){
  1505.     for(x=0; x<width; x++){
  1506.       if ( (data[(x + y*width)*3 + 2 ] - i)  == 0){
  1507.         u11 += (double)(x - cx)*(y - cy);
  1508.         u20 += (double)(x - cx)*(x - cx);
  1509.         u02 += (double)(y - cy)*(y - cy);
  1510.       }
  1511.     }
  1512.   }
  1513.   //angle = 0.5*atan(2.0*u11/(u20-u02));
  1514.   angle = 0.5*atan2(2.0*u11,(u20-u02));
  1515.   vx = 30*cos(angle);
  1516.   vy = 30*sin(angle);
  1517.   drawLine(cx-(int)vx, cy-(int)vy, cx+(int)vx, cy+(int)vy, 0, 128, 255);
  1518.   return angle;
  1519. }
  1520. /*
  1521.  * Purpose: Finds the center of gravity for all pixels whose value is greater
  1522.  *          than tol.
  1523.  */
  1524. int Image::markCenterOfGravityColor(int *cx, int *cy, unsigned char i){
  1525.   int x, y;
  1526.   int total_x, total_y, count;
  1527.   makeColor();
  1528.   total_x = total_y = count = 0;
  1529.   for(y=0; y<height; y++){
  1530.     for(x=0; x<width; x++){
  1531.       if ( (data[(x + y*width)*3 + 2 ] - i)  == 0){
  1532.         total_x += x;
  1533.         total_y += y;
  1534.         count ++;
  1535.       }
  1536.     }
  1537.   }
  1538.   *cx = (int)((double)total_x/(double)count);
  1539.   *cy = (int)((double)total_y/(double)count);
  1540.   drawLine(*cx-6, *cy, *cx+6, *cy, 255, 255, 255);
  1541.   drawLine(*cx, *cy-6, *cx, *cy+6, 255, 255, 255);
  1542.   return count;
  1543. }
  1544. /*
  1545.  * Purpose: draws a pixel line on the image from x0, y0, to x1, y1
  1546.  *          using the midpoint line algorithm.
  1547.  */
  1548. void Image::drawLine(int x0, int y0, int x1, int y1, unsigned char r, unsigned char g, unsigned char b) {
  1549.   int dy = y1 - y0;
  1550.   int dx = x1 - x0;
  1551.   int stepx, stepy;
  1552.   int fraction;
  1553.   if (dy < 0) { dy = -dy;  stepy = -1; } else { stepy = 1; }
  1554.   if (dx < 0) { dx = -dx;  stepx = -1; } else { stepx = 1; }
  1555.   dy <<= 1;                     // dy is now 2*dy
  1556.   dx <<= 1;                     // dx is now 2*dx
  1557.   if (x0 < width && x0 > 0 && y0 < height && y0 > 0){
  1558.     data[(x0+ y0*width)*3    ] = r;
  1559.     data[(x0+ y0*width)*3 + 1] = g;
  1560.     data[(x0+ y0*width)*3 + 2] = b;
  1561.   }
  1562.   if (dx > dy) {
  1563.       fraction = dy - (dx >> 1);          // same as 2*dy - dx
  1564.       while (x0 != x1) {
  1565.           if (fraction >= 0) {
  1566.               y0 += stepy;
  1567.               fraction -= dx;                // same as fraction -= 2*dx
  1568.           }
  1569.           x0 += stepx;
  1570.           fraction += dy;                    // same as fraction -= 2*dy
  1571.           if (x0 < width && x0 > 0 && y0 < height && y0 > 0){
  1572.             data[(x0+ y0*width)*3    ] = r;
  1573.             data[(x0+ y0*width)*3 + 1] = g;
  1574.             data[(x0+ y0*width)*3 + 2] = b;
  1575.           }
  1576.       }
  1577.   } else {
  1578.       fraction = dx - (dy >> 1);
  1579.       while (y0 != y1) {
  1580.           if (fraction >= 0) {
  1581.               x0 += stepx;
  1582.               fraction -= dy;
  1583.           }
  1584.           y0 += stepy;
  1585.           fraction += dx;
  1586.           if (x0 < width && x0 > 0 && y0 < height && y0 > 0){
  1587.             data[(x0+ y0*width)*3    ] = r;
  1588.             data[(x0+ y0*width)*3 + 1] = g;
  1589.             data[(x0+ y0*width)*3 + 2] = b;
  1590.           }
  1591.       }
  1592.   }
  1593. }
  1594. /*
  1595.  * Purpose: Finds the center of gravity for all pixels whose value is greater
  1596.  *          than tol.
  1597.  */
  1598. void Image::getCenterOfGravity(int *cx, int *cy, unsigned char tol){
  1599.   int x, y;
  1600.   int total_x, total_y, count;
  1601.   makeGray();
  1602.   total_x = total_y = count = 0;
  1603.   for(y=0; y<height; y++){
  1604.     for(x=0; x<width; x++){
  1605.       if (data[x + y*width] > tol){
  1606.         total_x += x;
  1607.         total_y += y;
  1608.         count ++;
  1609.       }
  1610.     }
  1611.   }
  1612.   *cx = (int)((double)total_x/(double)count);
  1613.   *cy = (int)((double)total_y/(double)count);
  1614. }
  1615. /*
  1616.  * Purpose: returns a pdf (probability density function) for the historgram of
  1617.  *          this image scaled to,  area = 1.0;
  1618.  */
  1619. void Image::getPDF(double *pdf, int n, double *max){
  1620.   double val, total;
  1621.   int x, y, i;
  1622.   if (2 != type && 5 != type){
  1623.     printf("Error: Image::getPDF(). must have gray image to get gray PDFn");
  1624.     return;
  1625.   }
  1626.   if (n != HIST_DIM){
  1627.     printf("Error: Image::getPDF(). Cannot handle a pdf of size other than %dn", HIST_DIM);
  1628.     return;
  1629.   }
  1630.   for (i=0; i<n; i++){
  1631.     pdf[i] = 0.0;
  1632.   }
  1633.   *max = 0.0;
  1634.   total = height*width;
  1635.   for (y=0; y<height; y++ ){
  1636.     for (x=0; x<width; x++){
  1637.       val = ++pdf[data[x + y*width]];
  1638.       if (val > *max) *max = val;
  1639.     }
  1640.   }
  1641.   //Normalize it afterwards for better floating point calculation
  1642.   //rather than adding little increments each step.
  1643.   for (i=0; i<n; i++){
  1644.     pdf[i] /= total;
  1645.   }
  1646. }
  1647. void Image::getItemBounds(unsigned char id, int *bot_x, int *bot_y,
  1648.                           int *top_x, int *top_y){
  1649.   int x, y, ind;
  1650.   makeColor();
  1651.   *bot_x = width;
  1652.   *bot_y = height;
  1653.   *top_x = 0;
  1654.   *top_y = 0;
  1655.   for(y=0; y<height; y++){
  1656.     for(x=0; x<width; x++){
  1657.       ind = (x+y*width)*3;
  1658.       if ( data[ind+2] == id){
  1659.         if (x > *top_x) *top_x = x;
  1660.         if (x < *bot_x) *bot_x = x;
  1661.         if (y > *top_y) *top_y = y;
  1662.         if (y < *bot_y) *bot_y = y;
  1663.       }
  1664.     }
  1665.   }
  1666. }
  1667. unsigned char Image::getBoundAverage(int bot_x, int bot_y,
  1668.                                     int top_x, int top_y){
  1669.   int x, y, ind, total, count;
  1670.   count = 0;
  1671.   total = 0;
  1672.   for(y=bot_y; y<top_y; y++){
  1673.     for(x=bot_x; x<top_x; x++){
  1674.       ind = x+y*width;
  1675.       total += (int)data[ind];
  1676.       count ++;
  1677.     }
  1678.   }
  1679.   return (int)((double)total/(double)count);
  1680. }
  1681. void Image::unionThreshRegion(Image *img, int bot_x, int bot_y, int top_x,
  1682.                               int top_y, int range, unsigned char t){
  1683.   int x, y, ind, l=0, m=0;
  1684.   double count, ave;
  1685.   makeGray();
  1686.   img->makeGray();
  1687.   for(y=bot_y+range; y<top_y-range; y++){
  1688.     for(x=bot_x+range; x<top_x-range; x++){
  1689.       //do averaging filter
  1690.       count = 0.0;
  1691.       ave = 0.0;//a[i + j*w];
  1692.       for(l=x-range; l<x+range; l++){
  1693.         for(m=y-range; m<y+range; m++){
  1694.           ave += (double)img->data[l + m*width];
  1695.           count ++;
  1696.         }
  1697.       }
  1698.       ave /= count;
  1699.       ind = l+m*width;
  1700.       if (img->data[ind] > t){
  1701.         data[ind] = 255;
  1702.       }
  1703.     }
  1704.   }
  1705. }
  1706. double *Image::getDouble(double offset){
  1707.   int i;
  1708.   double *d_data = new double[height*width];
  1709.   makeGray();
  1710.   for (i=0; i<height*width; i++){
  1711.     d_data[i] = (double)data[i] + offset;
  1712.   }
  1713.   return d_data;
  1714. }
  1715. void Image::invert(){
  1716.   int i;
  1717.   makeGray();
  1718.   for (i=0; i<height*width; i++){
  1719.     data[i] = 255 - data[i];
  1720.   }
  1721. }
  1722. /*
  1723.  * Purpose: We'll get the color for the shapes and set our own bg color
  1724.  */
  1725. void Image::getShapeColor(unsigned char *r, unsigned char *g, unsigned char *b,
  1726.                        unsigned char *br, unsigned char *bg, unsigned char *bb,
  1727.                        unsigned char tol){
  1728.   int x, y, ind;
  1729.   if (2 == type || 5 == type){
  1730.     printf("Error: Image::getShapeColor(). You must first make this a color image.n");
  1731.     return;
  1732.   }
  1733.   /*
  1734.   *br = data[0];
  1735.   *bg = data[1];
  1736.   *bb = data[2];
  1737.   */
  1738.   *br = 0;
  1739.   *bg = 0;
  1740.   *bb = 0;
  1741.   
  1742.   for (y=0; y<height; y++ ){
  1743.     for (x=0; x<width; x++){
  1744.       ind = (x + y*width)*3;
  1745.       if ( abs( (int)(*br)-(int)data[ind] )   > (int)tol ||
  1746.            abs( (int)(*bg)-(int)data[ind+1] ) > (int)tol ||
  1747.            abs( (int)(*bb)-(int)data[ind+2] ) > (int)tol){
  1748.         *r = data[ind];
  1749.         *g = data[ind+1];
  1750.         *b = data[ind+2];
  1751.         x = width;
  1752.         y = height;
  1753.       }
  1754.     }
  1755.   }
  1756.   //Now we'll recolor the back ground.
  1757.   for (y=0; y<height; y++ ){
  1758.     for (x=0; x<width; x++){
  1759.       ind = (x + y*width)*3;
  1760.       if ( abs( (int)(*br)-(int)data[ind] )   <= (int)tol ||
  1761.            abs( (int)(*bg)-(int)data[ind+1] ) <= (int)tol ||
  1762.            abs( (int)(*bb)-(int)data[ind+2] ) <= (int)tol){
  1763.         data[ind]   = *br;
  1764.         data[ind+1] = *bg;
  1765.         data[ind+2] = *bb;
  1766.       }
  1767.     }
  1768.   }
  1769. }
  1770. /*
  1771.  * Purpose: pads the image with black around the borders.
  1772.  */
  1773. void Image::pad(int new_width, int new_height){
  1774.   int i, x_off, y_off, nx, ny;
  1775.   unsigned char *new_data;
  1776.   if (NULL == data) return;
  1777.   if (width > new_width || height > new_height){
  1778.     printf("Error: size for pad not big enoughn");
  1779.     return;
  1780.   }
  1781.   makeGray();
  1782.   new_data = new unsigned char[new_height*new_width];
  1783.   for (i=0; i<new_height*new_width; i++){
  1784.     new_data[i] = 0;
  1785.   }
  1786.   x_off = (new_width-width)/2;
  1787.   y_off = (new_height-height)/2;
  1788.   for (ny = 0; ny<height; ny++){
  1789.     for (nx = 0; nx<width; nx++){
  1790.       new_data[nx+x_off + (ny+y_off)*new_width] = data[nx + ny*width];
  1791.     }
  1792.   }
  1793.   delete [] data;
  1794.   data = new_data;
  1795.   height = new_height;
  1796.   width = new_width;
  1797.   
  1798. }
  1799. /*
  1800.  * Purpose: rotates the image the specified amount.
  1801.  * Note:    The value passed in is assumed to be in degrees
  1802.  */
  1803. void Image::rotate(double theta){
  1804.   unsigned char *new_data;
  1805.   int nx, ny;
  1806.   double x1, x2, y1, y2; 
  1807.   double sint, cost;
  1808.   double xc, yc, tx, ty;
  1809.   double sample1, sample2;
  1810.   int ix2, iy2;
  1811.   makeGray();
  1812.   //convert to radians
  1813.   theta *= 6.28/360.0;
  1814.   sint = sin(theta);
  1815.   cost = cos(theta);
  1816.   new_data = new unsigned char [height*width];
  1817.   xc = (double)width/2.0;
  1818.   yc = (double)height/2.0;
  1819.   //To do the backwards mapping, we'll just take samples in the new image
  1820.   //and rotate back to the previous to get the value
  1821.   for (ny = 0; ny<height; ny++){
  1822.     for (nx = 0; nx<width; nx++){
  1823.       x1 = (double)nx;
  1824.       y1 = (double)ny;
  1825.       x2 = (x1-xc)*cost - (y1-yc)*sint + xc;
  1826.       y2 = (x1-xc)*sint + (y1-yc)*cost + yc;
  1827.       ix2 = (int)x2;
  1828.       iy2 = (int)y2;
  1829.       tx = x2 - (double)ix2;
  1830.       ty = y2 - (double)iy2;
  1831.       
  1832.       if (ix2+1<width && ix2>=0 && iy2+1<height && iy2>=0){
  1833.         sample1 = (1.0-tx)*(double)data[ix2 + iy2*width] +
  1834.                         tx*(double)data[ix2+1 + iy2*width];
  1835.         sample2 = (1.0-tx)*(double)data[ix2 + (iy2+1)*width] +
  1836.                         tx*(double)data[ix2+1 + (iy2+1)*width];
  1837.         new_data[nx + ny*width] = (unsigned char)((1-ty)*sample1 + ty*sample2);
  1838.       }
  1839.       else{
  1840.         new_data[nx + ny*width] = 128;
  1841.       }
  1842.     }
  1843.   }
  1844.   //Assign the new image
  1845.   delete [] data;
  1846.   data = new_data;
  1847. }
  1848. void Image::makeSinusoid(int size, double freq, double mag, int axis){
  1849.   int i, j;
  1850.   double N;
  1851.   double half = mag/2.0;
  1852.   unsigned char val;
  1853.   if (data != NULL) delete [] data;
  1854.   data = new unsigned char[size*size];
  1855.   height = width = size;
  1856.   type = 5;
  1857.   max = (int)mag;
  1858.   N = (double) size;
  1859.   if (axis == X_AXIS){
  1860.     for (i=0; i<size; i++){
  1861.       val = (unsigned char)(half*sin( (freq * (double)i)/N) + half);
  1862.       for (j=0; j<size; j++){
  1863.         data[i + j*size] = val;
  1864.       }
  1865.     }
  1866.   }
  1867.   else if (axis == Y_AXIS){
  1868.     for (j=0; j<size; j++){
  1869.       val = (unsigned char)(half*sin( (freq * (double)j)/N) + half);
  1870.       for (i=0; i<size; i++){
  1871.         data[i + j*size] = val;
  1872.       }
  1873.     }
  1874.   }
  1875. }
  1876. /*
  1877.  * Assigns the image the magnitued plot of the 2d signal
  1878.  */
  1879. void Image::makeMag(Complex *D_Data, int size){
  1880.   int i, j;
  1881.   double mag;
  1882.   double d_max, d_min;
  1883.   if (data != NULL) delete [] data;
  1884.   data = new unsigned char[size*size];
  1885.   height = width = size;
  1886.   type = 5;
  1887.   max = 255;
  1888.   //Set max and min
  1889.   d_max = d_min = D_Data[0].getMag();
  1890.   for (i=1; i<size*size; i++){
  1891.     mag = D_Data[i].getMag();
  1892.     if (d_max < mag) d_max = mag;
  1893.     else if (d_min > mag) d_min = mag;
  1894.   }
  1895.   for (i=0; i<size; i++){
  1896.     for (j=0; j<size; j++){
  1897.       data[i + j*size] = (unsigned char)( ((D_Data[i+j*size].getMag() - d_min)/(d_max-d_min))*(double)max);
  1898.     }
  1899.   }
  1900.   swapQuads();
  1901. }
  1902. /*
  1903.  * Assigns the image the magnitued plot of the 2d signal
  1904.  */
  1905. void Image::makePhase(Complex *D_Data, int size){
  1906.   int i, j;
  1907.   double mag;
  1908.   double d_max, d_min;
  1909.   if (data != NULL) delete [] data;
  1910.   data = new unsigned char[size*size];
  1911.   height = width = size;
  1912.   type = 5;
  1913.   d_max = 255;
  1914.   //Set max and min
  1915.   d_max = d_min = D_Data[0].getPhase();
  1916.   for (i=1; i<size*size; i++){
  1917.     mag = D_Data[i].getPhase();
  1918.     if (d_max < mag) d_max = mag;
  1919.     else if (d_min > mag) d_min = mag;
  1920.   }
  1921.   for (i=0; i<size; i++){
  1922.     for (j=0; j<size; j++){
  1923.       data[i + j*size] = (unsigned char)( ((D_Data[i+j*size].getPhase() - d_min)/(d_max-d_min))*(double)max);
  1924.     }
  1925.   }
  1926.   swapQuads();
  1927. }
  1928. void Image::mask(Image *m){
  1929.   int i, j;
  1930.   makeGray();
  1931.   m->makeGray();
  1932.   for (i=0; i<m->width; i++){
  1933.     for (j=0; j<m->height; j++){
  1934.       if (m->data[i + j*m->width] == 0){
  1935.         if (i<width && j<height){
  1936.           data[i + j*width] = 0;
  1937.         }
  1938.       }
  1939.     }
  1940.   }
  1941. }
  1942. void Image::signalToImage2D(double *d_data, int w, int h){
  1943.   int i, j;
  1944.   double mag;
  1945.   double d_max, d_min;
  1946.   if (data != NULL) delete [] data;
  1947.   data = new unsigned char[w*h];
  1948.   height = h;
  1949.   width = w;
  1950.   type = 5;
  1951.   max = 255;
  1952.   //Set max and min
  1953.   d_max = d_min = d_data[0];
  1954.   for (i=1; i<w*h; i++){
  1955.     mag = d_data[i];
  1956.     if (d_max < mag) d_max = mag;
  1957.     else if (d_min > mag) d_min = mag;
  1958.   }
  1959.   for (i=0; i<w; i++){
  1960.     for (j=0; j<h; j++){
  1961.       data[i + j*w] = (unsigned char)( ((d_data[i+j*w] - d_min)/(d_max-d_min))*(double)max);
  1962.     }
  1963.   }
  1964. }
  1965. void Image::swapQuads(){
  1966.   int i, j;
  1967.   int halfw, halfh;
  1968.   unsigned char *new_data;
  1969.   //now copy the image quadrants to where they're supposed to be.
  1970.   new_data = new unsigned char[width*height];
  1971.   halfw = width/2;
  1972.   halfh = height/2;
  1973.   for (i=0; i<width; i++){
  1974.     for (j=0; j<height; j++){
  1975.       if (i<halfw){
  1976.         if (j<halfh)
  1977.           new_data[i + j*width] = data[i + halfw + (j+halfh)*width];
  1978.         else
  1979.           new_data[i + j*width] = data[i + halfw + (j-halfh)*width];
  1980.       }
  1981.       else{
  1982.         if (j<halfh)
  1983.           new_data[i + j*width] = data[i - halfw + (j+halfh)*width];
  1984.         else
  1985.           new_data[i + j*width] = data[i - halfw + (j-halfh)*width];
  1986.       }
  1987.     }
  1988.   }
  1989.   delete [] data;
  1990.   data = new_data;
  1991. }
  1992. void Image::threshold(unsigned char thresh){
  1993.   int i;
  1994.   int dim = width*height;
  1995.   makeGray();
  1996.   for (i=0; i<dim; i++){
  1997.     if (data[i] < thresh) data[i] = 0;
  1998.     else data[i] = 255;
  1999.   }
  2000. }
  2001. /************************************************************************
  2002. writeImage - writes image out to the given filename in binary PPM format
  2003. note: forced_n being non negative will cause mica to only save out that
  2004.       particular image in the series.
  2005. ************************************************************************/
  2006. int Image::write(char *filename){
  2007.   int    i, j;
  2008.   FILE   *fpt;
  2009.   unsigned char *save_line;
  2010.   int count;
  2011.   
  2012.   printf("saving %s... ", filename);
  2013.   fpt = fopen( filename, "wb" );
  2014.   if ( !fpt ){
  2015.     printf("Error: Could not create file: %sn", filename);
  2016.     return FALSE;
  2017.   }
  2018.   if (2 == type){//gray scale
  2019.     fprintf( fpt, "P2n" );
  2020.     fprintf( fpt, "# Created by the illustrious Mike D. Smithn" );
  2021.     fprintf( fpt, "%d %dn", width, height);
  2022.     fprintf( fpt, "%dn", max);
  2023.     for( j=0, count=0; j<height; j++ ){
  2024.       for( i=0; i<width; i++, count++){
  2025.         if (count>10){//no line in the file should be longer than 70 chars
  2026.           fprintf( fpt, "n");
  2027.           count = 0;
  2028.         }
  2029.         fprintf( fpt, "%d ", data[i + j*width]);
  2030.       }
  2031.     }
  2032.     fclose( fpt );
  2033.     printf("donen");
  2034.     return TRUE;
  2035.   }
  2036.   else if (3 == type){
  2037.     fprintf( fpt, "P3n" );
  2038.     fprintf( fpt, "# Created by the illustrious Mike D. Smithn" );
  2039.     fprintf( fpt, "%d %dn", width, height);
  2040.     fprintf( fpt, "%dn", max);
  2041.     for( j=0, count=0; j<height; j++ ){
  2042.       for( i=0; i<width; i++, count++){
  2043.         if (count>4){//no line in the file should be longer than 70 chars
  2044.           fprintf( fpt, "n");
  2045.           count = 0;
  2046.         }
  2047.         fprintf( fpt, "%d %d %d ", data[(i + j*width)*3], data[(i + j*width)*3 + 1], data[(i + j*width)*3 + 2]);
  2048.       }
  2049.     }
  2050.     fclose( fpt );
  2051.     printf("donen");
  2052.     return TRUE;
  2053.   }
  2054.   else if (5 == type){
  2055.     fprintf( fpt, "P5n" );
  2056.     //if ( HeaderInfo!=NULL ) fprintf( F, "# %sn", HeaderInfo );
  2057.     //else fprintf( F, "# Grayscale Imagen" );
  2058.     fprintf( fpt, "%d %dn%dn", width, height, max );
  2059.     for ( i=0; i<height; i++ ) {
  2060.       save_line = &data[i*width];
  2061.       if ( fwrite( (void*)save_line, 1, (size_t)width, fpt ) !=
  2062.            (unsigned int)width ){
  2063.         printf("fail to writen");
  2064.         return FALSE;
  2065.       }
  2066.     }
  2067.     printf("donen");
  2068.     fclose( fpt );
  2069.     return TRUE;
  2070.   }
  2071.   else if (6 == type){
  2072.     fprintf( fpt, "P6n" );
  2073.     //if ( HeaderInfo!=NULL ) fprintf( F, "# %sn", HeaderInfo );
  2074.     //else fprintf( F, "# Grayscale Imagen" );
  2075.     fprintf( fpt, "%d %dn%dn", width, height, max );
  2076.     for ( i=0; i<height; i++ ) {
  2077.       save_line = &data[i*width*3];
  2078.       if ( fwrite( (void*)save_line, 1, (size_t)(width*3), fpt ) !=
  2079.            (unsigned int)(width*3)  ){
  2080.         puts("fail to write");
  2081.         return FALSE;
  2082.       }
  2083.     }
  2084.     printf("donen");
  2085.     fclose( fpt );
  2086.     return TRUE;
  2087.   }
  2088.   else{
  2089.     printf("Failure to write PPM type %d not supportedn", type);
  2090.     return FALSE;
  2091.   }
  2092. }
  2093. /*
  2094.  * Purpose: read - reads a PNM file into an unsigned byte array
  2095.  * 
  2096.  * filename = the image file you wish to open.
  2097.  * width  = the image dimension in the x direction.
  2098.  * height  = the image dimension in the y direction.
  2099.  */
  2100. int Image::read(char *filename) {
  2101.   FILE *fp;
  2102.   register int w, h;
  2103.   int ind;
  2104.   char buffer[250];
  2105.   fp = fopen(filename, "rb");
  2106.   if (fp == NULL) {
  2107.     printf("Could not open file %sn", filename);
  2108.     return FALSE;
  2109.   }
  2110.   printf("Reading %sn", filename);
  2111.     //get the ppm format type
  2112.   fscanf(fp,"P%dn", &type);
  2113.   //skip comments and empty spaces.
  2114.   buffer[0] = '#';
  2115.   while(buffer[0] == '#' || buffer[0] == ' ' || buffer[0] == 'n')
  2116.     fgets(buffer,250,fp);
  2117.     //then get the width and height
  2118.   sscanf(buffer,"%d%d",&w,&h);
  2119. //  printf("Image width: %d, height: %dn", w, h);
  2120.   width = w;
  2121.   height = h;
  2122.     //get the maximum value for the ppm file (should be 255)
  2123.   fgets(buffer,250,fp);
  2124.   sscanf(buffer,"%d",&max);
  2125. //  printf("max pixel value is: %d.n", max);
  2126.   if (type == 2){//Gray scale image
  2127.     data = new unsigned char [w * h];
  2128.     // code to handle P2 type ppms (ascII type)
  2129.     for (int y=0; y<h; y++ ){
  2130.       for (int x=0; x<w; x++){
  2131.         int c;
  2132.         fscanf(fp,"%d",&c);
  2133.         data[y*w +x] = (unsigned char)( ((double)c/(double)max)*255.0);
  2134.       }
  2135.     }
  2136.   }
  2137.   else if (type == 3){
  2138.     data = new unsigned char [w * h * 3];
  2139.     // code to handle P3 type ppms (ascII type)
  2140.     for (int y=0; y<h; y++ ){
  2141.       for (int x=0; x<w; x++){
  2142.         int r,g,b;
  2143.         fscanf(fp,"%d%d%d",&r,&g,&b);
  2144.         ind = (y*w +x) * 3;
  2145.         data[ind]     = (unsigned char)( ((double)r/(double)max)*255.0);
  2146.         data[ind + 1] = (unsigned char)( ((double)g/(double)max)*255.0);
  2147.         data[ind + 2] = (unsigned char)( ((double)b/(double)max)*255.0);
  2148.       }
  2149.     }
  2150.   }
  2151.   else if (type == 5){
  2152.     data = new unsigned char [w * h];
  2153.     fread(data, 1, h * w, fp );
  2154.   }
  2155.   else if (type == 6){
  2156.     data = new unsigned char [w * h * 3];
  2157.     fread(data, 1, h * w * 3, fp );
  2158.   }
  2159.   else{
  2160.     printf(" *.ppm File must be of type P2, P3, P5 or P6 to be read.n");
  2161.     return FALSE;
  2162.   }
  2163.   fclose(fp);
  2164. //  printf("Finished Reading %sn", filename);
  2165.   return TRUE;
  2166. }
  2167. /*
  2168.  *  Grey = 0.299 Red + 0.587 Green + 0.114 Blue
  2169.  */
  2170. void Image::makeGray() {
  2171.   int ind;
  2172.   unsigned char *gray_data;
  2173.   if (3 == type || 6 == type){
  2174.     //convert to gray.
  2175.     gray_data = new unsigned char [height*width];
  2176.     type = 5;
  2177.     for (int y=0; y<height; y++ ){
  2178.       for (int x=0; x<width; x++){
  2179.         ind = (y*width +x) * 3;
  2180.         gray_data[x + y*width] =  (unsigned char)( 0.299*(double)data[ind] +
  2181.                                                    0.587*(double)data[ind + 1]+
  2182.                                                    0.114*(double)data[ind + 2]);
  2183.       }
  2184.     }
  2185.     delete [] data;
  2186.     data = gray_data;
  2187.   }
  2188. }
  2189. /*
  2190.  * Don't do this.
  2191.  *  Red =   Gray/0.299 
  2192.  *  Green = Gray/0.587 
  2193.  *  Blue =  Gray/.114
  2194.  */
  2195. void Image::makeColor() {
  2196.   int ind;
  2197.   unsigned char *color_data;
  2198.   if (2 == type || 5 == type){
  2199.     //convert to color.
  2200.     color_data = new unsigned char [3*width*height];
  2201.     type = 6;
  2202.     for (int y=0; y<height; y++ ){
  2203.       for (int x=0; x<width; x++){
  2204.         ind = y*width +x;
  2205.         color_data[ind*3] =    data[ind];
  2206.         color_data[ind*3+1] =  data[ind];
  2207.         color_data[ind*3+2] =  data[ind];
  2208.       }
  2209.     }
  2210.     delete [] data;
  2211.     data = color_data;
  2212.   }
  2213. }
  2214. void Image::makePowerOf2(int new_size){
  2215.   makeGray();
  2216.   new_size = findPowerOf2(new_size);
  2217.   if (new_size == width && width == height) return;
  2218.   pad(new_size, new_size);
  2219. }
  2220. /*
  2221.  */
  2222. Image* Image::getEqualizeTableImage () {
  2223.   Image *table_img;
  2224.   double max;
  2225.   double *pdf = new double[HIST_DIM];
  2226.   double lup[HIST_DIM];
  2227.   int i, j, x, y;
  2228.   getPDF(pdf, HIST_DIM, &max);
  2229.   for (i=0; i<HIST_DIM; i++){
  2230.     lup[i] = 0.0;
  2231.     for (j=0; j<=i; j++){
  2232.       lup[i] += pdf[j];
  2233.     }
  2234.     lup[i] *= (HIST_DIM-1);
  2235.   }
  2236.   table_img = new Image(HIST_DIM, HIST_DIM, 5);
  2237.   for (x=0; x<HIST_DIM; x++){
  2238.     y=0;
  2239.     for (; y<(int)lup[x]; y++ ){
  2240.       //color this black
  2241.       table_img->data[x + (HIST_DIM-y-1)*HIST_DIM] = 127;
  2242.     }
  2243.     for (; y<HIST_DIM; y++ ){
  2244.       //color this gray
  2245.       table_img->data[x + (HIST_DIM-y-1)*HIST_DIM] = 255;
  2246.     }
  2247.   }
  2248.   delete [] pdf;
  2249.   return table_img;
  2250. }
  2251. /*
  2252.  *  Gray = 0.299 Red + 0.587 Green + 0.114 Blue
  2253.  */
  2254. Image* Image::getHistImage (int flag) {
  2255.   int x, y, i;
  2256.   int hist[HIST_DIM];
  2257.   Image *hist_img;
  2258.   int offset=0;
  2259.   int val, max=0;
  2260.   for (int i=0; i<HIST_DIM; i++ ){
  2261.     hist[i] = 0;
  2262.   }
  2263.   switch(flag){
  2264.     case GRY://gray
  2265.       if (2 != type && 5 != type){
  2266.         printf("Error: Image::getHistImage(). must have gray image to get gray histogramn");
  2267.         return NULL;
  2268.       }
  2269.       for (y=0; y<height; y++ ){
  2270.         for (x=0; x<width; x++){
  2271.           val = ++hist[data[x + y*width]];
  2272.           if (val > max) max = val;
  2273.         }
  2274.       }
  2275.       break;
  2276.     case BLU://blue
  2277.       offset++;
  2278.     case GRN://green
  2279.       offset++;
  2280.     case RED://red
  2281.       if (6 != type && 3 != type){
  2282.         printf("Error: Image::getHistImage(). must have color image to get gray histogramn");
  2283.         return NULL;
  2284.       }
  2285.       for (y=0; y<height; y++ ){
  2286.         for (x=0; x<width; x++){
  2287.           val = ++hist[data[(x+y*width)*3 + offset]];
  2288.           if (val > max) max = val;
  2289.         }
  2290.       }
  2291.       break;
  2292.     default:
  2293.         printf("Error: Image::getHistImage(). Histogram type not recognizedn");
  2294.         return NULL;
  2295.       break;
  2296.   }
  2297.   for (i=0; i<HIST_DIM; i++ ){
  2298.     hist[i] = (int) (((double)hist[i]/(double)(max+20))*255.0);
  2299.   }
  2300.   hist_img = new Image(HIST_DIM, HIST_DIM, 5);
  2301.   for (x=0; x<HIST_DIM; x++){
  2302.     y=0;
  2303.     for (; y<hist[x]; y++ ){
  2304.       //color this black
  2305.       hist_img->data[x + (HIST_DIM-y-1)*HIST_DIM] = 0;
  2306.     }
  2307.     for (; y<HIST_DIM; y++ ){
  2308.       //color this white
  2309.       hist_img->data[x + (HIST_DIM-y-1)*HIST_DIM] = 255;
  2310.     }
  2311.   }
  2312.   return hist_img;
  2313. }
  2314. /*
  2315.  * Purpose: Fills neighboring regions of the same color(r, g, b) within tol
  2316.  *          with (f_r, f_g, f_b);
  2317.  */
  2318. void Image::seedFill(int x, int y,
  2319.                      unsigned char r, unsigned char g, unsigned char b,
  2320.                      unsigned char f_r, unsigned char f_g, unsigned char f_b,
  2321.                      unsigned char tol, int *area){
  2322.   int ind;
  2323.   if (type != 6 && type != 3){
  2324.     printf("Error: Image::seedFill(). must have color image to color shapesn");
  2325.     return;
  2326.   }
  2327.   //Check edges of image
  2328.   if (x < 0 || x >= width || y < 0 || y >= height)
  2329.     return;
  2330.   ind = (x+y*width)*3;
  2331.   //Already this color
  2332.   if ( f_r == data[ind] && 
  2333.        f_g == data[ind+1] && 
  2334.        f_b == data[ind+2])
  2335.     return;
  2336.   //Eccedes tolerance for color that we're filling
  2337.   if ( abs( (int)r-(int)data[ind] )   > (int)tol ||
  2338.        abs( (int)g-(int)data[ind+1] ) > (int)tol ||
  2339.        abs( (int)b-(int)data[ind+2] ) > (int)tol)
  2340.     return;
  2341.   //Set the color
  2342.   data[ind]   = f_r;
  2343.   data[ind+1] = f_g;
  2344.   data[ind+2] = f_b;
  2345.   (*area) ++;
  2346.   //Now recurse in 4 directions.
  2347.   seedFill(x+1, y, r, g, b, f_r, f_g, f_b, tol, area);
  2348.   seedFill(x-1, y, r, g, b, f_r, f_g, f_b, tol, area);
  2349.   seedFill(x, y+1, r, g, b, f_r, f_g, f_b, tol, area);
  2350.   seedFill(x, y-1, r, g, b, f_r, f_g, f_b, tol, area);
  2351. }
  2352. //OPENGL FUNCTIONS
  2353. //OPENGL FUNCTIONS
  2354. //OPENGL FUNCTIONS
  2355. //OPENGL FUNCTIONS
  2356. //OPENGL FUNCTIONS
  2357. void drawPointList(){
  2358.   int i;
  2359.   int *x, *y;
  2360.   int *seeds;
  2361.   x = cutter->plist.x;
  2362.   y = cutter->plist.y;
  2363.   seeds = cutter->plist.seed;
  2364.   glColor3f(1.0, 1.0, 0.0);
  2365.   glBegin(GL_POINTS);
  2366.     for (i=0; i<cutter->plist.num_points; i++){
  2367.       glVertex2i(x[i], cut_yres - y[i]);
  2368.     }
  2369.   glEnd();
  2370.   glPointSize(5.0);
  2371.   glColor3f(0.0, 1.0, 0.0);
  2372.   glBegin(GL_POINTS);
  2373.     for (i=0; i<cutter->plist.num_seeds; i++){
  2374.       glVertex2i(x[seeds[i]], cut_yres - y[seeds[i]]);
  2375.     }
  2376.   glEnd();
  2377.   glPointSize(1.0);
  2378. }
  2379. void drawWaveFront(){
  2380.   int j;
  2381.   GraphNode *scan_wave;
  2382.   glColor3f(1.0, 1.0, 1.0);
  2383.   glBegin(GL_POINTS);
  2384.     for (j=0; j<WAVE_SIZE; j++){
  2385.       for (scan_wave = cutter->wavefront[j];
  2386.            NULL != scan_wave;
  2387.            scan_wave = scan_wave->next){
  2388.         glVertex2i(scan_wave->x, cut_yres - scan_wave->y);
  2389.       }
  2390.     }
  2391.   glEnd();
  2392. }
  2393. /******************************************************************************
  2394.  * Purpose:  used by OpenGL to display the image.
  2395.  *****************************************************************************/
  2396. void display (void) {
  2397.   int wait_time;
  2398.   /*
  2399.   glMatrixMode (GL_PROJECTION );
  2400.     glLoadIdentity();
  2401.     glScalef(1.0, -1.0, 0.0); glMatrixMode (GL_MODELVIEW);
  2402.    */
  2403.   glRasterPos2i( 0, 0 );
  2404.   glDrawPixels(cut_xres, cut_yres, GL_RGB, GL_UNSIGNED_BYTE, cut_image->data);
  2405.   drawPointList();
  2406.   drawWaveFront();
  2407.   glFlush();//finish rendering it all and dump it out
  2408.   glutSwapBuffers();//since we're using the GLUT_DOUBLE
  2409.   wait_time = clock() + g_delay;
  2410.   while(wait_time > clock()) cutter->updateWaveFront();
  2411.   glutPostRedisplay();
  2412. }
  2413. /******************************************************************************
  2414.  * Purpose:  initialize the opengl stuff
  2415.  *****************************************************************************/
  2416. void init () {
  2417.   glClearColor (0.0, 0.0, 0.0, 0.0);
  2418.   glMatrixMode (GL_PROJECTION);
  2419.   glLoadIdentity();
  2420. //  gluOrtho2D(0.0, cut_xres, cut_yres, 0.0);
  2421.   gluOrtho2D(0.0, cut_xres, 0.0, cut_yres);
  2422. //  glOrtho(0.0, xres, 0.0, yres, yres, -yres);
  2423.   glColor3f ( 1.0, 1.0, 1.0);
  2424.   glClear (GL_COLOR_BUFFER_BIT );//need to clear the buffer to do the z-buffer
  2425.   glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  2426.   g_delay = (int)( (double)CLOCKS_PER_SEC/30.0);
  2427. }
  2428. /******************************************************************************
  2429.  * Purpose:  called by opengl to handle mouse events
  2430.  *****************************************************************************/
  2431. void mouseButton(int button, int state, int x, int y) {
  2432.   if (x<0) x = 0; if (x>cut_xres-1) x = cut_xres-1;
  2433.   if (y<0) y=0; if (y>cut_yres-1) y = cut_yres-1;
  2434. //  y = cut_yres - y;//switch the coordiates
  2435.   switch (button) {
  2436.     case GLUT_LEFT_BUTTON:
  2437.       if (state == GLUT_DOWN) {
  2438.         button_down = TRUE;
  2439.         /*
  2440.         cutter->plist.x[cutter->plist.num_points] = x;
  2441.         cutter->plist.y[cutter->plist.num_points] = y;
  2442.         cutter->plist.num_points ++;
  2443.         cutter->addSeedPoint();
  2444.         cutter->initWaveFront(x, y);
  2445.         cutter->updatePointList(x, y);
  2446.         //TEMPORARY
  2447.         cutter->updateWaveFront();
  2448.         */
  2449. //        cutter->dumpWaveFront();
  2450. //        exit(0);
  2451.       }
  2452.       if (state == GLUT_UP) {
  2453.         button_down = FALSE;
  2454.       }
  2455.       break;
  2456.       /*
  2457.       if (state == GLUT_UP) {
  2458.         button_down = FALSE;
  2459.         to.x = (double)x;
  2460.         to.y = (double)y;
  2461.         warp();
  2462.       }
  2463.       break;
  2464.       */
  2465.     case GLUT_RIGHT_BUTTON:
  2466.       break;
  2467.   }
  2468. }
  2469. /******************************************************************************
  2470.  * Purpose:  called by opengl to handle mouse events
  2471.  *****************************************************************************/
  2472. void mouseMotion(int x, int y) {
  2473. //  y = cut_yres - y;//switch the coordiates
  2474.   if (x<0) x = 0; if (x>cut_xres-1) x = cut_xres-1;
  2475.   if (y<0) y=0; if (y>cut_yres-1) y = cut_yres-1;
  2476.   cutter->updatePointList(x, y);
  2477. //  printf("Motionn");
  2478.   /*
  2479.   if (button_down == TRUE) {
  2480.     y = yres - y;//switch the coordiates
  2481.     to.x = (double)x;
  2482.     to.y = (double)y;
  2483.     warp();
  2484.   }
  2485.   */
  2486. }
  2487. /******************************************************************************
  2488.  * Purpose:  called by opengl to handle keyboard events
  2489.  *****************************************************************************/
  2490. void keyboard(unsigned char key, int x, int y) {
  2491.   if (x<0) x = 0; if (x>cut_xres-1) x = cut_xres-1;
  2492.   if (y<0) y=0; if (y>cut_yres-1) y = cut_yres-1;
  2493.   switch ( key ) {
  2494.     case 'w':
  2495.       /*
  2496.       printf("warp the imagen");
  2497.       warp();
  2498.       */
  2499.       break;
  2500.     case 'a':
  2501.       if (!button_down) return;
  2502.       cutter->plist.x[cutter->plist.num_points] = x;
  2503.       cutter->plist.y[cutter->plist.num_points] = y;
  2504.       cutter->plist.num_points ++;
  2505.       cutter->addSeedPoint();
  2506.       cutter->initWaveFront(x, y);
  2507.       cutter->updateWaveFront();
  2508.       cutter->updatePointList(x, y);
  2509.       display();
  2510.       break;
  2511.     case 'c':
  2512.       if (!button_down) return;
  2513.       cutter->saveBoundary(x, y);
  2514.       break;
  2515.     case 'd':
  2516.       cutter->dumpCostGraph();
  2517.       break;
  2518.     case 's':
  2519.       cutter->removeSeedPoint();
  2520.       cutter->updateWaveFront();
  2521.       cutter->updatePointList(x, y);
  2522.       display();
  2523.       break;
  2524.     case 'p':
  2525.       cutter->dumpPaths();
  2526.       break;
  2527.     case 'q':
  2528.       exit (0);
  2529.     case 'h':
  2530.       printf(" KEYS:n");
  2531.       printf(" 'a' - Add a seed pointn");
  2532.       printf(" 'c' - Close the loop selection and save out in a file "bound.txt"n");
  2533.       printf(" 'd' - save out a cost graph image as "cost_graph.pgm"n");
  2534.       printf(" 's' - Subtract a seed pointn");
  2535.       printf(" 'p' - Dump bad path info (should be nothing)n");
  2536.       printf(" 'q' - EXITn");
  2537.       break;
  2538.     case 't':
  2539. //      ContourMaker *cont = new ContourMaker(cut_image);
  2540. //      delete cont;
  2541.       break;
  2542.   }
  2543. }
  2544. //FUNCTIONS
  2545. //FUNCTIONS
  2546. //FUNCTIONS
  2547. //FUNCTIONS
  2548. /*
  2549.  * Note   The Id of the shape is equvalent to the Red component value.
  2550.  */
  2551. void getColorFromID (int id, unsigned char *r,  unsigned char *g,  unsigned char *b ){
  2552.   *r = (unsigned char) ((77*id+67)%255);
  2553.   *g = (unsigned char)((id*21+90)%255);
  2554.   *b = (unsigned char)(id%255);//(unsigned char)((id*13+90)%255);
  2555.   /*
  2556.   c_r = 0;
  2557.   c_g = (unsigned char)(255.0 * .33333);
  2558.   c_b = (unsigned char)(255.0 * .66666);
  2559.   */
  2560.         /*
  2561.         c_r = (c_r + 41)%255;
  2562.         c_g = (unsigned char)((255.0 - (double)c_r)*0.33333);
  2563.         c_b = (unsigned char)((255.0 - (double)c_r)*0.66666);
  2564.         */
  2565. }
  2566. void getNeighbor8(int x, int y, int *nx, int *ny, int id){
  2567.   switch (id){
  2568.     case 0:
  2569.       *nx = x+1;
  2570.       *ny = y;
  2571.       break;
  2572.     case 1:
  2573.       *nx = x+1;
  2574.       *ny = y-1;
  2575.       break;
  2576.     case 2:
  2577.       *nx = x;
  2578.       *ny = y-1;
  2579.       break;
  2580.     case 3:
  2581.       *nx = x-1;
  2582.       *ny = y-1;
  2583.       break;
  2584.     case 4:
  2585.       *nx = x-1;
  2586.       *ny = y;
  2587.       break;
  2588.     case 5:
  2589.       *nx = x-1;
  2590.       *ny = y+1;
  2591.       break;
  2592.     case 6:
  2593.       *nx = x;
  2594.       *ny = y+1;
  2595.       break;
  2596.     case 7:
  2597.       *nx = x+1;
  2598.       *ny = y+1;
  2599.       break;
  2600.     default:
  2601.       printf("Error in Cutter::getNeighbor8() id: %dn", id);
  2602.       break;
  2603.   }
  2604. }
  2605. void getNeighbor4(int x, int y, int *nx, int *ny, int id){
  2606.   switch (id){
  2607.     case 0:
  2608.       *nx = x+1;
  2609.       *ny = y;
  2610.       break;
  2611.     case 1:
  2612.       *nx = x;
  2613.       *ny = y-1;
  2614.       break;
  2615.     case 2:
  2616.       *nx = x-1;
  2617.       *ny = y;
  2618.       break;
  2619.     case 3:
  2620.       *nx = x;
  2621.       *ny = y+1;
  2622.       break;
  2623.     default:
  2624.       printf("Error in Cutter::getNeighbor4() id: %dn", id);
  2625.       break;
  2626.   }
  2627. }
  2628. /*
  2629.  * Purpose: Convolves two signals and deposites their value in the third signal
  2630.  *          The ordering is to pass the second signal across the first.
  2631.  */
  2632. void convolve1D (double *s1, int size1,
  2633.                  double *s2, int size2,
  2634.                  double *s3, int *size3){
  2635.   int i, n;
  2636.   *size3 = size1 + size2 -1;
  2637.   
  2638.   for (n=0; n<*size3; n++){
  2639.     s3[n] = 0;
  2640.     for (i=0; i<size1; i++){ //first signal
  2641.       if (n-i >= 0 && n-i < size2)
  2642.         s3[n] += s2[n-i]*s1[i]; 
  2643.     }
  2644.   }
  2645.   
  2646. }
  2647. Image* sobelFilter (Image *img2){
  2648.   Image *img3;
  2649.   unsigned char *data2, *data3;
  2650.   int w1, h1, w2, h2, w3, h3;
  2651.   int i, j, n, m;
  2652.   double total;
  2653.   double *sobel;
  2654. //  int total_time;
  2655. //  clock_t total_time;
  2656.   img2->makeGray();
  2657.   w1 = 3;
  2658.   h1 = 3;
  2659.   w2 = img2->width;
  2660.   h2 = img2->height;
  2661.   w3 = w1 + w2 -1;
  2662.   h3 = h1 + h2 -1;
  2663.   
  2664.   //kernel = new double [w1*h1];
  2665.   img3 = new Image(w3, h3, 5);
  2666.   sobel = new double[9];
  2667.     /*
  2668.     sobel[6] = -1; sobel[7] = 0; sobel[8] = 1;
  2669.     sobel[3] = -2; sobel[4] = 0; sobel[5] = 2;
  2670.     sobel[0] = -1; sobel[1] = 0; sobel[2] = 1;
  2671.     */
  2672.     sobel[6] = -2; sobel[7] = -2; sobel[8] = 0;
  2673.     sobel[3] = -2; sobel[4] = 0; sobel[5] = 2;
  2674.     sobel[0] = 0; sobel[1] = 2; sobel[2] = 2;
  2675.     /*
  2676.     sobel[6] = -1; sobel[7] = -2; sobel[8] = -1;
  2677.     sobel[3] = 0; sobel[4] = 0; sobel[5] = 0;
  2678.     sobel[0] = 1; sobel[1] = 2; sobel[2] = 1;
  2679.     */
  2680.   data2 = img2->data;
  2681.   data3 = img3->data;
  2682. //  total_time = clock();
  2683. //  printf("clock: %dn", (int)clock());
  2684.   for (i=0; i<w3; i++){
  2685.     for (j=0; j<h3; j++){
  2686.       total = 0.0;
  2687.       for (m=0; m<w1; m++){
  2688.         for (n=0; n<h1; n++){
  2689.           if ( (i-m) >= 0 &&
  2690.                (i-m) < w2 &&
  2691.                (j-n) >= 0 &&
  2692.                (j-n) < h2 ){
  2693.             total += sobel[m + n*w1] * (double)data2[i-m + (j-n)*w2];
  2694.           }
  2695.         }
  2696.       }
  2697.       total/=2.0;
  2698.       total+= 128;
  2699.       if (total > 255.0) total = 255.0;
  2700.       if (total < 0.0) total = 0.0;
  2701.       data3[i + j*w3] = (unsigned char) total;
  2702.     }
  2703.   }
  2704. //  for (i=0; i<1000000; i++);
  2705. //  printf("clock: %dn", (int)clock());
  2706. //  total_time = clock() - total_time;
  2707. //  printf("Time to convolve sobel with image: %1.20f seconds.n", (double)total_time/(double)CLOCKS_PER_SEC);
  2708. //  printf("Time to convolve sobel with image: %d seconds.n", total_time/CLOCKS_PER_SEC);
  2709.   
  2710.   return img3;
  2711. }
  2712. /*
  2713.  * Purpose: Convolves two signals and deposites their value in the third signal
  2714.  *          The ordering is to pass the second signal across the first.
  2715.  */
  2716. Image* convolve2D (Image *img1, Image *img2){
  2717.   Image *img3;
  2718.   unsigned char *data1, *data2, *data3;
  2719.   int w1, h1, w2, h2, w3, h3;
  2720.   int i, j, n, m;
  2721.   int total_time;
  2722.   double total;
  2723.   double *kernel;
  2724.   img1->makeGray();
  2725.   img2->makeGray();
  2726.   w1 = img1->width;
  2727.   h1 = img1->height;
  2728.   w2 = img2->width;
  2729.   h2 = img2->height;
  2730.   w3 = w1 + w2 -1;
  2731.   h3 = h1 + h2 -1;
  2732.   
  2733.   kernel = new double [w1*h1];
  2734.   img3 = new Image(w3, h3, 5);
  2735.   data1 = img1->data;
  2736.   data2 = img2->data;
  2737.   data3 = img3->data;
  2738.   total = 0.0;
  2739.   for (i=0; i<w1*h1; i++)
  2740.     total += (double)data1[i];
  2741.   for (i=0; i<w1*h1; i++)
  2742.     kernel[i] = (double)(data1[i])/total;
  2743.   total_time = clock();
  2744.   for (i=0; i<w3; i++){
  2745.     for (j=0; j<h3; j++){
  2746.       total = 0.0;
  2747.       for (m=0; m<w1; m++){
  2748.         for (n=0; n<h1; n++){
  2749.           if ( (i-m) >= 0 &&
  2750.                (i-m) < w2 &&
  2751.                (j-n) >= 0 &&
  2752.                (j-n) < h2 ){
  2753.             total += kernel[m + n*w1] * (double)data2[i-m + (j-n)*w2];
  2754.           }
  2755.         }
  2756.       }
  2757.       if (total > 255.0) total = 255.0;
  2758.       data3[i + j*w3] = (unsigned char) total;
  2759.     }
  2760.   }
  2761.   total_time = clock() - total_time;
  2762.   printf("Time to convolve the 2 images: %f seconds.n", (double)total_time/(double)CLOCKS_PER_SEC);
  2763.   
  2764.   
  2765.   return img3;
  2766. }
  2767. void enhanceEdge1D (double *s1, int size1,
  2768.                     double *s2, int size2,
  2769.                     double *s3, int *size3){
  2770.   double s4[MAX_SIGNAL_SIZE];
  2771.     int size4;
  2772.   double s5[MAX_SIGNAL_SIZE];
  2773.     int size5;
  2774.   double delta[MAX_SIGNAL_SIZE];
  2775.     double total = 0.0;
  2776.     int size_delta;
  2777.   int i;
  2778.     total = 0.0;
  2779.     size_delta = size2;
  2780.     for (i=0; i<size_delta; i++){
  2781.         delta[i] = 0.0;
  2782.         total += s2[i];
  2783.     }
  2784.     delta[size2/2] = 2.0;
  2785.     //normalize the gaussian
  2786.     for (i=0; i<size2; i++){
  2787.         s2[i]/=total;
  2788.     }
  2789.     //Make a delta function the size of the gausian to make the signal have the 
  2790.     //same shift.
  2791.     //convolve s1 and s2
  2792.     convolve1D (s1, size1, s2, size2, s4, &size4);
  2793.     convolve1D (s1, size1, delta, size_delta, s5, &size5);
  2794.   
  2795.     /*
  2796.     //clear it out.
  2797.     for(i=0; i<*size3; i++){
  2798.       s3[i] = 0.0;
  2799.     }
  2800.     */
  2801.   
  2802.    *size3 = size5;
  2803.     for(i=0; i<size5; i++){
  2804.         printf("i: %d s3:%f =  s5:%f - s4:%fn",i, s3[i], s5[i], s4[i]);
  2805.       s3[i] = s5[i] - s4[i];
  2806.       //s3[i] = -s4[i];// - s4[i];
  2807.     }
  2808. }
  2809. /*
  2810.  * c1 is the image.
  2811.  * c2 is the filter.
  2812.  * c3 is  2*c1 - c2*c1;
  2813.  */
  2814. void enhanceEdge2D (Complex *c1, int size1,
  2815.                     Complex *c2, int size2,
  2816.                     Complex *c3, int *size3){
  2817.   int i;
  2818.   Complex m;
  2819.   *size3 = size1;
  2820.     
  2821.   for(i=0; i<size1*size1; i++){
  2822.     m = mult(c1[i], c2[i]);
  2823.     c3[i].r = 2*c1[i].r - m.r;
  2824.     c3[i].i = 2*c1[i].i - m.i;
  2825.   }
  2826. }
  2827. void medianFilter1D(int range, double *s1, int size1, double *s2, int *size2){
  2828.   int i, j, median;
  2829.   double *samples;
  2830.   
  2831.   median = range/2;
  2832.   samples = new double[range];
  2833.   *size2 = size1;
  2834.   for(i=0; i<*size2; i++)
  2835.     s2[i] = 0.0;
  2836.   for(i=median; i<size1-median; i++){
  2837.     for(j=-median; j<median; j++)
  2838.       samples[median+j] = s1[i+j];
  2839.     slowSort(samples, range);
  2840.     s2[i] = samples[median];
  2841.   }
  2842.   delete [] samples;
  2843. }
  2844. void medianFilter2D(int median, double *s1, int w, int h, double *s2){
  2845.   int i, j, k, l, range;
  2846.   double *samples;
  2847.   
  2848.   range = 2*median + 1;
  2849.   /*
  2850.   if (range%2 != 0) range ++;
  2851.   median = range/2;
  2852.   */
  2853.   samples = new double[range*range];
  2854.   for(i=0; i<w*h; i++)
  2855.     s2[i] = 0.0;
  2856.   for(i=median; i<w-median; i++){
  2857.     for(j=median; j<h-median; j++){
  2858.       //now get samples for the median filter
  2859.       for(k=-median; k<=median; k++)
  2860.         for(l=-median; l<=median; l++)
  2861.           samples[ k+median + (l+median)*range] = s1[ k+i + (l+j)*w] ;
  2862.       slowSort(samples, range*range);
  2863.       /*
  2864.       if (s1[i+j*w] > 0){
  2865.         for (q=0; q<range*range; q++){
  2866.           printf("sample[%d]: %fn", q, samples[q]);
  2867.         }
  2868.       }
  2869.       */
  2870.       s2[i+j*w] = samples[(range*range-1)/2];
  2871.      // s2[i+j*w] = s1[i+j*w];
  2872.     }
  2873.   }
  2874.   delete [] samples;
  2875. }
  2876. void normalize(double *d, int size){
  2877.   int i;
  2878.   double total;
  2879.   total = 0.0;
  2880.   for (i=0; i<size; i++){
  2881.     total += d[i];
  2882.   }
  2883.   for (i=0; i<size; i++){
  2884.     d[i]/=total;
  2885.   }
  2886. }
  2887. void slowSort(double *a, int size){
  2888.   int i, j;
  2889.     double temp;
  2890.     /*
  2891.   printf("pre sortn");
  2892.   for (i=0; i<size; i++){
  2893.     printf("a[%d]: %fn",i , a[i]);
  2894.   }
  2895.   */
  2896.   for (i=0; i<size; i++){
  2897.     for (j=i+1; j<size; j++){
  2898.             //swap
  2899.       if (a[j] < a[i]){
  2900.         temp = a[i];
  2901.         a[i] = a[j];
  2902.         a[j] = temp;
  2903.       }
  2904.     }
  2905.   }
  2906.   /*
  2907.   printf("post sortn");
  2908.   for (i=0; i<size; i++){
  2909.     printf("a[%d]: %fn",i , a[i]);
  2910.   }
  2911.   */
  2912. }
  2913. void threshold(double *a, int size, double val, double bot, double top){
  2914.   int i;
  2915.   for(i=0; i<size; i++){
  2916.     if (bot < a[i] && a[i] < top){
  2917.       a[i] = val;
  2918.     }
  2919.   }
  2920. }
  2921. /*
  2922.  * finds the average for a bunch of pixels and thresh holds the pixel based
  2923.  * on this average
  2924.  */
  2925. void threshBlock2D(double *a, int w, int h, int range, double val, double bot, double top, double *b){
  2926.   int i, j, l, m;
  2927.   double ave, count;
  2928.   for(i=0; i<w; i++) for(j=0; j<h; j++) b[i +j*w] = 0.0;
  2929. //  for(i=w-range; i<w; i++) for(j=h-range; j<h; j++) b[i +j*w] = 0.0;
  2930.   for(i=range; i<w-range; i++){
  2931.     for(j=range; j<h-range; j++){
  2932.       /*
  2933.       count = 1.0;
  2934.       ave = a[i + j*w];
  2935.       */
  2936.       count = 0.0;
  2937.       ave = 0.0;//a[i + j*w];
  2938.       for(l=i-range; l<i+range; l++){
  2939.         for(m=j-range; m<j+range; m++){
  2940.           ave += a[l + m*w];
  2941.           count ++;
  2942.         }
  2943.       }
  2944.       ave /= count;
  2945.       if (bot < ave && ave < top){ b[i + j*w] = val; }
  2946.       else{ b[i + j*w] = a[i + j*w]; }
  2947.     }
  2948.   }
  2949. }
  2950. void FFT1D(double *f, int size, Complex *F){
  2951.   int i;
  2952.   int offset;
  2953.   int new_size;
  2954.   new_size = findPowerOf2(size);
  2955.   if (new_size != size){
  2956.     offset = (new_size-size)/2;
  2957.     for (i=size-1; i>=0; i--){
  2958.       f[i+offset] = f[i];
  2959.     }
  2960.     //pad with zeros
  2961.     for (i=new_size-1; i>=size; i--){
  2962.       f[i] = 0.0;
  2963.       f[i-size] = 0.0;
  2964.     }
  2965.   }
  2966.   //now assign values into complex array
  2967.   for (i=0; i<new_size; i++){
  2968.     F[i].r = f[i];
  2969.     F[i].i = 0.0;
  2970.   }
  2971.   FFT1Drec(F, new_size, new_size/2, 1);
  2972. }
  2973. /*
  2974.  * Note: size refers to the size of both dimensions
  2975.  */
  2976. void FFT2D(double *f, int size, Complex *F){
  2977.   int i, r, c;
  2978.   int N;
  2979.   Complex *X;
  2980.   N = findPowerOf2(size);
  2981.   if (N != size){
  2982.     printf(" Error, image is not a power of 2n");
  2983.     return;
  2984.   }
  2985.   X = new Complex[N];
  2986.   //now assign values into complex array
  2987.  
  2988.   for (i=0; i<N*N; i++){
  2989.     F[i].r = f[i];
  2990.     F[i].i = 0.0;
  2991.   }
  2992.   //First the rows
  2993.   for (r=0; r<N; r++){
  2994.     for (i=0; i<N; i++){
  2995.       X[i].r = F[i + r*N].r; 
  2996.       X[i].i = F[i + r*N].i;
  2997.     }
  2998.     FFT1Drec(X, N, N/2, 1);
  2999.     for (i=0; i<N; i++){
  3000.       F[i + r*N].r = X[i].r;
  3001.       F[i + r*N].i = X[i].i;
  3002.     }
  3003.   }
  3004.   //now the columns
  3005.   for (c=0; c<N; c++){
  3006.     for (i=0; i<N; i++){
  3007.       X[i].r = F[c + i*N].r;
  3008.       X[i].i = F[c + i*N].i;
  3009.     }
  3010.     FFT1Drec(X, N, N/2, 1);
  3011.     for (i=0; i<N; i++){
  3012.       F[c + i*N].r = X[i].r;
  3013.       F[c + i*N].i = X[i].i;
  3014.     }
  3015.   }
  3016.   delete [] X;
  3017. }
  3018. /*
  3019.  * Note: size refers to the size of both dimensions
  3020.  */
  3021. void invFFT2D(Complex *F, int size, double *f){
  3022.   int i, r, c;
  3023.   int N;
  3024.   Complex *X;
  3025.   N = findPowerOf2(size);
  3026.   if (N != size){
  3027.     printf(" Error, image is not a power of 2n");
  3028.     return;
  3029.   }
  3030.   X = new Complex[N];
  3031.   //now assign values into complex array
  3032.  
  3033.   //First the rows
  3034.   for (r=0; r<N; r++){
  3035.     for (i=0; i<N; i++){
  3036.       X[i].r = F[i + r*N].r; 
  3037.       X[i].i = F[i + r*N].i;
  3038.     }
  3039.     FFT1Drec(X, N, N/2, -1);
  3040.     for (i=0; i<N; i++){
  3041.       F[i + r*N].r = X[i].r;
  3042.       F[i + r*N].i = X[i].i;
  3043.     }
  3044.   }
  3045.   //now the columns
  3046.   for (c=0; c<N; c++){
  3047.     for (i=0; i<N; i++){
  3048.       X[i].r = F[c + i*N].r;
  3049.       X[i].i = F[c + i*N].i;
  3050.     }
  3051.     FFT1Drec(X, N, N/2, -1);
  3052.     for (i=0; i<N; i++){
  3053.       F[c + i*N].r = X[i].r;
  3054.       F[c + i*N].i = X[i].i;
  3055.     }
  3056.   }
  3057.   for (i=0; i<N*N; i++){
  3058.     f[i] = F[i].r;
  3059.   }
  3060.   delete [] X;
  3061. }
  3062. void FFT1Drec(Complex *F, int size, int mid, int dir){
  3063.   int i, new_size;
  3064.   double tempr, tempi;
  3065.   Complex he, e;
  3066.   Complex *G = new Complex[size/2];
  3067.   Complex *H = new Complex[size/2];
  3068.   double g_total;
  3069.   double h_total;
  3070.   if (size == 2){
  3071.     tempr = F[0].r + F[1].r;
  3072.     tempi = F[0].i + F[1].i;
  3073.     F[1].r = F[0].r - F[1].r;
  3074.     F[1].i = F[0].i - F[1].i;
  3075.     F[0].r = tempr;
  3076.     F[0].i = tempi;
  3077.     delete [] G;
  3078.     delete [] H;
  3079.     return;
  3080.   }
  3081.   else{
  3082.     h_total = g_total = 0.0;
  3083.     for(i=0, new_size=0; i<size;){
  3084.       G[new_size].r = F[i].r;
  3085.       G[new_size].i = F[i].i;
  3086.       g_total += G[new_size].getMag();
  3087.       i++;
  3088.       H[new_size].r = F[i].r;
  3089.       H[new_size].i = F[i].i;
  3090.       h_total += H[new_size].getMag();
  3091.       i++;
  3092.       new_size ++;
  3093.     }
  3094.     if (g_total > DBL_TOL)
  3095.       FFT1Drec(G, new_size, mid, dir);
  3096.     if (h_total > DBL_TOL)
  3097.       FFT1Drec(H, new_size, mid, dir);
  3098.     for(i=0; i<size; i++){
  3099.       /*
  3100.       e.r = cos( (TWO_PIE*(double)(i-mid))/(double)size ); 
  3101.       e.i = -dir*sin( (TWO_PIE*(double)(i-mid))/(double)size ); 
  3102.       */
  3103.       e.r = cos( (TWO_PIE*(double)i)/(double)size ); 
  3104.       e.i = -dir*sin( (TWO_PIE*(double)i)/(double)size ); 
  3105.       he = mult(e, H[i%new_size]);
  3106.       F[i].r = G[i%new_size].r + he.r;
  3107.       F[i].i = G[i%new_size].i + he.i;
  3108.     }
  3109.   }
  3110.   delete [] G;
  3111.   delete [] H;
  3112. }
  3113. void invFFT1D(Complex *F, int size, double *f){
  3114.   int i;
  3115.   int offset;
  3116.   int new_size;
  3117.   double invN;
  3118.   new_size = findPowerOf2(size);
  3119.   offset = (new_size-size)/2;
  3120.   for (i=size-1; i>=0; i--){
  3121.     F[i+offset].r = F[i].r;
  3122.     F[i+offset].i = F[i].i;
  3123.   }
  3124.   //pad with zeros
  3125.   for (i=new_size-1; i>=size; i--){
  3126.     F[i].r = 0.0;
  3127.     F[i-size].r = 0.0;
  3128.     F[i].i = 0.0;
  3129.     F[i-size].i = 0.0;
  3130.   }
  3131.   FFT1Drec(F, new_size, new_size/2, -1);
  3132.   //now assign values into complex array
  3133.   invN = 1.0/(double)new_size;
  3134.   for (i=0; i<new_size; i++){
  3135.     f[i] = invN*F[i].r;
  3136.   }
  3137. }
  3138. void FourierTransform1D(double *f, int size, Complex *F){
  3139.   int u, k, half;
  3140.   double N;
  3141.   N = (double) size;
  3142.   half = size/2;
  3143.   for (u=0; u<size; u++){
  3144.     F[u].r = 0.0;
  3145.     F[u].i = 0.0;
  3146.     for (k=0; k<size; k++){
  3147.       //F[u].r += f[k] * cos( (TWO_PIE*(double)(u)*(double)(k))/N ); 
  3148.       F[u].r += f[k] * cos( (TWO_PIE*(double)(u-half)*(double)(k-half))/N ); 
  3149.       F[u].i -= f[k] * sin( (TWO_PIE*(double)(u-half)*(double)(k-half))/N ); 
  3150.     }
  3151.   }
  3152. }
  3153. void invFourierTransform1D(Complex *F, int size, double *f){
  3154.   int u, k, half;
  3155.   double N, invN;
  3156.   Complex f_c[size];
  3157.   Complex ex;
  3158.   N = (double) size;
  3159.   invN = 1.0/N;
  3160.   half = size/2;
  3161.   for (k=0; k<size; k++){
  3162.     f_c[k].r = 0.0;
  3163.     f_c[k].i = 0.0;
  3164.     for (u=0; u<size; u++){
  3165.       ex.r = cos( (TWO_PIE*(double)(u-half)*(double)(k-half))/N );
  3166.       ex.i = sin( (TWO_PIE*(double)(u-half)*(double)(k-half))/N );
  3167.       ex = mult(ex, F[u]);
  3168.       f_c[k] = add(f_c[k], ex);
  3169.     }
  3170.   }
  3171.   for (k=0; k<size; k++){
  3172.     f[k] = invN*f_c[k].r;
  3173.   }
  3174. }
  3175. int findPowerOf2(int size){
  3176.   int newdim=1;
  3177.   int power =0;
  3178.   int tmp_old = 0, test = 0;
  3179.   tmp_old = test = size;
  3180.   while(size != 0){
  3181.     ++power;
  3182.     size >>=1;
  3183.   }
  3184.   newdim <<= power;
  3185.   if ( (test <<= 1) == newdim)
  3186.     return tmp_old;//the original value is already a power of two
  3187.   return newdim;
  3188. }
  3189. Complex add(Complex c1, Complex c2){
  3190.   Complex result;
  3191.   result.r = c1.r + c2.r;
  3192.   result.i = c1.i + c2.i;
  3193.   return result;
  3194. }
  3195. Complex sub(Complex c1, Complex c2){
  3196.   Complex result;
  3197.   result.r = c1.r - c2.r;
  3198.   result.i = c1.i - c2.i;
  3199.   return result;
  3200. }
  3201. Complex mult(Complex c1, Complex c2){
  3202.   Complex result;
  3203.   /*
  3204.   double phase;
  3205.   double mag;
  3206.   mag = c1.getMag() * c2.getMag();
  3207.   phase = c1.getPhase() + c2.getPhase();
  3208.   result.r = mag*cos(phase);
  3209.   result.i = mag*sin(phase);
  3210.   */
  3211.   result.r = c1.r*c2.r - c1.i*c2.i;
  3212.   result.i = c1.r*c2.i + c1.i*c2.r;
  3213.   return result;
  3214. }
  3215. Complex div(Complex c1, Complex c2){
  3216.   double mag1, magsq;
  3217.   Complex result;
  3218.   /*
  3219.   double phase;
  3220.   double mag;
  3221.   
  3222.   mag = c2.getMag();
  3223.   if (mag < DBL_TOL){
  3224.     result.r = 0.0;
  3225.     result.i = 0.0;
  3226.   }
  3227.   else{
  3228.     mag = c1.getMag() / c2.getMag();
  3229.     phase = c1.getPhase() - c2.getPhase();
  3230.     result.r = mag*cos(phase);
  3231.     result.i = mag*sin(phase);
  3232.   }
  3233.   */
  3234.   mag1 = c1.getMag();
  3235.   magsq = c2.getMag();
  3236.   if (magsq < DBL_TOL || mag1 < DBL_TOL){
  3237.     result.r = 0.0;
  3238.     result.i = 0.0;
  3239.   }
  3240.   else{
  3241.     magsq*=magsq;
  3242.     result.r = (c1.r*c2.r + c1.i*c2.i)/magsq;
  3243.     result.i = (-c1.r*c2.i + c1.i*c2.r)/magsq;
  3244.   }
  3245.   return result;
  3246. }
  3247. void mult(double *s1, int size1, double *s2, int size2, double *s3, int *size3){
  3248.   int half_min, i;
  3249.   int half1, half2, half3;
  3250.   half1 = size1/2;
  3251.   half2 = size2/2;
  3252.   if (size1 > size2){ half3 = half1; *size3 = size1; half_min = half2; }
  3253.   else{ half3 = half2; *size3 = size2; half_min = half1; }
  3254.   //clean out the signal
  3255.   for (i=0; i<*size3; i++){ s3[i] = 0.0; }
  3256.   //multiply the two
  3257.   for (i=0; i<half_min; i++){
  3258.     s3[half3 + i] = s1[half1 + i]*s2[half2 + i];
  3259.     s3[half3 - i] = s1[half1 - i]*s2[half2 - i];
  3260.   }
  3261. }
  3262. void add(double *s1, int size1, double *s2, int size2, double *s3, int *size3){
  3263.   int half_min, i;
  3264.   int half1, half2, half3;
  3265.   half1 = size1/2;
  3266.   half2 = size2/2;
  3267.   if (size1 > size2){
  3268.     half3 = half1; *size3 = size1; half_min = half2;
  3269.     for (i=0; i<*size3; i++){ s3[i] = s1[i]; }
  3270.   }
  3271.   else{
  3272.     half3 = half2; *size3 = size2; half_min = half1;
  3273.     for (i=0; i<*size3; i++){ s3[i] = s2[i]; }
  3274.   }
  3275.   //add the two
  3276.   for (i=0; i<half_min; i++){
  3277.     s3[half3 + i] = s1[half1 + i] + s2[half2 + i];
  3278.     s3[half3 - i] = s1[half1 - i] + s2[half2 - i];
  3279.   }
  3280. }
  3281. void mult(Complex *s1, int size1, Complex *s2, int size2, Complex *s3, int *size3){
  3282.   int half_min, i;
  3283.   int half1, half2, half3;
  3284.   half1 = size1/2;
  3285.   half2 = size2/2;
  3286.   if (size1 > size2){ half3 = half1; *size3 = size1; half_min = half2; }
  3287.   else{ half3 = half2; *size3 = size2; half_min = half1; }
  3288.   //clean out the signal
  3289.   for (i=0; i<*size3; i++){ s3[i].r = 0.0; s3[i].i = 0.0; }
  3290.   //multiply the two
  3291.   for (i=0; i<half_min; i++){
  3292.     s3[half3 + i] = mult(s1[half1 + i], s2[half2 + i]);
  3293.     s3[half3 - i] = mult(s1[half1 - i], s2[half2 - i]);
  3294.   }
  3295. }
  3296. /*
  3297.  * s3  = s1/s2
  3298.  */
  3299. void div(Complex *s1, int size1, Complex *s2, int size2, Complex *s3, int *size3){
  3300.   int half_min, i;
  3301.   int half1, half2, half3;
  3302.   half1 = size1/2;
  3303.   half2 = size2/2;
  3304.   if (size1 > size2){ half3 = half1; *size3 = size1; half_min = half2; }
  3305.   else{ half3 = half2; *size3 = size2; half_min = half1; }
  3306.   //clean out the signal
  3307.   for (i=0; i<*size3; i++){ s3[i].r = 0.0; s3[i].i = 0.0; }
  3308.   //multiply the two
  3309.   /*
  3310.   for (i=0; i<half_min; i++){
  3311.     s3[half3 + i] = div(s1[half1 + i], s2[half2 + i]);
  3312.     s3[half3 - i] = div(s1[half1 - i], s2[half2 - i]);
  3313.   }
  3314.   */
  3315.   for (i=0; i<size1; i++){
  3316.     s3[i] = div(s1[i], s2[i]);
  3317.   }
  3318. }
  3319. /*
  3320.  * Purpose: Get bounding boxes for connected components that have enough area.
  3321.  *          if the  is big enough, do a local threshold and union the
  3322.  *          results into a final image.
  3323.  */
  3324. int  boxThreshold(Image *img, Image *img_out){
  3325.   double thresh;
  3326.   unsigned char t;
  3327.   double *img_dbl, *img_dbl2;
  3328.   Image img_temp;
  3329.   Image img2;
  3330.   int w, h;
  3331.   int bot_x, bot_y, top_x, top_y;
  3332.   int range = 3;
  3333.   unsigned char num_shapes, i;
  3334.   int j;
  3335.   //char name[255];
  3336.   
  3337.   img_temp.copy(img);
  3338.   w = img->width;
  3339.   h = img->height;
  3340.   img_dbl = img->getDouble(0);
  3341.   img_dbl2 = new double[w*h];
  3342.   printf("DBG07w: %d h: %dn", w, h);
  3343.   for (j=0; j<w*h; j++){
  3344.     img_dbl2[j] = 0.0;
  3345.   }
  3346.   thresh = getAverage(img_dbl, w*h, -1.0, 0);
  3347.   //make black
  3348.   threshold(img_dbl, w*h, 0, 0, 64);
  3349.   //make white
  3350.   threshold(img_dbl, w*h, 255, 255-64, 255);
  3351.   threshBlock2D(img_dbl, w, h, range, 255.0, thresh, 255.0, img_dbl2);
  3352.   threshold(img_dbl2, w*h, 0, 0, 250);
  3353.   img_temp.signalToImage2D(img_dbl2, w, h);
  3354.   img_temp.write("out_simple_thresh.pgm");
  3355.   //TEMP
  3356.   img2.copy(&img_temp);
  3357.   img2.erode(range);
  3358.   img2.erode(range);
  3359.   img2.dilate(range);
  3360.   img2.dilate(range);
  3361.   img2.write("out_simple_after_erode_dilate.pgm");
  3362.   num_shapes = img2.colorShapes(1, (w*h)>>9);
  3363.   img2.write("out_simple_thresh_color.pgm");
  3364.   num_shapes = img_temp.colorShapes(1, (w*h)>>9);
  3365. //  img_temp.write("out_simple_thresh_color.pgm");
  3366. //  num_shapes = 5;
  3367.   img_out->clear(0);
  3368.   //Now refine around each object
  3369.   for (i=1; i<=num_shapes; i++){
  3370.     img_temp.getItemBounds(i, &bot_x, &bot_y, &top_x, &top_y);
  3371.     bot_x -= (top_x- bot_x)/3;
  3372.     bot_y -= (top_y- bot_y)/3;
  3373.     top_x += (top_x- bot_x)/3;
  3374.     top_y += (top_y- bot_y)/3;
  3375.     if (top_x > img_temp.width) top_x = img_temp.width;
  3376.     if (top_y > img_temp.height) top_y = img_temp.height;
  3377.     if (bot_x < 0) bot_x = 0;
  3378.     if (bot_y < 0) bot_y = 0;
  3379.     t = img->getBoundAverage(bot_x, bot_y, top_x, top_y);
  3380.     img_out->unionThreshRegion(img, bot_x, bot_y, top_x, top_y, range, t);
  3381.     //img_out->write(name);
  3382.   }
  3383.   img_out->write("out_box_thresh.pgm");
  3384.   //Make back to black and white.
  3385.   img_out->threshold(2);
  3386.   /*
  3387.   delete [] img_dbl;
  3388.   img_dbl = img_out->getDouble(0);
  3389.   medianFilter2D(1, img_dbl, w, h, img_dbl2);
  3390.   img_out->signalToImage2D(img_dbl2, w, h); img_out->write("out_box_thresh_median.pgm");
  3391.   */
  3392.   img_out->erode(range);
  3393.   img_out->erode(range);
  3394.   img_out->write("out_after_erode.pgm");
  3395.   /*
  3396.   delete [] img_dbl;
  3397.   img_dbl = img_out->getDouble(0);
  3398.   medianFilter2D(1, img_dbl, w, h, img_dbl2);
  3399.   img_out->signalToImage2D(img_dbl2, w, h); img_out->write("out_box_thresh_median.pgm");
  3400.   */
  3401.   img_out->dilate(range);
  3402.   img_out->dilate(range);
  3403.   img_out->write("out_dilate.pgm");
  3404.   num_shapes = img_out->colorShapes(1, (w*h)>>9);
  3405.   img_out->write("out_after_color.pgm");
  3406.   delete [] img_dbl;
  3407.   delete [] img_dbl2;
  3408.   return num_shapes;
  3409. }
  3410. void copyDouble(double *dest, double *src, int size){
  3411.   int i;
  3412.   for (i=0; i<size; i++){
  3413.     dest[i] = src[i];
  3414.   }
  3415. }
  3416. void dilate(double *in, int w, int h, double *out, double *filter, int size, double tol){
  3417.   int i, j, l, m;
  3418.   int offset;
  3419.   offset = (size-1)/2;
  3420.   for(i=0; i<h*w; i++){
  3421.     out[i] = 0.0;
  3422.   }
  3423.   for(i=offset; i<w-offset; i++){
  3424.     for(j=offset; j<h-offset; j++){
  3425.       //If the pixel is set.
  3426.       if(in[i + j*w] > tol){
  3427.         for(l=0; l<size; l++){
  3428.           for(m=0; m<size; m++){
  3429.             if(filter[l + m*size] > tol)
  3430.               out[i-offset+l+(j-offset+m)*w] = filter[l + m*size];
  3431.           }
  3432.         }
  3433.       }
  3434.     }
  3435.   }
  3436. }
  3437. /*
  3438.  */
  3439. void erode(double *in, int w, int h, double *out, double *filter, int size, double tol){
  3440.   int i, j, l, m;
  3441.   int offset;
  3442.   int set;
  3443.   offset = (size-1)/2;
  3444.   for(i=0; i<h*w; i++){
  3445.     out[i] = 0.0;
  3446.   }
  3447.   for(i=offset; i<w-offset; i++){
  3448.     for(j=offset; j<h-offset; j++){
  3449.       set = 1;
  3450.       for(l=0; l<size; l++){
  3451.         for(m=0; m<size; m++){
  3452.           if(filter[l + m*size] > tol && in[i-offset+l+(j-offset+m)*w]<=tol){
  3453.               set = 0;
  3454.               m = l = size;//break out
  3455.           }
  3456.         }
  3457.       }
  3458.       //Now set the pixel in the outgoing image if needs be.
  3459.       if (1 == set)
  3460.         out[i + j*w] = in[i + j*w];
  3461.       else
  3462.         out[i + j*w] = 0.0;
  3463.     }
  3464.   }
  3465. }
  3466. /*
  3467.  * purpose: returns the average value of the all non zero pixels in the scene
  3468.  */
  3469. double getAverage(double *s, int size, double zero_intensity, double tol){
  3470.   int i;
  3471.   double total;
  3472.   double count;
  3473.   total = 0.0;
  3474.   count = 0.0;
  3475.   for (i=0; i<size; i++){
  3476.     if (s[i]+tol > zero_intensity && s[i]-tol < zero_intensity)
  3477.       continue;
  3478.     total += s[i];
  3479.     count ++;
  3480.   }
  3481.   return total/count;
  3482. }
  3483. /*
  3484.  * returns the percentage of this value (within a tol) in the array
  3485.  */
  3486. double getPercentValue(double *s, int size, double val, double tol){
  3487.   int i;
  3488.   double total;
  3489.   total = 0.0;
  3490.   for (i=0; i<size; i++){
  3491.     if ( s[i] + tol >= val && s[i] - tol <= val){
  3492.       total ++;
  3493.     }
  3494.   }
  3495.   return total/(double)size;
  3496. }
  3497. /*
  3498.  * get the difference between two arrays
  3499.  */
  3500. double getDiff(double *s1, double *s2, int size){
  3501.   double total;
  3502.   int i;
  3503.   total = 0.0;
  3504.   for (i=0; i<size; i++){
  3505.     total += fabs(s1[i] - s2[i]);
  3506.   }
  3507.   return total;
  3508. }
  3509. /*
  3510.  * get the average difference between two arrays
  3511.  */
  3512. double getAveDiff(double *s1, double *s2, int size){
  3513.   return getDiff(s1, s2, size)/(double)size;
  3514. }
  3515. void invert(Complex *s1, int size1){
  3516.   double magsq; 
  3517.   for (int i=0; i<size1; i++){
  3518.     magsq = s1[i].getMag()*s1[i].getMag();
  3519.     if (magsq < DBL_TOL){
  3520.       s1[i].r = 0.0;
  3521.       s1[i].i = 0.0; 
  3522.     }
  3523.     else{
  3524.       s1[i].r /= magsq;
  3525.       s1[i].i /= -magsq;
  3526.     }
  3527.   }
  3528. }
  3529. /*
  3530.  * Purpose: rotates the image the specified amount.
  3531.  * Note:    The value passed in is assumed to be in degrees
  3532.  */
  3533. void rotate2D(double *s, int n, double *d, double theta){
  3534.   int nx, ny;
  3535.   double x1, x2, y1, y2; 
  3536.   double sint, cost;
  3537.   double xc, yc, tx, ty;
  3538.   double sample1, sample2;
  3539.   int ix2, iy2;
  3540.   //convert to radians
  3541.   theta *= PIE/180.0;
  3542.   sint = sin(theta);
  3543.   cost = cos(theta);
  3544.   xc = (double)n/2.0;
  3545.   yc = (double)n/2.0;
  3546.   //To do the backwards mapping, we'll just take samples in the new image
  3547.   //and rotate back to the previous to get the value
  3548.   for (ny = 0; ny<n; ny++){
  3549.     for (nx = 0; nx<n; nx++){
  3550.       x1 = (double)nx;
  3551.       y1 = (double)ny;
  3552.       x2 = (x1-xc)*cost - (y1-yc)*sint + xc;
  3553.       y2 = (x1-xc)*sint + (y1-yc)*cost + yc;
  3554.       ix2 = (int)x2;
  3555.       iy2 = (int)y2;
  3556.       tx = x2 - (double)ix2;
  3557.       ty = y2 - (double)iy2;
  3558.       
  3559.       if (ix2+1<n && ix2>=0 && iy2+1<n && iy2>=0){
  3560.         sample1 = (1.0-tx)*s[ix2 + iy2*n] + tx*s[ix2+1 + iy2*n];
  3561.         sample2 = (1.0-tx)*s[ix2 + (iy2+1)*n] + tx*s[ix2+1 + (iy2+1)*n];
  3562.         d[nx + ny*n] = (1-ty)*sample1 + ty*sample2;
  3563.       }
  3564.       else{
  3565.         d[nx + ny*n] = 0.0;
  3566.       }
  3567.     }
  3568.   }
  3569. }
  3570. /*
  3571.  * Swap the mag and phase of the two
  3572.  */
  3573. void swapMagPhase(Complex *c1, Complex *c2, int size){
  3574.   int i;
  3575.   Complex c_temp;
  3576.   for (i=0; i<size; i++){
  3577.     c_temp.r = c1[i].r;
  3578.     c_temp.i = c1[i].i;
  3579.     if (c1[i].getMag() > DBL_TOL && c2[i].getMag() > DBL_TOL){
  3580.       c1[i].r = c1[i].getMag() * c2[i].r/c2[i].getMag();
  3581.       c1[i].i = c1[i].getMag() * c2[i].i/c2[i].getMag();
  3582.     }
  3583.     if ( c2[i].getMag() > DBL_TOL && c_temp.getMag() > DBL_TOL){
  3584.       c2[i].r = c2[i].getMag() * c_temp.r/c_temp.getMag();
  3585.       c2[i].i = c2[i].getMag() * c_temp.i/c_temp.getMag();
  3586.     }
  3587.   }
  3588. }
  3589. /*
  3590.  * Purpose: Go through the image, and find the offending frequency and remove it
  3591.  */
  3592. void removeInterference(Complex *D_Data, int size){
  3593.   int x, y, i;
  3594.   double ave, val;
  3595.   ave = 0.0;
  3596.   for (i=0; i<(size*size); i++){
  3597.     ave += D_Data[i].getMag()/(double)(size*size);
  3598.   }
  3599.   printf("ave: %fn",ave);
  3600.   for (x=1; x<size-1; x++){
  3601.     for (y=1; y<size-1; y++){
  3602.       val = D_Data[x+y*size].getMag();
  3603.       val -= ave*ave*0.1;
  3604.       if (val >  D_Data[x+1+y*size].getMag() &&
  3605.           val >  D_Data[x-1+y*size].getMag() &&
  3606.           val >  D_Data[x+(y+1)*size].getMag() &&
  3607.           val >  D_Data[x+(y-1)*size].getMag() ){
  3608.         D_Data[x+y*size].r = 0.0;
  3609.         D_Data[x+y*size].i = 0.0;
  3610.       }
  3611.     }
  3612.   }
  3613. }
  3614. void createSinusoid(int size, double w, double *s, double mag, int flag){
  3615.   int i;
  3616.   double N;
  3617.   int half = size/2;
  3618.   N = (double) size;
  3619.   if (flag == SINE){
  3620.     for (i=0; i<size; i++){
  3621.       s[i] = mag*sin( (w * (double)(i-half))/N);
  3622.     }
  3623.   }
  3624.   else if (flag == COSINE){
  3625.     for (i=0; i<size; i++){
  3626.       s[i] = mag*cos( (w * (double)(i-half))/N);
  3627.     }
  3628.   }
  3629. }
  3630. void createSobel (int s, double *sobel){
  3631.   int i;
  3632.   int h = s/2;
  3633.   for (i=0; i<s*s; i++){
  3634.     sobel[i] = 0.0;
  3635.   }
  3636.   sobel[(h-1)+(h+1)*s]=-2; sobel[(h)+(h+1)*s]=-2; sobel[(h+1)+(h+1)*s]= 0;
  3637.   sobel[(h-1)+(h)*s]=  -2; sobel[(h)+(h)*s]=   0; sobel[(h+1)+(h)*s]=   2;
  3638.   sobel[(h-1)+(h-1)*s]= 0; sobel[(h)+(h-1)*s]= 2; sobel[(h+1)+(h-1)*s]= 2;
  3639. }
  3640. void createSinusoid2D(int size, double freq, double *s, double mag, int axis, int flag){
  3641.   int i, j;
  3642.   double N;
  3643.   double val = 0.0;
  3644.   N = (double) size;
  3645.   if (axis == X_AXIS){
  3646.     for (i=0; i<size; i++){
  3647.       if (SINE == flag) val = mag*sin( (freq * (double)i)/N);
  3648.       else if (COSINE == flag) val = mag*cos( (freq * (double)i)/N);
  3649.       for (j=0; j<size; j++){
  3650.         s[i + j*size] = val;
  3651.       }
  3652.     }
  3653.   }
  3654.   else if (axis == Y_AXIS){
  3655.     for (j=0; j<size; j++){
  3656.       if(SINE == flag) val = mag*sin( (freq * (double)j)/N);
  3657.       else if(COSINE == flag) val = mag*cos( (freq * (double)j)/N);
  3658.       for (i=0; i<size; i++){
  3659.         s[i + j*size] = val;
  3660.       }
  3661.     }
  3662.   }
  3663. }
  3664. /*
  3665.  * Purpose: Implements a butterworth filter.
  3666.  */
  3667. void butterFilter(Complex *d, int size, double cut_off, double n){
  3668.   int i;
  3669.   Complex c;
  3670.   c.i = 0;
  3671.   c.r = 0;
  3672.   cut_off *= cut_off;
  3673.   for (i=0; i<size/2; i++){
  3674.     c.r = 1.0/(1.0 + pow( ((double)i*i)/ cut_off, n) );
  3675.     d[i] = mult(d[i], c);
  3676.     d[size - i] = mult(d[size-i], c);
  3677.   }
  3678. }
  3679. void plot1D(double *a, int size, char *plot_name, char *filename){
  3680.   FILE   *fpt;
  3681.   char plotfname[255];
  3682.   char buff[255];
  3683.   char fname[255];
  3684.   char ticks[255];
  3685.   sprintf(ticks, "set xtics %d", size/20);
  3686.   printf("saving plots %sn", filename);
  3687.   sprintf(plotfname, "plot_%s.dat", filename);
  3688.   sprintf(fname, "plot_real_%s.plot", filename);
  3689.   
  3690.   fpt = fopen( fname, "wb" );
  3691.   fprintf(fpt, "set title "%s"n", plot_name);
  3692.   fprintf(fpt, "set xrange [0:%d]n", size);
  3693.   fprintf(fpt, "%sn", ticks);
  3694.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3695.   fclose(fpt);
  3696.   fpt = fopen( plotfname, "wb" );
  3697.   for (int i=0; i<size; i++){
  3698.     fprintf(fpt, "%fn", a[i]);
  3699.   }
  3700.   fclose(fpt);
  3701.   sprintf(buff, "gnuplot -persist %sn", fname);
  3702.   system(buff);
  3703. }
  3704. void plotComplex1D(Complex *a, int size, char *plot_name, char *filename){
  3705.   FILE   *fpt, *fpt_data;
  3706.   char plotfname[255];
  3707.   char plot_title[255];
  3708.   char fname[255];
  3709.   char buff[255];
  3710.   char ticks[255];
  3711.   sprintf(ticks, "set xtics %d", size/20);
  3712.   
  3713.   printf("saving complex plots %sn", filename);
  3714.   /*
  3715.   //Real plot
  3716.   sprintf(plot_title, "%s (Real part)", plot_name);
  3717.   sprintf(plotfname, "plot_real_%s.dat", filename);
  3718.   sprintf(fname, "plot_real_%s.plot", filename);
  3719.   fpt = fopen( fname, "wb" );
  3720.   fprintf(fpt, "set title "%s"n", plot_title);
  3721.   fprintf(fpt, "set xrange [0:%d]n", size);
  3722.   fprintf(fpt, "%sn", ticks);
  3723.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3724.   fclose(fpt);
  3725.   fpt_data = fopen( plotfname, "wb" );
  3726.   for (int i=0; i<size; i++){
  3727.     fprintf(fpt_data, "%fn", a[i].r);
  3728.   }
  3729.   fclose(fpt_data);
  3730.   sprintf(buff, "gnuplot -persist %sn", fname);
  3731.   system(buff);
  3732.   //Imaginary plot
  3733.   sprintf(plot_title, "%s (Imaginary part)", plot_name);
  3734.   sprintf(plotfname, "plot_img_%s.dat", filename);
  3735.   sprintf(fname, "plot_img_%s.plot", filename);
  3736.   fpt = fopen( fname, "wb" );
  3737.   fprintf(fpt, "set title "%s"n", plot_title);
  3738.   fprintf(fpt, "set xrange [0:%d]n", size);
  3739.   fprintf(fpt, "%sn", ticks);
  3740.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3741.   fclose(fpt);
  3742.   fpt_data = fopen( plotfname, "wb" );
  3743.   for (int i=0; i<size; i++){
  3744.     fprintf(fpt_data, "%fn", a[i].i);
  3745.   }
  3746.   fclose(fpt_data);
  3747.   sprintf(buff, "gnuplot -persist %sn", fname);
  3748.   system(buff);
  3749.   //Magnitude plot
  3750.   sprintf(plot_title, "%s (Magnitude)", plot_name);
  3751.   sprintf(plotfname, "plot_mag_%s.dat", filename);
  3752.   sprintf(fname, "plot_mag_%s.plot", filename);
  3753.   fpt = fopen( fname, "wb" );
  3754.   fprintf(fpt, "set title "%s"n", plot_title);
  3755.   fprintf(fpt, "set xrange [0:%d]n", size);
  3756.   fprintf(fpt, "%sn", ticks);
  3757.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3758.   fclose(fpt);
  3759.   fpt_data = fopen( plotfname, "wb" );
  3760.   for (int i=0; i<size; i++){
  3761.     fprintf(fpt_data, "%fn", a[i].getMag());
  3762.   }
  3763.   fclose(fpt_data);
  3764.   sprintf(buff, "gnuplot -persist %sn", fname);
  3765.   system(buff);
  3766.   //Phase plot
  3767.   sprintf(plot_title, "%s (Phase)", plot_name);
  3768.   sprintf(plotfname, "plot_phase_%s.dat", filename);
  3769.   sprintf(fname, "plot_phase_%s.plot", filename);
  3770.   fpt = fopen( fname, "wb" );
  3771.   fprintf(fpt, "set title "%s"n", plot_title);
  3772.   fprintf(fpt, "set xrange [0:%d]n", size);
  3773.   fprintf(fpt, "%sn", ticks);
  3774.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3775.   fclose(fpt);
  3776.   fpt_data = fopen( plotfname, "wb" );
  3777.   for (int i=0; i<size; i++){
  3778.     fprintf(fpt_data, "%fn", a[i].getPhase());
  3779.   }
  3780.   fclose(fpt_data);
  3781.   sprintf(buff, "gnuplot -persist %sn", fname);
  3782.   system(buff);
  3783.   */
  3784.   //Power plot
  3785.   sprintf(plot_title, "%s (Power Spectrum)", plot_name);
  3786.   sprintf(plotfname, "plot_power_%s.dat", filename);
  3787.   sprintf(fname, "plot_power_%s.plot", filename);
  3788.   fpt = fopen( fname, "wb" );
  3789.   fprintf(fpt, "set title "%s"n", plot_title);
  3790.   fprintf(fpt, "set xrange [0:%d]n", size);
  3791.   //fprintf(fpt, "set xtics 16n");
  3792.   fprintf(fpt, "%sn", ticks);
  3793.   fprintf(fpt, "plot "%s" with linesn", plotfname);
  3794.   fclose(fpt);
  3795.   fpt_data = fopen( plotfname, "wb" );
  3796.   for (int i=0; i<size; i++){
  3797.     fprintf(fpt_data, "%fn", a[i].getMag()*a[i].getMag() );
  3798.   }
  3799.   fclose(fpt_data);
  3800.   sprintf(buff, "gnuplot -persist %sn", fname);
  3801.   system(buff);
  3802. }
  3803. void writeUB(double *s, int size, char *filename){
  3804.   int i;
  3805.   unsigned char uc;
  3806.   FILE *fptr;
  3807.   printf("saving ub file %sn", filename);
  3808.   fptr = fopen(filename, "wb");
  3809.   for (i=0; i<size; i++){
  3810.     uc = (unsigned char)(s[i] + 128);
  3811.     fputc(uc, fptr);
  3812.   }
  3813.   fclose(fptr);
  3814.   /*
  3815.   if ( fwrite( uc, 1, (size_t)1, fptr ) !=
  3816.        (unsigned int)1 ){
  3817.     puts("fail to write");
  3818.     return FALSE;
  3819.   }
  3820.   */
  3821. }
  3822. /*
  3823.   else if (5 == type){
  3824.     fprintf( fpt, "P5n" );
  3825.     //if ( HeaderInfo!=NULL ) fprintf( F, "# %sn", HeaderInfo );
  3826.     //else fprintf( F, "# Grayscale Imagen" );
  3827.     fprintf( fpt, "%d %dn%dn", width, height, max );
  3828.     for ( i=0; i<height; i++ ) {
  3829.       save_line = &data[i*width];
  3830.       if ( fwrite( (void*)save_line, 1, (size_t)width, fpt ) !=
  3831.            (unsigned int)width ){
  3832.         puts("fail to write");
  3833.         return FALSE;
  3834.       }
  3835.     }
  3836.     fclose( fpt );
  3837.     return TRUE;
  3838.   }
  3839.   */
  3840. void filterFourier2D(Complex *c, int size, int begin, int end){
  3841.   int i,j;
  3842.   for (i=begin; i<end; i++)
  3843.     for (j=begin; j<end; j++){
  3844.       c[i + j*size].r = 0.0; 
  3845.       c[i + j*size].i = 0.0; 
  3846.     }
  3847. }
  3848. //MAIN
  3849. //MAIN
  3850. //MAIN
  3851. //MAIN
  3852. //MAIN
  3853. //MAIN
  3854. /******************************************************************************
  3855.  * Purpose:  main... The head hancho function.
  3856.  *****************************************************************************/
  3857. int main(int argc, char** argv) {
  3858.   Image image;
  3859.   double *img_dbl = NULL;
  3860.   double *img_dbl2 = NULL;
  3861.   if (argc < 3 || argc >4) {
  3862.     printf("  Error:nIncorrect command line arguments.  The format is:n");
  3863.     printf(">dig [-flag] [image.ppm or image.pgm] n");
  3864.     printf("Example:ndig -p1 parrots.ppmnn");
  3865.     printf("flags are as follows:n");
  3866.     printf("-p2    :Get the difference of two images.n");
  3867.     printf("-p3    :Segment out objects in the image using thresholding.n");
  3868.     printf("-cont  :Find the boundaries of the shapes using contour following.n");
  3869.     printf("-cut   :Use Cutter on the image (Intelligent Scissors).n");
  3870.     exit (1);
  3871.   }
  3872.   if (strcmp( "-p3", argv[1]) == 0 ){
  3873.     Image orig_img, im2, info_image;
  3874.     int num_shapes;
  3875.     
  3876.     if ( !image.read(argv[2]) ){
  3877.       printf("Error in reading the imagen");
  3878.       exit (1);
  3879.     }
  3880.     image.makeGray();
  3881.     orig_img.copy(&image);
  3882.     im2.copy(&image);
  3883.     num_shapes = boxThreshold(&image, &im2);
  3884.     info_image.copy(&im2);
  3885.     info_image.createInfoImage(num_shapes);
  3886.     info_image.write("out_info.ppm");
  3887.     orig_img.mask(&im2);
  3888.     orig_img.write("out_segmented.pgm");
  3889.   }
  3890.   else if (strcmp( "-p2", argv[1]) == 0 ){
  3891.     Image image2;
  3892.     if ( !image.read(argv[2]) ){
  3893.       printf("Error in reading the imagen");
  3894.       exit (1);
  3895.     }
  3896.     if ( !image2.read(argv[3]) ){
  3897.       printf("Error in reading the imagen");
  3898.       exit (1);
  3899.     }
  3900.     image.makeGray();
  3901.     image2.makeGray();
  3902.     image.diff(&image2);
  3903.     image.write("out_diff.pgm");
  3904.   }
  3905.   else if (strcmp( "-cont", argv[1]) == 0 ){
  3906.     ContourMaker *cont = NULL;
  3907.     if ( !image.read(argv[2]) ){
  3908.       printf("Error in reading the imagen");
  3909.       exit (1);
  3910.     }
  3911.     cont = new ContourMaker(&image);
  3912.     /*
  3913.     //TEMP
  3914.     int size, i;
  3915.     PointList pl2;
  3916.     FILE *fp;
  3917.     fp = fopen(argv[3], "r");
  3918.     if (fp == NULL) { printf("Could not open file %sn", argv[3]); exit (1); }
  3919.     fscanf(fp,"%d", &size);
  3920.     pl2.num_points = size;
  3921.     for (i=0; i<size; i++){
  3922.       fscanf(fp,"%d", &pl2.x[i]);
  3923.       fscanf(fp,"%d", &pl2.y[i]);
  3924.     }
  3925.     fclose(fp);
  3926.     image.clear(0);
  3927.     cont->drawBoundaries(&image, &pl2);
  3928.     image.write("out_myst.ppm");
  3929.     //END TEMP
  3930.     */
  3931.     delete cont;
  3932.   }
  3933.   else if (strcmp( "-match", argv[1]) == 0 ){
  3934.     FILE *fp; 
  3935.     int i, size;
  3936.     double dist;
  3937.     PointList pl1, pl2;
  3938.     fp = fopen(argv[2], "r");
  3939.     if (fp == NULL) { printf("Could not open file %sn", argv[2]); exit (1); }
  3940.     fscanf(fp,"%d", &size);
  3941.     pl1.num_points = size;
  3942.     for (i=0; i<size; i++){
  3943.       fscanf(fp,"%d", &pl1.x[i]);
  3944.       fscanf(fp,"%d", &pl1.y[i]);
  3945.     }
  3946.     fclose(fp);
  3947.     fp = fopen(argv[3], "r");
  3948.     if (fp == NULL) { printf("Could not open file %sn", argv[3]); exit (1); }
  3949.     fscanf(fp,"%d", &size);
  3950.     pl2.num_points = size;
  3951.     for (i=0; i<size; i++){
  3952.       fscanf(fp,"%d", &pl2.x[i]);
  3953.       fscanf(fp,"%d", &pl2.y[i]);
  3954.     }
  3955.     fclose(fp);
  3956.     pl1.createShapeRep();
  3957.     pl2.createShapeRep();
  3958.     dist = pl1.getDist(&pl2);
  3959.     pl1.writeShapeRep("out_boundary");
  3960.     pl2.writeShapeRep("out_bound_shape");
  3961.     printf("dist between pl1 & pl2: %fn", dist);
  3962.     /*
  3963.     for (i=0; i<size; i++){
  3964.       printf("%d ", pl1.x[i]);
  3965.       printf("%dn", pl1.y[i]);
  3966.     }
  3967.     */
  3968.     //size1 = 0; EOF != fscanf(fp,"%lf", &signal1[size1]); size1++)
  3969.     //    ;
  3970.   }
  3971.   else if (strcmp( "-cut", argv[1]) == 0 ){
  3972.     char name[255];
  3973.     cut_image = new Image();
  3974.     if ( !cut_image->read(argv[2]) ){
  3975.       printf("Error in reading image %sn", argv[2]);
  3976.       return 1;
  3977.     }
  3978.     image.copy(cut_image);
  3979.     cut_xres = cut_image->width;
  3980.     cut_yres = cut_image->height;
  3981.     //FOR GL RENDERING
  3982.     cut_image->flip();
  3983.     cut_image->makeColor();
  3984.     cutter = new Cutter(&image);
  3985. //    cutter->dumpCostGraph();
  3986.     glutInit(&argc, argv);
  3987.       //GLUT_DOUBLE used for double buffering
  3988.     glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE );
  3989.     glutInitWindowSize (cut_xres, cut_yres);
  3990.     glutInitWindowPosition (0, 0);
  3991.     strcpy(name, argv[0]);//name of the program
  3992.     strcat(name, " - Intelligent Scissors");
  3993.     glutCreateWindow (name);
  3994.     init();
  3995.       // The names of the functions that you use for displaying and events go
  3996.       // as the parameters of this functions below
  3997.     glutDisplayFunc(display);
  3998.     glutMouseFunc(mouseButton);
  3999.     glutMotionFunc(mouseMotion);
  4000.     glutKeyboardFunc(keyboard);
  4001.     glutMainLoop();
  4002.   }
  4003.   else{
  4004.     printf("Flag '%s' not supportedn", argv[1]);
  4005.   }
  4006.   if (NULL != img_dbl) delete [] img_dbl;
  4007.   if (NULL != img_dbl2) delete [] img_dbl2;
  4008.   return 0;
  4009. }
  4010. /*
  4011.  * ADDITIONAL CODE FOR PROJECTS. CODE NOT WRITTEN BY MIKE SMITH
  4012.  *
  4013.  *
  4014.  *
  4015.  */
  4016. /*==================================================*
  4017.  * Code framework for BYU CS 450 Assignment 4, audio part.
  4018.  *
  4019.  * This source code is covered under the GNU Public License, Version 2
  4020.  * ( see http://www.gnu.org/ ).
  4021. *====================================================*/
  4022. /* Write a long out to a stream in little-endian order */
  4023. void fwrite_long(unsigned long number, FILE *stream) {
  4024.   fputc((number & 0xff), stream);            fputc((number & 0xff00) >> 8, stream);
  4025.   fputc((number & 0xff0000) >> 16, stream);  fputc((number & 0xff000000) >> 24, stream);
  4026. }
  4027. /* Write a short out to a stream in little-endian order */
  4028. void fwrite_short(unsigned short number, FILE *stream) {
  4029.   fputc((number & 0x00ff), stream);
  4030.   fputc((number & 0xff00) >> 8, stream);
  4031. }
  4032. /* Write a WAV sound file (1-channel, 16-bit only).
  4033.  * Works on both big/little endian machines.
  4034.  */
  4035. void write_wav(char *filename, double *sample, unsigned long sample_len, unsigned long freq) {
  4036.   FILE *out = fopen(filename, "wb");
  4037.   double min, max, offset, scale;
  4038.   unsigned int i;
  4039.   /* Check output file opened OK */
  4040.   if (!out) {
  4041.     fprintf(stderr, "nCould not open .wav filenn");
  4042.     exit(1);
  4043.   }
  4044.   /* Write the .wav header */
  4045.   fprintf(out, "RIFF");
  4046.   fwrite_long(2 * sample_len + 36, out);
  4047.   fprintf(out, "WAVEfmt ");
  4048.   fputc(16, out);  fputc(0, out);  fputc(0, out);   fputc(0, out);
  4049.   fputc(1, out);   fputc(0, out);  fputc(1, out);   fputc(0, out);
  4050.   fwrite_long(freq, out);  fwrite_long(freq * 2, out);
  4051.   fputc(2, out);   fputc(0, out);  fputc(16, out);  fputc(0, out);
  4052.   fprintf(out, "data");
  4053.   fwrite_long(sample_len * 2, out);
  4054.   
  4055.   /* Find the range of values in the sample */
  4056.   min = max = sample[0];
  4057.   for (i = 1; i < sample_len; i++) {
  4058.     if (sample[i] > max)  max = sample[i];
  4059.     if (sample[i] < min)  min = sample[i];
  4060.   }
  4061.   if (min == max) {
  4062.     max = 1.0;
  4063.     min = -1.0;    /* Avoid division by zero */
  4064.   }
  4065.   offset = (max + min) / 2.0;
  4066.   scale = ((double) SHRT_MAX * 2.0) / (max - min);
  4067.   /* Write the scaled samples to the file */
  4068.   for (i = 0; i < sample_len; i++) {
  4069.     /* Center the sample about 0, and scale it so its maximum value is +/- SHRT_MAX */
  4070.     short sample_val = (short)((sample[i] - offset) * scale +.5);
  4071.     // Uncomment the next line to output the first 250 sample values for plotting:
  4072.     //if (i < 250) printf("%dn", (int) sample_val);
  4073.     fwrite_short((unsigned short) sample_val, out);  /* Convert to unsigned short for output */
  4074.   }
  4075.   fclose(out);
  4076. }
  4077. /*-----------------------------------------------------------------------*
  4078.  * Create a sine wave sound sample with a fraction of random noise
  4079. *----------------------------------------------------------*/
  4080. double *create_wav(int sample_len, double freq_scale, double noise_frac) {
  4081.   double *sample = (double *) malloc(sample_len * sizeof(double));
  4082.   int i;
  4083.   /* Create the sample */
  4084.   for (i = 0; i < sample_len; i++) {
  4085.     double a = (double) i * freq_scale;
  4086.     double s = sin(a);
  4087.     double r = noise_frac * (2.0 * (double) random() / ((double) RAND_MAX + 1.0) - 1.0);
  4088.     sample[i] = s + r;
  4089.   }
  4090.   return sample;
  4091. }
  4092. //EOF