segmentation
msImageProcessor.cpp
1 /*******************************************************
2 
3  Mean Shift Analysis Library
4  =============================================
5 
6 
7  The mean shift library is a collection of routines
8  that use the mean shift algorithm. Using this algorithm,
9  the necessary output will be generated needed
10  to analyze a given input set of data.
11 
12  Mean Shift Image Processor Class:
13  ================================
14 
15  The following class inherits from the mean shift library
16  in order to perform the specialized tasks of image
17  segmentation and filtering.
18 
19  The definition of the Mean Shift Image Processor Class
20  is provided below. Its prototype is provided in
21  'msImageProcessor.h'.
22 
23 The theory is described in the papers:
24 
25  D. Comaniciu, P. Meer: Mean Shift: A robust approach toward feature
26  space analysis.
27 
28  C. Christoudias, B. Georgescu, P. Meer: Synergism in low level vision.
29 
30 and they are is available at:
31  http://www.caip.rutgers.edu/riul/research/papers/
32 
33 Implemented by Chris M. Christoudias, Bogdan Georgescu
34 ********************************************************/
35 
36 //include image processor class prototype
37 #include "segm/msImageProcessor.h"
38 
39 //include needed libraries
40 #include <math.h>
41 #include <stdio.h>
42 #include <assert.h>
43 #include <string.h>
44 #include <stdlib.h>
45 
46 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
47 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
48 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ PUBLIC METHODS @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
49 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
50 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
51 
52  /*/\/\/\/\/\/\/\/\/\/\/\/\*/
53  /* Constructor/Destructor */
54  /*\/\/\/\/\/\/\/\/\/\/\/\/*/
55 
56 /*******************************************************/
57 /*Class Constructor */
58 /*******************************************************/
59 /*Post: */
60 /* The msImageProcessor class has been properly */
61 /* initialized. */
62 /*******************************************************/
63 
64 msImageProcessor::msImageProcessor( void )
65 {
66 
67  //intialize basin of attraction structure
68  //used by the filtering algorithm
69  modeTable = NULL;
70  pointList = NULL;
71  pointCount = 0;
72 
73  //initialize region list
74  regionList = NULL;
75 
76  //initialize output structures...
77  msRawData = NULL;
78  labels = NULL;
79  modes = NULL;
80  modePointCounts = NULL;
81  regionCount = 0;
82 
83  //intialize temporary buffers used for
84  //performing connected components
85  indexTable = NULL;
86  LUV_data = NULL;
87 
88  //initialize region adjacency matrix
89  raList = NULL;
90  freeRAList = NULL;
91  raPool = NULL;
92 
93  //intialize visit table to having NULL entries
94  visitTable = NULL;
95 
96  //initialize epsilon such that transitive closure
97  //does not take edge strength into consideration when
98  //fusing regions of similar color
99  epsilon = 1.0;
100 
101  //initialize class state to indicate that
102  //an output data structure has not yet been
103  //created...
104  class_state.OUTPUT_DEFINED = false;
105 
106 //Changed by Sushil from 1.0 to 0.1, 11/11/2008
107  LUV_treshold = 0.1F;
108 }
109 
110 /*******************************************************/
111 /*Class Destructor */
112 /*******************************************************/
113 /*Post: */
114 /* The msImageProcessor class has been properly */
115 /* destroyed. */
116 /*******************************************************/
117 
118 msImageProcessor::~msImageProcessor( void )
119 {
120 
121  //de-allocate memory
122  if(class_state.OUTPUT_DEFINED) DestroyOutput();
123  if(regionList) delete regionList;
124  regionList = NULL;
125 
126  //done.
127 
128 }
129 
130  /*/\/\/\/\/\/\/\/\/\/\/\/\/\*/
131  /* Input Image Declaration */
132  /*\/\/\/\/\/\/\/\/\/\/\/\/\/*/
133 
134 /*******************************************************/
135 /*Define Image */
136 /*******************************************************/
137 /*Uploads an image into the image segmenter class to */
138 /*be segmented. */
139 /*******************************************************/
140 /*Pre: */
141 /* - data_ is a one dimensional array of unsigned */
142 /* char RGB vectors */
143 /* - type is the type of the image: COLOR or */
144 /* GREYSCALE */
145 /* - height_ and width_ define the dimension of */
146 /* the image */
147 /* - if the image is of type GREYSCALE then */
148 /* data containes only one number per pixel */
149 /* location, where a pixel location is defined */
150 /* by the index into the data array */
151 /*Post: */
152 /* - the image specified has been uploaded into */
153 /* the image segmenter class to be segmented. */
154 /*******************************************************/
155 
156 void msImageProcessor::DefineImage(byte *data_, imageType type, int height_, int width_)
157 {
158 
159  //obtain image dimension from image type
160  int dim;
161  if(type == COLOR)
162  dim = 3;
163  else
164  dim = 1;
165 
166  //perfor rgb to luv conversion
167  int i;
168  float *luv = new float [height_*width_*dim];
169  if(dim == 1)
170  {
171  for(i = 0; i < height_*width_; i++)
172  luv[i] = (float)(data_[i]);
173  }
174  else
175  {
176  for(i = 0; i < height_*width_; i++)
177  {
178  RGBtoLUV(&data_[dim*i], &luv[dim*i]);
179  }
180  }
181 
182  //define input defined on a lattice using mean shift base class
183  DefineLInput(luv, height_, width_, dim);
184 
185  //Define a default kernel if it has not been already
186  //defined by user
187  if(!h)
188  {
189  //define default kernel paramerters...
190  kernelType k[2] = {Uniform, Uniform};
191  int P[2] = {2, N};
192  float tempH[2] = {1.0 , 1.0};
193 
194  //define default kernel in mean shift base class
195  DefineKernel(k, tempH, P, 2);
196  }
197 
198  //de-allocate memory
199  delete [] luv;
200 
201  //done.
202  return;
203 
204 }
205 
206 void msImageProcessor::DefineBgImage(byte* data_, imageType type, int height_, int width_)
207 {
208 
209  //obtain image dimension from image type
210  int dim;
211  if(type == COLOR)
212  dim = 3;
213  else
214  dim = 1;
215 
216  //perform texton classification
217  int i;
218 
219  float *luv = new float [height_*width_*dim];
220 
221  if(dim == 1)
222  {
223  for(i = 0; i < height_*width_; i++)
224  luv[i] = (float)(data_[i]);
225  }
226  else
227  {
228  for(i = 0; i < height_*width_; i++)
229  RGBtoLUV(&data_[dim*i], &luv[dim*i]);
230 
231  }
232 
233  //define input defined on a lattice using mean shift base class
234  DefineLInput(luv, height_, width_, dim);
235 
236 
237  //Define a default kernel if it has not been already
238  //defined by user
239  if(!h)
240  {
241  //define default kernel paramerters...
242  kernelType k[2] = {Uniform, Uniform};
243  int P[2] = {2, N};
244  float tempH[2] = {1.0 , 1.0};
245 
246  //define default kernel in mean shift base class
247  DefineKernel(k, tempH, P, 2);
248  }
249 
250  //de-allocate memory
251  delete [] luv;
252 
253  //done.
254  return;
255 
256 }
257 
258  /*/\/\/\/\/\/\/\/\*/
259  /* Weight Map */
260  /*\/\/\/\/\/\/\/\/*/
261 
262 /*******************************************************/
263 /*Set Weight Map */
264 /*******************************************************/
265 /*Populates the weight map with specified edge */
266 /*strengths. */
267 /*******************************************************/
268 /*Pre: */
269 /* - wm is a floating point array of size */
270 /* (height x width) specifying for each pixel */
271 /* edge strength. */
272 /* - eps is a threshold used to fuse similar */
273 /* regions during transitive closure. */
274 /*Post: */
275 /* - wm has been used to populate the weight */
276 /* map. */
277 /* - the threshold used during transitive closure */
278 /* is taken as eps. */
279 /*******************************************************/
280 
281 void msImageProcessor::SetWeightMap(float *wm, float eps)
282 {
283 
284  //initlaize confmap using wm
285  SetLatticeWeightMap(wm);
286 
287  //set threshold value
288  if((epsilon = eps) < 0)
289  ErrorHandler("msImageProcessor", "SetWeightMap", "Threshold is negative.");
290 
291  //done.
292  return;
293 
294 }
295 
296 /*******************************************************/
297 /*Remove Weight Map */
298 /*******************************************************/
299 /*Removes the weight map. */
300 /*******************************************************/
301 /*Post: */
302 /* - the weight map has been removed. */
303 /* - if a weight map did not exist NO error */
304 /* is flagged. */
305 /*******************************************************/
306 
307 void msImageProcessor::RemoveWeightMap( void )
308 {
309 
310  //remove confmap
311  RemoveLatticeWeightMap();
312 
313  //set threshold value to zero
314  epsilon = 0;
315 
316  //done.
317  return;
318 
319 }
320 
321  /*/\/\/\/\/\/\/\/\/\*/
322  /* Image Filtering */
323  /*\/\/\/\/\/\/\/\/\/*/
324 
325 /*******************************************************/
326 /*Filter */
327 /*******************************************************/
328 /*Performs mean shift filtering on the specified input */
329 /*image using a user defined kernel. */
330 /*******************************************************/
331 /*Pre: */
332 /* - the user defined kernel used to apply mean */
333 /* shift filtering to the defined input image */
334 /* has spatial bandwidth sigmaS and range band- */
335 /* width sigmaR */
336 /* - speedUpLevel determines whether or not the */
337 /* filtering should be optimized for faster */
338 /* execution: a value of NO_SPEEDUP turns this */
339 /* optimization off and a value SPEEDUP turns */
340 /* this optimization on */
341 /* - a data set has been defined */
342 /* - the height and width of the lattice has been */
343 /* specified using method DefineLattice() */
344 /*Post: */
345 /* - mean shift filtering has been applied to the */
346 /* input image using a user defined kernel */
347 /* - the filtered image is stored in the private */
348 /* data members of the msImageProcessor class. */
349 /*******************************************************/
350 
351 void msImageProcessor::Filter(int sigmaS, float sigmaR, SpeedUpLevel speedUpLevel)
352 {
353 
354  //Check Class consistency...
355 
356  //check:
357  // (1) if this operation is consistent
358  // (2) if kernel was created
359  // (3) if data set is defined
360  // (4) if the dimension of the kernel agrees with that
361  // of the defined data set
362  // if not ... flag an error!
363  classConsistencyCheck(N+2, true);
364  if(ErrorStatus == EL_ERROR)
365  return;
366 
367  //If the algorithm has been halted, then exit
368  if((ErrorStatus = msSys.Progress((float)(0.0))) == EL_HALT)
369  {
370  return;
371  }
372 
373  //If the image has just been read then allocate memory
374  //for and initialize output data structure used to store
375  //image modes and their corresponding regions...
376  if(class_state.OUTPUT_DEFINED == false)
377  {
378  InitializeOutput();
379 
380  //check for errors...
381  if(ErrorStatus == EL_ERROR)
382  return;
383  }
384 
385  //****************** Allocate Memory ******************
386 
387  //Allocate memory for basin of attraction mode structure...
388  if((!(modeTable = new unsigned char [L]))||(!(pointList = new int [L])))
389  {
390  ErrorHandler("msImageProcessor", "Allocate", "Not enough memory.");
391  return;
392  }
393 
394  //start timer
395 #ifdef PROMPT
396  double timer;
397  msSys.StartTimer();
398 #endif
399 
400  //*****************************************************
401 
402  //filter image according to speedup level...
403  switch(speedUpLevel)
404  {
405  //no speedup...
406  case NO_SPEEDUP:
407  //NonOptimizedFilter((float)(sigmaS), sigmaR); break;
408  NewNonOptimizedFilter((float)(sigmaS), sigmaR);
409  break;
410  //medium speedup
411  case MED_SPEEDUP:
412  //OptimizedFilter1((float)(sigmaS), sigmaR); break;
413  NewOptimizedFilter1((float)(sigmaS), sigmaR);
414  break;
415  //high speedup
416  case HIGH_SPEEDUP:
417  //OptimizedFilter2((float)(sigmaS), sigmaR); break;
418  NewOptimizedFilter2((float)(sigmaS), sigmaR);
419  break;
420  // new speedup
421  }
422 
423  //****************** Deallocate Memory ******************
424 
425  //de-allocate memory used by basin of attraction mode structure
426  delete [] modeTable;
427  delete [] pointList;
428 
429  //re-initialize structure
430  modeTable = NULL;
431  pointList = NULL;
432  pointCount = 0;
433 
434  //*******************************************************
435 
436  //If the algorithm has been halted, then de-allocate the output
437  //and exit
438  if((ErrorStatus = msSys.Progress((float)(0.8))) == EL_HALT)
439  {
440  DestroyOutput();
441  return;
442  }
443 
444  //Label image regions, also if segmentation is not to be
445  //performed use the resulting classification structure to
446  //calculate the image boundaries...
447 
448  /*
449  //copy msRawData into LUV_data, rounding each component of each
450  //LUV value stored by msRawData to the nearest integer
451  int i;
452  for(i = 0; i < L*N; i++)
453  {
454  if(msRawData[i] < 0)
455  LUV_data[i] = (int)(msRawData[i] - 0.5);
456  else
457  LUV_data[i] = (int)(msRawData[i] + 0.5);
458  }
459  */
460  int i;
461  for (i=0; i<L*N; i++)
462  {
463  LUV_data[i] = msRawData[i];
464  }
465 
466 
467 #ifdef PROMPT
468  timer = msSys.ElapsedTime();
469  msSys.Prompt("(%6.2f sec)\nConnecting regions ...", timer);
470  msSys.StartTimer();
471 #endif
472 
473  //Perform connecting (label image regions) using LUV_data
474  Connect();
475 
476 #ifdef PROMPT
477  timer = msSys.ElapsedTime();
478  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount);
479  msSys.StartTimer();
480 #endif
481 
482  //done.
483  return;
484 
485 }
486 
487  /*/\/\/\/\/\/\/\/\/\/\/\*/
488  /* Image Region Fusing */
489  /*\/\/\/\/\/\/\/\/\/\/\/*/
490 
491 /*******************************************************/
492 /*Fuse Regions */
493 /*******************************************************/
494 /*Fuses the regions of a filtered image. */
495 /*******************************************************/
496 /*Pre: */
497 /* - the range radius is specified by sigmaR */
498 /* - minRegion is the minimum point density that */
499 /* a region may have in the resulting segment- */
500 /* ed image */
501 /* - a data set has been defined */
502 /* - the height and width of the lattice has been */
503 /* specified using method DefineLattice() */
504 /*Post: */
505 /* - the image regions have been fused. */
506 /* - if an result is stored by this class then */
507 /* this result is used as input to this method. */
508 /* - if no result is stored by this class, */
509 /* the input image defined by calling the */
510 /* method DefineImage is used. */
511 /*******************************************************/
512 
513 void msImageProcessor::FuseRegions(float sigmaS, int minRegion)
514 {
515 
516  //Check Class consistency...
517 
518  //check:
519  // (1) if this operation is consistent
520  // (2) if kernel was created
521  // (3) if data set is defined
522  // (4) if the dimension of the kernel agrees with that
523  // of the defined data set
524  // if not ... flag an error!
525  classConsistencyCheck(N+2, true);
526  if(ErrorStatus == EL_ERROR)
527  return;
528 
529  //Check to see if the algorithm is to be halted, if so then
530  //destroy output and exit
531  if((ErrorStatus = msSys.Progress((float)(0.8))) == EL_HALT)
532  {
533  if(class_state.OUTPUT_DEFINED) DestroyOutput();
534  return;
535  }
536 
537  //obtain sigmaS (make sure it is not zero or negative, if not
538  //flag an error)
539  if((h[1] = sigmaS) <= 0)
540  {
541  ErrorHandler("msImageProcessor", "FuseRegions", "The feature radius must be greater than or equal to zero.");
542  return;
543  }
544 
545  //if output has not yet been generated then classify the input
546  //image regions to be fused...
547  if(!(class_state.OUTPUT_DEFINED))
548  {
549 
550  //Initialize output data structure used to store
551  //image modes and their corresponding regions...
552  InitializeOutput();
553 
554  //check for errors...
555  if(ErrorStatus == EL_ERROR)
556  return;
557 
558  //copy data into LUV_data used to classify
559  //image regions
560  /*
561  int i;
562  for(i = 0; i < L*N; i++)
563  {
564  if(data[i] < 0)
565  LUV_data[i] = (int)(data[i] - 0.5);
566  else
567  LUV_data[i] = (int)(data[i] + 0.5);
568  }
569  */
570  int i;
571  for (i=0; i<L*N; i++)
572  {
573  LUV_data[i] = data[i];
574  }
575 
576 #ifdef PROMPT
577  msSys.Prompt("Connecting regions ...");
578  msSys.StartTimer();
579 #endif
580 
581  //Perform connecting (label image regions) using LUV_data
582  Connect();
583 
584  //check for errors
585  if(ErrorStatus == EL_ERROR)
586  return;
587 
588 #ifdef PROMPT
589  double timer = msSys.ElapsedTime();
590  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount);
591 #endif
592 
593  }
594 
595  //Check to see if the algorithm is to be halted, if so then
596  //destroy output and exit
597  if((ErrorStatus = msSys.Progress((float)(0.85))) == EL_HALT)
598  {
599  DestroyOutput();
600  return;
601  }
602 
603 #ifdef PROMPT
604  msSys.Prompt("Applying transitive closure...");
605  msSys.StartTimer();
606 #endif
607 
608  //allocate memory visit table
609  visitTable = new unsigned char [L];
610 
611  //Apply transitive closure iteratively to the regions classified
612  //by the RAM updating labels and modes until the color of each neighboring
613  //region is within sqrt(rR2) of one another.
614  rR2 = (float)(h[1]*h[1]*0.25);
615  TransitiveClosure();
616  int oldRC = regionCount;
617  int deltaRC, counter = 0;
618  do {
619  TransitiveClosure();
620  deltaRC = oldRC-regionCount;
621  oldRC = regionCount;
622  counter++;
623  } while ((deltaRC <= 0)&&(counter < 10));
624 
625  //de-allocate memory for visit table
626  delete [] visitTable;
627  visitTable = NULL;
628 
629  //Check to see if the algorithm is to be halted, if so then
630  //destroy output and region adjacency matrix and exit
631  if((ErrorStatus = msSys.Progress((float)(1.0))) == EL_HALT)
632  {
633  DestroyRAM();
634  DestroyOutput();
635  return;
636  }
637 
638 #ifdef PROMPT
639  double timer = msSys.ElapsedTime();
640  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d)\nPruning spurious regions ...", timer, regionCount);
641  msSys.StartTimer();
642 #endif
643 
644  //Prune spurious regions (regions whose area is under
645  //minRegion) using RAM
646  Prune(minRegion);
647 
648 #ifdef PROMPT
649  timer = msSys.ElapsedTime();
650  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount);
651  msSys.StartTimer();
652 #endif
653 
654  //Check to see if the algorithm is to be halted, if so then
655  //destroy output and region adjacency matrix and exit
656  if((ErrorStatus = msSys.Progress((float)(1.0))) == EL_HALT)
657  {
658  DestroyRAM();
659  DestroyOutput();
660  return;
661  }
662 
663  //de-allocate memory for region adjacency matrix
664  DestroyRAM();
665 
666  //output to msRawData
667  int i, j, label;
668  for(i = 0; i < L; i++)
669  {
670  label = labels[i];
671  for(j = 0; j < N; j++)
672  {
673  msRawData[N*i+j] = modes[N*label+j];
674  }
675  }
676 
677  //done.
678  return;
679 
680 }
681 
682  /*/\/\/\/\/\/\/\/\/\/\*/
683  /* Image Segmentation */
684  /*\/\/\/\/\/\/\/\/\/\/*/
685 
686 /*******************************************************/
687 /*Segment */
688 /*******************************************************/
689 /*Segments the defined image. */
690 /*******************************************************/
691 /*Pre: */
692 /* - sigmaS and sigmaR are the spatial and range */
693 /* radii of the search window respectively */
694 /* - minRegion is the minimum point density that */
695 /* a region may have in the resulting segment- */
696 /* ed image */
697 /* - speedUpLevel determines whether or not the */
698 /* filtering should be optimized for faster */
699 /* execution: a value of NO_SPEEDUP turns this */
700 /* optimization off and a value SPEEDUP turns */
701 /* this optimization on */
702 /*Post: */
703 /* - the defined image is segmented and the */
704 /* resulting segmented image is stored in the */
705 /* private data members of the image segmenter */
706 /* class. */
707 /* - any regions whose point densities are less */
708 /* than or equal to minRegion have been pruned */
709 /* from the segmented image. */
710 /*******************************************************/
711 
712 void msImageProcessor::Segment(int sigmaS, float sigmaR, int minRegion, SpeedUpLevel speedUpLevel)
713 {
714 
715  //make sure kernel is properly defined...
716  if((!h)||(kp < 2))
717  {
718  ErrorHandler("msImageProcessor", "Segment", "Kernel corrupt or undefined.");
719  return;
720  }
721 
722  //Apply mean shift to data set using sigmaS and sigmaR...
723  Filter(sigmaS, sigmaR, speedUpLevel);
724 
725  //check for errors
726  if(ErrorStatus == EL_ERROR)
727  return;
728 
729  //check to see if the system has been halted, if so exit
730  if(ErrorStatus == EL_HALT)
731  return;
732 
733  //Check to see if the algorithm is to be halted, if so then
734  //destroy output and exit
735  if((ErrorStatus = msSys.Progress((float)(0.85))) == EL_HALT)
736  {
737  DestroyOutput();
738  return;
739  }
740 
741 #ifdef PROMPT
742  msSys.Prompt("Applying transitive closure...");
743  msSys.StartTimer();
744 #endif
745 
746  //allocate memory visit table
747  visitTable = new unsigned char [L];
748 
749  //Apply transitive closure iteratively to the regions classified
750  //by the RAM updating labels and modes until the color of each neighboring
751  //region is within sqrt(rR2) of one another.
752  rR2 = (float)(h[1]*h[1]*0.25);
753  TransitiveClosure();
754  int oldRC = regionCount;
755  int deltaRC, counter = 0;
756  do {
757  TransitiveClosure();
758  deltaRC = oldRC-regionCount;
759  oldRC = regionCount;
760  counter++;
761  } while ((deltaRC <= 0)&&(counter < 10));
762 
763  //de-allocate memory for visit table
764  delete [] visitTable;
765  visitTable = NULL;
766 
767  //Check to see if the algorithm is to be halted, if so then
768  //destroy output and regions adjacency matrix and exit
769  if((ErrorStatus = msSys.Progress((float)(0.95))) == EL_HALT)
770  {
771  DestroyRAM();
772  DestroyOutput();
773  return;
774  }
775 
776 #ifdef PROMPT
777  double timer = msSys.ElapsedTime();
778  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d).\nPruning spurious regions\t... ", timer, regionCount);
779  msSys.StartTimer();
780 #endif
781 
782  //Prune spurious regions (regions whose area is under
783  //minRegion) using RAM
784  Prune(minRegion);
785 
786 #ifdef PROMPT
787  timer = msSys.ElapsedTime();
788  msSys.Prompt("done. (%6.2f seconds, numRegions = %6d)\nPruning spurious regions ...", timer, regionCount);
789  msSys.StartTimer();
790 #endif
791 
792  //Check to see if the algorithm is to be halted, if so then
793  //destroy output and regions adjacency matrix and exit
794  if((ErrorStatus = msSys.Progress(1.0)) == EL_HALT)
795  {
796  DestroyRAM();
797  DestroyOutput();
798  return;
799  }
800 
801  //de-allocate memory for region adjacency matrix
802  DestroyRAM();
803 
804  //output to msRawData
805  int j, i, label;
806  for(i = 0; i < L; i++)
807  {
808  label = labels[i];
809  for(j = 0; j < N; j++)
810  {
811  msRawData[N*i+j] = modes[N*label+j];
812  }
813  }
814 
815  //done.
816  return;
817 
818 }
819 
820  /*/\/\/\/\/\/\/\/\/\/\/\/\*/
821  /* Data Space Conversion */
822  /*\/\/\/\/\/\/\/\/\/\/\/\/*/
823 
824 /*******************************************************/
825 /*RGB To LUV */
826 /*******************************************************/
827 /*Converts an RGB vector to LUV. */
828 /* */
829 /*See: */
830 /* G. Wyszecki and W.S. Stiles: Color Science: */
831 /* Concepts and Methods, Quantitative Data and */
832 /* Formulae, Wiley, New York, 1982. */
833 /*******************************************************/
834 /*Pre: */
835 /* - rgbVal is an unsigned char array containing */
836 /* the RGB vector */
837 /* - luvVal is a floating point array containing */
838 /* the resulting LUV vector */
839 /*Post: */
840 /* - rgbVal has been converted to LUV and the */
841 /* result has been stored in luvVal. */
842 /*******************************************************/
843 
844 void msImageProcessor::RGBtoLUV(byte *rgbVal, float *luvVal)
845 {
846 
847  //delcare variables
848  double x, y, z, L0, u_prime, v_prime, constant;
849 
850  //convert RGB to XYZ...
851  x = XYZ[0][0]*rgbVal[0] + XYZ[0][1]*rgbVal[1] + XYZ[0][2]*rgbVal[2];
852  y = XYZ[1][0]*rgbVal[0] + XYZ[1][1]*rgbVal[1] + XYZ[1][2]*rgbVal[2];
853  z = XYZ[2][0]*rgbVal[0] + XYZ[2][1]*rgbVal[1] + XYZ[2][2]*rgbVal[2];
854 
855  //convert XYZ to LUV...
856 
857  //compute L*
858  L0 = y / (255.0 * Yn);
859  if(L0 > Lt)
860  luvVal[0] = (float)(116.0 * (pow(L0, 1.0/3.0)) - 16.0);
861  else
862  luvVal[0] = (float)(903.3 * L0);
863 
864  //compute u_prime and v_prime
865  constant = x + 15 * y + 3 * z;
866  if(constant != 0)
867  {
868  u_prime = (4 * x) / constant;
869  v_prime = (9 * y) / constant;
870  }
871  else
872  {
873  u_prime = 4.0;
874  v_prime = 9.0/15.0;
875  }
876 
877  //compute u* and v*
878  luvVal[1] = (float) (13 * luvVal[0] * (u_prime - Un_prime));
879  luvVal[2] = (float) (13 * luvVal[0] * (v_prime - Vn_prime));
880 
881  //done.
882  return;
883 
884 }
885 
886 /*******************************************************/
887 /*LUV To RGB */
888 /*******************************************************/
889 /*Converts an LUV vector to RGB. */
890 /*******************************************************/
891 /*Pre: */
892 /* - luvVal is a floating point array containing */
893 /* the LUV vector */
894 /* - rgbVal is an unsigned char array containing */
895 /* the resulting RGB vector */
896 /*Post: */
897 /* - luvVal has been converted to RGB and the */
898 /* result has been stored in rgbVal. */
899 /*******************************************************/
900 
901 //define inline rounding function...
902 inline int my_round(double in_x)
903 {
904  if (in_x < 0)
905  return (int)(in_x - 0.5);
906  else
907  return (int)(in_x + 0.5);
908 }
909 
910 void msImageProcessor::LUVtoRGB(float *luvVal, byte *rgbVal)
911 {
912 
913  //declare variables...
914  int r, g, b;
915  double x, y, z, u_prime, v_prime;
916 
917  //perform conversion
918  if(luvVal[0] < 0.1)
919  r = g = b = 0;
920  else
921  {
922  //convert luv to xyz...
923  if(luvVal[0] < 8.0)
924  y = Yn * luvVal[0] / 903.3;
925  else
926  {
927  y = (luvVal[0] + 16.0) / 116.0;
928  y *= Yn * y * y;
929  }
930 
931  u_prime = luvVal[1] / (13 * luvVal[0]) + Un_prime;
932  v_prime = luvVal[2] / (13 * luvVal[0]) + Vn_prime;
933 
934  x = 9 * u_prime * y / (4 * v_prime);
935  z = (12 - 3 * u_prime - 20 * v_prime) * y / (4 * v_prime);
936 
937  //convert xyz to rgb...
938  //[r, g, b] = RGB*[x, y, z]*255.0
939  r = my_round((RGB[0][0]*x + RGB[0][1]*y + RGB[0][2]*z)*255.0);
940  g = my_round((RGB[1][0]*x + RGB[1][1]*y + RGB[1][2]*z)*255.0);
941  b = my_round((RGB[2][0]*x + RGB[2][1]*y + RGB[2][2]*z)*255.0);
942 
943  //check bounds...
944  if(r < 0) r = 0; if(r > 255) r = 255;
945  if(g < 0) g = 0; if(g > 255) g = 255;
946  if(b < 0) b = 0; if(b > 255) b = 255;
947 
948  }
949 
950  //assign rgb values to rgb vector rgbVal
951  rgbVal[0] = r;
952  rgbVal[1] = g;
953  rgbVal[2] = b;
954 
955  //done.
956  return;
957 
958 }
959 
960  /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/
961  /* Filtered and Segmented Image Output */
962  /*\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/*/
963 
964 /*******************************************************/
965 /*Get Raw Data */
966 /*******************************************************/
967 /*The output image data is returned. */
968 /*******************************************************/
969 /*Pre: */
970 /* - outputImageData is a pre-allocated floating */
971 /* point array used to store the filtered or */
972 /* segmented image pixels. */
973 /*Post: */
974 /* - the filtered or segmented image data is */
975 /* stored by outputImageData. */
976 /*******************************************************/
977 
978 void msImageProcessor::GetRawData(float *outputImageData)
979 {
980  //make sure that outputImageData is not NULL
981  if(!outputImageData)
982  {
983  ErrorHandler("msImageProcessor", "GetRawData", "Output image data buffer is NULL.");
984  return;
985  }
986 
987  //copy msRawData to outputImageData
988  int i;
989  for(i = 0; i < L*N; i++)
990  outputImageData[i] = msRawData[i];
991 
992  //done.
993  return;
994 }
995 
996 /*******************************************************/
997 /*Get Results */
998 /*******************************************************/
999 /*The output image is returned. */
1000 /*******************************************************/
1001 /*Pre: */
1002 /* - outputImage is a pre-allocated unsinged char */
1003 /* array used to store the filtered or segment- */
1004 /* ed image pixels */
1005 /*Post: */
1006 /* - the filtered or segmented image is stored by */
1007 /* outputImage. */
1008 /*******************************************************/
1009 
1010 void msImageProcessor::GetResults(byte *outputImage)
1011 {
1012 
1013  //make sure that outpuImage is not NULL
1014  if(!outputImage)
1015  {
1016  ErrorHandler("msImageProcessor", "GetResults", "Output image buffer is NULL.");
1017  return;
1018  }
1019 
1020  //if the image type is GREYSCALE simply
1021  //copy it over to the segmentedImage
1022  if(N == 1)
1023  {
1024  //copy over msRawData to segmentedImage checking
1025  //bounds
1026  int i, pxValue;
1027  for(i = 0; i < L; i++)
1028  {
1029 
1030  //get value
1031  pxValue = (int)(msRawData[i]+0.5);
1032 
1033  //store into segmented image checking bounds...
1034  if(pxValue < 0)
1035  outputImage[i] = (byte)(0);
1036  else if(pxValue > 255)
1037  outputImage[i] = (byte)(255);
1038  else
1039  outputImage[i] = (byte)(pxValue);
1040 
1041  }
1042 
1043  }
1044  else if (N == 3)
1045  {
1046 
1047  //otherwise convert msRawData from LUV to RGB
1048  //storing the result in segmentedImage
1049  int i;
1050  for(i = 0; i < L; i++)
1051  LUVtoRGB(&msRawData[N*i], &outputImage[N*i]);
1052 
1053  }
1054  else
1055  //Unknown image type: should use MeanShift::GetRawData()...
1056  ErrorHandler("msImageProcessor", "GetResults", "Unknown image type. Try using MeanShift::GetRawData().");
1057 
1058  //done.
1059  return;
1060 
1061 }
1062 
1063 /*******************************************************/
1064 /*Get Boundaries */
1065 /*******************************************************/
1066 /*A region list containing the boundary locations for */
1067 /*each region is returned. */
1068 /*******************************************************/
1069 /*Post: */
1070 /* - a region list object containing the boundary */
1071 /* locations for each region is constructed */
1072 /* - the region list is returned */
1073 /* - NULL is returned if the image has not been */
1074 /* filtered or segmented */
1075 /*******************************************************/
1076 
1077 RegionList *msImageProcessor::GetBoundaries( void )
1078 {
1079 
1080  //define bounds using label information
1081  if(class_state.OUTPUT_DEFINED)
1082  DefineBoundaries();
1083 
1084  //return region list structure
1085  return regionList;
1086 
1087 }
1088 
1089 /*******************************************************/
1090 /*Get Regions */
1091 /*******************************************************/
1092 /*Returns the regions of the processed image. */
1093 /*******************************************************/
1094 /*Pre: */
1095 /* - labels_out is an integer array of size */
1096 /* height*width that stores for each pixel a */
1097 /* label relating that pixel to a corresponding */
1098 /* region in the image */
1099 /* - modes_out is floating point array of size */
1100 /* regionCount*N storing the feature component */
1101 /* of each region, and indexed by region label */
1102 /* - modePointCounts is an integer array of size */
1103 /* regionCount, indexed by region label, that */
1104 /* stores the area of each region in pixels. */
1105 /*Post: */
1106 /* If an input image was defined and processed, */
1107 /* - memory has been allocated for labels_out, */
1108 /* modes_out and MPC_out. */
1109 /* - labels_out, modes_out, and MPC_out have been */
1110 /* populated. */
1111 /* - the number of regions contained by the segm- */
1112 /* ented image has been returned. */
1113 /* If the image has not been defined or processed */
1114 /* or if there is in-sufficient memory, */
1115 /* - no memory has been allocated for labels_out, */
1116 /* modes_out, and MPC_out. */
1117 /* - -1 is returned for regionCount. */
1118 /*******************************************************/
1119 
1120 int msImageProcessor::GetRegions(int **labels_out, float **modes_out, int **MPC_out)
1121 {
1122  //check to see if output has been defined for the given input image...
1123  if(class_state.OUTPUT_DEFINED == false)
1124  return -1;
1125 
1126  //allocate memory for labels_out, modes_out and MPC_out based
1127  //on output storage structure
1128 
1129  //ALEX BERNARDINO 29/07/09 - I had to change this in order to get the outputs
1130  //int *labels_ = *labels_out, *MPC_out_ = *MPC_out;
1131  //float *modes_ = *modes_out;
1132  int *labels_;
1133  int *MPC_out_;
1134  float *modes_;
1135  //END MODIFICATION
1136 
1137 
1138  if(!(labels_ = new int [L]))
1139  {
1140  ErrorHandler("msImageProcessor", "GetRegions", "Not enough memory.");
1141  return -1;
1142  }
1143  if(!(modes_ = new float [regionCount*N]))
1144  {
1145  ErrorHandler("msImageProcessor", "GetRegions", "Not enough memory.");
1146  return -1;
1147  }
1148  if(!(MPC_out_ = new int [regionCount]))
1149  {
1150  ErrorHandler("msImageProcessor", "GetRegions", "Not enough memory.");
1151  return -1;
1152  }
1153 
1154  //populate labels_out with image labels
1155  int i;
1156  for(i = 0; i < L; i++)
1157  labels_[i] = labels[i];
1158 
1159  //populate modes_out and MPC_out with the color and point
1160  //count of each region
1161  for(i = 0; i < regionCount*N; i++)
1162  modes_[i] = modes[i];
1163  for(i = 0; i < regionCount; i++)
1164  MPC_out_[i] = modePointCounts[i];
1165 
1166  //done. Return the number of regions resulting from filtering or segmentation.
1167 
1168  //ALEX BERNARDINO 29/07/09 - I had to change this in order to get the outputs
1169  if(labels_out)
1170  *labels_out = labels;
1171  if(modes_out)
1172  *modes_out = modes_;
1173  if(MPC_out)
1174  *MPC_out = MPC_out_;
1175  //END MODIFICATION
1176 
1177 
1178 
1179  return regionCount;
1180 }
1181 
1182 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
1183 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
1184 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ PRIVATE METHODS @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
1185 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
1186 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
1187 
1188  /*/\/\/\/\/\/\/\/\/\*/
1189  /* Image Filtering */
1190  /*\/\/\/\/\/\/\/\/\/*/
1191 
1192 /*******************************************************/
1193 /*Non Optimized Filter */
1194 /*******************************************************/
1195 /*Performs mean shift filtering on the specified input */
1196 /*image using a user defined kernel. */
1197 /*******************************************************/
1198 /*Pre: */
1199 /* - the user defined kernel used to apply mean */
1200 /* shift filtering to the defined input image */
1201 /* has spatial bandwidth sigmaS and range band- */
1202 /* width sigmaR */
1203 /* - a data set has been defined */
1204 /* - the height and width of the lattice has been */
1205 /* specified using method DefineLattice() */
1206 /*Post: */
1207 /* - mean shift filtering has been applied to the */
1208 /* input image using a user defined kernel */
1209 /* - the filtered image is stored in the private */
1210 /* data members of the msImageProcessor class. */
1211 /*******************************************************/
1212 
1213 void msImageProcessor::NonOptimizedFilter(float sigmaS, float sigmaR)
1214 {
1215 
1216  // Declare Variables
1217  int iterationCount, i, j;
1218  double mvAbs;
1219 
1220  //make sure that a lattice height and width have
1221  //been defined...
1222  if(!height)
1223  {
1224  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
1225  return;
1226  }
1227 
1228  //re-assign bandwidths to sigmaS and sigmaR
1229  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
1230  {
1231  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
1232  return;
1233  }
1234 
1235  //define input data dimension with lattice
1236  int lN = N + 2;
1237 
1238  // Traverse each data point applying mean shift
1239  // to each data point
1240 
1241  // Allcocate memory for yk
1242  double *yk = new double [lN];
1243 
1244  // Allocate memory for Mh
1245  double *Mh = new double [lN];
1246 
1247  // proceed ...
1248 #ifdef PROMPT
1249  msSys.Prompt("done.\nApplying mean shift (Using Lattice)... ");
1250 #ifdef SHOW_PROGRESS
1251  msSys.Prompt("\n 0%%");
1252 #endif
1253 #endif
1254 
1255  for(i = 0; i < L; i++)
1256  {
1257 
1258  // Assign window center (window centers are
1259  // initialized by createLattice to be the point
1260  // data[i])
1261  yk[0] = i%width;
1262  yk[1] = i/width;
1263  for(j = 0; j < N; j++)
1264  yk[j+2] = data[N*i+j];
1265 
1266  // Calculate the mean shift vector using the lattice
1267  LatticeMSVector(Mh, yk);
1268 
1269  // Calculate its magnitude squared
1270  mvAbs = 0;
1271  for(j = 0; j < lN; j++)
1272  mvAbs += Mh[j]*Mh[j];
1273 
1274  // Keep shifting window center until the magnitude squared of the
1275  // mean shift vector calculated at the window center location is
1276  // under a specified threshold (Epsilon)
1277 
1278  // NOTE: iteration count is for speed up purposes only - it
1279  // does not have any theoretical importance
1280  iterationCount = 1;
1281  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
1282  {
1283 
1284  // Shift window location
1285  for(j = 0; j < lN; j++)
1286  yk[j] += Mh[j];
1287 
1288  // Calculate the mean shift vector at the new
1289  // window location using lattice
1290  LatticeMSVector(Mh, yk);
1291 
1292  // Calculate its magnitude squared
1293  mvAbs = 0;
1294  for(j = 0; j < lN; j++)
1295  mvAbs += Mh[j]*Mh[j];
1296 
1297  // Increment interation count
1298  iterationCount++;
1299 
1300  }
1301 
1302  // Shift window location
1303  for(j = 0; j < lN; j++)
1304  yk[j] += Mh[j];
1305 
1306  //store result into msRawData...
1307  for(j = 0; j < N; j++)
1308  msRawData[N*i+j] = (float)(yk[j+2]);
1309 
1310  // Prompt user on progress
1311 #ifdef SHOW_PROGRESS
1312  percent_complete = (float)(i/(float)(L))*100;
1313  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
1314 #endif
1315 
1316  // Check to see if the algorithm has been halted
1317  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
1318  break;
1319  }
1320 
1321  // Prompt user that filtering is completed
1322 #ifdef PROMPT
1323 #ifdef SHOW_PROGRESS
1324  msSys.Prompt("\r");
1325 #endif
1326  msSys.Prompt("done.");
1327 #endif
1328 
1329  // de-allocate memory
1330  delete [] yk;
1331  delete [] Mh;
1332 
1333  // done.
1334  return;
1335 
1336 }
1337 
1338 
1339 
1340 /*******************************************************/
1341 /*Optimized Filter 1 */
1342 /*******************************************************/
1343 /*Performs mean shift filtering on the specified input */
1344 /*image using a user defined kernel. Previous mode */
1345 /*information is used to avoid re-applying mean shift */
1346 /*on certain data points to improve performance. */
1347 /*******************************************************/
1348 /*Pre: */
1349 /* - the user defined kernel used to apply mean */
1350 /* shift filtering to the defined input image */
1351 /* has spatial bandwidth sigmaS and range band- */
1352 /* width sigmaR */
1353 /* - a data set has been defined */
1354 /* - the height and width of the lattice has been */
1355 /* specified using method DefineLattice() */
1356 /*Post: */
1357 /* - mean shift filtering has been applied to the */
1358 /* input image using a user defined kernel */
1359 /* - the filtered image is stored in the private */
1360 /* data members of the msImageProcessor class. */
1361 /*******************************************************/
1362 
1363 void msImageProcessor::OptimizedFilter1(float sigmaS, float sigmaR)
1364 {
1365 
1366  // Declare Variables
1367  int iterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;
1368  float *modeCandidatePoint;
1369  double mvAbs, diff, el;
1370 
1371  //make sure that a lattice height and width have
1372  //been defined...
1373  if(!height)
1374  {
1375  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
1376  return;
1377  }
1378 
1379  //re-assign bandwidths to sigmaS and sigmaR
1380  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
1381  {
1382  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
1383  return;
1384  }
1385 
1386  //define input data dimension with lattice
1387  int lN = N + 2;
1388 
1389  // Traverse each data point applying mean shift
1390  // to each data point
1391 
1392  // Allcocate memory for yk
1393  double *yk = new double [lN];
1394 
1395  // Allocate memory for Mh
1396  double *Mh = new double [lN];
1397 
1398  // Initialize mode table used for basin of attraction
1399  memset(modeTable, 0, width*height);
1400 
1401  // Allocate memory mode candidate data point...
1402  //floating point version
1403  modeCandidatePoint = new float [N];
1404 
1405  // proceed ...
1406 #ifdef PROMPT
1407  msSys.Prompt("done.\nApplying mean shift (Using Lattice) ... ");
1408 #ifdef SHOW_PROGRESS
1409  msSys.Prompt("\n 0%%");
1410 #endif
1411 #endif
1412 
1413 
1414  for(i = 0; i < L; i++)
1415  {
1416  // if a mode was already assigned to this data point
1417  // then skip this point, otherwise proceed to
1418  // find its mode by applying mean shift...
1419  if (modeTable[i] == 1)
1420  continue;
1421 
1422  // initialize point list...
1423  pointCount = 0;
1424 
1425  // Assign window center (window centers are
1426  // initialized by createLattice to be the point
1427  // data[i])
1428  yk[0] = i%width;
1429  yk[1] = i/width;
1430  for(j = 0; j < N; j++)
1431  yk[j+2] = data[N*i+j];
1432 
1433  // Calculate the mean shift vector using the lattice
1434  LatticeMSVector(Mh, yk);
1435 
1436  // Calculate its magnitude squared
1437  mvAbs = 0;
1438  for(j = 0; j < lN; j++)
1439  mvAbs += Mh[j]*Mh[j];
1440 
1441  // Keep shifting window center until the magnitude squared of the
1442  // mean shift vector calculated at the window center location is
1443  // under a specified threshold (Epsilon)
1444 
1445  // NOTE: iteration count is for speed up purposes only - it
1446  // does not have any theoretical importance
1447  iterationCount = 1;
1448  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
1449  {
1450 
1451  // Shift window location
1452  for(j = 0; j < lN; j++)
1453  yk[j] += Mh[j];
1454 
1455  // check to see if the current mode location is in the
1456  // basin of attraction...
1457 
1458  // calculate the location of yk on the lattice
1459  modeCandidateX = (int) (yk[0]+0.5);
1460  modeCandidateY = (int) (yk[1]+0.5);
1461  modeCandidate_i = modeCandidateY*width + modeCandidateX;
1462 
1463  // if mvAbs != 0 (yk did indeed move) then check
1464  // location basin_i in the mode table to see if
1465  // this data point either:
1466 
1467  // (1) has not been associated with a mode yet
1468  // (modeTable[basin_i] = 0), so associate
1469  // it with this one
1470  //
1471  // (2) it has been associated with a mode other
1472  // than the one that this data point is converging
1473  // to (modeTable[basin_i] = 1), so assign to
1474  // this data point the same mode as that of basin_i
1475 
1476  if ((modeTable[modeCandidate_i] != 2) && (modeCandidate_i != i))
1477  {
1478  // obtain the data point at basin_i to
1479  // see if it is within h*TC_DIST_FACTOR of
1480  // of yk
1481  for (j = 0; j < N; j++)
1482  modeCandidatePoint[j] = data[N*modeCandidate_i + j];
1483 
1484  // check basin on non-spatial data spaces only
1485  k = 1;
1486  s = 0;
1487  diff = 0;
1488  while ((diff < TC_DIST_FACTOR) && (k<kp))
1489  {
1490  diff = 0;
1491  for (p=0; p<P[k]; p++)
1492  {
1493  el = (modeCandidatePoint[p+s]-yk[p+s+2])/h[k];
1494  diff += el*el;
1495  }
1496  s+=P[k];
1497  k++;
1498  }
1499 
1500  // if the data point at basin_i is within
1501  // a distance of h*TC_DIST_FACTOR of yk
1502  // then depending on modeTable[basin_i] perform
1503  // either (1) or (2)
1504  if (diff < TC_DIST_FACTOR)
1505  {
1506  // if the data point at basin_i has not
1507  // been associated to a mode then associate
1508  // it with the mode that this one will converge
1509  // to
1510  if (modeTable[modeCandidate_i] == 0)
1511  {
1512  // no mode associated yet so associate
1513  // it with this one...
1514  pointList[pointCount++] = modeCandidate_i;
1515  modeTable[modeCandidate_i] = 2;
1516 
1517  } else
1518  {
1519 
1520  // the mode has already been associated with
1521  // another mode, thererfore associate this one
1522  // mode and the modes in the point list with
1523  // the mode associated with data[basin_i]...
1524 
1525  // store the mode info into yk using msRawData...
1526  for (j = 0; j < N; j++)
1527  yk[j+2] = msRawData[modeCandidate_i*N+j];
1528 
1529  // update mode table for this data point
1530  // indicating that a mode has been associated
1531  // with it
1532  modeTable[i] = 1;
1533 
1534  // indicate that a mode has been associated
1535  // to this data point (data[i])
1536  mvAbs = -1;
1537 
1538  // stop mean shift calculation...
1539  break;
1540  }
1541  }
1542  }
1543 
1544  // Calculate the mean shift vector at the new
1545  // window location using lattice
1546  LatticeMSVector(Mh, yk);
1547 
1548  // Calculate its magnitude squared
1549  mvAbs = 0;
1550  for(j = 0; j < lN; j++)
1551  mvAbs += Mh[j]*Mh[j];
1552 
1553  // Increment iteration count
1554  iterationCount++;
1555 
1556  }
1557 
1558  // if a mode was not associated with this data point
1559  // yet associate it with yk...
1560  if (mvAbs >= 0)
1561  {
1562  // Shift window location
1563  for(j = 0; j < lN; j++)
1564  yk[j] += Mh[j];
1565 
1566  // update mode table for this data point
1567  // indicating that a mode has been associated
1568  // with it
1569  modeTable[i] = 1;
1570  }
1571 
1572  // associate the data point indexed by
1573  // the point list with the mode stored
1574  // by yk
1575  for (j = 0; j < pointCount; j++)
1576  {
1577  // obtain the point location from the
1578  // point list
1579  modeCandidate_i = pointList[j];
1580 
1581  // update the mode table for this point
1582  modeTable[modeCandidate_i] = 1;
1583 
1584  //store result into msRawData...
1585  for(k = 0; k < N; k++)
1586  msRawData[N*modeCandidate_i+k] = (float)(yk[k+2]);
1587  }
1588 
1589 
1590  //store result into msRawData...
1591  for(j = 0; j < N; j++)
1592  msRawData[N*i+j] = (float)(yk[j+2]);
1593 
1594  // Prompt user on progress
1595 #ifdef SHOW_PROGRESS
1596  percent_complete = (float)(i/(float)(L))*100;
1597  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
1598 #endif
1599 
1600  // Check to see if the algorithm has been halted
1601  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
1602  break;
1603  }
1604 
1605  // Prompt user that filtering is completed
1606 #ifdef PROMPT
1607 #ifdef SHOW_PROGRESS
1608  msSys.Prompt("\r");
1609 #endif
1610  msSys.Prompt("done.");
1611 #endif
1612 
1613  // de-allocate memory
1614  delete [] modeCandidatePoint;
1615  delete [] yk;
1616  delete [] Mh;
1617 
1618  // done.
1619  return;
1620 
1621 }
1622 
1623 /*******************************************************/
1624 /*Optimized Filter 2 */
1625 /*******************************************************/
1626 /*Performs mean shift filtering on the specified input */
1627 /*image using a user defined kernel. Previous mode */
1628 /*information is used to avoid re-applying mean shift */
1629 /*on certain data points to improve performance. To */
1630 /*further improve perfmance (during segmentation) poi- */
1631 /*nts within h of a window center during the window */
1632 /*center's traversal to a mode are associated with the */
1633 /*mode that the window converges to. */
1634 /*******************************************************/
1635 /*Pre: */
1636 /* - the user defined kernel used to apply mean */
1637 /* shift filtering to the defined input image */
1638 /* has spatial bandwidth sigmaS and range band- */
1639 /* width sigmaR */
1640 /* - a data set has been defined */
1641 /* - the height and width of the lattice has been */
1642 /* specified using method DefineLattice() */
1643 /*Post: */
1644 /* - mean shift filtering has been applied to the */
1645 /* input image using a user defined kernel */
1646 /* - the filtered image is stored in the private */
1647 /* data members of the msImageProcessor class. */
1648 /*******************************************************/
1649 
1650 void msImageProcessor::OptimizedFilter2(float sigmaS, float sigmaR)
1651 {
1652 
1653  //if confidence map is null set it to zero
1654  if(!weightMap)
1655  {
1656  weightMap = new float [L];
1657  int i;
1658  for(i = 0; i < L; i++)
1659  weightMap[i] = 0;
1660  }
1661 
1662  // Declare Variables
1663  int iterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;
1664  float *modeCandidatePoint;
1665  double mvAbs, diff, el;
1666 
1667  //make sure that a lattice height and width have
1668  //been defined...
1669  if(!height)
1670  {
1671  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
1672  return;
1673  }
1674 
1675  //re-assign bandwidths to sigmaS and sigmaR
1676  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
1677  {
1678  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
1679  return;
1680  }
1681 
1682  //define input data dimension with lattice
1683  int lN = N + 2;
1684 
1685  // Traverse each data point applying mean shift
1686  // to each data point
1687 
1688  // Allcocate memory for yk
1689  double *yk = new double [lN];
1690 
1691  // Allocate memory for Mh
1692  double *Mh = new double [lN];
1693 
1694  // Initialize mode table used for basin of attraction
1695  memset(modeTable, 0, width*height);
1696 
1697  // Allocate memory mode candidate data point...
1698  //floating point version
1699  modeCandidatePoint = new float [N];
1700 
1701  // proceed ...
1702 #ifdef PROMPT
1703  msSys.Prompt("done.\nApplying mean shift (Using Lattice)... ");
1704 #ifdef SHOW_PROGRESS
1705  msSys.Prompt("\n 0%%");
1706 #endif
1707 #endif
1708 
1709  for(i = 0; i < L; i++)
1710  {
1711  // if a mode was already assigned to this data point
1712  // then skip this point, otherwise proceed to
1713  // find its mode by applying mean shift...
1714  if (modeTable[i] == 1)
1715  continue;
1716 
1717  // initialize point list...
1718  pointCount = 0;
1719 
1720  // Assign window center (window centers are
1721  // initialized by createLattice to be the point
1722  // data[i])
1723  yk[0] = i%width;
1724  yk[1] = i/width;
1725  for(j = 0; j < N; j++)
1726  yk[j+2] = data[N*i+j];
1727 
1728  // Calculate the mean shift vector using the lattice
1729  OptLatticeMSVector(Mh, yk);
1730 
1731  // Calculate its magnitude squared
1732  mvAbs = 0;
1733  for(j = 0; j < lN; j++)
1734  mvAbs += Mh[j]*Mh[j];
1735 
1736  // Keep shifting window center until the magnitude squared of the
1737  // mean shift vector calculated at the window center location is
1738  // under a specified threshold (Epsilon)
1739 
1740  // NOTE: iteration count is for speed up purposes only - it
1741  // does not have any theoretical importance
1742  iterationCount = 1;
1743  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
1744  {
1745 
1746  // Shift window location
1747  for(j = 0; j < lN; j++)
1748  yk[j] += Mh[j];
1749 
1750  // check to see if the current mode location is in the
1751  // basin of attraction...
1752 
1753  // calculate the location of yk on the lattice
1754  modeCandidateX = (int) (yk[0]+0.5);
1755  modeCandidateY = (int) (yk[1]+0.5);
1756  modeCandidate_i = modeCandidateY*width + modeCandidateX;
1757 
1758  // if mvAbs != 0 (yk did indeed move) then check
1759  // location basin_i in the mode table to see if
1760  // this data point either:
1761 
1762  // (1) has not been associated with a mode yet
1763  // (modeTable[basin_i] = 0), so associate
1764  // it with this one
1765  //
1766  // (2) it has been associated with a mode other
1767  // than the one that this data point is converging
1768  // to (modeTable[basin_i] = 1), so assign to
1769  // this data point the same mode as that of basin_i
1770 
1771  if ((modeTable[modeCandidate_i] != 2) && (modeCandidate_i != i))
1772  {
1773  // obtain the data point at basin_i to
1774  // see if it is within h*TC_DIST_FACTOR of
1775  // of yk
1776  for (j = 0; j < N; j++)
1777  modeCandidatePoint[j] = data[N*modeCandidate_i + j];
1778 
1779  // check basin on non-spatial data spaces only
1780  k = 1;
1781  s = 0;
1782  diff = 0;
1783  while ((diff < TC_DIST_FACTOR) && (k<kp))
1784  {
1785  diff = 0;
1786  for (p=0; p<P[k]; p++)
1787  {
1788  el = (modeCandidatePoint[p+s]-yk[p+s+2])/h[k];
1789  diff += el*el;
1790  }
1791  s+=P[k];
1792  k++;
1793  }
1794 
1795  // if the data point at basin_i is within
1796  // a distance of h*TC_DIST_FACTOR of yk
1797  // then depending on modeTable[basin_i] perform
1798  // either (1) or (2)
1799  if (diff < TC_DIST_FACTOR)
1800  {
1801  // if the data point at basin_i has not
1802  // been associated to a mode then associate
1803  // it with the mode that this one will converge
1804  // to
1805  if (modeTable[modeCandidate_i] == 0)
1806  {
1807  // no mode associated yet so associate
1808  // it with this one...
1809  pointList[pointCount++] = modeCandidate_i;
1810  modeTable[modeCandidate_i] = 2;
1811 
1812  } else
1813  {
1814 
1815  // the mode has already been associated with
1816  // another mode, thererfore associate this one
1817  // mode and the modes in the point list with
1818  // the mode associated with data[basin_i]...
1819 
1820  // store the mode infor int yk using msRawData...
1821  for (j = 0; j < N; j++)
1822  yk[j+2] = msRawData[modeCandidate_i*N+j];
1823 
1824  // update mode table for this data point
1825  // indicating that a mode has been associated
1826  // with it
1827  modeTable[i] = 1;
1828 
1829  // indicate that a mode has been associated
1830  // to this data point (data[i])
1831  mvAbs = -1;
1832 
1833  // stop mean shift calculation...
1834  break;
1835  }
1836  }
1837  }
1838 
1839  // Calculate the mean shift vector at the new
1840  // window location using lattice
1841  OptLatticeMSVector(Mh, yk);
1842 
1843  // Calculate its magnitude squared
1844  mvAbs = 0;
1845  for(j = 0; j < lN; j++)
1846  mvAbs += Mh[j]*Mh[j];
1847 
1848  // Increment interation count
1849  iterationCount++;
1850 
1851  }
1852 
1853  // if a mode was not associated with this data point
1854  // yet then perform a shift the window center yk one
1855  // last time using the mean shift vector...
1856  if (mvAbs >= 0)
1857  {
1858  // Shift window location
1859  for(j = 0; j < lN; j++)
1860  yk[j] += Mh[j];
1861 
1862  // update mode table for this data point
1863  // indicating that a mode has been associated
1864  // with it
1865  modeTable[i] = 1;
1866  }
1867 
1868  // associate the data point indexed by
1869  // the point list with the mode stored
1870  // by yk
1871  for (j = 0; j < pointCount; j++)
1872  {
1873  // obtain the point location from the
1874  // point list
1875  modeCandidate_i = pointList[j];
1876 
1877  // update the mode table for this point
1878  modeTable[modeCandidate_i] = 1;
1879 
1880  //store result into msRawData...
1881  for(k = 0; k < N; k++)
1882  msRawData[N*modeCandidate_i+k] = (float)(yk[k+2]);
1883  }
1884 
1885 
1886  //store result into msRawData...
1887  for(j = 0; j < N; j++)
1888  msRawData[N*i+j] = (float)(yk[j+2]);
1889 
1890  // Prompt user on progress
1891 #ifdef SHOW_PROGRESS
1892  percent_complete = (float)(i/(float)(L))*100;
1893  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
1894 #endif
1895 
1896  // Check to see if the algorithm has been halted
1897  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
1898  break;
1899 
1900  }
1901 
1902  // Prompt user that filtering is completed
1903 #ifdef PROMPT
1904 #ifdef SHOW_PROGRESS
1905  msSys.Prompt("\r");
1906 #endif
1907  msSys.Prompt("done.");
1908 #endif
1909 
1910  // de-allocate memory
1911  delete [] modeCandidatePoint;
1912  delete [] yk;
1913  delete [] Mh;
1914 
1915  // done.
1916  return;
1917 
1918 }
1919 
1920  /*/\/\/\/\/\/\/\/\/\/\/\*/
1921  /* Image Classification */
1922  /*\/\/\/\/\/\/\/\/\/\/\/*/
1923 
1924 /*******************************************************/
1925 /*Connect */
1926 /*******************************************************/
1927 /*Classifies the regions of the mean shift filtered */
1928 /*image. */
1929 /*******************************************************/
1930 /*Post: */
1931 /* - the regions of the mean shift image have been*/
1932 /* classified using the private classification */
1933 /* structure of the msImageProcessor Class. */
1934 /* Namely, each region uniquely identified by */
1935 /* its LUV color (stored by LUV_data) and loc- */
1936 /* ation has been labeled and its area computed */
1937 /* via an eight-connected fill. */
1938 /*******************************************************/
1939 
1940 void msImageProcessor::Connect( void )
1941 {
1942 
1943  //define eight connected neighbors
1944  neigh[0] = 1;
1945  neigh[1] = 1-width;
1946  neigh[2] = -width;
1947  neigh[3] = -(1+width);
1948  neigh[4] = -1;
1949  neigh[5] = width-1;
1950  neigh[6] = width;
1951  neigh[7] = width+1;
1952 
1953  //initialize labels and modePointCounts
1954  int i;
1955  for(i = 0; i < width*height; i++)
1956  {
1957  labels[i] = -1;
1958  modePointCounts[i] = 0;
1959  }
1960 
1961  //Traverse the image labeling each new region encountered
1962  int k, label = -1;
1963  for(i = 0; i < height*width; i++)
1964  {
1965  //if this region has not yet been labeled - label it
1966  if(labels[i] < 0)
1967  {
1968  //assign new label to this region
1969  labels[i] = ++label;
1970 
1971  //copy region color into modes
1972  for(k = 0; k < N; k++)
1973  modes[(N*label)+k] = LUV_data[(N*i)+k];
1974 // modes[(N*label)+k] = (float)(LUV_data[(N*i)+k]);
1975 
1976  //populate labels with label for this specified region
1977  //calculating modePointCounts[label]...
1978  Fill(i, label);
1979  }
1980  }
1981 
1982  //calculate region count using label
1983  regionCount = label+1;
1984 
1985  //done.
1986  return;
1987 }
1988 
1989 /*******************************************************/
1990 /*Fill */
1991 /*******************************************************/
1992 /*Given a region seed and a region label, Fill uses */
1993 /*the region seed to perform an eight-connected fill */
1994 /*for the specified region, labeling all pixels con- */
1995 /*tained by the region with the specified label: */
1996 /*label. */
1997 /*******************************************************/
1998 /*Pre: */
1999 /* - regionLoc is a region seed - a pixel that is */
2000 /* identified as being part of the region */
2001 /* labled using the label, label. */
2002 /*Post: */
2003 /* - all pixels belonging to the region specified */
2004 /* by regionLoc (having the same integer LUV */
2005 /* value specified by LUV_data) are classified */
2006 /* as one region by labeling each pixel in the */
2007 /* image clasification structure using label */
2008 /* via an eight-connected fill. */
2009 /*******************************************************/
2010 
2011 void msImageProcessor::Fill(int regionLoc, int label)
2012 {
2013 
2014  //declare variables
2015  int i, k, neighLoc, neighborsFound, imageSize = width*height;
2016 
2017  //Fill region starting at region location
2018  //using labels...
2019 
2020  //initialzie indexTable
2021  int index = 0;
2022  indexTable[0] = regionLoc;
2023 
2024  //increment mode point counts for this region to
2025  //indicate that one pixel belongs to this region
2026  modePointCounts[label]++;
2027 
2028  while(true)
2029  {
2030 
2031  //assume no neighbors will be found
2032  neighborsFound = 0;
2033 
2034  //check the eight connected neighbors at regionLoc -
2035  //if a pixel has similar color to that located at
2036  //regionLoc then declare it as part of this region
2037  for(i = 0; i < 8; i++)
2038  {
2039  // no need
2040  /*
2041  //if at boundary do not check certain neighbors because
2042  //they do not exist...
2043  if((regionLoc%width == 0)&&((i == 3)||(i == 4)||(i == 5)))
2044  continue;
2045  if((regionLoc%(width-1) == 0)&&((i == 0)||(i == 1)||(i == 7)))
2046  continue;
2047  */
2048 
2049  //check bounds and if neighbor has been already labeled
2050  neighLoc = regionLoc + neigh[i];
2051  if((neighLoc >= 0)&&(neighLoc < imageSize)&&(labels[neighLoc] < 0))
2052  {
2053  for(k = 0; k < N; k++)
2054  {
2055 // if(LUV_data[(regionLoc*N)+k] != LUV_data[(neighLoc*N)+k])
2056  if (fabs(LUV_data[(regionLoc*N)+k]-LUV_data[(neighLoc*N)+k])>=LUV_treshold)
2057  break;
2058  }
2059 
2060  //neighbor i belongs to this region so label it and
2061  //place it onto the index table buffer for further
2062  //processing
2063  if(k == N)
2064  {
2065  //assign label to neighbor i
2066  labels[neighLoc] = label;
2067 
2068  //increment region point count
2069  modePointCounts[label]++;
2070 
2071  //place index of neighbor i onto the index tabel buffer
2072  indexTable[++index] = neighLoc;
2073 
2074  //indicate that a neighboring region pixel was
2075  //identified
2076  neighborsFound = 1;
2077  }
2078  }
2079  }
2080 
2081  //check the indexTable to see if there are any more
2082  //entries to be explored - if so explore them, otherwise
2083  //exit the loop - we are finished
2084  if(neighborsFound)
2085  regionLoc = indexTable[index];
2086  else if (index > 1)
2087  regionLoc = indexTable[--index];
2088  else
2089  break; //fill complete
2090  }
2091 
2092  //done.
2093  return;
2094 
2095 }
2096 
2097  /*/\/\/\/\/\/\/\/\*/
2098  /* Image Pruning */
2099  /*\/\/\/\/\/\/\/\/*/
2100 
2101 /*******************************************************/
2102 /*Build Region Adjacency Matrix */
2103 /*******************************************************/
2104 /*Constructs a region adjacency matrix. */
2105 /*******************************************************/
2106 /*Pre: */
2107 /* - the classification data structure has been */
2108 /* constructed. */
2109 /*Post: */
2110 /* - a region adjacency matrix has been built */
2111 /* using the classification data structure. */
2112 /*******************************************************/
2113 
2114 void msImageProcessor::BuildRAM( void )
2115 {
2116 
2117  //Allocate memory for region adjacency matrix if it hasn't already been allocated
2118  if((!raList)&&((!(raList = new RAList [regionCount]))||(!(raPool = new RAList [NODE_MULTIPLE*regionCount]))))
2119  {
2120  ErrorHandler("msImageProcessor", "Allocate", "Not enough memory.");
2121  return;
2122  }
2123 
2124  //initialize the region adjacency list
2125  int i;
2126  for(i = 0; i < regionCount; i++)
2127  {
2128  raList[i].edgeStrength = 0;
2129  raList[i].edgePixelCount = 0;
2130  raList[i].label = i;
2131  raList[i].next = NULL;
2132  }
2133 
2134  //initialize RAM free list
2135  freeRAList = raPool;
2136  for(i = 0; i < NODE_MULTIPLE*regionCount-1; i++)
2137  {
2138  raPool[i].edgeStrength = 0;
2139  raPool[i].edgePixelCount = 0;
2140  raPool[i].next = &raPool[i+1];
2141  }
2142  raPool[NODE_MULTIPLE*regionCount-1].next = NULL;
2143 
2144  //traverse the labeled image building
2145  //the RAM by looking to the right of
2146  //and below the current pixel location thus
2147  //determining if a given region is adjacent
2148  //to another
2149  int j, curLabel, rightLabel, bottomLabel, exists;
2150  RAList *raNode1, *raNode2, *oldRAFreeList;
2151  for(i = 0; i < height - 1; i++)
2152  {
2153  //check the right and below neighbors
2154  //for pixel locations whose x < width - 1
2155  for(j = 0; j < width - 1; j++)
2156  {
2157  //calculate pixel labels
2158  curLabel = labels[i*width+j ]; //current pixel
2159  rightLabel = labels[i*width+j+1 ]; //right pixel
2160  bottomLabel = labels[(i+1)*width+j]; //bottom pixel
2161 
2162  //check to the right, if the label of
2163  //the right pixel is not the same as that
2164  //of the current one then region[j] and region[j+1]
2165  //are adjacent to one another - update the RAM
2166  if(curLabel != rightLabel)
2167  {
2168  //obtain RAList object from region adjacency free
2169  //list
2170  raNode1 = freeRAList;
2171  raNode2 = freeRAList->next;
2172 
2173  //keep a pointer to the old region adj. free
2174  //list just in case nodes already exist in respective
2175  //region lists
2176  oldRAFreeList = freeRAList;
2177 
2178  //update region adjacency free list
2179  freeRAList = freeRAList->next->next;
2180 
2181  //populate RAList nodes
2182  raNode1->label = curLabel;
2183  raNode2->label = rightLabel;
2184 
2185  //insert nodes into the RAM
2186  exists = 0;
2187  raList[curLabel ].Insert(raNode2);
2188  exists = raList[rightLabel].Insert(raNode1);
2189 
2190  //if the node already exists then place
2191  //nodes back onto the region adjacency
2192  //free list
2193  if(exists)
2194  freeRAList = oldRAFreeList;
2195 
2196  }
2197 
2198  //check below, if the label of
2199  //the bottom pixel is not the same as that
2200  //of the current one then region[j] and region[j+width]
2201  //are adjacent to one another - update the RAM
2202  if(curLabel != bottomLabel)
2203  {
2204  //obtain RAList object from region adjacency free
2205  //list
2206  raNode1 = freeRAList;
2207  raNode2 = freeRAList->next;
2208 
2209  //keep a pointer to the old region adj. free
2210  //list just in case nodes already exist in respective
2211  //region lists
2212  oldRAFreeList = freeRAList;
2213 
2214  //update region adjacency free list
2215  freeRAList = freeRAList->next->next;
2216 
2217  //populate RAList nodes
2218  raNode1->label = curLabel;
2219  raNode2->label = bottomLabel;
2220 
2221  //insert nodes into the RAM
2222  exists = 0;
2223  raList[curLabel ].Insert(raNode2);
2224  exists = raList[bottomLabel].Insert(raNode1);
2225 
2226  //if the node already exists then place
2227  //nodes back onto the region adjacency
2228  //free list
2229  if(exists)
2230  freeRAList = oldRAFreeList;
2231 
2232  }
2233 
2234  }
2235 
2236  //check only to the bottom neighbors of the right boundary
2237  //pixels...
2238 
2239  //calculate pixel locations (j = width-1)
2240  curLabel = labels[i*width+j ]; //current pixel
2241  bottomLabel = labels[(i+1)*width+j]; //bottom pixel
2242 
2243  //check below, if the label of
2244  //the bottom pixel is not the same as that
2245  //of the current one then region[j] and region[j+width]
2246  //are adjacent to one another - update the RAM
2247  if(curLabel != bottomLabel)
2248  {
2249  //obtain RAList object from region adjacency free
2250  //list
2251  raNode1 = freeRAList;
2252  raNode2 = freeRAList->next;
2253 
2254  //keep a pointer to the old region adj. free
2255  //list just in case nodes already exist in respective
2256  //region lists
2257  oldRAFreeList = freeRAList;
2258 
2259  //update region adjacency free list
2260  freeRAList = freeRAList->next->next;
2261 
2262  //populate RAList nodes
2263  raNode1->label = curLabel;
2264  raNode2->label = bottomLabel;
2265 
2266  //insert nodes into the RAM
2267  exists = 0;
2268  raList[curLabel ].Insert(raNode2);
2269  exists = raList[bottomLabel].Insert(raNode1);
2270 
2271  //if the node already exists then place
2272  //nodes back onto the region adjacency
2273  //free list
2274  if(exists)
2275  freeRAList = oldRAFreeList;
2276 
2277  }
2278  }
2279 
2280  //check only to the right neighbors of the bottom boundary
2281  //pixels...
2282 
2283  //check the right for pixel locations whose x < width - 1
2284  for(j = 0; j < width - 1; j++)
2285  {
2286  //calculate pixel labels (i = height-1)
2287  curLabel = labels[i*width+j ]; //current pixel
2288  rightLabel = labels[i*width+j+1 ]; //right pixel
2289 
2290  //check to the right, if the label of
2291  //the right pixel is not the same as that
2292  //of the current one then region[j] and region[j+1]
2293  //are adjacent to one another - update the RAM
2294  if(curLabel != rightLabel)
2295  {
2296  //obtain RAList object from region adjacency free
2297  //list
2298  raNode1 = freeRAList;
2299  raNode2 = freeRAList->next;
2300 
2301  //keep a pointer to the old region adj. free
2302  //list just in case nodes already exist in respective
2303  //region lists
2304  oldRAFreeList = freeRAList;
2305 
2306  //update region adjacency free list
2307  freeRAList = freeRAList->next->next;
2308 
2309  //populate RAList nodes
2310  raNode1->label = curLabel;
2311  raNode2->label = rightLabel;
2312 
2313  //insert nodes into the RAM
2314  exists = 0;
2315  raList[curLabel ].Insert(raNode2);
2316  exists = raList[rightLabel].Insert(raNode1);
2317 
2318  //if the node already exists then place
2319  //nodes back onto the region adjacency
2320  //free list
2321  if(exists)
2322  freeRAList = oldRAFreeList;
2323 
2324  }
2325 
2326  }
2327 
2328  //done.
2329  return;
2330 
2331 }
2332 
2333 /*******************************************************/
2334 /*Destroy Region Adjacency Matrix */
2335 /*******************************************************/
2336 /*Destroy a region adjacency matrix. */
2337 /*******************************************************/
2338 /*Post: */
2339 /* - the region adjacency matrix has been destr- */
2340 /* oyed: (1) its memory has been de-allocated, */
2341 /* (2) the RAM structure has been initialize */
2342 /* for re-use. */
2343 /*******************************************************/
2344 
2345 void msImageProcessor::DestroyRAM( void )
2346 {
2347 
2348  //de-allocate memory for region adjaceny list
2349  if (raList) delete [] raList;
2350  if (raPool) delete [] raPool;
2351 
2352  //initialize region adjacency matrix
2353  raList = NULL;
2354  freeRAList = NULL;
2355  raPool = NULL;
2356 
2357  //done.
2358  return;
2359 
2360 }
2361 
2362 /*******************************************************/
2363 /*Transitive Closure */
2364 /*******************************************************/
2365 /*Applies transitive closure to the RAM updating */
2366 /*labels, modes and modePointCounts to reflect the new */
2367 /*set of merged regions resulting from transitive clo- */
2368 /*sure. */
2369 /*******************************************************/
2370 /*Post: */
2371 /* - transitive closure has been applied to the */
2372 /* regions classified by the RAM and labels, */
2373 /* modes and modePointCounts have been updated */
2374 /* to reflect the new set of mergd regions res- */
2375 /* ulting from transitive closure. */
2376 /*******************************************************/
2377 
2378 void msImageProcessor::TransitiveClosure( void )
2379 {
2380 
2381  //Step (1):
2382 
2383  // Build RAM using classifiction structure originally
2384  // generated by the method GridTable::Connect()
2385  BuildRAM();
2386 
2387  //Step (1a):
2388  //Compute weights of weight graph using confidence map
2389  //(if defined)
2390  if(weightMapDefined) ComputeEdgeStrengths();
2391 
2392  //Step (2):
2393 
2394  //Treat each region Ri as a disjoint set:
2395 
2396  // - attempt to join Ri and Rj for all i != j that are neighbors and
2397  // whose associated modes are a normalized distance of < 0.5 from one
2398  // another
2399 
2400  // - the label of each region in the raList is treated as a pointer to the
2401  // canonical element of that region (e.g. raList[i], initially has raList[i].label = i,
2402  // namely each region is initialized to have itself as its canonical element).
2403 
2404  //Traverse RAM attempting to join raList[i] with its neighbors...
2405  int i, iCanEl, neighCanEl;
2406  float threshold;
2407  RAList *neighbor;
2408  for(i = 0; i < regionCount; i++)
2409  {
2410  //aquire first neighbor in region adjacency list pointed to
2411  //by raList[i]
2412  neighbor = raList[i].next;
2413 
2414  //compute edge strenght threshold using global and local
2415  //epsilon
2416  if(epsilon > raList[i].edgeStrength)
2417  threshold = epsilon;
2418  else
2419  threshold = raList[i].edgeStrength;
2420 
2421  //traverse region adjacency list of region i, attempting to join
2422  //it with regions whose mode is a normalized distance < 0.5 from
2423  //that of region i...
2424  while(neighbor)
2425  {
2426  //attempt to join region and neighbor...
2427  if((InWindow(i, neighbor->label))&&(neighbor->edgeStrength < epsilon))
2428  {
2429  //region i and neighbor belong together so join them
2430  //by:
2431 
2432  // (1) find the canonical element of region i
2433  iCanEl = i;
2434  while(raList[iCanEl].label != iCanEl)
2435  iCanEl = raList[iCanEl].label;
2436 
2437  // (2) find the canonical element of neighboring region
2438  neighCanEl = neighbor->label;
2439  while(raList[neighCanEl].label != neighCanEl)
2440  neighCanEl = raList[neighCanEl].label;
2441 
2442  // if the canonical elements of are not the same then assign
2443  // the canonical element having the smaller label to be the parent
2444  // of the other region...
2445  if(iCanEl < neighCanEl)
2446  raList[neighCanEl].label = iCanEl;
2447  else
2448  {
2449  //must replace the canonical element of previous
2450  //parent as well
2451  raList[raList[iCanEl].label].label = neighCanEl;
2452 
2453  //re-assign canonical element
2454  raList[iCanEl].label = neighCanEl;
2455  }
2456  }
2457 
2458  //check the next neighbor...
2459  neighbor = neighbor->next;
2460 
2461  }
2462  }
2463 
2464  // Step (3):
2465 
2466  // Level binary trees formed by canonical elements
2467  for(i = 0; i < regionCount; i++)
2468  {
2469  iCanEl = i;
2470  while(raList[iCanEl].label != iCanEl)
2471  iCanEl = raList[iCanEl].label;
2472  raList[i].label = iCanEl;
2473  }
2474 
2475  // Step (4):
2476 
2477  //Traverse joint sets, relabeling image.
2478 
2479  // (a)
2480 
2481  // Accumulate modes and re-compute point counts using canonical
2482  // elements generated by step 2.
2483 
2484  //allocate memory for mode and point count temporary buffers...
2485  float *modes_buffer = new float [N*regionCount];
2486  int *MPC_buffer = new int [regionCount];
2487 
2488  //initialize buffers to zero
2489  for(i = 0; i < regionCount; i++)
2490  MPC_buffer[i] = 0;
2491  for(i = 0; i < N*regionCount; i++)
2492  modes_buffer[i] = 0;
2493 
2494  //traverse raList accumulating modes and point counts
2495  //using canoncial element information...
2496  int k, iMPC;
2497  for(i = 0; i < regionCount; i++)
2498  {
2499 
2500  //obtain canonical element of region i
2501  iCanEl = raList[i].label;
2502 
2503  //obtain mode point count of region i
2504  iMPC = modePointCounts[i];
2505 
2506  //accumulate modes_buffer[iCanEl]
2507  for(k = 0; k < N; k++)
2508  modes_buffer[(N*iCanEl)+k] += iMPC*modes[(N*i)+k];
2509 
2510  //accumulate MPC_buffer[iCanEl]
2511  MPC_buffer[iCanEl] += iMPC;
2512 
2513  }
2514 
2515  // (b)
2516 
2517  // Re-label new regions of the image using the canonical
2518  // element information generated by step (2)
2519 
2520  // Also use this information to compute the modes of the newly
2521  // defined regions, and to assign new region point counts in
2522  // a consecute manner to the modePointCounts array
2523 
2524  //allocate memory for label buffer
2525  int *label_buffer = new int [regionCount];
2526 
2527  //initialize label buffer to -1
2528  for(i = 0; i < regionCount; i++)
2529  label_buffer[i] = -1;
2530 
2531  //traverse raList re-labeling the regions
2532  int label = -1;
2533  for(i = 0; i < regionCount; i++)
2534  {
2535  //obtain canonical element of region i
2536  iCanEl = raList[i].label;
2537  if(label_buffer[iCanEl] < 0)
2538  {
2539  //assign a label to the new region indicated by canonical
2540  //element of i
2541  label_buffer[iCanEl] = ++label;
2542 
2543  //recompute mode storing the result in modes[label]...
2544  iMPC = MPC_buffer[iCanEl];
2545  for(k = 0; k < N; k++)
2546  modes[(N*label)+k] = (modes_buffer[(N*iCanEl)+k])/(iMPC);
2547 
2548  //assign a corresponding mode point count for this region into
2549  //the mode point counts array using the MPC buffer...
2550  modePointCounts[label] = MPC_buffer[iCanEl];
2551  }
2552  }
2553 
2554  //re-assign region count using label counter
2555  int oldRegionCount = regionCount;
2556  regionCount = label+1;
2557 
2558  // (c)
2559 
2560  // Use the label buffer to reconstruct the label map, which specified
2561  // the new image given its new regions calculated above
2562 
2563  for(i = 0; i < height*width; i++)
2564  labels[i] = label_buffer[raList[labels[i]].label];
2565 
2566  //de-allocate memory
2567  delete [] modes_buffer;
2568  delete [] MPC_buffer;
2569  delete [] label_buffer;
2570 
2571  //done.
2572  return;
2573 
2574 }
2575 
2576 /*******************************************************/
2577 /*Compute Edge Strengths */
2578 /*******************************************************/
2579 /*Computes the a weight for each link in the region */
2580 /*graph maintined by the RAM, resulting in a weighted */
2581 /*graph in which the weights consist of a confidence */
2582 /*between zero and one indicating if the regions are */
2583 /*separated by a strong or weak edge. */
2584 /*******************************************************/
2585 /*Post: */
2586 /* - an edge strength has been computed between */
2587 /* each region of the image and placed as a */
2588 /* weight in the RAM to be used during transi- */
2589 /* tive closure. */
2590 /*******************************************************/
2591 
2592 void msImageProcessor::ComputeEdgeStrengths( void )
2593 {
2594 
2595  //initialize visit table - used to keep track
2596  //of which pixels have already been visited such
2597  //as not to contribute their strength value to
2598  //a boundary sum multiple times...
2599  memset(visitTable, 0, L*sizeof(unsigned char));
2600 
2601  //traverse labeled image computing edge strengths
2602  //(excluding image boundary)...
2603  int x, y, dp, curLabel, rightLabel, bottomLabel;
2604  RAList *curRegion;
2605  for(y = 1; y < height-1; y++)
2606  {
2607  for(x = 1; x < width-1; x++)
2608  {
2609  //compute data point location using x and y
2610  dp = y*width + x;
2611 
2612  //obtain labels at different pixel locations
2613  curLabel = labels[dp ]; //current pixel
2614  rightLabel = labels[dp+1 ]; //right pixel
2615  bottomLabel = labels[dp+width]; //bottom pixel
2616 
2617  //check right and bottom neighbor to see if there is a
2618  //change in label then we are at an edge therefore record
2619  //the edge strength at this edge accumulating its value
2620  //in the RAM...
2621  if(curLabel != rightLabel)
2622  {
2623  //traverse into RAM...
2624  curRegion = &raList[curLabel];
2625  while((curRegion)&&(curRegion->label != rightLabel))
2626  curRegion = curRegion->next;
2627 
2628  //this should not occur...
2629  assert(curRegion);
2630 
2631  //accumulate edge strength
2632  curRegion->edgeStrength += weightMap[dp] + weightMap[dp+1];
2633  curRegion->edgePixelCount += 2;
2634  }
2635 
2636  if(curLabel != bottomLabel)
2637  {
2638  //traverse into RAM...
2639  curRegion = &raList[curLabel];
2640  while((curRegion)&&(curRegion->label != bottomLabel))
2641  curRegion = curRegion->next;
2642 
2643  //this should not occur...
2644  assert(curRegion);
2645 
2646  //accumulate edge strength
2647  if(curLabel == rightLabel)
2648  {
2649  curRegion->edgeStrength += weightMap[dp] + weightMap[dp+width];
2650  curRegion->edgePixelCount += 2;
2651  }
2652  else
2653  {
2654  curRegion->edgeStrength += weightMap[dp+width];
2655  curRegion->edgePixelCount += 1;
2656  }
2657 
2658  }
2659  }
2660  }
2661 
2662  //compute strengths using accumulated strengths obtained above...
2663  RAList *neighborRegion;
2664  float edgeStrength;
2665  int edgePixelCount;
2666  for(x = 0; x < regionCount; x++)
2667  {
2668  //traverse the region list of the current region
2669  curRegion = &raList[x];
2670  curRegion = curRegion->next;
2671  while(curRegion)
2672  {
2673  //with the assumption that regions having a smaller
2674  //label in the current region list have already
2675  //had their edge strengths computed, only compute
2676  //edge strengths for the regions whose label is greater
2677  //than x, the current region (region list) under
2678  //consideration...
2679  curLabel = curRegion->label;
2680  if(curLabel > x)
2681  {
2682  //obtain pointer to the element identifying the
2683  //current region in the neighbors region list...
2684  neighborRegion = &raList[curLabel];
2685  while((neighborRegion)&&(neighborRegion->label != x))
2686  neighborRegion = neighborRegion->next;
2687 
2688  //this should not occur...
2689  assert(neighborRegion);
2690 
2691  //compute edge strengths using accumulated confidence
2692  //value and pixel count
2693  if((edgePixelCount = curRegion->edgePixelCount + neighborRegion->edgePixelCount) != 0)
2694  {
2695  //compute edge strength
2696  edgeStrength = curRegion->edgeStrength + neighborRegion->edgeStrength;
2697  edgeStrength /= edgePixelCount;
2698 
2699  //store edge strength and pixel count for corresponding regions
2700  curRegion->edgeStrength = neighborRegion->edgeStrength = edgeStrength;
2701  curRegion->edgePixelCount = neighborRegion->edgePixelCount = edgePixelCount;
2702  }
2703  }
2704 
2705  //traverse to the next region in the region adjacency list
2706  //of the current region x
2707  curRegion = curRegion->next;
2708 
2709  }
2710  }
2711 
2712  //compute average edge strength amongst the edges connecting
2713  //it to each of its neighbors
2714  int numNeighbors;
2715  for(x = 0; x < regionCount; x++)
2716  {
2717  //traverse the region list of the current region
2718  //accumulating weights
2719  curRegion = &raList[x];
2720  curRegion = curRegion->next;
2721  edgeStrength = 0;
2722  numNeighbors = 0;
2723  while(curRegion)
2724  {
2725  numNeighbors++;
2726  edgeStrength += curRegion->edgeStrength;
2727  curRegion = curRegion->next;
2728  }
2729 
2730  //divide by the number of regions connected
2731  //to the current region
2732  if(numNeighbors) edgeStrength /= numNeighbors;
2733 
2734  //store the result in the raList for region
2735  //x
2736  raList[x].edgeStrength = edgeStrength;
2737  }
2738 
2739  //traverse raList and output the resulting list
2740  //to a file
2741 
2742  //done.
2743  return;
2744 
2745 }
2746 
2747 /*******************************************************/
2748 /*Prune */
2749 /*******************************************************/
2750 /*Prunes regions from the image whose pixel density */
2751 /*is less than a specified threshold. */
2752 /*******************************************************/
2753 /*Pre: */
2754 /* - minRegion is the minimum allowable pixel de- */
2755 /* nsity a region may have without being pruned */
2756 /* from the image */
2757 /*Post: */
2758 /* - regions whose pixel density is less than */
2759 /* or equal to minRegion have been pruned from */
2760 /* the image. */
2761 /*******************************************************/
2762 
2763 void msImageProcessor::Prune(int minRegion)
2764 {
2765 
2766  //Allocate Memory for temporary buffers...
2767 
2768  //allocate memory for mode and point count temporary buffers...
2769  float *modes_buffer = new float [N*regionCount];
2770  int *MPC_buffer = new int [regionCount];
2771 
2772  //allocate memory for label buffer
2773  int *label_buffer = new int [regionCount];
2774 
2775  //Declare variables
2776  int i, k, candidate, iCanEl, neighCanEl, iMPC, label, oldRegionCount, minRegionCount;
2777  double minSqDistance, neighborDistance;
2778  RAList *neighbor;
2779 
2780  //Apply pruning algorithm to classification structure, removing all regions whose area
2781  //is under the threshold area minRegion (pixels)
2782  do
2783  {
2784  //Assume that no region has area under threshold area of
2785  minRegionCount = 0;
2786 
2787  //Step (1):
2788 
2789  // Build RAM using classifiction structure originally
2790  // generated by the method GridTable::Connect()
2791  BuildRAM();
2792 
2793  // Step (2):
2794 
2795  // Traverse the RAM joining regions whose area is less than minRegion (pixels)
2796  // with its respective candidate region.
2797 
2798  // A candidate region is a region that displays the following properties:
2799 
2800  // - it is adjacent to the region being pruned
2801 
2802  // - the distance of its mode is a minimum to that of the region being pruned
2803  // such that or it is the only adjacent region having an area greater than
2804  // minRegion
2805 
2806  for(i = 0; i < regionCount; i++)
2807  {
2808  //if the area of the ith region is less than minRegion
2809  //join it with its candidate region...
2810 
2811  //*******************************************************************************
2812 
2813  //Note: Adjust this if statement if a more sophisticated pruning criterion
2814  // is desired. Basically in this step a region whose area is less than
2815  // minRegion is pruned by joining it with its "closest" neighbor (in color).
2816  // Therefore, by placing a different criterion for fusing a region the
2817  // pruning method may be altered to implement a more sophisticated algorithm.
2818 
2819  //*******************************************************************************
2820 
2821  if(modePointCounts[i] < minRegion)
2822  {
2823  //update minRegionCount to indicate that a region
2824  //having area less than minRegion was found
2825  minRegionCount++;
2826 
2827  //obtain a pointer to the first region in the
2828  //region adjacency list of the ith region...
2829  neighbor = raList[i].next;
2830 
2831  //calculate the distance between the mode of the ith
2832  //region and that of the neighboring region...
2833  candidate = neighbor->label;
2834  minSqDistance = SqDistance(i, candidate);
2835 
2836  //traverse region adjacency list of region i and select
2837  //a candidate region
2838  neighbor = neighbor->next;
2839  while(neighbor)
2840  {
2841 
2842  //calculate the square distance between region i
2843  //and current neighbor...
2844  neighborDistance = SqDistance(i, neighbor->label);
2845 
2846  //if this neighbors square distance to region i is less
2847  //than minSqDistance, then select this neighbor as the
2848  //candidate region for region i
2849  if(neighborDistance < minSqDistance)
2850  {
2851  minSqDistance = neighborDistance;
2852  candidate = neighbor->label;
2853  }
2854 
2855  //traverse region list of region i
2856  neighbor = neighbor->next;
2857 
2858  }
2859 
2860  //join region i with its candidate region:
2861 
2862  // (1) find the canonical element of region i
2863  iCanEl = i;
2864  while(raList[iCanEl].label != iCanEl)
2865  iCanEl = raList[iCanEl].label;
2866 
2867  // (2) find the canonical element of neighboring region
2868  neighCanEl = candidate;
2869  while(raList[neighCanEl].label != neighCanEl)
2870  neighCanEl = raList[neighCanEl].label;
2871 
2872  // if the canonical elements of are not the same then assign
2873  // the canonical element having the smaller label to be the parent
2874  // of the other region...
2875  if(iCanEl < neighCanEl)
2876  raList[neighCanEl].label = iCanEl;
2877  else
2878  {
2879  //must replace the canonical element of previous
2880  //parent as well
2881  raList[raList[iCanEl].label].label = neighCanEl;
2882 
2883  //re-assign canonical element
2884  raList[iCanEl].label = neighCanEl;
2885  }
2886  }
2887  }
2888 
2889  // Step (3):
2890 
2891  // Level binary trees formed by canonical elements
2892  for(i = 0; i < regionCount; i++)
2893  {
2894  iCanEl = i;
2895  while(raList[iCanEl].label != iCanEl)
2896  iCanEl = raList[iCanEl].label;
2897  raList[i].label = iCanEl;
2898  }
2899 
2900  // Step (4):
2901 
2902  //Traverse joint sets, relabeling image.
2903 
2904  // Accumulate modes and re-compute point counts using canonical
2905  // elements generated by step 2.
2906 
2907  //initialize buffers to zero
2908  for(i = 0; i < regionCount; i++)
2909  MPC_buffer[i] = 0;
2910  for(i = 0; i < N*regionCount; i++)
2911  modes_buffer[i] = 0;
2912 
2913  //traverse raList accumulating modes and point counts
2914  //using canoncial element information...
2915  for(i = 0; i < regionCount; i++)
2916  {
2917 
2918  //obtain canonical element of region i
2919  iCanEl = raList[i].label;
2920 
2921  //obtain mode point count of region i
2922  iMPC = modePointCounts[i];
2923 
2924  //accumulate modes_buffer[iCanEl]
2925  for(k = 0; k < N; k++)
2926  modes_buffer[(N*iCanEl)+k] += iMPC*modes[(N*i)+k];
2927 
2928  //accumulate MPC_buffer[iCanEl]
2929  MPC_buffer[iCanEl] += iMPC;
2930 
2931  }
2932 
2933  // (b)
2934 
2935  // Re-label new regions of the image using the canonical
2936  // element information generated by step (2)
2937 
2938  // Also use this information to compute the modes of the newly
2939  // defined regions, and to assign new region point counts in
2940  // a consecute manner to the modePointCounts array
2941 
2942  //initialize label buffer to -1
2943  for(i = 0; i < regionCount; i++)
2944  label_buffer[i] = -1;
2945 
2946  //traverse raList re-labeling the regions
2947  label = -1;
2948  for(i = 0; i < regionCount; i++)
2949  {
2950  //obtain canonical element of region i
2951  iCanEl = raList[i].label;
2952  if(label_buffer[iCanEl] < 0)
2953  {
2954  //assign a label to the new region indicated by canonical
2955  //element of i
2956  label_buffer[iCanEl] = ++label;
2957 
2958  //recompute mode storing the result in modes[label]...
2959  iMPC = MPC_buffer[iCanEl];
2960  for(k = 0; k < N; k++)
2961  modes[(N*label)+k] = (modes_buffer[(N*iCanEl)+k])/(iMPC);
2962 
2963  //assign a corresponding mode point count for this region into
2964  //the mode point counts array using the MPC buffer...
2965  modePointCounts[label] = MPC_buffer[iCanEl];
2966  }
2967  }
2968 
2969  //re-assign region count using label counter
2970  oldRegionCount = regionCount;
2971  regionCount = label+1;
2972 
2973  // (c)
2974 
2975  // Use the label buffer to reconstruct the label map, which specified
2976  // the new image given its new regions calculated above
2977 
2978  for(i = 0; i < height*width; i++)
2979  labels[i] = label_buffer[raList[labels[i]].label];
2980 
2981 
2982  } while(minRegionCount > 0);
2983 
2984  //de-allocate memory
2985  delete [] modes_buffer;
2986  delete [] MPC_buffer;
2987  delete [] label_buffer;
2988 
2989  //done.
2990  return;
2991 
2992 }
2993 
2994 /*******************************************************/
2995 /*Define Boundaries */
2996 /*******************************************************/
2997 /*Defines the boundaries for each region of the segm- */
2998 /*ented image storing the result into a region list */
2999 /*object. */
3000 /*******************************************************/
3001 /*Pre: */
3002 /* - the image has been segmented and a classifi- */
3003 /* cation structure has been created for this */
3004 /* image */
3005 /*Post: */
3006 /* - the boundaries of the segmented image have */
3007 /* been defined and the boundaries of each reg- */
3008 /* ion has been stored into a region list obj- */
3009 /* ect. */
3010 /*******************************************************/
3011 
3012 void msImageProcessor::DefineBoundaries( void )
3013 {
3014 
3015  //declare and allocate memory for boundary map and count
3016  int *boundaryMap, *boundaryCount;
3017  if((!(boundaryMap = new int [L]))||(!(boundaryCount = new int [regionCount])))
3018  ErrorHandler("msImageProcessor", "DefineBoundaries", "Not enough memory.");
3019 
3020  //initialize boundary map and count
3021  int i;
3022  for(i = 0; i < L; i++)
3023  boundaryMap[i] = -1;
3024  for(i = 0; i < regionCount; i++)
3025  boundaryCount[i] = 0;
3026 
3027  //initialize and declare total boundary count -
3028  //the total number of boundary pixels present in
3029  //the segmented image
3030  int totalBoundaryCount = 0;
3031 
3032  //traverse the image checking the right and bottom
3033  //four connected neighbors of each pixel marking
3034  //boundary map with the boundaries of each region and
3035  //incrementing boundaryCount using the label information
3036 
3037  //***********************************************************************
3038  //***********************************************************************
3039 
3040  int j, label, dataPoint;
3041 
3042  //first row (every pixel is a boundary pixel)
3043  for(i = 0; i < width; i++)
3044  {
3045  boundaryMap[i] = label = labels[i];
3046  boundaryCount[label]++;
3047  totalBoundaryCount++;
3048  }
3049 
3050  //define boundaries for all rows except for the first
3051  //and last one...
3052  for(i = 1; i < height - 1; i++)
3053  {
3054  //mark the first pixel in an image row as an image boundary...
3055  dataPoint = i*width;
3056  boundaryMap[dataPoint] = label = labels[dataPoint];
3057  boundaryCount[label]++;
3058  totalBoundaryCount++;
3059 
3060  for(j = 1; j < width - 1; j++)
3061  {
3062  //define datapoint and its right and bottom
3063  //four connected neighbors
3064  dataPoint = i*width+j;
3065 
3066  //check four connected neighbors if they are
3067  //different this pixel is a boundary pixel
3068  label = labels[dataPoint];
3069  if((label != labels[dataPoint-1]) ||(label != labels[dataPoint+1])||
3070  (label != labels[dataPoint-width])||(label != labels[dataPoint+width]))
3071  {
3072  boundaryMap[dataPoint] = label = labels[dataPoint];
3073  boundaryCount[label]++;
3074  totalBoundaryCount++;
3075  }
3076  }
3077 
3078  //mark the last pixel in an image row as an image boundary...
3079  dataPoint = (i+1)*width-1;
3080  boundaryMap[dataPoint] = label = labels[dataPoint];
3081  boundaryCount[label]++;
3082  totalBoundaryCount++;
3083 
3084  }
3085 
3086  //last row (every pixel is a boundary pixel) (i = height-1)
3087  register int start = (height-1)*width, stop = height*width;
3088  for(i = start; i < stop; i++)
3089  {
3090  boundaryMap[i] = label = labels[i];
3091  boundaryCount[label]++;
3092  totalBoundaryCount++;
3093  }
3094 
3095  //***********************************************************************
3096  //***********************************************************************
3097 
3098  //store boundary locations into a boundary buffer using
3099  //boundary map and count
3100 
3101  //***********************************************************************
3102  //***********************************************************************
3103 
3104  int *boundaryBuffer = new int [totalBoundaryCount], *boundaryIndex = new int [regionCount];
3105 
3106  //use boundary count to initialize boundary index...
3107  int counter = 0;
3108  for(i = 0; i < regionCount; i++)
3109  {
3110  boundaryIndex[i] = counter;
3111  counter += boundaryCount[i];
3112  }
3113 
3114  //traverse boundary map placing the boundary pixel
3115  //locations into the boundaryBuffer
3116  for(i = 0; i < L; i++)
3117  {
3118  //if its a boundary pixel store it into
3119  //the boundary buffer
3120  if((label = boundaryMap[i]) >= 0)
3121  {
3122  boundaryBuffer[boundaryIndex[label]] = i;
3123  boundaryIndex[label]++;
3124  }
3125  }
3126 
3127  //***********************************************************************
3128  //***********************************************************************
3129 
3130  //store the boundary locations stored by boundaryBuffer into
3131  //the region list for each region
3132 
3133  //***********************************************************************
3134  //***********************************************************************
3135 
3136  //destroy the old region list
3137  if(regionList) delete regionList;
3138 
3139  //create a new region list
3140  if(!(regionList = new RegionList(regionCount, totalBoundaryCount, N)))
3141  ErrorHandler("msImageProcessor", "DefineBoundaries", "Not enough memory.");
3142 
3143  //add boundary locations for each region using the boundary
3144  //buffer and boundary counts
3145  counter = 0;
3146  for(i = 0; i < regionCount; i++)
3147  {
3148  regionList->AddRegion(i, boundaryCount[i], &boundaryBuffer[counter]);
3149  counter += boundaryCount[i];
3150  }
3151 
3152  //***********************************************************************
3153  //***********************************************************************
3154 
3155  // dealocate local used memory
3156  delete [] boundaryMap;
3157  delete [] boundaryCount;
3158  delete [] boundaryBuffer;
3159  delete [] boundaryIndex;
3160 
3161  //done.
3162  return;
3163 
3164 }
3165 
3166  /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/
3167  /* Image Data Searching/Distance Calculation */
3168  /*\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/*/
3169 
3170 /*******************************************************/
3171 /*In Window */
3172 /*******************************************************/
3173 /*Returns true if the two specified data points are */
3174 /*within rR of each other. */
3175 /*******************************************************/
3176 /*Pre: */
3177 /* - mode1 and mode2 are indeces into msRawData */
3178 /* specifying the modes of the pixels having */
3179 /* these indeces. */
3180 /*Post: */
3181 /* - true is returned if mode1 and mode2 are wi- */
3182 /* thin rR of one another, false is returned */
3183 /* otherwise. */
3184 /*******************************************************/
3185 
3186 bool msImageProcessor::InWindow(int mode1, int mode2)
3187 {
3188  int k = 1, s = 0, p;
3189  double diff = 0, el;
3190  while((diff < 0.25)&&(k != kp)) // Partial Distortion Search
3191  {
3192  //Calculate distance squared of sub-space s
3193  diff = 0;
3194  for(p = 0; p < P[k]; p++)
3195  {
3196  el = (modes[mode1*N+p+s]-modes[mode2*N+p+s])/(h[k]*offset[k]);
3197  if((!p)&&(k == 1)&&(modes[mode1*N] > 80))
3198  diff += 4*el*el;
3199  else
3200  diff += el*el;
3201  }
3202 
3203  //next subspace
3204  s += P[k];
3205  k++;
3206  }
3207  return (bool)(diff < 0.25);
3208 }
3209 
3210 /*******************************************************/
3211 /*Square Distance */
3212 /*******************************************************/
3213 /*Computs the normalized square distance between two */
3214 /*modes. */
3215 /*******************************************************/
3216 /*Pre: */
3217 /* - mode1 and mode2 are indeces into the modes */
3218 /* array specifying two modes of the image */
3219 /*Post: */
3220 /* - the normalized square distance between modes */
3221 /* indexed by mode1 and mode2 has been calc- */
3222 /* ulated and the result has been returned. */
3223 /*******************************************************/
3224 
3225 float msImageProcessor::SqDistance(int mode1, int mode2)
3226 {
3227 
3228  int k = 1, s = 0, p;
3229  float dist = 0, el;
3230  for(k = 1; k < kp; k++)
3231  {
3232  //Calculate distance squared of sub-space s
3233  for(p = 0; p < P[k]; p++)
3234  {
3235  el = (modes[mode1*N+p+s]-modes[mode2*N+p+s])/(h[k]*offset[k]);
3236  dist += el*el;
3237  }
3238 
3239  //next subspace
3240  s += P[k];
3241  k++;
3242  }
3243 
3244  //return normalized square distance between modes
3245  //1 and 2
3246  return dist;
3247 
3248 }
3249 
3250  /*/\/\/\/\/\/\/\/\/\/\*/
3251  /* Memory Management */
3252  /*\/\/\/\/\/\/\/\/\/\/*/
3253 
3254 /*******************************************************/
3255 /*Initialize Output */
3256 /*******************************************************/
3257 /*Allocates memory needed by the mean shift image pro- */
3258 /*cessor class output storage data structure. */
3259 /*******************************************************/
3260 /*Post: */
3261 /* - the memory needed by the output storage */
3262 /* structure of this class has been (re-)allo- */
3263 /* cated. */
3264 /*******************************************************/
3265 
3266 void msImageProcessor::InitializeOutput( void )
3267 {
3268 
3269  //De-allocate memory if output was defined for previous image
3270  DestroyOutput();
3271 
3272  //Allocate memory for msRawData (filtered image output)
3273  if(!(msRawData = new float [L*N]))
3274  {
3275  ErrorHandler("msImageProcessor", "Allocate", "Not enough memory.");
3276  return;
3277  }
3278 
3279  //Allocate memory used to store image modes and their corresponding regions...
3280  if((!(modes = new float [L*(N+2)]))||(!(labels = new int [L]))||(!(modePointCounts = new int [L]))||(!(indexTable = new int [L])))
3281  {
3282  ErrorHandler("msImageProcessor", "Allocate", "Not enough memory");
3283  return;
3284  }
3285 
3286  //Allocate memory for integer modes used to perform connected components
3287  //(image labeling)...
3288 // if(!(LUV_data = new int [N*L]))
3289  if (!(LUV_data = new float[N*L]))
3290  {
3291  ErrorHandler("msImageProcessor", "Allocate", "Not enough memory");
3292  return;
3293  }
3294 
3295  //indicate that the class output storage structure has been defined
3296  class_state.OUTPUT_DEFINED = true;
3297 
3298 }
3299 
3300 /*******************************************************/
3301 /*Destroy Output */
3302 /*******************************************************/
3303 /*De-allocates memory needed by the mean shift image */
3304 /*processor class output storage data structure. */
3305 /*******************************************************/
3306 /*Post: */
3307 /* - the memory needed by the output storage */
3308 /* structure of this class has been de-alloc- */
3309 /* ated. */
3310 /* - the output storage structure has been init- */
3311 /* ialized for re-use. */
3312 /*******************************************************/
3313 
3314 void msImageProcessor::DestroyOutput( void )
3315 {
3316 
3317  //de-allocate memory for msRawData (filtered image output)
3318  if (msRawData) delete [] msRawData;
3319 
3320  //de-allocate memory used by output storage and image
3321  //classification structure
3322  if (modes) delete [] modes;
3323  if (labels) delete [] labels;
3324  if (modePointCounts) delete [] modePointCounts;
3325  if (indexTable) delete [] indexTable;
3326 
3327  //de-allocate memory for LUV_data
3328  if (LUV_data) delete [] LUV_data;
3329 
3330  //initialize data members for re-use...
3331 
3332  //initialize output structures...
3333  msRawData = NULL;
3334 
3335  //re-initialize classification structure
3336  modes = NULL;
3337  labels = NULL;
3338  modePointCounts = NULL;
3339  regionCount = 0;
3340 
3341  //indicate that the output has been destroyed
3342  class_state.OUTPUT_DEFINED = false;
3343 
3344  //done.
3345  return;
3346 
3347 }
3348 
3349 // NEW
3350 void msImageProcessor::NewOptimizedFilter1(float sigmaS, float sigmaR)
3351 {
3352  // Declare Variables
3353  int iterationCount, i, j, k, modeCandidateX, modeCandidateY, modeCandidate_i;
3354  double mvAbs, diff, el;
3355 
3356  //make sure that a lattice height and width have
3357  //been defined...
3358  if(!height)
3359  {
3360  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
3361  return;
3362  }
3363 
3364  //re-assign bandwidths to sigmaS and sigmaR
3365  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
3366  {
3367  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
3368  return;
3369  }
3370 
3371  //define input data dimension with lattice
3372  int lN = N + 2;
3373 
3374  // Traverse each data point applying mean shift
3375  // to each data point
3376 
3377  // Allcocate memory for yk
3378  double *yk = new double [lN];
3379 
3380  // Allocate memory for Mh
3381  double *Mh = new double [lN];
3382 
3383  // let's use some temporary data
3384  float* sdata;
3385  sdata = new float[lN*L];
3386 
3387  // copy the scaled data
3388  int idxs, idxd;
3389  idxs = idxd = 0;
3390  if (N==3)
3391  {
3392  for(i=0; i<L; i++)
3393  {
3394  sdata[idxs++] = (i%width)/sigmaS;
3395  sdata[idxs++] = (i/width)/sigmaS;
3396  sdata[idxs++] = data[idxd++]/sigmaR;
3397  sdata[idxs++] = data[idxd++]/sigmaR;
3398  sdata[idxs++] = data[idxd++]/sigmaR;
3399  }
3400  } else if (N==1)
3401  {
3402  for(i=0; i<L; i++)
3403  {
3404  sdata[idxs++] = (i%width)/sigmaS;
3405  sdata[idxs++] = (i/width)/sigmaS;
3406  sdata[idxs++] = data[idxd++]/sigmaR;
3407  }
3408  } else
3409  {
3410  for(i=0; i<L; i++)
3411  {
3412  sdata[idxs++] = (i%width)/sigmaS;
3413  sdata[idxs++] = (i/width)/sigmaS;
3414  for (j=0; j<N; j++)
3415  sdata[idxs++] = data[idxd++]/sigmaR;
3416  }
3417  }
3418  // index the data in the 3d buckets (x, y, L)
3419  int* buckets;
3420  int* slist;
3421  slist = new int[L];
3422  int bucNeigh[27];
3423 
3424  float sMins; // just for L
3425  float sMaxs[3]; // for all
3426  sMaxs[0] = width/sigmaS;
3427  sMaxs[1] = height/sigmaS;
3428  sMins = sMaxs[2] = sdata[2];
3429  idxs = 2;
3430  float cval;
3431  for(i=0; i<L; i++)
3432  {
3433  cval = sdata[idxs];
3434  if (cval < sMins)
3435  sMins = cval;
3436  else if (cval > sMaxs[2])
3437  sMaxs[2] = cval;
3438 
3439  idxs += lN;
3440  }
3441 
3442  int nBuck1, nBuck2, nBuck3;
3443  int cBuck1, cBuck2, cBuck3, cBuck;
3444  nBuck1 = (int) (sMaxs[0] + 3);
3445  nBuck2 = (int) (sMaxs[1] + 3);
3446  nBuck3 = (int) (sMaxs[2] - sMins + 3);
3447  buckets = new int[nBuck1*nBuck2*nBuck3];
3448  for(i=0; i<(nBuck1*nBuck2*nBuck3); i++)
3449  buckets[i] = -1;
3450 
3451  idxs = 0;
3452  for(i=0; i<L; i++)
3453  {
3454  // find bucket for current data and add it to the list
3455  cBuck1 = (int) sdata[idxs] + 1;
3456  cBuck2 = (int) sdata[idxs+1] + 1;
3457  cBuck3 = (int) (sdata[idxs+2] - sMins) + 1;
3458  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3459 
3460  slist[i] = buckets[cBuck];
3461  buckets[cBuck] = i;
3462 
3463  idxs += lN;
3464  }
3465  // init bucNeigh
3466  idxd = 0;
3467  for (cBuck1=-1; cBuck1<=1; cBuck1++)
3468  {
3469  for (cBuck2=-1; cBuck2<=1; cBuck2++)
3470  {
3471  for (cBuck3=-1; cBuck3<=1; cBuck3++)
3472  {
3473  bucNeigh[idxd++] = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3474  }
3475  }
3476  }
3477  double wsuml, weight;
3478  double hiLTr = 80.0/sigmaR;
3479  // done indexing/hashing
3480 
3481 
3482  // Initialize mode table used for basin of attraction
3483  memset(modeTable, 0, width*height);
3484 
3485  // proceed ...
3486 #ifdef PROMPT
3487  msSys.Prompt("done.\nApplying mean shift (Using Lattice) ... ");
3488 #ifdef SHOW_PROGRESS
3489  msSys.Prompt("\n 0%%");
3490 #endif
3491 #endif
3492 
3493 
3494  for(i = 0; i < L; i++)
3495  {
3496  // if a mode was already assigned to this data point
3497  // then skip this point, otherwise proceed to
3498  // find its mode by applying mean shift...
3499  if (modeTable[i] == 1)
3500  continue;
3501 
3502  // initialize point list...
3503  pointCount = 0;
3504 
3505  // Assign window center (window centers are
3506  // initialized by createLattice to be the point
3507  // data[i])
3508  idxs = i*lN;
3509  for (j=0; j<lN; j++)
3510  yk[j] = sdata[idxs+j];
3511 
3512  // Calculate the mean shift vector using the lattice
3513  // LatticeMSVector(Mh, yk); // modify to new
3514  /*****************************************************/
3515  // Initialize mean shift vector
3516  for(j = 0; j < lN; j++)
3517  Mh[j] = 0;
3518  wsuml = 0;
3519  // uniformLSearch(Mh, yk_ptr); // modify to new
3520  // find bucket of yk
3521  cBuck1 = (int) yk[0] + 1;
3522  cBuck2 = (int) yk[1] + 1;
3523  cBuck3 = (int) (yk[2] - sMins) + 1;
3524  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3525  for (j=0; j<27; j++)
3526  {
3527  idxd = buckets[cBuck+bucNeigh[j]];
3528  // list parse, crt point is cHeadList
3529  while (idxd>=0)
3530  {
3531  idxs = lN*idxd;
3532  // determine if inside search window
3533  el = sdata[idxs+0]-yk[0];
3534  diff = el*el;
3535  el = sdata[idxs+1]-yk[1];
3536  diff += el*el;
3537 
3538  if (diff < 1.0)
3539  {
3540  el = sdata[idxs+2]-yk[2];
3541  if (yk[2] > hiLTr)
3542  diff = 4*el*el;
3543  else
3544  diff = el*el;
3545 
3546  if (N>1)
3547  {
3548  el = sdata[idxs+3]-yk[3];
3549  diff += el*el;
3550  el = sdata[idxs+4]-yk[4];
3551  diff += el*el;
3552  }
3553 
3554  if (diff < 1.0)
3555  {
3556  weight = 1-weightMap[idxd];
3557  for (k=0; k<lN; k++)
3558  Mh[k] += weight*sdata[idxs+k];
3559  wsuml += weight;
3560  }
3561  }
3562  idxd = slist[idxd];
3563  }
3564  }
3565  if (wsuml > 0)
3566  {
3567  for(j = 0; j < lN; j++)
3568  Mh[j] = Mh[j]/wsuml - yk[j];
3569  }
3570  else
3571  {
3572  for(j = 0; j < lN; j++)
3573  Mh[j] = 0;
3574  }
3575  /*****************************************************/
3576  // Calculate its magnitude squared
3577  //mvAbs = 0;
3578  //for(j = 0; j < lN; j++)
3579  // mvAbs += Mh[j]*Mh[j];
3580  mvAbs = (Mh[0]*Mh[0]+Mh[1]*Mh[1])*sigmaS*sigmaS;
3581  if (N==3)
3582  mvAbs += (Mh[2]*Mh[2]+Mh[3]*Mh[3]+Mh[4]*Mh[4])*sigmaR*sigmaR;
3583  else
3584  mvAbs += Mh[2]*Mh[2]*sigmaR*sigmaR;
3585 
3586 
3587  // Keep shifting window center until the magnitude squared of the
3588  // mean shift vector calculated at the window center location is
3589  // under a specified threshold (Epsilon)
3590 
3591  // NOTE: iteration count is for speed up purposes only - it
3592  // does not have any theoretical importance
3593  iterationCount = 1;
3594  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
3595  {
3596 
3597  // Shift window location
3598  for(j = 0; j < lN; j++)
3599  yk[j] += Mh[j];
3600 
3601  // check to see if the current mode location is in the
3602  // basin of attraction...
3603 
3604  // calculate the location of yk on the lattice
3605  modeCandidateX = (int) (sigmaS*yk[0]+0.5);
3606  modeCandidateY = (int) (sigmaS*yk[1]+0.5);
3607  modeCandidate_i = modeCandidateY*width + modeCandidateX;
3608 
3609  // if mvAbs != 0 (yk did indeed move) then check
3610  // location basin_i in the mode table to see if
3611  // this data point either:
3612 
3613  // (1) has not been associated with a mode yet
3614  // (modeTable[basin_i] = 0), so associate
3615  // it with this one
3616  //
3617  // (2) it has been associated with a mode other
3618  // than the one that this data point is converging
3619  // to (modeTable[basin_i] = 1), so assign to
3620  // this data point the same mode as that of basin_i
3621 
3622  if ((modeTable[modeCandidate_i] != 2) && (modeCandidate_i != i))
3623  {
3624  // obtain the data point at basin_i to
3625  // see if it is within h*TC_DIST_FACTOR of
3626  // of yk
3627  diff = 0;
3628  idxs = lN*modeCandidate_i;
3629  for (k=2; k<lN; k++)
3630  {
3631  el = sdata[idxs+k] - yk[k];
3632  diff += el*el;
3633  }
3634 
3635  // if the data point at basin_i is within
3636  // a distance of h*TC_DIST_FACTOR of yk
3637  // then depending on modeTable[basin_i] perform
3638  // either (1) or (2)
3639  if (diff < TC_DIST_FACTOR)
3640  {
3641  // if the data point at basin_i has not
3642  // been associated to a mode then associate
3643  // it with the mode that this one will converge
3644  // to
3645  if (modeTable[modeCandidate_i] == 0)
3646  {
3647  // no mode associated yet so associate
3648  // it with this one...
3649  pointList[pointCount++] = modeCandidate_i;
3650  modeTable[modeCandidate_i] = 2;
3651 
3652  } else
3653  {
3654 
3655  // the mode has already been associated with
3656  // another mode, thererfore associate this one
3657  // mode and the modes in the point list with
3658  // the mode associated with data[basin_i]...
3659 
3660  // store the mode info into yk using msRawData...
3661  for (j = 0; j < N; j++)
3662  yk[j+2] = msRawData[modeCandidate_i*N+j]/sigmaR;
3663 
3664  // update mode table for this data point
3665  // indicating that a mode has been associated
3666  // with it
3667  modeTable[i] = 1;
3668 
3669  // indicate that a mode has been associated
3670  // to this data point (data[i])
3671  mvAbs = -1;
3672 
3673  // stop mean shift calculation...
3674  break;
3675  }
3676  }
3677  }
3678 
3679  // Calculate the mean shift vector at the new
3680  // window location using lattice
3681  // Calculate the mean shift vector using the lattice
3682  // LatticeMSVector(Mh, yk); // modify to new
3683  /*****************************************************/
3684  // Initialize mean shift vector
3685  for(j = 0; j < lN; j++)
3686  Mh[j] = 0;
3687  wsuml = 0;
3688  // uniformLSearch(Mh, yk_ptr); // modify to new
3689  // find bucket of yk
3690  cBuck1 = (int) yk[0] + 1;
3691  cBuck2 = (int) yk[1] + 1;
3692  cBuck3 = (int) (yk[2] - sMins) + 1;
3693  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3694  for (j=0; j<27; j++)
3695  {
3696  idxd = buckets[cBuck+bucNeigh[j]];
3697  // list parse, crt point is cHeadList
3698  while (idxd>=0)
3699  {
3700  idxs = lN*idxd;
3701  // determine if inside search window
3702  el = sdata[idxs+0]-yk[0];
3703  diff = el*el;
3704  el = sdata[idxs+1]-yk[1];
3705  diff += el*el;
3706 
3707  if (diff < 1.0)
3708  {
3709  el = sdata[idxs+2]-yk[2];
3710  if (yk[2] > hiLTr)
3711  diff = 4*el*el;
3712  else
3713  diff = el*el;
3714 
3715  if (N>1)
3716  {
3717  el = sdata[idxs+3]-yk[3];
3718  diff += el*el;
3719  el = sdata[idxs+4]-yk[4];
3720  diff += el*el;
3721  }
3722 
3723  if (diff < 1.0)
3724  {
3725  weight = 1-weightMap[idxd];
3726  for (k=0; k<lN; k++)
3727  Mh[k] += weight*sdata[idxs+k];
3728  wsuml += weight;
3729  }
3730  }
3731  idxd = slist[idxd];
3732  }
3733  }
3734  if (wsuml > 0)
3735  {
3736  for(j = 0; j < lN; j++)
3737  Mh[j] = Mh[j]/wsuml - yk[j];
3738  }
3739  else
3740  {
3741  for(j = 0; j < lN; j++)
3742  Mh[j] = 0;
3743  }
3744  /*****************************************************/
3745 
3746  // Calculate its magnitude squared
3747  //mvAbs = 0;
3748  //for(j = 0; j < lN; j++)
3749  // mvAbs += Mh[j]*Mh[j];
3750  mvAbs = (Mh[0]*Mh[0]+Mh[1]*Mh[1])*sigmaS*sigmaS;
3751  if (N==3)
3752  mvAbs += (Mh[2]*Mh[2]+Mh[3]*Mh[3]+Mh[4]*Mh[4])*sigmaR*sigmaR;
3753  else
3754  mvAbs += Mh[2]*Mh[2]*sigmaR*sigmaR;
3755 
3756  // Increment iteration count
3757  iterationCount++;
3758 
3759  }
3760 
3761  // if a mode was not associated with this data point
3762  // yet associate it with yk...
3763  if (mvAbs >= 0)
3764  {
3765  // Shift window location
3766  for(j = 0; j < lN; j++)
3767  yk[j] += Mh[j];
3768 
3769  // update mode table for this data point
3770  // indicating that a mode has been associated
3771  // with it
3772  modeTable[i] = 1;
3773 
3774  }
3775 
3776  for (k=0; k<N; k++)
3777  yk[k+2] *= sigmaR;
3778 
3779  // associate the data point indexed by
3780  // the point list with the mode stored
3781  // by yk
3782  for (j = 0; j < pointCount; j++)
3783  {
3784  // obtain the point location from the
3785  // point list
3786  modeCandidate_i = pointList[j];
3787 
3788  // update the mode table for this point
3789  modeTable[modeCandidate_i] = 1;
3790 
3791  //store result into msRawData...
3792  for(k = 0; k < N; k++)
3793  msRawData[N*modeCandidate_i+k] = (float)(yk[k+2]);
3794  }
3795 
3796  //store result into msRawData...
3797  for(j = 0; j < N; j++)
3798  msRawData[N*i+j] = (float)(yk[j+2]);
3799 
3800  // Prompt user on progress
3801 #ifdef SHOW_PROGRESS
3802  percent_complete = (float)(i/(float)(L))*100;
3803  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
3804 #endif
3805 
3806  // Check to see if the algorithm has been halted
3807  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
3808  break;
3809  }
3810 
3811  // Prompt user that filtering is completed
3812 #ifdef PROMPT
3813 #ifdef SHOW_PROGRESS
3814  msSys.Prompt("\r");
3815 #endif
3816  msSys.Prompt("done.");
3817 #endif
3818  // de-allocate memory
3819  delete [] buckets;
3820  delete [] slist;
3821  delete [] sdata;
3822 
3823  delete [] yk;
3824  delete [] Mh;
3825 
3826  // done.
3827  return;
3828 
3829 }
3830 
3831 // NEW
3832 void msImageProcessor::NewOptimizedFilter2(float sigmaS, float sigmaR)
3833 {
3834  // Declare Variables
3835  int iterationCount, i, j, k, modeCandidateX, modeCandidateY, modeCandidate_i;
3836  double mvAbs, diff, el;
3837 
3838  //make sure that a lattice height and width have
3839  //been defined...
3840  if(!height)
3841  {
3842  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
3843  return;
3844  }
3845 
3846  //re-assign bandwidths to sigmaS and sigmaR
3847  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
3848  {
3849  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
3850  return;
3851  }
3852 
3853  //define input data dimension with lattice
3854  int lN = N + 2;
3855 
3856  // Traverse each data point applying mean shift
3857  // to each data point
3858 
3859  // Allcocate memory for yk
3860  double *yk = new double [lN];
3861 
3862  // Allocate memory for Mh
3863  double *Mh = new double [lN];
3864 
3865  // let's use some temporary data
3866  float* sdata;
3867  sdata = new float[lN*L];
3868 
3869  // copy the scaled data
3870  int idxs, idxd;
3871  idxs = idxd = 0;
3872  if (N==3)
3873  {
3874  for(i=0; i<L; i++)
3875  {
3876  sdata[idxs++] = (i%width)/sigmaS;
3877  sdata[idxs++] = (i/width)/sigmaS;
3878  sdata[idxs++] = data[idxd++]/sigmaR;
3879  sdata[idxs++] = data[idxd++]/sigmaR;
3880  sdata[idxs++] = data[idxd++]/sigmaR;
3881  }
3882  } else if (N==1)
3883  {
3884  for(i=0; i<L; i++)
3885  {
3886  sdata[idxs++] = (i%width)/sigmaS;
3887  sdata[idxs++] = (i/width)/sigmaS;
3888  sdata[idxs++] = data[idxd++]/sigmaR;
3889  }
3890  } else
3891  {
3892  for(i=0; i<L; i++)
3893  {
3894  sdata[idxs++] = (i%width)/sigmaS;
3895  sdata[idxs++] = (i/width)/sigmaS;
3896  for (j=0; j<N; j++)
3897  sdata[idxs++] = data[idxd++]/sigmaR;
3898  }
3899  }
3900  // index the data in the 3d buckets (x, y, L)
3901  int* buckets;
3902  int* slist;
3903  slist = new int[L];
3904  int bucNeigh[27];
3905 
3906  float sMins; // just for L
3907  float sMaxs[3]; // for all
3908  sMaxs[0] = width/sigmaS;
3909  sMaxs[1] = height/sigmaS;
3910  sMins = sMaxs[2] = sdata[2];
3911  idxs = 2;
3912  float cval;
3913  for(i=0; i<L; i++)
3914  {
3915  cval = sdata[idxs];
3916  if (cval < sMins)
3917  sMins = cval;
3918  else if (cval > sMaxs[2])
3919  sMaxs[2] = cval;
3920 
3921  idxs += lN;
3922  }
3923 
3924  int nBuck1, nBuck2, nBuck3;
3925  int cBuck1, cBuck2, cBuck3, cBuck;
3926  nBuck1 = (int) (sMaxs[0] + 3);
3927  nBuck2 = (int) (sMaxs[1] + 3);
3928  nBuck3 = (int) (sMaxs[2] - sMins + 3);
3929  buckets = new int[nBuck1*nBuck2*nBuck3];
3930  for(i=0; i<(nBuck1*nBuck2*nBuck3); i++)
3931  buckets[i] = -1;
3932 
3933  idxs = 0;
3934  for(i=0; i<L; i++)
3935  {
3936  // find bucket for current data and add it to the list
3937  cBuck1 = (int) sdata[idxs] + 1;
3938  cBuck2 = (int) sdata[idxs+1] + 1;
3939  cBuck3 = (int) (sdata[idxs+2] - sMins) + 1;
3940  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3941 
3942  slist[i] = buckets[cBuck];
3943  buckets[cBuck] = i;
3944 
3945  idxs += lN;
3946  }
3947  // init bucNeigh
3948  idxd = 0;
3949  for (cBuck1=-1; cBuck1<=1; cBuck1++)
3950  {
3951  for (cBuck2=-1; cBuck2<=1; cBuck2++)
3952  {
3953  for (cBuck3=-1; cBuck3<=1; cBuck3++)
3954  {
3955  bucNeigh[idxd++] = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
3956  }
3957  }
3958  }
3959  double wsuml, weight;
3960  double hiLTr = 80.0/sigmaR;
3961  // done indexing/hashing
3962 
3963 
3964  // Initialize mode table used for basin of attraction
3965  memset(modeTable, 0, width*height);
3966 
3967  // proceed ...
3968 #ifdef PROMPT
3969  msSys.Prompt("done.\nApplying mean shift (Using Lattice) ... ");
3970 #ifdef SHOW_PROGRESS
3971  msSys.Prompt("\n 0%%");
3972 #endif
3973 #endif
3974 
3975 
3976  for(i = 0; i < L; i++)
3977  {
3978  // if a mode was already assigned to this data point
3979  // then skip this point, otherwise proceed to
3980  // find its mode by applying mean shift...
3981  if (modeTable[i] == 1)
3982  continue;
3983 
3984  // initialize point list...
3985  pointCount = 0;
3986 
3987  // Assign window center (window centers are
3988  // initialized by createLattice to be the point
3989  // data[i])
3990  idxs = i*lN;
3991  for (j=0; j<lN; j++)
3992  yk[j] = sdata[idxs+j];
3993 
3994  // Calculate the mean shift vector using the lattice
3995  // LatticeMSVector(Mh, yk); // modify to new
3996  /*****************************************************/
3997  // Initialize mean shift vector
3998  for(j = 0; j < lN; j++)
3999  Mh[j] = 0;
4000  wsuml = 0;
4001  // uniformLSearch(Mh, yk_ptr); // modify to new
4002  // find bucket of yk
4003  cBuck1 = (int) yk[0] + 1;
4004  cBuck2 = (int) yk[1] + 1;
4005  cBuck3 = (int) (yk[2] - sMins) + 1;
4006  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4007  for (j=0; j<27; j++)
4008  {
4009  idxd = buckets[cBuck+bucNeigh[j]];
4010  // list parse, crt point is cHeadList
4011  while (idxd>=0)
4012  {
4013  idxs = lN*idxd;
4014  // determine if inside search window
4015  el = sdata[idxs+0]-yk[0];
4016  diff = el*el;
4017  el = sdata[idxs+1]-yk[1];
4018  diff += el*el;
4019 
4020  if (diff < 1.0)
4021  {
4022  el = sdata[idxs+2]-yk[2];
4023  if (yk[2] > hiLTr)
4024  diff = 4*el*el;
4025  else
4026  diff = el*el;
4027 
4028  if (N>1)
4029  {
4030  el = sdata[idxs+3]-yk[3];
4031  diff += el*el;
4032  el = sdata[idxs+4]-yk[4];
4033  diff += el*el;
4034  }
4035 
4036  if (diff < 1.0)
4037  {
4038  weight = 1-weightMap[idxd];
4039  for (k=0; k<lN; k++)
4040  Mh[k] += weight*sdata[idxs+k];
4041  wsuml += weight;
4042 
4043  //set basin of attraction mode table
4044  if (diff < speedThreshold)
4045  {
4046  if(modeTable[idxd] == 0)
4047  {
4048  pointList[pointCount++] = idxd;
4049  modeTable[idxd] = 2;
4050  }
4051  }
4052  }
4053  }
4054  idxd = slist[idxd];
4055  }
4056  }
4057  if (wsuml > 0)
4058  {
4059  for(j = 0; j < lN; j++)
4060  Mh[j] = Mh[j]/wsuml - yk[j];
4061  }
4062  else
4063  {
4064  for(j = 0; j < lN; j++)
4065  Mh[j] = 0;
4066  }
4067  /*****************************************************/
4068  // Calculate its magnitude squared
4069  //mvAbs = 0;
4070  //for(j = 0; j < lN; j++)
4071  // mvAbs += Mh[j]*Mh[j];
4072  mvAbs = (Mh[0]*Mh[0]+Mh[1]*Mh[1])*sigmaS*sigmaS;
4073  if (N==3)
4074  mvAbs += (Mh[2]*Mh[2]+Mh[3]*Mh[3]+Mh[4]*Mh[4])*sigmaR*sigmaR;
4075  else
4076  mvAbs += Mh[2]*Mh[2]*sigmaR*sigmaR;
4077 
4078 
4079  // Keep shifting window center until the magnitude squared of the
4080  // mean shift vector calculated at the window center location is
4081  // under a specified threshold (Epsilon)
4082 
4083  // NOTE: iteration count is for speed up purposes only - it
4084  // does not have any theoretical importance
4085  iterationCount = 1;
4086  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
4087  {
4088 
4089  // Shift window location
4090  for(j = 0; j < lN; j++)
4091  yk[j] += Mh[j];
4092 
4093  // check to see if the current mode location is in the
4094  // basin of attraction...
4095 
4096  // calculate the location of yk on the lattice
4097  modeCandidateX = (int) (sigmaS*yk[0]+0.5);
4098  modeCandidateY = (int) (sigmaS*yk[1]+0.5);
4099  modeCandidate_i = modeCandidateY*width + modeCandidateX;
4100 
4101  // if mvAbs != 0 (yk did indeed move) then check
4102  // location basin_i in the mode table to see if
4103  // this data point either:
4104 
4105  // (1) has not been associated with a mode yet
4106  // (modeTable[basin_i] = 0), so associate
4107  // it with this one
4108  //
4109  // (2) it has been associated with a mode other
4110  // than the one that this data point is converging
4111  // to (modeTable[basin_i] = 1), so assign to
4112  // this data point the same mode as that of basin_i
4113 
4114  if ((modeTable[modeCandidate_i] != 2) && (modeCandidate_i != i))
4115  {
4116  // obtain the data point at basin_i to
4117  // see if it is within h*TC_DIST_FACTOR of
4118  // of yk
4119  diff = 0;
4120  idxs = lN*modeCandidate_i;
4121  for (k=2; k<lN; k++)
4122  {
4123  el = sdata[idxs+k] - yk[k];
4124  diff += el*el;
4125  }
4126 
4127  // if the data point at basin_i is within
4128  // a distance of h*TC_DIST_FACTOR of yk
4129  // then depending on modeTable[basin_i] perform
4130  // either (1) or (2)
4131  if (diff < speedThreshold)
4132  {
4133  // if the data point at basin_i has not
4134  // been associated to a mode then associate
4135  // it with the mode that this one will converge
4136  // to
4137  if (modeTable[modeCandidate_i] == 0)
4138  {
4139  // no mode associated yet so associate
4140  // it with this one...
4141  pointList[pointCount++] = modeCandidate_i;
4142  modeTable[modeCandidate_i] = 2;
4143 
4144  } else
4145  {
4146 
4147  // the mode has already been associated with
4148  // another mode, thererfore associate this one
4149  // mode and the modes in the point list with
4150  // the mode associated with data[basin_i]...
4151 
4152  // store the mode info into yk using msRawData...
4153  for (j = 0; j < N; j++)
4154  yk[j+2] = msRawData[modeCandidate_i*N+j]/sigmaR;
4155 
4156  // update mode table for this data point
4157  // indicating that a mode has been associated
4158  // with it
4159  modeTable[i] = 1;
4160 
4161  // indicate that a mode has been associated
4162  // to this data point (data[i])
4163  mvAbs = -1;
4164 
4165  // stop mean shift calculation...
4166  break;
4167  }
4168  }
4169  }
4170 
4171  // Calculate the mean shift vector at the new
4172  // window location using lattice
4173  // Calculate the mean shift vector using the lattice
4174  // LatticeMSVector(Mh, yk); // modify to new
4175  /*****************************************************/
4176  // Initialize mean shift vector
4177  for(j = 0; j < lN; j++)
4178  Mh[j] = 0;
4179  wsuml = 0;
4180  // uniformLSearch(Mh, yk_ptr); // modify to new
4181  // find bucket of yk
4182  cBuck1 = (int) yk[0] + 1;
4183  cBuck2 = (int) yk[1] + 1;
4184  cBuck3 = (int) (yk[2] - sMins) + 1;
4185  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4186  for (j=0; j<27; j++)
4187  {
4188  idxd = buckets[cBuck+bucNeigh[j]];
4189  // list parse, crt point is cHeadList
4190  while (idxd>=0)
4191  {
4192  idxs = lN*idxd;
4193  // determine if inside search window
4194  el = sdata[idxs+0]-yk[0];
4195  diff = el*el;
4196  el = sdata[idxs+1]-yk[1];
4197  diff += el*el;
4198 
4199  if (diff < 1.0)
4200  {
4201  el = sdata[idxs+2]-yk[2];
4202  if (yk[2] > hiLTr)
4203  diff = 4*el*el;
4204  else
4205  diff = el*el;
4206 
4207  if (N>1)
4208  {
4209  el = sdata[idxs+3]-yk[3];
4210  diff += el*el;
4211  el = sdata[idxs+4]-yk[4];
4212  diff += el*el;
4213  }
4214 
4215  if (diff < 1.0)
4216  {
4217  weight = 1-weightMap[idxd];
4218  for (k=0; k<lN; k++)
4219  Mh[k] += weight*sdata[idxs+k];
4220  wsuml += weight;
4221 
4222  //set basin of attraction mode table
4223  if (diff < speedThreshold)
4224  {
4225  if(modeTable[idxd] == 0)
4226  {
4227  pointList[pointCount++] = idxd;
4228  modeTable[idxd] = 2;
4229  }
4230  }
4231 
4232  }
4233  }
4234  idxd = slist[idxd];
4235  }
4236  }
4237  if (wsuml > 0)
4238  {
4239  for(j = 0; j < lN; j++)
4240  Mh[j] = Mh[j]/wsuml - yk[j];
4241  }
4242  else
4243  {
4244  for(j = 0; j < lN; j++)
4245  Mh[j] = 0;
4246  }
4247  /*****************************************************/
4248 
4249  // Calculate its magnitude squared
4250  //mvAbs = 0;
4251  //for(j = 0; j < lN; j++)
4252  // mvAbs += Mh[j]*Mh[j];
4253  mvAbs = (Mh[0]*Mh[0]+Mh[1]*Mh[1])*sigmaS*sigmaS;
4254  if (N==3)
4255  mvAbs += (Mh[2]*Mh[2]+Mh[3]*Mh[3]+Mh[4]*Mh[4])*sigmaR*sigmaR;
4256  else
4257  mvAbs += Mh[2]*Mh[2]*sigmaR*sigmaR;
4258 
4259  // Increment iteration count
4260  iterationCount++;
4261 
4262  }
4263 
4264  // if a mode was not associated with this data point
4265  // yet associate it with yk...
4266  if (mvAbs >= 0)
4267  {
4268  // Shift window location
4269  for(j = 0; j < lN; j++)
4270  yk[j] += Mh[j];
4271 
4272  // update mode table for this data point
4273  // indicating that a mode has been associated
4274  // with it
4275  modeTable[i] = 1;
4276 
4277  }
4278 
4279  for (k=0; k<N; k++)
4280  yk[k+2] *= sigmaR;
4281 
4282  // associate the data point indexed by
4283  // the point list with the mode stored
4284  // by yk
4285  for (j = 0; j < pointCount; j++)
4286  {
4287  // obtain the point location from the
4288  // point list
4289  modeCandidate_i = pointList[j];
4290 
4291  // update the mode table for this point
4292  modeTable[modeCandidate_i] = 1;
4293 
4294  //store result into msRawData...
4295  for(k = 0; k < N; k++)
4296  msRawData[N*modeCandidate_i+k] = (float)(yk[k+2]);
4297  }
4298 
4299  //store result into msRawData...
4300  for(j = 0; j < N; j++)
4301  msRawData[N*i+j] = (float)(yk[j+2]);
4302 
4303  // Prompt user on progress
4304 #ifdef SHOW_PROGRESS
4305  percent_complete = (float)(i/(float)(L))*100;
4306  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
4307 #endif
4308 
4309  // Check to see if the algorithm has been halted
4310  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
4311  break;
4312  }
4313 
4314  // Prompt user that filtering is completed
4315 #ifdef PROMPT
4316 #ifdef SHOW_PROGRESS
4317  msSys.Prompt("\r");
4318 #endif
4319  msSys.Prompt("done.");
4320 #endif
4321  // de-allocate memory
4322  delete [] buckets;
4323  delete [] slist;
4324  delete [] sdata;
4325 
4326  delete [] yk;
4327  delete [] Mh;
4328 
4329  // done.
4330  return;
4331 
4332 }
4333 
4334 void msImageProcessor::NewNonOptimizedFilter(float sigmaS, float sigmaR)
4335 {
4336 
4337  // Declare Variables
4338  int iterationCount, i, j, k;
4339  double mvAbs, diff, el;
4340 
4341  //make sure that a lattice height and width have
4342  //been defined...
4343  if(!height)
4344  {
4345  ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
4346  return;
4347  }
4348 
4349  //re-assign bandwidths to sigmaS and sigmaR
4350  if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
4351  {
4352  ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
4353  return;
4354  }
4355 
4356  //define input data dimension with lattice
4357  int lN = N + 2;
4358 
4359  // Traverse each data point applying mean shift
4360  // to each data point
4361 
4362  // Allcocate memory for yk
4363  double *yk = new double [lN];
4364 
4365  // Allocate memory for Mh
4366  double *Mh = new double [lN];
4367 
4368  // let's use some temporary data
4369  double* sdata;
4370  sdata = new double[lN*L];
4371 
4372  // copy the scaled data
4373  int idxs, idxd;
4374  idxs = idxd = 0;
4375  if (N==3)
4376  {
4377  for(i=0; i<L; i++)
4378  {
4379  sdata[idxs++] = (i%width)/sigmaS;
4380  sdata[idxs++] = (i/width)/sigmaS;
4381  sdata[idxs++] = data[idxd++]/sigmaR;
4382  sdata[idxs++] = data[idxd++]/sigmaR;
4383  sdata[idxs++] = data[idxd++]/sigmaR;
4384  }
4385  } else if (N==1)
4386  {
4387  for(i=0; i<L; i++)
4388  {
4389  sdata[idxs++] = (i%width)/sigmaS;
4390  sdata[idxs++] = (i/width)/sigmaS;
4391  sdata[idxs++] = data[idxd++]/sigmaR;
4392  }
4393  } else
4394  {
4395  for(i=0; i<L; i++)
4396  {
4397  sdata[idxs++] = (i%width)/sigmaS;
4398  sdata[idxs++] = (i%width)/sigmaS;
4399  for (j=0; j<N; j++)
4400  sdata[idxs++] = data[idxd++]/sigmaR;
4401  }
4402  }
4403  // index the data in the 3d buckets (x, y, L)
4404  int* buckets;
4405  int* slist;
4406  slist = new int[L];
4407  int bucNeigh[27];
4408 
4409  double sMins; // just for L
4410  double sMaxs[3]; // for all
4411  sMaxs[0] = width/sigmaS;
4412  sMaxs[1] = height/sigmaS;
4413  sMins = sMaxs[2] = sdata[2];
4414  idxs = 2;
4415  double cval;
4416  for(i=0; i<L; i++)
4417  {
4418  cval = sdata[idxs];
4419  if (cval < sMins)
4420  sMins = cval;
4421  else if (cval > sMaxs[2])
4422  sMaxs[2] = cval;
4423 
4424  idxs += lN;
4425  }
4426 
4427  int nBuck1, nBuck2, nBuck3;
4428  int cBuck1, cBuck2, cBuck3, cBuck;
4429  nBuck1 = (int) (sMaxs[0] + 3);
4430  nBuck2 = (int) (sMaxs[1] + 3);
4431  nBuck3 = (int) (sMaxs[2] - sMins + 3);
4432  buckets = new int[nBuck1*nBuck2*nBuck3];
4433  for(i=0; i<(nBuck1*nBuck2*nBuck3); i++)
4434  buckets[i] = -1;
4435 
4436  idxs = 0;
4437  for(i=0; i<L; i++)
4438  {
4439  // find bucket for current data and add it to the list
4440  cBuck1 = (int) sdata[idxs] + 1;
4441  cBuck2 = (int) sdata[idxs+1] + 1;
4442  cBuck3 = (int) (sdata[idxs+2] - sMins) + 1;
4443  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4444 
4445  slist[i] = buckets[cBuck];
4446  buckets[cBuck] = i;
4447 
4448  idxs += lN;
4449  }
4450  // init bucNeigh
4451  idxd = 0;
4452  for (cBuck1=-1; cBuck1<=1; cBuck1++)
4453  {
4454  for (cBuck2=-1; cBuck2<=1; cBuck2++)
4455  {
4456  for (cBuck3=-1; cBuck3<=1; cBuck3++)
4457  {
4458  bucNeigh[idxd++] = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4459  }
4460  }
4461  }
4462  double wsuml, weight;
4463  double hiLTr = 80.0/sigmaR;
4464  // done indexing/hashing
4465 
4466  // proceed ...
4467 #ifdef PROMPT
4468  msSys.Prompt("done.\nApplying mean shift (Using Lattice)... ");
4469 #ifdef SHOW_PROGRESS
4470  msSys.Prompt("\n 0%%");
4471 #endif
4472 #endif
4473 
4474  for(i = 0; i < L; i++)
4475  {
4476 
4477  // Assign window center (window centers are
4478  // initialized by createLattice to be the point
4479  // data[i])
4480  idxs = i*lN;
4481  for (j=0; j<lN; j++)
4482  yk[j] = sdata[idxs+j];
4483 
4484  // Calculate the mean shift vector using the lattice
4485  // LatticeMSVector(Mh, yk);
4486  /*****************************************************/
4487  // Initialize mean shift vector
4488  for(j = 0; j < lN; j++)
4489  Mh[j] = 0;
4490  wsuml = 0;
4491  // uniformLSearch(Mh, yk_ptr); // modify to new
4492  // find bucket of yk
4493  cBuck1 = (int) yk[0] + 1;
4494  cBuck2 = (int) yk[1] + 1;
4495  cBuck3 = (int) (yk[2] - sMins) + 1;
4496  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4497  for (j=0; j<27; j++)
4498  {
4499  idxd = buckets[cBuck+bucNeigh[j]];
4500  // list parse, crt point is cHeadList
4501  while (idxd>=0)
4502  {
4503  idxs = lN*idxd;
4504  // determine if inside search window
4505  el = sdata[idxs+0]-yk[0];
4506  diff = el*el;
4507  el = sdata[idxs+1]-yk[1];
4508  diff += el*el;
4509 
4510  if (diff < 1.0)
4511  {
4512  el = sdata[idxs+2]-yk[2];
4513  if (yk[2] > hiLTr)
4514  diff = 4*el*el;
4515  else
4516  diff = el*el;
4517 
4518  if (N>1)
4519  {
4520  el = sdata[idxs+3]-yk[3];
4521  diff += el*el;
4522  el = sdata[idxs+4]-yk[4];
4523  diff += el*el;
4524  }
4525 
4526  if (diff < 1.0)
4527  {
4528  weight = 1-weightMap[idxd];
4529  for (k=0; k<lN; k++)
4530  Mh[k] += weight*sdata[idxs+k];
4531  wsuml += weight;
4532  }
4533  }
4534  idxd = slist[idxd];
4535  }
4536  }
4537  if (wsuml > 0)
4538  {
4539  for(j = 0; j < lN; j++)
4540  Mh[j] = Mh[j]/wsuml - yk[j];
4541  }
4542  else
4543  {
4544  for(j = 0; j < lN; j++)
4545  Mh[j] = 0;
4546  }
4547  /*****************************************************/
4548 
4549  // Calculate its magnitude squared
4550  mvAbs = 0;
4551  for(j = 0; j < lN; j++)
4552  mvAbs += Mh[j]*Mh[j];
4553 
4554  // Keep shifting window center until the magnitude squared of the
4555  // mean shift vector calculated at the window center location is
4556  // under a specified threshold (Epsilon)
4557 
4558  // NOTE: iteration count is for speed up purposes only - it
4559  // does not have any theoretical importance
4560  iterationCount = 1;
4561  while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
4562  {
4563 
4564  // Shift window location
4565  for(j = 0; j < lN; j++)
4566  yk[j] += Mh[j];
4567 
4568  // Calculate the mean shift vector at the new
4569  // window location using lattice
4570  // LatticeMSVector(Mh, yk);
4571  /*****************************************************/
4572  // Initialize mean shift vector
4573  for(j = 0; j < lN; j++)
4574  Mh[j] = 0;
4575  wsuml = 0;
4576  // uniformLSearch(Mh, yk_ptr); // modify to new
4577  // find bucket of yk
4578  cBuck1 = (int) yk[0] + 1;
4579  cBuck2 = (int) yk[1] + 1;
4580  cBuck3 = (int) (yk[2] - sMins) + 1;
4581  cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);
4582  for (j=0; j<27; j++)
4583  {
4584  idxd = buckets[cBuck+bucNeigh[j]];
4585  // list parse, crt point is cHeadList
4586  while (idxd>=0)
4587  {
4588  idxs = lN*idxd;
4589  // determine if inside search window
4590  el = sdata[idxs+0]-yk[0];
4591  diff = el*el;
4592  el = sdata[idxs+1]-yk[1];
4593  diff += el*el;
4594 
4595  if (diff < 1.0)
4596  {
4597  el = sdata[idxs+2]-yk[2];
4598  if (yk[2] > hiLTr)
4599  diff = 4*el*el;
4600  else
4601  diff = el*el;
4602 
4603  if (N>1)
4604  {
4605  el = sdata[idxs+3]-yk[3];
4606  diff += el*el;
4607  el = sdata[idxs+4]-yk[4];
4608  diff += el*el;
4609  }
4610 
4611  if (diff < 1.0)
4612  {
4613  weight = 1-weightMap[idxd];
4614  for (k=0; k<lN; k++)
4615  Mh[k] += weight*sdata[idxs+k];
4616  wsuml += weight;
4617  }
4618  }
4619  idxd = slist[idxd];
4620  }
4621  }
4622  if (wsuml > 0)
4623  {
4624  for(j = 0; j < lN; j++)
4625  Mh[j] = Mh[j]/wsuml - yk[j];
4626  }
4627  else
4628  {
4629  for(j = 0; j < lN; j++)
4630  Mh[j] = 0;
4631  }
4632  /*****************************************************/
4633 
4634  // Calculate its magnitude squared
4635  //mvAbs = 0;
4636  //for(j = 0; j < lN; j++)
4637  // mvAbs += Mh[j]*Mh[j];
4638  mvAbs = (Mh[0]*Mh[0]+Mh[1]*Mh[1])*sigmaS*sigmaS;
4639  if (N==3)
4640  mvAbs += (Mh[2]*Mh[2]+Mh[3]*Mh[3]+Mh[4]*Mh[4])*sigmaR*sigmaR;
4641  else
4642  mvAbs += Mh[2]*Mh[2]*sigmaR*sigmaR;
4643 
4644  // Increment interation count
4645  iterationCount++;
4646  }
4647 
4648  // Shift window location
4649  for(j = 0; j < lN; j++)
4650  yk[j] += Mh[j];
4651 
4652  //store result into msRawData...
4653  for(j = 0; j < N; j++)
4654  msRawData[N*i+j] = (float)(yk[j+2]*sigmaR);
4655 
4656  // Prompt user on progress
4657 #ifdef SHOW_PROGRESS
4658  percent_complete = (float)(i/(float)(L))*100;
4659  msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
4660 #endif
4661 
4662  // Check to see if the algorithm has been halted
4663  if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
4664  break;
4665  }
4666 
4667  // Prompt user that filtering is completed
4668 #ifdef PROMPT
4669 #ifdef SHOW_PROGRESS
4670  msSys.Prompt("\r");
4671 #endif
4672  msSys.Prompt("done.");
4673 #endif
4674 
4675  // de-allocate memory
4676  delete [] buckets;
4677  delete [] slist;
4678  delete [] sdata;
4679 
4680  delete [] yk;
4681  delete [] Mh;
4682 
4683  // done.
4684  return;
4685 
4686 }
4687 
4688 void msImageProcessor::SetSpeedThreshold(float speedUpThreshold)
4689 {
4690  speedThreshold = speedUpThreshold;
4691 }
4692 
4693 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
4694 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
4695 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ END OF CLASS DEFINITION @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
4696 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
4697 /*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
Generated by doxygen 1.9.1