NICE_Segmentation/edisonSegm/msImageProcessor.cpp

/*******************************************************

                 Mean Shift Analysis Library
 =============================================


 The mean shift library is a collection of routines
 that use the mean shift algorithm. Using this algorithm,
 the necessary output will be generated needed
 to analyze a given input set of data.

  Mean Shift Image Processor Class:
  ================================

 The following class inherits from the mean shift library
 in order to perform the specialized tasks of image
 segmentation and filtering.

 The definition of the Mean Shift Image Processor Class
 is provided below. Its prototype is provided in
 'msImageProcessor.h'.

The theory is described in the papers:

  D. Comaniciu, P. Meer: Mean Shift: A robust approach toward feature
          space analysis.

  C. Christoudias, B. Georgescu, P. Meer: Synergism in low level vision.

and they are is available at:
  http://www.caip.rutgers.edu/riul/research/papers/

Implemented by Chris M. Christoudias, Bogdan Georgescu
********************************************************/
#ifdef NICE_USELIB_OPENMP
#include <omp.h>
#endif

//include image processor class prototype
#include "msImageProcessor.h"

//include needed libraries
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <iostream>
using namespace std;

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@      PUBLIC METHODS     @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

/*/\/\/\/\/\/\/\/\/\/\/\/\*/
/* Constructor/Destructor */
/*\/\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Class Constructor                                    */
/*******************************************************/
/*Post:                                                */
/*      The msImageProcessor class has been properly   */
/*      initialized.                                   */
/*******************************************************/

msImageProcessor::msImageProcessor ( void )
{
  clog << "[log] msImageProcessor::msImageProcessor: use Edison ";
#ifdef NICE_USELIB_OPENMP
  clog << "parallel!" << endl;
  omp_set_dynamic ( 0 );
#else
  clog << "seriell!" << endl;
#endif

  //intialize basin of attraction structure
  //used by the filtering algorithm
  modeTable   = NULL;
  pointList   = NULL;
  pointCount   = 0;

  //initialize region list
  regionList   = NULL;

  //initialize output structures...
  msRawData   = NULL;
  labels    = NULL;
  modes    = NULL;
  modePointCounts   = NULL;
  regionCount   = 0;

  //intialize temporary buffers used for
  //performing connected components
  indexTable   = NULL;
  LUV_data   = NULL;

  //initialize region adjacency matrix
  raList    = NULL;
  freeRAList   = NULL;
  raPool    = NULL;

  //intialize visit table to having NULL entries
  visitTable   = NULL;

  //initialize epsilon such that transitive closure
  //does not take edge strength into consideration when
  //fusing regions of similar color
  epsilon    = 1.0;

  //initialize class state to indicate that
  //an output data structure has not yet been
  //created...
  class_state.OUTPUT_DEFINED = false;

//Changed by Sushil from 1.0 to 0.1, 11/11/2008
  LUV_treshold = 0.1;
}

/*******************************************************/
/*Class Destructor                                     */
/*******************************************************/
/*Post:                                                */
/*      The msImageProcessor class has been properly   */
/*      destroyed.                                     */
/*******************************************************/

msImageProcessor::~msImageProcessor ( void )
{

  //de-allocate memory
  if ( class_state.OUTPUT_DEFINED ) DestroyOutput();
  if ( regionList )     delete regionList;
  regionList = NULL;

  //done.

}

/*/\/\/\/\/\/\/\/\/\/\/\/\/\*/
/*  Input Image Declaration */
/*\/\/\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Define Image                                         */
/*******************************************************/
/*Uploads an image into the image segmenter class to   */
/*be segmented.                                        */
/*******************************************************/
/*Pre:                                                 */
/*      - data_ is a one dimensional array of unsigned */
/*        char RGB vectors                             */
/*      - type is the type of the image: COLOR or      */
/*        GREYSCALE                                    */
/*      - height_ and width_ define the dimension of   */
/*        the image                                    */
/*      - if the image is of type GREYSCALE then       */
/*        data containes only one number per pixel     */
/*        location, where a pixel location is defined  */
/*        by the index into the data array             */
/*Post:                                                */
/*      - the image specified has been uploaded into   */
/*        the image segmenter class to be segmented.   */
/*******************************************************/

void msImageProcessor::DefineImage ( byte *data_, imageType type, int height_, int width_ )
{
  /* Ein neuer LUV-Vektor wird angelegt. In diesen wird der Inhalt von data_ gespeichert. Kann data_ auch mit 'const' Qualifier definiert werden ?

  */

  //obtain image dimension from image type
  int dim;
  if ( type == COLOR )
    dim = 3;
  else
    dim = 1;

  //perfor rgb to luv conversion
  int  i;
  float *luv = new float [height_*width_*dim];
  if ( dim == 1 )
  {
    for ( i = 0; i < height_*width_; i++ )
      luv[i] = ( float ) ( data_[i] );
  }
  else
  {
    for ( i = 0; i < height_*width_; i++ )
    {
      RGBtoLUV ( &data_[dim*i], &luv[dim*i] );
    }
  }

  //define input defined on a lattice using mean shift base class
  DefineLInput ( luv, height_, width_, dim );

  //Define a default kernel if it has not been already
  //defined by user
  if ( !h )
  {
    //define default kernel paramerters...
    kernelType k[2]  = {Uniform, Uniform};
    int   P[2]  = {2, N};
    float  tempH[2] = {1.0 , 1.0};

    //define default kernel in mean shift base class
    DefineKernel ( k, tempH, P, 2 );
  }

  //de-allocate memory
  delete [] luv;

  //done.
  return;

}

void msImageProcessor::DefineBgImage ( byte* data_, imageType type, int height_, int width_ )
{

  //obtain image dimension from image type
  int dim;
  if ( type == COLOR )
    dim = 3;
  else
    dim = 1;

  //perform texton classification
  int  i;

  float *luv = new float [height_*width_*dim];

  if ( dim == 1 )
  {
    for ( i = 0; i < height_*width_; i++ )
      luv[i] = ( float ) ( data_[i] );
  }
  else
  {
    for ( i = 0; i < height_*width_; i++ )
      RGBtoLUV ( &data_[dim*i], &luv[dim*i] );

  }

  //define input defined on a lattice using mean shift base class
  DefineLInput ( luv, height_, width_, dim );


  //Define a default kernel if it has not been already
  //defined by user
  if ( !h )
  {
    //define default kernel paramerters...
    kernelType k[2]  = {Uniform, Uniform};
    int   P[2]  = {2, N};
    float  tempH[2] = {1.0 , 1.0};

    //define default kernel in mean shift base class
    DefineKernel ( k, tempH, P, 2 );
  }

  //de-allocate memory
  delete [] luv;

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\*/
/*   Weight Map   */
/*\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Set Weight Map                                       */
/*******************************************************/
/*Populates the weight map with specified edge         */
/*strengths.                                           */
/*******************************************************/
/*Pre:                                                 */
/*      - wm is a floating point array of size         */
/*        (height x width) specifying for each pixel   */
/*        edge strength.                               */
/*      - eps is a threshold used to fuse similar      */
/*        regions during transitive closure.           */
/*Post:                                                */
/*      - wm has been used to populate the weight      */
/*        map.                                         */
/*      - the threshold used during transitive closure */
/*        is taken as eps.                             */
/*******************************************************/

void msImageProcessor::SetWeightMap ( float *wm, float eps )
{

  //initlaize confmap using wm
  SetLatticeWeightMap ( wm );

  //set threshold value
  if ( ( epsilon = eps ) < 0 )
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "SetWeightMap", ( char* ) "Threshold is negative." );

  //done.
  return;

}

/*******************************************************/
/*Remove Weight Map                                    */
/*******************************************************/
/*Removes the weight map.                              */
/*******************************************************/
/*Post:                                                */
/*      - the weight map has been removed.             */
/*      - if a weight map did not exist NO error       */
/*        is flagged.                                  */
/*******************************************************/

void msImageProcessor::RemoveWeightMap ( void )
{

  //remove confmap
  RemoveLatticeWeightMap();

  //set threshold value to zero
  epsilon = 0;

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\*/
/* Image Filtering  */
/*\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Filter                                               */
/*******************************************************/
/*Performs mean shift filtering on the specified input */
/*image using a user defined kernel.                   */
/*******************************************************/
/*Pre:                                                 */
/*      - the user defined kernel used to apply mean   */
/*        shift filtering to the defined input image   */
/*        has spatial bandwidth sigmaS and range band- */
/*        width sigmaR                                 */
/*      - speedUpLevel determines whether or not the   */
/*        filtering should be optimized for faster     */
/*        execution: a value of NO_SPEEDUP turns this  */
/*        optimization off and a value SPEEDUP turns   */
/*        this optimization on                         */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - mean shift filtering has been applied to the */
/*        input image using a user defined kernel      */
/*      - the filtered image is stored in the private  */
/*        data members of the msImageProcessor class.  */
/*******************************************************/

void msImageProcessor::Filter ( int sigmaS, float sigmaR, SpeedUpLevel speedUpLevel )
{

  //Check Class consistency...

  //check:
  // (1) if this operation is consistent
  // (2) if kernel was created
  // (3) if data set is defined
  // (4) if the dimension of the kernel agrees with that
  //     of the defined data set
  // if not ... flag an error!
  classConsistencyCheck ( N + 2, true );
  if ( ErrorStatus == EL_ERROR )
    return;

  //If the algorithm has been halted, then exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.0 ) ) ) == EL_HALT )
  {
    return;
  }

  //If the image has just been read then allocate memory
  //for and initialize output data structure used to store
  //image modes and their corresponding regions...
  if ( class_state.OUTPUT_DEFINED == false )
  {
    InitializeOutput();

    //check for errors...
    if ( ErrorStatus == EL_ERROR )
      return;
  }

  //****************** Allocate Memory ******************

  //Allocate memory for basin of attraction mode structure...
  if ( ( ! ( modeTable = new unsigned char [L] ) ) || ( ! ( pointList = new int [L] ) ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Allocate", ( char* ) "Not enough memory." );
    return;
  }

  //start timer
#ifdef PROMPT
  double timer;
  msSys.StartTimer();
#endif

  //*****************************************************

// Parallelisieren !!!

  //filter image according to speedup level...
  switch ( speedUpLevel )
  {
      //no speedup...
    case NO_SPEEDUP:
      //NonOptimizedFilter((float)(sigmaS), sigmaR); break;
      NewNonOptimizedFilter ( ( float ) ( sigmaS ), sigmaR );
      break;
      //medium speedup
    case MED_SPEEDUP:
      //OptimizedFilter1((float)(sigmaS), sigmaR);  break;
      NewOptimizedFilter1 ( ( float ) ( sigmaS ), sigmaR );
      break;
      //high speedup
    case HIGH_SPEEDUP:
      //OptimizedFilter2((float)(sigmaS), sigmaR);  break;
      NewOptimizedFilter2 ( ( float ) ( sigmaS ), sigmaR );
      break;
      // new speedup
  }

  //****************** Deallocate Memory ******************

  //de-allocate memory used by basin of attraction mode structure
  delete [] modeTable;
  delete [] pointList;

  //re-initialize structure
  modeTable = NULL;
  pointList = NULL;
  pointCount = 0;

  //*******************************************************

  //If the algorithm has been halted, then de-allocate the output
  //and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.8 ) ) ) == EL_HALT )
  {
    DestroyOutput();
    return;
  }

  //Label image regions, also if segmentation is not to be
  //performed use the resulting classification structure to
  //calculate the image boundaries...

  /*
  //copy msRawData into LUV_data, rounding each component of each
  //LUV value stored by msRawData to the nearest integer
  int i;
  for(i = 0; i < L*N; i++)
  {
  if(msRawData[i] < 0)
  LUV_data[i] = (int)(msRawData[i] - 0.5);
  else
  LUV_data[i] = (int)(msRawData[i] + 0.5);
  }
  */
// Parallelisieren !!!
  int i;
  for ( i = 0; i < L*N; i++ )
  {
    LUV_data[i] = msRawData[i];
  }


#ifdef PROMPT
  timer = msSys.ElapsedTime();
  printf ( ( char* ) "(%6.2f sec)\nConnecting regions         ...", timer );
  msSys.StartTimer();
#endif

  //Perform connecting (label image regions) using LUV_data
  Connect();

#ifdef PROMPT
  timer = msSys.ElapsedTime();
  printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount );
  msSys.StartTimer();
#endif

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\/\*/
/* Image Region Fusing  */
/*\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Fuse Regions                                         */
/*******************************************************/
/*Fuses the regions of a filtered image.               */
/*******************************************************/
/*Pre:                                                 */
/*      - the range radius is specified by sigmaR      */
/*      - minRegion is the minimum point density that  */
/*        a region may have in the resulting segment-  */
/*        ed image                                     */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - the image regions have been fused.           */
/*      - if an result is stored by this class then    */
/*        this result is used as input to this method. */
/*      - if no result is stored by this class,        */
/*        the input image defined by calling the       */
/*        method DefineImage is used.                  */
/*******************************************************/

void msImageProcessor::FuseRegions ( float sigmaS, int minRegion )
{

  //Check Class consistency...

  //check:
  // (1) if this operation is consistent
  // (2) if kernel was created
  // (3) if data set is defined
  // (4) if the dimension of the kernel agrees with that
  //     of the defined data set
  // if not ... flag an error!
  classConsistencyCheck ( N + 2, true );
  if ( ErrorStatus == EL_ERROR )
    return;

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.8 ) ) ) == EL_HALT )
  {
    if ( class_state.OUTPUT_DEFINED ) DestroyOutput();
    return;
  }

  //obtain sigmaS (make sure it is not zero or negative, if not
  //flag an error)
  if ( ( h[1] = sigmaS ) <= 0 )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "FuseRegions", ( char* ) "The feature radius must be greater than or equal to zero." );
    return;
  }

  //if output has not yet been generated then classify the input
  //image regions to be fused...
  if ( ! ( class_state.OUTPUT_DEFINED ) )
  {

    //Initialize output data structure used to store
    //image modes and their corresponding regions...
    InitializeOutput();

    //check for errors...
    if ( ErrorStatus == EL_ERROR )
      return;

    //copy data into LUV_data used to classify
    //image regions
    /*
    int i;
    for(i = 0; i < L*N; i++)
    {
    if(data[i] < 0)
    LUV_data[i] = (int)(data[i] - 0.5);
    else
    LUV_data[i] = (int)(data[i] + 0.5);
    }
    */
    int i;
    for ( i = 0; i < L*N; i++ )
    {
      LUV_data[i] = data[i];
    }

#ifdef PROMPT
    printf ( ( char* ) "Connecting regions         ..." );
    msSys.StartTimer();
#endif

    //Perform connecting (label image regions) using LUV_data
    Connect();

    //check for errors
    if ( ErrorStatus == EL_ERROR )
      return;

#ifdef PROMPT
    double timer = msSys.ElapsedTime();
    printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount );
#endif

  }

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.85 ) ) ) == EL_HALT )
  {
    DestroyOutput();
    return;
  }

#ifdef PROMPT
  printf ( ( char* ) "Applying transitive closure..." );
  msSys.StartTimer();
#endif

  //allocate memory visit table
  visitTable = new unsigned char [L];

  //Apply transitive closure iteratively to the regions classified
  //by the RAM updating labels and modes until the color of each neighboring
  //region is within sqrt(rR2) of one another.
  rR2 = ( float ) ( h[1] * h[1] * 0.25 );
  TransitiveClosure();
  int oldRC = regionCount;
  int deltaRC, counter = 0;
  do {
    TransitiveClosure();
    deltaRC = oldRC - regionCount;
    oldRC = regionCount;
    counter++;
  } while ( ( deltaRC <= 0 ) && ( counter < 10 ) );

  //de-allocate memory for visit table
  delete [] visitTable;
  visitTable = NULL;

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and region adjacency matrix and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 1.0 ) ) ) == EL_HALT )
  {
    DestroyRAM();
    DestroyOutput();
    return;
  }

#ifdef PROMPT
  double timer = msSys.ElapsedTime();
  printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d)\nPruning spurious regions   ...", timer, regionCount );
  msSys.StartTimer();
#endif

  //Prune spurious regions (regions whose area is under
  //minRegion) using RAM
  Prune ( minRegion );

#ifdef PROMPT
  timer = msSys.ElapsedTime();
  printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d)\n", timer, regionCount );
  msSys.StartTimer();
#endif

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and region adjacency matrix and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 1.0 ) ) ) == EL_HALT )
  {
    DestroyRAM();
    DestroyOutput();
    return;
  }

  //de-allocate memory for region adjacency matrix
  DestroyRAM();

  //output to msRawData
  int i, j, label;
  for ( i = 0; i < L; i++ )
  {
    label = labels[i];
    for ( j = 0; j < N; j++ )
    {
      msRawData[N*i+j] = modes[N*label+j];
    }
  }

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\*/
/* Image Segmentation */
/*\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Segment                                              */
/*******************************************************/
/*Segments the defined image.                          */
/*******************************************************/
/*Pre:                                                 */
/*      - sigmaS and sigmaR are the spatial and range  */
/*        radii of the search window respectively      */
/*      - minRegion is the minimum point density that  */
/*        a region may have in the resulting segment-  */
/*        ed image                                     */
/*      - speedUpLevel determines whether or not the   */
/*        filtering should be optimized for faster     */
/*        execution: a value of NO_SPEEDUP turns this  */
/*        optimization off and a value SPEEDUP turns   */
/*        this optimization on                         */
/*Post:                                                */
/*      - the defined image is segmented and the       */
/*        resulting segmented image is stored in the   */
/*        private data members of the image segmenter  */
/*        class.                                       */
/*      - any regions whose point densities are less   */
/*        than or equal to minRegion have been pruned  */
/*        from the segmented image.                    */
/*******************************************************/

void msImageProcessor::Segment ( int sigmaS, float sigmaR, int minRegion, SpeedUpLevel speedUpLevel )
{

  //make sure kernel is properly defined...
  if ( ( !h ) || ( kp < 2 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "Kernel corrupt or undefined." );
    return;
  }

  //Apply mean shift to data set using sigmaS and sigmaR...
  Filter ( sigmaS, sigmaR, speedUpLevel );

  //check for errors
  if ( ErrorStatus == EL_ERROR )
    return;

  //check to see if the system has been halted, if so exit
  if ( ErrorStatus == EL_HALT )
    return;

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.85 ) ) ) == EL_HALT )
  {
    DestroyOutput();
    return;
  }

#ifdef PROMPT
  printf ( ( char* ) "Applying transitive closure..." );
  msSys.StartTimer();
#endif

  //allocate memory visit table
  visitTable = new unsigned char [L];

  //Apply transitive closure iteratively to the regions classified
  //by the RAM updating labels and modes until the color of each neighboring
  //region is within sqrt(rR2) of one another.
  rR2 = ( float ) ( h[1] * h[1] * 0.25 );
  TransitiveClosure();
  int oldRC = regionCount;
  int deltaRC, counter = 0;
  do {
    TransitiveClosure();
    deltaRC = oldRC - regionCount;
    oldRC = regionCount;
    counter++;
  } while ( ( deltaRC <= 0 ) && ( counter < 10 ) );

  //de-allocate memory for visit table
  delete [] visitTable;
  visitTable = NULL;

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and regions adjacency matrix and exit
  if ( ( ErrorStatus = msSys.Progress ( ( float ) ( 0.95 ) ) ) == EL_HALT )
  {
    DestroyRAM();
    DestroyOutput();
    return;
  }

#ifdef PROMPT
  double timer = msSys.ElapsedTime();
  printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d).\nPruning spurious regions\t... ", timer, regionCount );
  msSys.StartTimer();
#endif

  //Prune spurious regions (regions whose area is under
  //minRegion) using RAM
  Prune ( minRegion );

#ifdef PROMPT
  timer = msSys.ElapsedTime();
  printf ( ( char* ) "done. (%6.2f seconds, numRegions = %6d)\nPruning spurious regions    ...", timer, regionCount );
  msSys.StartTimer();
#endif

  //Check to see if the algorithm is to be halted, if so then
  //destroy output and regions adjacency matrix and exit
  if ( ( ErrorStatus = msSys.Progress ( 1.0 ) ) == EL_HALT )
  {
    DestroyRAM();
    DestroyOutput();
    return;
  }

  //de-allocate memory for region adjacency matrix
  DestroyRAM();

// Parallelisieren !!!

  //output to msRawData
  int j, i, label;

  for ( i = 0; i < L; i++ )
  {
    label = labels[i];
    for ( j = 0; j < N; j++ )
    {
      msRawData[N*i+j] = modes[N*label+j];
    }
  }

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\/\/\*/
/*  Data Space Conversion */
/*\/\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*RGB To LUV                                           */
/*******************************************************/
/*Converts an RGB vector to LUV.                       */
/*                                                     */
/*See:                                                 */
/*   G. Wyszecki and W.S. Stiles: Color Science:       */
/*   Concepts and Methods, Quantitative Data and       */
/*   Formulae, Wiley, New York, 1982.                  */
/*******************************************************/
/*Pre:                                                 */
/*      - rgbVal is an unsigned char array containing  */
/*        the RGB vector                               */
/*      - luvVal is a floating point array containing  */
/*        the resulting LUV vector                     */
/*Post:                                                */
/*      - rgbVal has been converted to LUV and the     */
/*        result has been stored in luvVal.            */
/*******************************************************/

void msImageProcessor::RGBtoLUV ( byte *rgbVal, float *luvVal )
{

  //delcare variables
  double x, y, z, L0, u_prime, v_prime, constant;

  //convert RGB to XYZ...
  x  = XYZ[0][0] * rgbVal[0] + XYZ[0][1] * rgbVal[1] + XYZ[0][2] * rgbVal[2];
  y  = XYZ[1][0] * rgbVal[0] + XYZ[1][1] * rgbVal[1] + XYZ[1][2] * rgbVal[2];
  z  = XYZ[2][0] * rgbVal[0] + XYZ[2][1] * rgbVal[1] + XYZ[2][2] * rgbVal[2];

  //convert XYZ to LUV...

  //compute L*
  L0  = y / ( 255.0 * Yn );
  if ( L0 > Lt )
    luvVal[0] = ( float ) ( 116.0 * ( pow ( L0, 1.0 / 3.0 ) ) - 16.0 );
  else
    luvVal[0] = ( float ) ( 903.3 * L0 );

  //compute u_prime and v_prime
  constant = x + 15 * y + 3 * z;
  if ( constant != 0 )
  {
    u_prime = ( 4 * x ) / constant;
    v_prime = ( 9 * y ) / constant;
  }
  else
  {
    u_prime = 4.0;
    v_prime = 9.0 / 15.0;
  }

  //compute u* and v*
  luvVal[1] = ( float ) ( 13 * luvVal[0] * ( u_prime - Un_prime ) );
  luvVal[2] = ( float ) ( 13 * luvVal[0] * ( v_prime - Vn_prime ) );

  //done.
  return;

}

/*******************************************************/
/*LUV To RGB                                           */
/*******************************************************/
/*Converts an LUV vector to RGB.                       */
/*******************************************************/
/*Pre:                                                 */
/*      - luvVal is a floating point array containing  */
/*        the LUV vector                               */
/*      - rgbVal is an unsigned char array containing  */
/*        the resulting RGB vector                     */
/*Post:                                                */
/*      - luvVal has been converted to RGB and the     */
/*        result has been stored in rgbVal.            */
/*******************************************************/

//define inline rounding function...
inline int my_round ( double in_x )
{
  if ( in_x < 0 )
    return ( int ) ( in_x - 0.5 );
  else
    return ( int ) ( in_x + 0.5 );
}

void msImageProcessor::LUVtoRGB ( float *luvVal, byte *rgbVal )
{

  //declare variables...
  int  r, g, b;
  double x, y, z, u_prime, v_prime;

  //perform conversion
  if ( luvVal[0] < 0.1 )
    r = g = b = 0;
  else
  {
    //convert luv to xyz...
    if ( luvVal[0] < 8.0 )
      y = Yn * luvVal[0] / 903.3;
    else
    {
      y = ( luvVal[0] + 16.0 ) / 116.0;
      y  *= Yn * y * y;
    }

    u_prime = luvVal[1] / ( 13 * luvVal[0] ) + Un_prime;
    v_prime = luvVal[2] / ( 13 * luvVal[0] ) + Vn_prime;

    x  = 9 * u_prime * y / ( 4 * v_prime );
    z  = ( 12 - 3 * u_prime - 20 * v_prime ) * y / ( 4 * v_prime );

    //convert xyz to rgb...
    //[r, g, b] = RGB*[x, y, z]*255.0
    r  = my_round ( ( RGB[0][0] * x + RGB[0][1] * y + RGB[0][2] * z ) * 255.0 );
    g  = my_round ( ( RGB[1][0] * x + RGB[1][1] * y + RGB[1][2] * z ) * 255.0 );
    b  = my_round ( ( RGB[2][0] * x + RGB[2][1] * y + RGB[2][2] * z ) * 255.0 );

    //check bounds...
    if ( r < 0 ) r = 0;
    if ( r > 255 ) r = 255;
    if ( g < 0 ) g = 0;
    if ( g > 255 ) g = 255;
    if ( b < 0 ) b = 0;
    if ( b > 255 ) b = 255;

  }

  //assign rgb values to rgb vector rgbVal
  rgbVal[0] = r;
  rgbVal[1] = g;
  rgbVal[2] = b;

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/
/*  Filtered and Segmented Image Output */
/*\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Get Raw Data                                         */
/*******************************************************/
/*The output image data is returned.                   */
/*******************************************************/
/*Pre:                                                 */
/*      - outputImageData is a pre-allocated floating  */
/*        point array used to store the filtered or    */
/*        segmented image pixels.                      */
/*Post:                                                */
/*      - the filtered or segmented image data is      */
/*        stored by outputImageData.                   */
/*******************************************************/

void msImageProcessor::GetRawData ( float *outputImageData )
{
  //make sure that outputImageData is not NULL
  if ( !outputImageData )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetRawData", ( char* ) "Output image data buffer is NULL." );
    return;
  }

  //copy msRawData to outputImageData
  int i;
  for ( i = 0; i < L*N; i++ )
    outputImageData[i] = msRawData[i];

  //done.
  return;
}

/*******************************************************/
/*Get Results                                          */
/*******************************************************/
/*The output image is returned.                        */
/*******************************************************/
/*Pre:                                                 */
/*      - outputImage is a pre-allocated unsinged char */
/*        array used to store the filtered or segment- */
/*        ed image pixels                              */
/*Post:                                                */
/*      - the filtered or segmented image is stored by */
/*        outputImage.                                 */
/*******************************************************/

void msImageProcessor::GetResults ( byte *outputImage )
{

  //make sure that outpuImage is not NULL
  if ( !outputImage )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetResults", ( char* ) "Output image buffer is NULL." );
    return;
  }

  //if the image type is GREYSCALE simply
  //copy it over to the segmentedImage
  if ( N == 1 )
  {
    //copy over msRawData to segmentedImage checking
    //bounds
    int i, pxValue;
    for ( i = 0; i < L; i++ )
    {

      //get value
      pxValue = ( int ) ( msRawData[i] + 0.5 );

      //store into segmented image checking bounds...
      if ( pxValue < 0 )
        outputImage[i] = ( byte ) ( 0 );
      else if ( pxValue > 255 )
        outputImage[i] = ( byte ) ( 255 );
      else
        outputImage[i] = ( byte ) ( pxValue );

    }

  }
  else if ( N == 3 )
  {

    //otherwise convert msRawData from LUV to RGB
    //storing the result in segmentedImage
    int i;
    for ( i = 0; i < L; i++ )
      LUVtoRGB ( &msRawData[N*i], &outputImage[N*i] );

  }
  else
    //Unknown image type: should use MeanShift::GetRawData()...
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetResults", ( char* ) "Unknown image type. Try using MeanShift::GetRawData()." );

  //done.
  return;

}

/*******************************************************/
/*Get Boundaries                                       */
/*******************************************************/
/*A region list containing the boundary locations for  */
/*each region is returned.                             */
/*******************************************************/
/*Post:                                                */
/*      - a region list object containing the boundary */
/*        locations for each region is constructed     */
/*      - the region list is returned                  */
/*      - NULL is returned if the image has not been   */
/*        filtered or segmented                        */
/*******************************************************/

RegionList *msImageProcessor::GetBoundaries ( void )
{

  //define bounds using label information
  if ( class_state.OUTPUT_DEFINED )
    DefineBoundaries();

  //return region list structure
  return regionList;

}

/*******************************************************/
/*Get Regions                                          */
/*******************************************************/
/*Returns the regions of the processed image.          */
/*******************************************************/
/*Pre:                                                 */
/*      - labels_out is an integer array of size       */
/*        height*width that stores for each pixel a    */
/*        label relating that pixel to a corresponding */
/*        region in the image                          */
/*      - modes_out is floating point array of size    */
/*        regionCount*N storing the feature component  */
/*        of each region, and indexed by region label  */
/*      - modePointCounts is an integer array of size  */
/*        regionCount, indexed by region label, that   */
/*        stores the area of each region in pixels.    */
/*Post:                                                */
/*      If an input image was defined and processed,   */
/*      - memory has been allocated for labels_out,    */
/*        modes_out and MPC_out.                       */
/*      - labels_out, modes_out, and MPC_out have been */
/*        populated.                                   */
/*      - the number of regions contained by the segm- */
/*        ented image has been returned.               */
/*      If the image has not been defined or processed */
/*      or if there is in-sufficient memory,           */
/*      - no memory has been allocated for labels_out, */
/*        modes_out, and MPC_out.                      */
/*      - -1 is returned for regionCount.              */
/*******************************************************/

int msImageProcessor::GetRegions ( int **labels_out, float **modes_out, int **MPC_out )
{
  //check to see if output has been defined for the given input image...
  if ( class_state.OUTPUT_DEFINED == false )
    return -1;

  //allocate memory for labels_out, modes_out and MPC_out based
  //on output storage structure
  int  *labels_ = *labels_out, *MPC_out_ = *MPC_out;
  float *modes_  = *modes_out;
  if ( ! ( labels_ = new int [L] ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetRegions", ( char* ) "Not enough memory." );
    return -1;
  }
  if ( ! ( modes_ = new float [regionCount*N] ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetRegions", ( char* ) "Not enough memory." );
    return -1;
  }
  if ( ! ( MPC_out_ = new int [regionCount] ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "GetRegions", ( char* ) "Not enough memory." );
    return -1;
  }

  //populate labels_out with image labels
  int i;
  for ( i = 0; i < L; i++ )
    labels_[i] = labels[i];

  //populate modes_out and MPC_out with the color and point
  //count of each region
  for ( i = 0; i < regionCount*N; i++ )
    modes_[i] = modes[i];
  for ( i = 0; i < regionCount; i++ )
    MPC_out_[i] = modePointCounts[i];

  //done. Return the number of regions resulting from filtering or segmentation.
  return regionCount;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@     PRIVATE METHODS     @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

/*/\/\/\/\/\/\/\/\/\*/
/*  Image Filtering */
/*\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Non Optimized Filter                                 */
/*******************************************************/
/*Performs mean shift filtering on the specified input */
/*image using a user defined kernel.                   */
/*******************************************************/
/*Pre:                                                 */
/*      - the user defined kernel used to apply mean   */
/*        shift filtering to the defined input image   */
/*        has spatial bandwidth sigmaS and range band- */
/*        width sigmaR                                 */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - mean shift filtering has been applied to the */
/*        input image using a user defined kernel      */
/*      - the filtered image is stored in the private  */
/*        data members of the msImageProcessor class.  */
/*******************************************************/

void msImageProcessor::NonOptimizedFilter ( float sigmaS, float sigmaR )
{

  // Declare Variables
  int   iterationCount, i, j;
  double mvAbs;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice)... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif

// Parallelisieren !!!

  for ( i = 0; i < L; i++ )
  {

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    yk[0] = i % width;
    yk[1] = i / width;
    for ( j = 0; j < N; j++ )
      yk[j+2] = data[N*i+j];

    // Calculate the mean shift vector using the lattice
    LatticeMSVector ( Mh, yk );

    // Calculate its magnitude squared
    mvAbs = 0;
    for ( j = 0; j < lN; j++ )
      mvAbs += Mh[j] * Mh[j];

    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // Calculate the mean shift vector at the new
      // window location using lattice
      LatticeMSVector ( Mh, yk );

      // Calculate its magnitude squared
      mvAbs = 0;
      for ( j = 0; j < lN; j++ )
        mvAbs += Mh[j] * Mh[j];

      // Increment interation count
      iterationCount++;

    }

    // Shift window location
    for ( j = 0; j < lN; j++ )
      yk[j] += Mh[j];

    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;
  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif

  // de-allocate memory
  delete [] yk;
  delete [] Mh;

  // done.
  return;

}



/*******************************************************/
/*Optimized Filter 1                                   */
/*******************************************************/
/*Performs mean shift filtering on the specified input */
/*image using a user defined kernel. Previous mode     */
/*information is used to avoid re-applying mean shift  */
/*on certain data points to improve performance.       */
/*******************************************************/
/*Pre:                                                 */
/*      - the user defined kernel used to apply mean   */
/*        shift filtering to the defined input image   */
/*        has spatial bandwidth sigmaS and range band- */
/*        width sigmaR                                 */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - mean shift filtering has been applied to the */
/*        input image using a user defined kernel      */
/*      - the filtered image is stored in the private  */
/*        data members of the msImageProcessor class.  */
/*******************************************************/

void msImageProcessor::OptimizedFilter1 ( float sigmaS, float sigmaR )
{

  // Declare Variables
  int  iterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;
  float *modeCandidatePoint;
  double mvAbs, diff, el;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // Initialize mode table used for basin of attraction
  memset ( modeTable, 0, width*height );

  // Allocate memory mode candidate data point...
  //floating point version
  modeCandidatePoint = new float [N];

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice) ... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif

// Parallelisieren !!!

  for ( i = 0; i < L; i++ )
  {
    // if a mode was already assigned to this data point
    // then skip this point, otherwise proceed to
    // find its mode by applying mean shift...
    if ( modeTable[i] == 1 )
      continue;

    // initialize point list...
    pointCount = 0;

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    yk[0] = i % width;
    yk[1] = i / width;
    for ( j = 0; j < N; j++ )
      yk[j+2] = data[N*i+j];

    // Calculate the mean shift vector using the lattice
    LatticeMSVector ( Mh, yk );

    // Calculate its magnitude squared
    mvAbs = 0;
    for ( j = 0; j < lN; j++ )
      mvAbs += Mh[j] * Mh[j];

    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // check to see if the current mode location is in the
      // basin of attraction...

      // calculate the location of yk on the lattice
      modeCandidateX = ( int ) ( yk[0] + 0.5 );
      modeCandidateY = ( int ) ( yk[1] + 0.5 );
      modeCandidate_i = modeCandidateY * width + modeCandidateX;

      // if mvAbs != 0 (yk did indeed move) then check
      // location basin_i in the mode table to see if
      // this data point either:

      // (1) has not been associated with a mode yet
      //     (modeTable[basin_i] = 0), so associate
      //     it with this one
      //
      // (2) it has been associated with a mode other
      //     than the one that this data point is converging
      //     to (modeTable[basin_i] = 1), so assign to
      //     this data point the same mode as that of basin_i

      if ( ( modeTable[modeCandidate_i] != 2 ) && ( modeCandidate_i != i ) )
      {
        // obtain the data point at basin_i to
        // see if it is within h*TC_DIST_FACTOR of
        // of yk
        for ( j = 0; j < N; j++ )
          modeCandidatePoint[j] = data[N*modeCandidate_i + j];

        // check basin on non-spatial data spaces only
        k = 1;
        s = 0;
        diff = 0;
        while ( ( diff < TC_DIST_FACTOR ) && ( k < kp ) )
        {
          diff = 0;
          for ( p = 0; p < P[k]; p++ )
          {
            el = ( modeCandidatePoint[p+s] - yk[p+s+2] ) / h[k];
            diff += el * el;
          }
          s += P[k];
          k++;
        }

        // if the data point at basin_i is within
        // a distance of h*TC_DIST_FACTOR of yk
        // then depending on modeTable[basin_i] perform
        // either (1) or (2)
        if ( diff < TC_DIST_FACTOR )
        {
          // if the data point at basin_i has not
          // been associated to a mode then associate
          // it with the mode that this one will converge
          // to
          if ( modeTable[modeCandidate_i] == 0 )
          {
            // no mode associated yet so associate
            // it with this one...
            pointList[pointCount++]  = modeCandidate_i;
            modeTable[modeCandidate_i] = 2;

          } else
          {

            // the mode has already been associated with
            // another mode, thererfore associate this one
            // mode and the modes in the point list with
            // the mode associated with data[basin_i]...

            // store the mode info into yk using msRawData...
            for ( j = 0; j < N; j++ )
              yk[j+2] = msRawData[modeCandidate_i*N+j];

            // update mode table for this data point
            // indicating that a mode has been associated
            // with it
            modeTable[i] = 1;

            // indicate that a mode has been associated
            // to this data point (data[i])
            mvAbs = -1;

            // stop mean shift calculation...
            break;
          }
        }
      }

      // Calculate the mean shift vector at the new
      // window location using lattice
      LatticeMSVector ( Mh, yk );

      // Calculate its magnitude squared
      mvAbs = 0;
      for ( j = 0; j < lN; j++ )
        mvAbs += Mh[j] * Mh[j];

      // Increment iteration count
      iterationCount++;

    }

    // if a mode was not associated with this data point
    // yet associate it with yk...
    if ( mvAbs >= 0 )
    {
      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // update mode table for this data point
      // indicating that a mode has been associated
      // with it
      modeTable[i] = 1;
    }

    // associate the data point indexed by
    // the point list with the mode stored
    // by yk
    for ( j = 0; j < pointCount; j++ )
    {
      // obtain the point location from the
      // point list
      modeCandidate_i = pointList[j];

      // update the mode table for this point
      modeTable[modeCandidate_i] = 1;

      //store result into msRawData...
      for ( k = 0; k < N; k++ )
        msRawData[N*modeCandidate_i+k] = ( float ) ( yk[k+2] );
    }


    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;
  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif

  // de-allocate memory
  delete [] modeCandidatePoint;
  delete [] yk;
  delete [] Mh;

  // done.
  return;

}

/*******************************************************/
/*Optimized Filter 2                                   */
/*******************************************************/
/*Performs mean shift filtering on the specified input */
/*image using a user defined kernel. Previous mode     */
/*information is used to avoid re-applying mean shift  */
/*on certain data points to improve performance. To    */
/*further improve perfmance (during segmentation) poi- */
/*nts within h of a window center during the window    */
/*center's traversal to a mode are associated with the */
/*mode that the window converges to.                   */
/*******************************************************/
/*Pre:                                                 */
/*      - the user defined kernel used to apply mean   */
/*        shift filtering to the defined input image   */
/*        has spatial bandwidth sigmaS and range band- */
/*        width sigmaR                                 */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - mean shift filtering has been applied to the */
/*        input image using a user defined kernel      */
/*      - the filtered image is stored in the private  */
/*        data members of the msImageProcessor class.  */
/*******************************************************/

void msImageProcessor::OptimizedFilter2 ( float sigmaS, float sigmaR )
{

  //if confidence map is null set it to zero
  if ( !weightMap )
  {
    weightMap = new float [L];
    int i;
    for ( i = 0; i < L; i++ )
      weightMap[i] = 0;
  }

  // Declare Variables
  int  iterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;
  float *modeCandidatePoint;
  double mvAbs, diff, el;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // Initialize mode table used for basin of attraction
  memset ( modeTable, 0, width*height );

  // Allocate memory mode candidate data point...
  //floating point version
  modeCandidatePoint = new float [N];

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice)... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif

// Parallelisieren !!!

  for ( i = 0; i < L; i++ )
  {
    // if a mode was already assigned to this data point
    // then skip this point, otherwise proceed to
    // find its mode by applying mean shift...
    if ( modeTable[i] == 1 )
      continue;

    // initialize point list...
    pointCount = 0;

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    yk[0] = i % width;
    yk[1] = i / width;
    for ( j = 0; j < N; j++ )
      yk[j+2] = data[N*i+j];

    // Calculate the mean shift vector using the lattice
    OptLatticeMSVector ( Mh, yk );

    // Calculate its magnitude squared
    mvAbs = 0;
    for ( j = 0; j < lN; j++ )
      mvAbs += Mh[j] * Mh[j];

    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // check to see if the current mode location is in the
      // basin of attraction...

      // calculate the location of yk on the lattice
      modeCandidateX = ( int ) ( yk[0] + 0.5 );
      modeCandidateY = ( int ) ( yk[1] + 0.5 );
      modeCandidate_i = modeCandidateY * width + modeCandidateX;

      // if mvAbs != 0 (yk did indeed move) then check
      // location basin_i in the mode table to see if
      // this data point either:

      // (1) has not been associated with a mode yet
      //     (modeTable[basin_i] = 0), so associate
      //     it with this one
      //
      // (2) it has been associated with a mode other
      //     than the one that this data point is converging
      //     to (modeTable[basin_i] = 1), so assign to
      //     this data point the same mode as that of basin_i

      if ( ( modeTable[modeCandidate_i] != 2 ) && ( modeCandidate_i != i ) )
      {
        // obtain the data point at basin_i to
        // see if it is within h*TC_DIST_FACTOR of
        // of yk
        for ( j = 0; j < N; j++ )
          modeCandidatePoint[j] = data[N*modeCandidate_i + j];

        // check basin on non-spatial data spaces only
        k = 1;
        s = 0;
        diff = 0;
        while ( ( diff < TC_DIST_FACTOR ) && ( k < kp ) )
        {
          diff = 0;
          for ( p = 0; p < P[k]; p++ )
          {
            el = ( modeCandidatePoint[p+s] - yk[p+s+2] ) / h[k];
            diff += el * el;
          }
          s += P[k];
          k++;
        }

        // if the data point at basin_i is within
        // a distance of h*TC_DIST_FACTOR of yk
        // then depending on modeTable[basin_i] perform
        // either (1) or (2)
        if ( diff < TC_DIST_FACTOR )
        {
          // if the data point at basin_i has not
          // been associated to a mode then associate
          // it with the mode that this one will converge
          // to
          if ( modeTable[modeCandidate_i] == 0 )
          {
            // no mode associated yet so associate
            // it with this one...
            pointList[pointCount++]  = modeCandidate_i;
            modeTable[modeCandidate_i] = 2;

          } else
          {

            // the mode has already been associated with
            // another mode, thererfore associate this one
            // mode and the modes in the point list with
            // the mode associated with data[basin_i]...

            // store the mode infor int yk using msRawData...
            for ( j = 0; j < N; j++ )
              yk[j+2] = msRawData[modeCandidate_i*N+j];

            // update mode table for this data point
            // indicating that a mode has been associated
            // with it
            modeTable[i] = 1;

            // indicate that a mode has been associated
            // to this data point (data[i])
            mvAbs = -1;

            // stop mean shift calculation...
            break;
          }
        }
      }

      // Calculate the mean shift vector at the new
      // window location using lattice
      OptLatticeMSVector ( Mh, yk );

      // Calculate its magnitude squared
      mvAbs = 0;
      for ( j = 0; j < lN; j++ )
        mvAbs += Mh[j] * Mh[j];

      // Increment interation count
      iterationCount++;

    }

    // if a mode was not associated with this data point
    // yet then perform a shift the window center yk one
    // last time using the mean shift vector...
    if ( mvAbs >= 0 )
    {
      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // update mode table for this data point
      // indicating that a mode has been associated
      // with it
      modeTable[i] = 1;
    }

    // associate the data point indexed by
    // the point list with the mode stored
    // by yk
    for ( j = 0; j < pointCount; j++ )
    {
      // obtain the point location from the
      // point list
      modeCandidate_i = pointList[j];

      // update the mode table for this point
      modeTable[modeCandidate_i] = 1;

      //store result into msRawData...
      for ( k = 0; k < N; k++ )
        msRawData[N*modeCandidate_i+k] = ( float ) ( yk[k+2] );
    }


    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;

  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif

  // de-allocate memory
  delete [] modeCandidatePoint;
  delete [] yk;
  delete [] Mh;

  // done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\/\*/
/* Image Classification */
/*\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Connect                                              */
/*******************************************************/
/*Classifies the regions of the mean shift filtered    */
/*image.                                               */
/*******************************************************/
/*Post:                                                */
/*      - the regions of the mean shift image have been*/
/*        classified using the private classification  */
/*        structure of the msImageProcessor Class.     */
/*        Namely, each region uniquely identified by   */
/*        its LUV color  (stored by LUV_data) and loc- */
/*        ation has been labeled and its area computed */
/*        via an eight-connected fill.                 */
/*******************************************************/

void msImageProcessor::Connect ( void )
{

  //define eight connected neighbors
  neigh[0] = 1;
  neigh[1] = 1 - width;
  neigh[2] = -width;
  neigh[3] = - ( 1 + width );
  neigh[4] = -1;
  neigh[5] = width - 1;
  neigh[6] = width;
  neigh[7] = width + 1;

  //initialize labels and modePointCounts
  int i;
  for ( i = 0; i < width*height; i++ )
  {
    labels[i]   = -1;
    modePointCounts[i] =  0;
  }

  //Traverse the image labeling each new region encountered
  int k, label = -1;
  for ( i = 0; i < height*width; i++ )
  {
    //if this region has not yet been labeled - label it
    if ( labels[i] < 0 )
    {
      //assign new label to this region
      labels[i] = ++label;

      //copy region color into modes
      for ( k = 0; k < N; k++ )
        modes[ ( N*label ) +k] = LUV_data[ ( N*i ) +k];
//    modes[(N*label)+k] = (float)(LUV_data[(N*i)+k]);

      //populate labels with label for this specified region
      //calculating modePointCounts[label]...
      Fill ( i, label );
    }
  }

  //calculate region count using label
  regionCount = label + 1;

  //done.
  return;
}

/*******************************************************/
/*Fill                                                 */
/*******************************************************/
/*Given a region seed and a region label, Fill uses    */
/*the region seed to perform an eight-connected fill   */
/*for the specified region, labeling all pixels con-   */
/*tained by the region with the specified label:       */
/*label.                                               */
/*******************************************************/
/*Pre:                                                 */
/*      - regionLoc is a region seed - a pixel that is */
/*        identified as being part of the region       */
/*        labled using the label, label.               */
/*Post:                                                */
/*      - all pixels belonging to the region specified */
/*        by regionLoc (having the same integer LUV    */
/*        value specified by LUV_data) are classified  */
/*        as one region by labeling each pixel in the  */
/*        image clasification structure using label    */
/*        via an eight-connected fill.                 */
/*******************************************************/

void msImageProcessor::Fill ( int regionLoc, int label )
{

  //declare variables
  int i, k, neighLoc, neighborsFound, imageSize = width * height;

  //Fill region starting at region location
  //using labels...

  //initialzie indexTable
  int index  = 0;
  indexTable[0] = regionLoc;

  //increment mode point counts for this region to
  //indicate that one pixel belongs to this region
  modePointCounts[label]++;

  while ( true )
  {

    //assume no neighbors will be found
    neighborsFound = 0;

    //check the eight connected neighbors at regionLoc -
    //if a pixel has similar color to that located at
    //regionLoc then declare it as part of this region
    for ( i = 0; i < 8; i++ )
    {
      // no need
      /*
      //if at boundary do not check certain neighbors because
      //they do not exist...
      if((regionLoc%width == 0)&&((i == 3)||(i == 4)||(i == 5)))
      continue;
      if((regionLoc%(width-1) == 0)&&((i == 0)||(i == 1)||(i == 7)))
      continue;
      */

      //check bounds and if neighbor has been already labeled
      neighLoc   = regionLoc + neigh[i];
      if ( ( neighLoc >= 0 ) && ( neighLoc < imageSize ) && ( labels[neighLoc] < 0 ) )
      {
        for ( k = 0; k < N; k++ )
        {
//     if(LUV_data[(regionLoc*N)+k] != LUV_data[(neighLoc*N)+k])
          if ( fabs ( LUV_data[ ( regionLoc*N ) +k] - LUV_data[ ( neighLoc*N ) +k] ) >= LUV_treshold )
            break;
        }

        //neighbor i belongs to this region so label it and
        //place it onto the index table buffer for further
        //processing
        if ( k == N )
        {
          //assign label to neighbor i
          labels[neighLoc] = label;

          //increment region point count
          modePointCounts[label]++;

          //place index of neighbor i onto the index tabel buffer
          indexTable[++index] = neighLoc;

          //indicate that a neighboring region pixel was
          //identified
          neighborsFound = 1;
        }
      }
    }

    //check the indexTable to see if there are any more
    //entries to be explored - if so explore them, otherwise
    //exit the loop - we are finished
    if ( neighborsFound )
      regionLoc = indexTable[index];
    else if ( index > 1 )
      regionLoc = indexTable[--index];
    else
      break; //fill complete
  }

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\*/
/*  Image Pruning */
/*\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Build Region Adjacency Matrix                        */
/*******************************************************/
/*Constructs a region adjacency matrix.                */
/*******************************************************/
/*Pre:                                                 */
/*      - the classification data structure has been   */
/*        constructed.                                 */
/*Post:                                                */
/*      - a region adjacency matrix has been built     */
/*        using the classification data structure.     */
/*******************************************************/

void msImageProcessor::BuildRAM ( void )
{

  //Allocate memory for region adjacency matrix if it hasn't already been allocated
  if ( ( !raList ) && ( ( ! ( raList = new RAList [regionCount] ) ) || ( ! ( raPool = new RAList [NODE_MULTIPLE*regionCount] ) ) ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Allocate", ( char* ) "Not enough memory." );
    return;
  }

  //initialize the region adjacency list
  int i;
  for ( i = 0; i < regionCount; i++ )
  {
    raList[i].edgeStrength  = 0;
    raList[i].edgePixelCount = 0;
    raList[i].label    = i;
    raList[i].next    = NULL;
  }

  //initialize RAM free list
  freeRAList = raPool;
  for ( i = 0; i < NODE_MULTIPLE*regionCount - 1; i++ )
  {
    raPool[i].edgeStrength  = 0;
    raPool[i].edgePixelCount = 0;
    raPool[i].next = &raPool[i+1];
  }
  raPool[NODE_MULTIPLE*regionCount-1].next = NULL;

  //traverse the labeled image building
  //the RAM by looking to the right of
  //and below the current pixel location thus
  //determining if a given region is adjacent
  //to another
  int  j, curLabel, rightLabel, bottomLabel, exists;
  RAList *raNode1, *raNode2, *oldRAFreeList;
  for ( i = 0; i < height - 1; i++ )
  {
    //check the right and below neighbors
    //for pixel locations whose x < width - 1
    for ( j = 0; j < width - 1; j++ )
    {
      //calculate pixel labels
      curLabel = labels[i*width+j    ]; //current pixel
      rightLabel = labels[i*width+j+1  ]; //right   pixel
      bottomLabel = labels[ ( i+1 ) *width+j]; //bottom  pixel

      //check to the right, if the label of
      //the right pixel is not the same as that
      //of the current one then region[j] and region[j+1]
      //are adjacent to one another - update the RAM
      if ( curLabel != rightLabel )
      {
        //obtain RAList object from region adjacency free
        //list
        raNode1   = freeRAList;
        raNode2   = freeRAList->next;

        //keep a pointer to the old region adj. free
        //list just in case nodes already exist in respective
        //region lists
        oldRAFreeList = freeRAList;

        //update region adjacency free list
        freeRAList  = freeRAList->next->next;

        //populate RAList nodes
        raNode1->label = curLabel;
        raNode2->label = rightLabel;

        //insert nodes into the RAM
        exists   = 0;
        raList[curLabel  ].Insert ( raNode2 );
        exists   = raList[rightLabel].Insert ( raNode1 );

        //if the node already exists then place
        //nodes back onto the region adjacency
        //free list
        if ( exists )
          freeRAList = oldRAFreeList;

      }

      //check below, if the label of
      //the bottom pixel is not the same as that
      //of the current one then region[j] and region[j+width]
      //are adjacent to one another - update the RAM
      if ( curLabel != bottomLabel )
      {
        //obtain RAList object from region adjacency free
        //list
        raNode1   = freeRAList;
        raNode2   = freeRAList->next;

        //keep a pointer to the old region adj. free
        //list just in case nodes already exist in respective
        //region lists
        oldRAFreeList = freeRAList;

        //update region adjacency free list
        freeRAList  = freeRAList->next->next;

        //populate RAList nodes
        raNode1->label = curLabel;
        raNode2->label = bottomLabel;

        //insert nodes into the RAM
        exists   = 0;
        raList[curLabel  ].Insert ( raNode2 );
        exists   = raList[bottomLabel].Insert ( raNode1 );

        //if the node already exists then place
        //nodes back onto the region adjacency
        //free list
        if ( exists )
          freeRAList = oldRAFreeList;

      }

    }

    //check only to the bottom neighbors of the right boundary
    //pixels...

    //calculate pixel locations (j = width-1)
    curLabel = labels[i*width+j    ]; //current pixel
    bottomLabel = labels[ ( i+1 ) *width+j]; //bottom  pixel

    //check below, if the label of
    //the bottom pixel is not the same as that
    //of the current one then region[j] and region[j+width]
    //are adjacent to one another - update the RAM
    if ( curLabel != bottomLabel )
    {
      //obtain RAList object from region adjacency free
      //list
      raNode1   = freeRAList;
      raNode2   = freeRAList->next;

      //keep a pointer to the old region adj. free
      //list just in case nodes already exist in respective
      //region lists
      oldRAFreeList = freeRAList;

      //update region adjacency free list
      freeRAList  = freeRAList->next->next;

      //populate RAList nodes
      raNode1->label = curLabel;
      raNode2->label = bottomLabel;

      //insert nodes into the RAM
      exists   = 0;
      raList[curLabel  ].Insert ( raNode2 );
      exists   = raList[bottomLabel].Insert ( raNode1 );

      //if the node already exists then place
      //nodes back onto the region adjacency
      //free list
      if ( exists )
        freeRAList = oldRAFreeList;

    }
  }

  //check only to the right neighbors of the bottom boundary
  //pixels...

  //check the right for pixel locations whose x < width - 1
  for ( j = 0; j < width - 1; j++ )
  {
    //calculate pixel labels (i = height-1)
    curLabel = labels[i*width+j    ]; //current pixel
    rightLabel = labels[i*width+j+1  ]; //right   pixel

    //check to the right, if the label of
    //the right pixel is not the same as that
    //of the current one then region[j] and region[j+1]
    //are adjacent to one another - update the RAM
    if ( curLabel != rightLabel )
    {
      //obtain RAList object from region adjacency free
      //list
      raNode1   = freeRAList;
      raNode2   = freeRAList->next;

      //keep a pointer to the old region adj. free
      //list just in case nodes already exist in respective
      //region lists
      oldRAFreeList = freeRAList;

      //update region adjacency free list
      freeRAList  = freeRAList->next->next;

      //populate RAList nodes
      raNode1->label = curLabel;
      raNode2->label = rightLabel;

      //insert nodes into the RAM
      exists   = 0;
      raList[curLabel  ].Insert ( raNode2 );
      exists   = raList[rightLabel].Insert ( raNode1 );

      //if the node already exists then place
      //nodes back onto the region adjacency
      //free list
      if ( exists )
        freeRAList = oldRAFreeList;

    }

  }

  //done.
  return;

}

/*******************************************************/
/*Destroy Region Adjacency Matrix                      */
/*******************************************************/
/*Destroy a region adjacency matrix.                   */
/*******************************************************/
/*Post:                                                */
/*      - the region adjacency matrix has been destr-  */
/*        oyed: (1) its memory has been de-allocated,  */
/*        (2) the RAM structure has been initialize    */
/*        for re-use.                                  */
/*******************************************************/

void msImageProcessor::DestroyRAM ( void )
{

  //de-allocate memory for region adjaceny list
  if ( raList )    delete [] raList;
  if ( raPool )    delete [] raPool;

  //initialize region adjacency matrix
  raList    = NULL;
  freeRAList   = NULL;
  raPool    = NULL;

  //done.
  return;

}

/*******************************************************/
/*Transitive Closure                                   */
/*******************************************************/
/*Applies transitive closure to the RAM updating       */
/*labels, modes and modePointCounts to reflect the new */
/*set of merged regions resulting from transitive clo- */
/*sure.                                                */
/*******************************************************/
/*Post:                                                */
/*      - transitive closure has been applied to the   */
/*        regions classified by the RAM and labels,    */
/*        modes and modePointCounts have been updated  */
/*        to reflect the new set of mergd regions res- */
/*        ulting from transitive closure.              */
/*******************************************************/

void msImageProcessor::TransitiveClosure ( void )
{

  //Step (1):

  // Build RAM using classifiction structure originally
  // generated by the method GridTable::Connect()
  BuildRAM();

  //Step (1a):
  //Compute weights of weight graph using confidence map
  //(if defined)
  if ( weightMapDefined ) ComputeEdgeStrengths();

  //Step (2):

  //Treat each region Ri as a disjoint set:

  // - attempt to join Ri and Rj for all i != j that are neighbors and
  //   whose associated modes are a normalized distance of < 0.5 from one
  //   another

  // - the label of each region in the raList is treated as a pointer to the
  //   canonical element of that region (e.g. raList[i], initially has raList[i].label = i,
  //   namely each region is initialized to have itself as its canonical element).

  //Traverse RAM attempting to join raList[i] with its neighbors...
  int  i, iCanEl, neighCanEl;
  float threshold;
  RAList *neighbor;
  for ( i = 0; i < regionCount; i++ )
  {
    //aquire first neighbor in region adjacency list pointed to
    //by raList[i]
    neighbor = raList[i].next;

    //compute edge strenght threshold using global and local
    //epsilon
    if ( epsilon > raList[i].edgeStrength )
      threshold   = epsilon;
    else
      threshold   = raList[i].edgeStrength;

    //traverse region adjacency list of region i, attempting to join
    //it with regions whose mode is a normalized distance < 0.5 from
    //that of region i...
    while ( neighbor )
    {
      //attempt to join region and neighbor...
      if ( ( InWindow ( i, neighbor->label ) ) && ( neighbor->edgeStrength < epsilon ) )
      {
        //region i and neighbor belong together so join them
        //by:

        // (1) find the canonical element of region i
        iCanEl  = i;
        while ( raList[iCanEl].label != iCanEl )
          iCanEl  = raList[iCanEl].label;

        // (2) find the canonical element of neighboring region
        neighCanEl = neighbor->label;
        while ( raList[neighCanEl].label != neighCanEl )
          neighCanEl = raList[neighCanEl].label;

        // if the canonical elements of are not the same then assign
        // the canonical element having the smaller label to be the parent
        // of the other region...
        if ( iCanEl < neighCanEl )
          raList[neighCanEl].label = iCanEl;
        else
        {
          //must replace the canonical element of previous
          //parent as well
          raList[raList[iCanEl].label].label = neighCanEl;

          //re-assign canonical element
          raList[iCanEl].label    = neighCanEl;
        }
      }

      //check the next neighbor...
      neighbor = neighbor->next;

    }
  }

  // Step (3):

  // Level binary trees formed by canonical elements
  for ( i = 0; i < regionCount; i++ )
  {
    iCanEl = i;
    while ( raList[iCanEl].label != iCanEl )
      iCanEl = raList[iCanEl].label;
    raList[i].label = iCanEl;
  }

  // Step (4):

  //Traverse joint sets, relabeling image.

  // (a)

  // Accumulate modes and re-compute point counts using canonical
  // elements generated by step 2.

  //allocate memory for mode and point count temporary buffers...
  float *modes_buffer = new float [N*regionCount];
  int  *MPC_buffer  = new int [regionCount];

  //initialize buffers to zero
  for ( i = 0; i < regionCount; i++ )
    MPC_buffer[i] = 0;
  for ( i = 0; i < N*regionCount; i++ )
    modes_buffer[i] = 0;

  //traverse raList accumulating modes and point counts
  //using canoncial element information...
  int k, iMPC;
  for ( i = 0; i < regionCount; i++ )
  {

    //obtain canonical element of region i
    iCanEl = raList[i].label;

    //obtain mode point count of region i
    iMPC = modePointCounts[i];

    //accumulate modes_buffer[iCanEl]
    for ( k = 0; k < N; k++ )
      modes_buffer[ ( N*iCanEl ) +k] += iMPC * modes[ ( N*i ) +k];

    //accumulate MPC_buffer[iCanEl]
    MPC_buffer[iCanEl] += iMPC;

  }

  // (b)

  // Re-label new regions of the image using the canonical
  // element information generated by step (2)

  // Also use this information to compute the modes of the newly
  // defined regions, and to assign new region point counts in
  // a consecute manner to the modePointCounts array

  //allocate memory for label buffer
  int *label_buffer = new int [regionCount];

  //initialize label buffer to -1
  for ( i = 0; i < regionCount; i++ )
    label_buffer[i] = -1;

  //traverse raList re-labeling the regions
  int label = -1;
  for ( i = 0; i < regionCount; i++ )
  {
    //obtain canonical element of region i
    iCanEl = raList[i].label;
    if ( label_buffer[iCanEl] < 0 )
    {
      //assign a label to the new region indicated by canonical
      //element of i
      label_buffer[iCanEl] = ++label;

      //recompute mode storing the result in modes[label]...
      iMPC = MPC_buffer[iCanEl];
      for ( k = 0; k < N; k++ )
        modes[ ( N*label ) +k] = ( modes_buffer[ ( N*iCanEl ) +k] ) / ( iMPC );

      //assign a corresponding mode point count for this region into
      //the mode point counts array using the MPC buffer...
      modePointCounts[label] = MPC_buffer[iCanEl];
    }
  }

  //re-assign region count using label counter
  //int oldRegionCount = regionCount;
  regionCount = label + 1;

  // (c)

  // Use the label buffer to reconstruct the label map, which specified
  // the new image given its new regions calculated above

  for ( i = 0; i < height*width; i++ )
    labels[i] = label_buffer[raList[labels[i]].label];

  //de-allocate memory
  delete [] modes_buffer;
  delete [] MPC_buffer;
  delete [] label_buffer;

  //done.
  return;

}

/*******************************************************/
/*Compute Edge Strengths                               */
/*******************************************************/
/*Computes the a weight for each link in the region    */
/*graph maintined by the RAM, resulting in a weighted  */
/*graph in which the weights consist of a confidence   */
/*between zero and one indicating if the regions are   */
/*separated by a strong or weak edge.                  */
/*******************************************************/
/*Post:                                                */
/*      - an edge strength has been computed between   */
/*        each region of the image and placed as a     */
/*        weight in the RAM to be used during transi-  */
/*        tive closure.                                */
/*******************************************************/

void msImageProcessor::ComputeEdgeStrengths ( void )
{

  //initialize visit table - used to keep track
  //of which pixels have already been visited such
  //as not to contribute their strength value to
  //a boundary sum multiple times...
  memset ( visitTable, 0, L*sizeof ( unsigned char ) );

  //traverse labeled image computing edge strengths
  //(excluding image boundary)...
  int    x, y, dp, curLabel, rightLabel, bottomLabel;
  RAList *curRegion;
  for ( y = 1; y < height - 1; y++ )
  {
    for ( x = 1; x < width - 1; x++ )
    {
      //compute data point location using x and y
      dp = y * width + x;

      //obtain labels at different pixel locations
      curLabel = labels[dp      ]; //current pixel
      rightLabel = labels[dp+1    ]; //right   pixel
      bottomLabel = labels[dp+width]; //bottom  pixel

      //check right and bottom neighbor to see if there is a
      //change in label then we are at an edge therefore record
      //the edge strength at this edge accumulating its value
      //in the RAM...
      if ( curLabel != rightLabel )
      {
        //traverse into RAM...
        curRegion = &raList[curLabel];
        while ( ( curRegion ) && ( curRegion->label != rightLabel ) )
          curRegion = curRegion->next;

        //this should not occur...
        assert ( curRegion );

        //accumulate edge strength
        curRegion->edgeStrength   += weightMap[dp] + weightMap[dp+1];
        curRegion->edgePixelCount += 2;
      }

      if ( curLabel != bottomLabel )
      {
        //traverse into RAM...
        curRegion = &raList[curLabel];
        while ( ( curRegion ) && ( curRegion->label != bottomLabel ) )
          curRegion = curRegion->next;

        //this should not occur...
        assert ( curRegion );

        //accumulate edge strength
        if ( curLabel == rightLabel )
        {
          curRegion->edgeStrength   += weightMap[dp] + weightMap[dp+width];
          curRegion->edgePixelCount += 2;
        }
        else
        {
          curRegion->edgeStrength   += weightMap[dp+width];
          curRegion->edgePixelCount += 1;
        }

      }
    }
  }

  //compute strengths using accumulated strengths obtained above...
  RAList *neighborRegion;
  float edgeStrength;
  int  edgePixelCount;
  for ( x = 0; x < regionCount; x++ )
  {
    //traverse the region list of the current region
    curRegion = &raList[x];
    curRegion = curRegion->next;
    while ( curRegion )
    {
      //with the assumption that regions having a smaller
      //label in the current region list have already
      //had their edge strengths computed, only compute
      //edge strengths for the regions whose label is greater
      //than x, the current region (region list) under
      //consideration...
      curLabel = curRegion->label;
      if ( curLabel > x )
      {
        //obtain pointer to the element identifying the
        //current region in the neighbors region list...
        neighborRegion = &raList[curLabel];
        while ( ( neighborRegion ) && ( neighborRegion->label != x ) )
          neighborRegion = neighborRegion->next;

        //this should not occur...
        assert ( neighborRegion );

        //compute edge strengths using accumulated confidence
        //value and pixel count
        if ( ( edgePixelCount = curRegion->edgePixelCount + neighborRegion->edgePixelCount ) != 0 )
        {
          //compute edge strength
          edgeStrength = curRegion->edgeStrength + neighborRegion->edgeStrength;
          edgeStrength /= edgePixelCount;

          //store edge strength and pixel count for corresponding regions
          curRegion->edgeStrength  = neighborRegion->edgeStrength  = edgeStrength;
          curRegion->edgePixelCount = neighborRegion->edgePixelCount = edgePixelCount;
        }
      }

      //traverse to the next region in the region adjacency list
      //of the current region x
      curRegion = curRegion->next;

    }
  }

  //compute average edge strength amongst the edges connecting
  //it to each of its neighbors
  int numNeighbors;
  for ( x = 0; x < regionCount; x++ )
  {
    //traverse the region list of the current region
    //accumulating weights
    curRegion  = &raList[x];
    curRegion  = curRegion->next;
    edgeStrength = 0;
    numNeighbors = 0;
    while ( curRegion )
    {
      numNeighbors++;
      edgeStrength   += curRegion->edgeStrength;
      curRegion  = curRegion->next;
    }

    //divide by the number of regions connected
    //to the current region
    if ( numNeighbors ) edgeStrength /= numNeighbors;

    //store the result in the raList for region
    //x
    raList[x].edgeStrength = edgeStrength;
  }

  //traverse raList and output the resulting list
  //to a file

  //done.
  return;

}

/*******************************************************/
/*Prune                                                */
/*******************************************************/
/*Prunes regions from the image whose pixel density    */
/*is less than a specified threshold.                  */
/*******************************************************/
/*Pre:                                                 */
/*      - minRegion is the minimum allowable pixel de- */
/*        nsity a region may have without being pruned */
/*        from the image                               */
/*Post:                                                */
/*      - regions whose pixel density is less than     */
/*        or equal to minRegion have been pruned from  */
/*        the image.                                   */
/*******************************************************/

void msImageProcessor::Prune ( int minRegion )
{

  //Allocate Memory for temporary buffers...

  //allocate memory for mode and point count temporary buffers...
  float *modes_buffer = new float [N*regionCount];
  int  *MPC_buffer  = new int [regionCount];

  //allocate memory for label buffer
  int *label_buffer  = new int [regionCount];

  //Declare variables
  int  i, k, candidate, iCanEl, neighCanEl, iMPC, label, oldRegionCount, minRegionCount;
  double minSqDistance, neighborDistance;
  RAList *neighbor;

  //Apply pruning algorithm to classification structure, removing all regions whose area
  //is under the threshold area minRegion (pixels)
  do
  {
    //Assume that no region has area under threshold area  of
    minRegionCount = 0;

    //Step (1):

    // Build RAM using classifiction structure originally
    // generated by the method GridTable::Connect()
    BuildRAM();

    // Step (2):

    // Traverse the RAM joining regions whose area is less than minRegion (pixels)
    // with its respective candidate region.

    // A candidate region is a region that displays the following properties:

    // - it is adjacent to the region being pruned

    //  - the distance of its mode is a minimum to that of the region being pruned
    //    such that or it is the only adjacent region having an area greater than
    //    minRegion

    for ( i = 0; i < regionCount; i++ )
    {
      //if the area of the ith region is less than minRegion
      //join it with its candidate region...

      //*******************************************************************************

      //Note: Adjust this if statement if a more sophisticated pruning criterion
      //      is desired. Basically in this step a region whose area is less than
      //      minRegion is pruned by joining it with its "closest" neighbor (in color).
      //      Therefore, by placing a different criterion for fusing a region the
      //      pruning method may be altered to implement a more sophisticated algorithm.

      //*******************************************************************************

      if ( modePointCounts[i] < minRegion )
      {
        //update minRegionCount to indicate that a region
        //having area less than minRegion was found
        minRegionCount++;

        //obtain a pointer to the first region in the
        //region adjacency list of the ith region...
        neighbor = raList[i].next;

        //calculate the distance between the mode of the ith
        //region and that of the neighboring region...
        candidate  = neighbor->label;
        minSqDistance = SqDistance ( i, candidate );

        //traverse region adjacency list of region i and select
        //a candidate region
        neighbor = neighbor->next;
        while ( neighbor )
        {

          //calculate the square distance between region i
          //and current neighbor...
          neighborDistance = SqDistance ( i, neighbor->label );

          //if this neighbors square distance to region i is less
          //than minSqDistance, then select this neighbor as the
          //candidate region for region i
          if ( neighborDistance < minSqDistance )
          {
            minSqDistance = neighborDistance;
            candidate  = neighbor->label;
          }

          //traverse region list of region i
          neighbor = neighbor->next;

        }

        //join region i with its candidate region:

        // (1) find the canonical element of region i
        iCanEl  = i;
        while ( raList[iCanEl].label != iCanEl )
          iCanEl  = raList[iCanEl].label;

        // (2) find the canonical element of neighboring region
        neighCanEl = candidate;
        while ( raList[neighCanEl].label != neighCanEl )
          neighCanEl = raList[neighCanEl].label;

        // if the canonical elements of are not the same then assign
        // the canonical element having the smaller label to be the parent
        // of the other region...
        if ( iCanEl < neighCanEl )
          raList[neighCanEl].label = iCanEl;
        else
        {
          //must replace the canonical element of previous
          //parent as well
          raList[raList[iCanEl].label].label = neighCanEl;

          //re-assign canonical element
          raList[iCanEl].label    = neighCanEl;
        }
      }
    }

    // Step (3):

    // Level binary trees formed by canonical elements
    for ( i = 0; i < regionCount; i++ )
    {
      iCanEl = i;
      while ( raList[iCanEl].label != iCanEl )
        iCanEl = raList[iCanEl].label;
      raList[i].label = iCanEl;
    }

    // Step (4):

    //Traverse joint sets, relabeling image.

    // Accumulate modes and re-compute point counts using canonical
    // elements generated by step 2.

    //initialize buffers to zero
    for ( i = 0; i < regionCount; i++ )
      MPC_buffer[i] = 0;
    for ( i = 0; i < N*regionCount; i++ )
      modes_buffer[i] = 0;

    //traverse raList accumulating modes and point counts
    //using canoncial element information...
    for ( i = 0; i < regionCount; i++ )
    {

      //obtain canonical element of region i
      iCanEl = raList[i].label;

      //obtain mode point count of region i
      iMPC = modePointCounts[i];

      //accumulate modes_buffer[iCanEl]
      for ( k = 0; k < N; k++ )
        modes_buffer[ ( N*iCanEl ) +k] += iMPC * modes[ ( N*i ) +k];

      //accumulate MPC_buffer[iCanEl]
      MPC_buffer[iCanEl] += iMPC;

    }

    // (b)

    // Re-label new regions of the image using the canonical
    // element information generated by step (2)

    // Also use this information to compute the modes of the newly
    // defined regions, and to assign new region point counts in
    // a consecute manner to the modePointCounts array

    //initialize label buffer to -1
    for ( i = 0; i < regionCount; i++ )
      label_buffer[i] = -1;

    //traverse raList re-labeling the regions
    label = -1;
    for ( i = 0; i < regionCount; i++ )
    {
      //obtain canonical element of region i
      iCanEl = raList[i].label;
      if ( label_buffer[iCanEl] < 0 )
      {
        //assign a label to the new region indicated by canonical
        //element of i
        label_buffer[iCanEl] = ++label;

        //recompute mode storing the result in modes[label]...
        iMPC = MPC_buffer[iCanEl];
        for ( k = 0; k < N; k++ )
          modes[ ( N*label ) +k] = ( modes_buffer[ ( N*iCanEl ) +k] ) / ( iMPC );

        //assign a corresponding mode point count for this region into
        //the mode point counts array using the MPC buffer...
        modePointCounts[label] = MPC_buffer[iCanEl];
      }
    }

    //re-assign region count using label counter
    oldRegionCount = regionCount;
    regionCount  = label + 1;

    // (c)

    // Use the label buffer to reconstruct the label map, which specified
    // the new image given its new regions calculated above

    for ( i = 0; i < height*width; i++ )
      labels[i] = label_buffer[raList[labels[i]].label];


  } while ( minRegionCount > 0 );

  //de-allocate memory
  delete [] modes_buffer;
  delete [] MPC_buffer;
  delete [] label_buffer;

  //done.
  return;

}

/*******************************************************/
/*Define Boundaries                                    */
/*******************************************************/
/*Defines the boundaries for each region of the segm-  */
/*ented image storing the result into a region list    */
/*object.                                              */
/*******************************************************/
/*Pre:                                                 */
/*      - the image has been segmented and a classifi- */
/*        cation structure has been created for this   */
/*        image                                        */
/*Post:                                                */
/*      - the boundaries of the segmented image have   */
/*        been defined and the boundaries of each reg- */
/*        ion has been stored into a region list obj-  */
/*        ect.                                         */
/*******************************************************/

void msImageProcessor::DefineBoundaries ( void )
{

  //declare and allocate memory for boundary map and count
  int *boundaryMap, *boundaryCount;
  if ( ( ! ( boundaryMap = new int [L] ) ) || ( ! ( boundaryCount = new int [regionCount] ) ) )
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "DefineBoundaries", ( char* ) "Not enough memory." );

  //initialize boundary map and count
  int i;
  for ( i = 0; i < L; i++ )
    boundaryMap[i]  = -1;
  for ( i = 0; i < regionCount; i++ )
    boundaryCount[i] =  0;

  //initialize and declare total boundary count -
  //the total number of boundary pixels present in
  //the segmented image
  int totalBoundaryCount = 0;

  //traverse the image checking the right and bottom
  //four connected neighbors of each pixel marking
  //boundary map with the boundaries of each region and
  //incrementing boundaryCount using the label information

  //***********************************************************************
  //***********************************************************************

  int  j, label, dataPoint;

  //first row (every pixel is a boundary pixel)
  for ( i = 0; i < width; i++ )
  {
    boundaryMap[i]  = label = labels[i];
    boundaryCount[label]++;
    totalBoundaryCount++;
  }

  //define boundaries for all rows except for the first
  //and last one...
  for ( i = 1; i < height - 1; i++ )
  {
    //mark the first pixel in an image row as an image boundary...
    dataPoint    = i * width;
    boundaryMap[dataPoint] = label = labels[dataPoint];
    boundaryCount[label]++;
    totalBoundaryCount++;

    for ( j = 1; j < width - 1; j++ )
    {
      //define datapoint and its right and bottom
      //four connected neighbors
      dataPoint  = i * width + j;

      //check four connected neighbors if they are
      //different this pixel is a boundary pixel
      label = labels[dataPoint];
      if ( ( label != labels[dataPoint-1] )    || ( label != labels[dataPoint+1] ) ||
           ( label != labels[dataPoint-width] ) || ( label != labels[dataPoint+width] ) )
      {
        boundaryMap[dataPoint]  = label = labels[dataPoint];
        boundaryCount[label]++;
        totalBoundaryCount++;
      }
    }

    //mark the last pixel in an image row as an image boundary...
    dataPoint    = ( i + 1 ) * width - 1;
    boundaryMap[dataPoint] = label = labels[dataPoint];
    boundaryCount[label]++;
    totalBoundaryCount++;

  }

  //last row (every pixel is a boundary pixel) (i = height-1)
  register int start = ( height - 1 ) * width, stop = height * width;
  for ( i = start; i < stop; i++ )
  {
    boundaryMap[i]  = label = labels[i];
    boundaryCount[label]++;
    totalBoundaryCount++;
  }

  //***********************************************************************
  //***********************************************************************

  //store boundary locations into a boundary buffer using
  //boundary map and count

  //***********************************************************************
  //***********************************************************************

  int *boundaryBuffer = new int [totalBoundaryCount], *boundaryIndex = new int [regionCount];

  //use boundary count to initialize boundary index...
  int counter = 0;
  for ( i = 0; i < regionCount; i++ )
  {
    boundaryIndex[i] = counter;
    counter      += boundaryCount[i];
  }

  //traverse boundary map placing the boundary pixel
  //locations into the boundaryBuffer
  for ( i = 0; i < L; i++ )
  {
    //if its a boundary pixel store it into
    //the boundary buffer
    if ( ( label = boundaryMap[i] ) >= 0 )
    {
      boundaryBuffer[boundaryIndex[label]] = i;
      boundaryIndex[label]++;
    }
  }

  //***********************************************************************
  //***********************************************************************

  //store the boundary locations stored by boundaryBuffer into
  //the region list for each region

  //***********************************************************************
  //***********************************************************************

  //destroy the old region list
  if ( regionList ) delete regionList;

  //create a new region list
  if ( ! ( regionList = new RegionList ( regionCount, totalBoundaryCount, N ) ) )
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "DefineBoundaries", ( char* ) "Not enough memory." );

  //add boundary locations for each region using the boundary
  //buffer and boundary counts
  counter = 0;
  for ( i = 0; i < regionCount; i++ )
  {
    regionList->AddRegion ( i, boundaryCount[i], &boundaryBuffer[counter] );
    counter += boundaryCount[i];
  }

  //***********************************************************************
  //***********************************************************************

  // dealocate local used memory
  delete [] boundaryMap;
  delete [] boundaryCount;
  delete [] boundaryBuffer;
  delete [] boundaryIndex;

  //done.
  return;

}

/*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/
/*  Image Data Searching/Distance Calculation */
/*\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*In Window                                            */
/*******************************************************/
/*Returns true if the two specified data points are    */
/*within rR of each other.                             */
/*******************************************************/
/*Pre:                                                 */
/*      - mode1 and mode2 are indeces into msRawData   */
/*        specifying the modes of the pixels having    */
/*        these indeces.                               */
/*Post:                                                */
/*      - true is returned if mode1 and mode2 are wi-  */
/*        thin rR of one another, false is returned    */
/*        otherwise.                                   */
/*******************************************************/

bool msImageProcessor::InWindow ( int mode1, int mode2 )
{
  int  k  = 1, s = 0, p;
  double diff = 0, el;
  while ( ( diff < 0.25 ) && ( k != kp ) ) // Partial Distortion Search
  {
    //Calculate distance squared of sub-space s
    diff = 0;
    for ( p = 0; p < P[k]; p++ )
    {
      el    = ( modes[mode1*N+p+s] - modes[mode2*N+p+s] ) / ( h[k] * offset[k] );
      if ( ( !p ) && ( k == 1 ) && ( modes[mode1*N] > 80 ) )
        diff += 4 * el * el;
      else
        diff += el * el;
    }

    //next subspace
    s += P[k];
    k++;
  }
  return ( bool ) ( diff < 0.25 );
}

/*******************************************************/
/*Square Distance                                      */
/*******************************************************/
/*Computs the normalized square distance between two   */
/*modes.                                               */
/*******************************************************/
/*Pre:                                                 */
/*      - mode1 and mode2 are indeces into the modes   */
/*        array specifying two modes of the image      */
/*Post:                                                */
/*      - the normalized square distance between modes */
/*        indexed by mode1 and mode2 has been calc-    */
/*        ulated and the result has been returned.     */
/*******************************************************/

float msImageProcessor::SqDistance ( int mode1, int mode2 )
{

  int  k  = 1, s = 0, p;
  float dist = 0, el;
  for ( k = 1; k < kp; k++ )
  {
    //Calculate distance squared of sub-space s
    for ( p = 0; p < P[k]; p++ )
    {
      el    = ( modes[mode1*N+p+s] - modes[mode2*N+p+s] ) / ( h[k] * offset[k] );
      dist += el * el;
    }

    //next subspace
    s += P[k];
    k++;
  }

  //return normalized square distance between modes
  //1 and 2
  return dist;

}

/*/\/\/\/\/\/\/\/\/\/\*/
/*  Memory Management */
/*\/\/\/\/\/\/\/\/\/\/*/

/*******************************************************/
/*Initialize Output                                    */
/*******************************************************/
/*Allocates memory needed by the mean shift image pro- */
/*cessor class output storage data structure.          */
/*******************************************************/
/*Post:                                                */
/*      - the memory needed by the output storage      */
/*        structure of this class has been (re-)allo-  */
/*        cated.                                       */
/*******************************************************/

void msImageProcessor::InitializeOutput ( void )
{

  //De-allocate memory if output was defined for previous image
  DestroyOutput();

  //Allocate memory for msRawData (filtered image output)
  if ( ! ( msRawData = new float [L*N] ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Allocate", ( char* ) "Not enough memory." );
    return;
  }

  //Allocate memory used to store image modes and their corresponding regions...
  if ( ( ! ( modes = new float [L* ( N+2 ) ] ) ) || ( ! ( labels = new int [L] ) ) || ( ! ( modePointCounts = new int [L] ) ) || ( ! ( indexTable = new int [L] ) ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Allocate", ( char* ) "Not enough memory" );
    return;
  }

  //Allocate memory for integer modes used to perform connected components
  //(image labeling)...
// if(!(LUV_data = new int [N*L]))
  if ( ! ( LUV_data = new float[N*L] ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Allocate", ( char* ) "Not enough memory" );
    return;
  }

  //indicate that the class output storage structure has been defined
  class_state.OUTPUT_DEFINED = true;

}

/*******************************************************/
/*Destroy Output                                       */
/*******************************************************/
/*De-allocates memory needed by the mean shift image   */
/*processor class output storage data structure.       */
/*******************************************************/
/*Post:                                                */
/*      - the memory needed by the output storage      */
/*        structure of this class has been de-alloc-   */
/*        ated.                                        */
/*      - the output storage structure has been init-  */
/*        ialized for re-use.                          */
/*******************************************************/

void msImageProcessor::DestroyOutput ( void )
{

  //de-allocate memory for msRawData (filtered image output)
  if ( msRawData )   delete [] msRawData;

  //de-allocate memory used by output storage and image
  //classification structure
  if ( modes )    delete [] modes;
  if ( labels )    delete [] labels;
  if ( modePointCounts ) delete [] modePointCounts;
  if ( indexTable )   delete [] indexTable;

  //de-allocate memory for LUV_data
  if ( LUV_data )   delete [] LUV_data;

  //initialize data members for re-use...

  //initialize output structures...
  msRawData   = NULL;

  //re-initialize classification structure
  modes      = NULL;
  labels      = NULL;
  modePointCounts    = NULL;
  regionCount     = 0;

  //indicate that the output has been destroyed
  class_state.OUTPUT_DEFINED = false;

  //done.
  return;

}

// NEW
void msImageProcessor::NewOptimizedFilter1 ( float sigmaS, float sigmaR )
{
  // Declare Variables
  int  iterationCount, i, j, k, modeCandidateX, modeCandidateY, modeCandidate_i;
  double mvAbs, diff, el;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // let's use some temporary data
  float* sdata;
  sdata = new float[lN*L];

  // copy the scaled data
  int idxs, idxd;
  idxs = idxd = 0;
  if ( N == 3 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else if ( N == 1 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      for ( j = 0; j < N; j++ )
        sdata[idxs++] = data[idxd++] / sigmaR;
    }
  }
  // index the data in the 3d buckets (x, y, L)
  int* buckets;
  int* slist;
  slist = new int[L];
  int bucNeigh[27];

  float sMins; // just for L
  float sMaxs[3]; // for all
  sMaxs[0] = width / sigmaS;
  sMaxs[1] = height / sigmaS;
  sMins = sMaxs[2] = sdata[2];
  idxs = 2;
  float cval;
  for ( i = 0; i < L; i++ )
  {
    cval = sdata[idxs];
    if ( cval < sMins )
      sMins = cval;
    else if ( cval > sMaxs[2] )
      sMaxs[2] = cval;

    idxs += lN;
  }

  int nBuck1, nBuck2, nBuck3;
  int cBuck1, cBuck2, cBuck3, cBuck;
  nBuck1 = ( int ) ( sMaxs[0] + 3 );
  nBuck2 = ( int ) ( sMaxs[1] + 3 );
  nBuck3 = ( int ) ( sMaxs[2] - sMins + 3 );
  buckets = new int[nBuck1*nBuck2*nBuck3];

  for ( i = 0; i < ( nBuck1*nBuck2*nBuck3 ); i++ )
    buckets[i] = -1;

  idxs = 0;

  for ( i = 0; i < L; i++ )
  {
    // find bucket for current data and add it to the list
    cBuck1 = ( int )  sdata[idxs] + 1;
    cBuck2 = ( int )  sdata[idxs+1] + 1;
    cBuck3 = ( int ) ( sdata[idxs+2] - sMins ) + 1;

    idxs += lN;

    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );

    slist[i] = buckets[cBuck];
    buckets[cBuck] = i;
  }
  // init bucNeigh

  idxd = 0;
  for ( cBuck1 = -1; cBuck1 <= 1; cBuck1++ )
  {
    for ( cBuck2 = -1; cBuck2 <= 1; cBuck2++ )
    {
      for ( cBuck3 = -1; cBuck3 <= 1; cBuck3++ )
      {
        bucNeigh[idxd++] = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      }
    }
  }

  double wsuml, weight;
  double hiLTr = 80.0 / sigmaR;
  // done indexing/hashing


  // Initialize mode table used for basin of attraction
  memset ( modeTable, 0, width*height );

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice) ... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif


  for ( i = 0; i < L; i++ )
  {
    // if a mode was already assigned to this data point
    // then skip this point, otherwise proceed to
    // find its mode by applying mean shift...
    if ( modeTable[i] == 1 )
      continue;

    // initialize point list...
    pointCount = 0;

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    idxs = i * lN;
    for ( j = 0; j < lN; j++ )
      yk[j] = sdata[idxs+j];

    // Calculate the mean shift vector using the lattice
    // LatticeMSVector(Mh, yk); // modify to new
    /*****************************************************/
    // Initialize mean shift vector

    for ( j = 0; j < lN; j++ )
      Mh[j] = 0;
    wsuml = 0;
    // uniformLSearch(Mh, yk_ptr); // modify to new
    // find bucket of yk
    cBuck1 = ( int ) yk[0] + 1;
    cBuck2 = ( int ) yk[1] + 1;
    cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
    for ( j = 0; j < 27; j++ )
    {
      idxd = buckets[cBuck+bucNeigh[j]];
      // list parse, crt point is cHeadList
      while ( idxd >= 0 )
      {
        idxs = lN * idxd;
        // determine if inside search window
        el = sdata[idxs+0] - yk[0];
        diff = el * el;
        el = sdata[idxs+1] - yk[1];
        diff += el * el;

        if ( diff < 1.0 )
        {
          el = sdata[idxs+2] - yk[2];
          if ( yk[2] > hiLTr )
            diff = 4 * el * el;
          else
            diff = el * el;

          if ( N > 1 )
          {
            el = sdata[idxs+3] - yk[3];
            diff += el * el;
            el = sdata[idxs+4] - yk[4];
            diff += el * el;
          }

          if ( diff < 1.0 )
          {
            weight = 1 - weightMap[idxd];
            for ( k = 0; k < lN; k++ )
              Mh[k] += weight * sdata[idxs+k];
            wsuml += weight;
          }
        }
        idxd = slist[idxd];
      }
    }
    if ( wsuml > 0 )
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = Mh[j] / wsuml - yk[j];
    }
    else
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
    }
    /*****************************************************/
    // Calculate its magnitude squared
    //mvAbs = 0;
    //for(j = 0; j < lN; j++)
    // mvAbs += Mh[j]*Mh[j];
    mvAbs = ( Mh[0] * Mh[0] + Mh[1] * Mh[1] ) * sigmaS * sigmaS;
    if ( N == 3 )
      mvAbs += ( Mh[2] * Mh[2] + Mh[3] * Mh[3] + Mh[4] * Mh[4] ) * sigmaR * sigmaR;
    else
      mvAbs += Mh[2] * Mh[2] * sigmaR * sigmaR;


    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // check to see if the current mode location is in the
      // basin of attraction...

      // calculate the location of yk on the lattice
      modeCandidateX = ( int ) ( sigmaS * yk[0] + 0.5 );
      modeCandidateY = ( int ) ( sigmaS * yk[1] + 0.5 );
      modeCandidate_i = modeCandidateY * width + modeCandidateX;

      // if mvAbs != 0 (yk did indeed move) then check
      // location basin_i in the mode table to see if
      // this data point either:

      // (1) has not been associated with a mode yet
      //     (modeTable[basin_i] = 0), so associate
      //     it with this one
      //
      // (2) it has been associated with a mode other
      //     than the one that this data point is converging
      //     to (modeTable[basin_i] = 1), so assign to
      //     this data point the same mode as that of basin_i

      if ( ( modeTable[modeCandidate_i] != 2 ) && ( modeCandidate_i != i ) )
      {
        // obtain the data point at basin_i to
        // see if it is within h*TC_DIST_FACTOR of
        // of yk
        diff = 0;
        idxs = lN * modeCandidate_i;
        for ( k = 2; k < lN; k++ )
        {
          el = sdata[idxs+k] - yk[k];
          diff += el * el;
        }

        // if the data point at basin_i is within
        // a distance of h*TC_DIST_FACTOR of yk
        // then depending on modeTable[basin_i] perform
        // either (1) or (2)
        if ( diff < TC_DIST_FACTOR )
        {
          // if the data point at basin_i has not
          // been associated to a mode then associate
          // it with the mode that this one will converge
          // to
          if ( modeTable[modeCandidate_i] == 0 )
          {
            // no mode associated yet so associate
            // it with this one...
            pointList[pointCount++]  = modeCandidate_i;
            modeTable[modeCandidate_i] = 2;

          } else
          {

            // the mode has already been associated with
            // another mode, thererfore associate this one
            // mode and the modes in the point list with
            // the mode associated with data[basin_i]...

            // store the mode info into yk using msRawData...
            for ( j = 0; j < N; j++ )
              yk[j+2] = msRawData[modeCandidate_i*N+j] / sigmaR;

            // update mode table for this data point
            // indicating that a mode has been associated
            // with it
            modeTable[i] = 1;

            // indicate that a mode has been associated
            // to this data point (data[i])
            mvAbs = -1;

            // stop mean shift calculation...
            break;
          }
        }
      }

      // Calculate the mean shift vector at the new
      // window location using lattice
      // Calculate the mean shift vector using the lattice
      // LatticeMSVector(Mh, yk); // modify to new
      /*****************************************************/
      // Initialize mean shift vector
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
      wsuml = 0;
      // uniformLSearch(Mh, yk_ptr); // modify to new
      // find bucket of yk
      cBuck1 = ( int ) yk[0] + 1;
      cBuck2 = ( int ) yk[1] + 1;
      cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
      cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      for ( j = 0; j < 27; j++ )
      {
        idxd = buckets[cBuck+bucNeigh[j]];
        // list parse, crt point is cHeadList
        while ( idxd >= 0 )
        {
          idxs = lN * idxd;
          // determine if inside search window
          el = sdata[idxs+0] - yk[0];
          diff = el * el;
          el = sdata[idxs+1] - yk[1];
          diff += el * el;

          if ( diff < 1.0 )
          {
            el = sdata[idxs+2] - yk[2];
            if ( yk[2] > hiLTr )
              diff = 4 * el * el;
            else
              diff = el * el;

            if ( N > 1 )
            {
              el = sdata[idxs+3] - yk[3];
              diff += el * el;
              el = sdata[idxs+4] - yk[4];
              diff += el * el;
            }

            if ( diff < 1.0 )
            {
              weight = 1 - weightMap[idxd];
              for ( k = 0; k < lN; k++ )
                Mh[k] += weight * sdata[idxs+k];
              wsuml += weight;
            }
          }
          idxd = slist[idxd];
        }
      }
      if ( wsuml > 0 )
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = Mh[j] / wsuml - yk[j];
      }
      else
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = 0;
      }
      /*****************************************************/

      // Calculate its magnitude squared
      //mvAbs = 0;
      //for(j = 0; j < lN; j++)
      // mvAbs += Mh[j]*Mh[j];
      mvAbs = ( Mh[0] * Mh[0] + Mh[1] * Mh[1] ) * sigmaS * sigmaS;
      if ( N == 3 )
        mvAbs += ( Mh[2] * Mh[2] + Mh[3] * Mh[3] + Mh[4] * Mh[4] ) * sigmaR * sigmaR;
      else
        mvAbs += Mh[2] * Mh[2] * sigmaR * sigmaR;

      // Increment iteration count
      iterationCount++;

    }

    // if a mode was not associated with this data point
    // yet associate it with yk...
    if ( mvAbs >= 0 )
    {
      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // update mode table for this data point
      // indicating that a mode has been associated
      // with it
      modeTable[i] = 1;

    }

    for ( k = 0; k < N; k++ )
      yk[k+2] *= sigmaR;

    // associate the data point indexed by
    // the point list with the mode stored
    // by yk
    for ( j = 0; j < pointCount; j++ )
    {
      // obtain the point location from the
      // point list
      modeCandidate_i = pointList[j];

      // update the mode table for this point
      modeTable[modeCandidate_i] = 1;

      //store result into msRawData...
      for ( k = 0; k < N; k++ )
        msRawData[N*modeCandidate_i+k] = ( float ) ( yk[k+2] );
    }

    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;
  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif
  // de-allocate memory
  delete [] buckets;
  delete [] slist;
  delete [] sdata;

  delete [] yk;
  delete [] Mh;
  // done.
  return;

}

// NEW
void msImageProcessor::NewOptimizedFilter2 ( float sigmaS, float sigmaR )
{
  // Declare Variables
  int  iterationCount, i, j, k, modeCandidateX, modeCandidateY, modeCandidate_i;
  double mvAbs, diff, el;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // let's use some temporary data
  float* sdata;
  sdata = new float[lN*L];

  // copy the scaled data
  int idxs, idxd;
  idxs = idxd = 0;
  if ( N == 3 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else if ( N == 1 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      for ( j = 0; j < N; j++ )
        sdata[idxs++] = data[idxd++] / sigmaR;
    }
  }
  // index the data in the 3d buckets (x, y, L)
  int* buckets;
  int* slist;
  slist = new int[L];
  int bucNeigh[27];

  float sMins; // just for L
  float sMaxs[3]; // for all
  sMaxs[0] = width / sigmaS;
  sMaxs[1] = height / sigmaS;
  sMins = sMaxs[2] = sdata[2];
  idxs = 2;
  float cval;
  for ( i = 0; i < L; i++ )
  {
    cval = sdata[idxs];
    if ( cval < sMins )
      sMins = cval;
    else if ( cval > sMaxs[2] )
      sMaxs[2] = cval;

    idxs += lN;
  }

  int nBuck1, nBuck2, nBuck3;
  int cBuck1, cBuck2, cBuck3, cBuck;
  nBuck1 = ( int ) ( sMaxs[0] + 3 );
  nBuck2 = ( int ) ( sMaxs[1] + 3 );
  nBuck3 = ( int ) ( sMaxs[2] - sMins + 3 );
  buckets = new int[nBuck1*nBuck2*nBuck3];
  for ( i = 0; i < ( nBuck1*nBuck2*nBuck3 ); i++ )
    buckets[i] = -1;

  idxs = 0;
  for ( i = 0; i < L; i++ )
  {
    // find bucket for current data and add it to the list
    cBuck1 = ( int ) sdata[idxs] + 1;
    cBuck2 = ( int ) sdata[idxs+1] + 1;
    cBuck3 = ( int ) ( sdata[idxs+2] - sMins ) + 1;
    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );

    slist[i] = buckets[cBuck];
    buckets[cBuck] = i;

    idxs += lN;
  }
  // init bucNeigh
  idxd = 0;
  for ( cBuck1 = -1; cBuck1 <= 1; cBuck1++ )
  {
    for ( cBuck2 = -1; cBuck2 <= 1; cBuck2++ )
    {
      for ( cBuck3 = -1; cBuck3 <= 1; cBuck3++ )
      {
        bucNeigh[idxd++] = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      }
    }
  }
  double wsuml, weight;
  double hiLTr = 80.0 / sigmaR;
  // done indexing/hashing


  // Initialize mode table used for basin of attraction
  memset ( modeTable, 0, width*height );

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice) ... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif


  for ( i = 0; i < L; i++ )
  {
    // if a mode was already assigned to this data point
    // then skip this point, otherwise proceed to
    // find its mode by applying mean shift...
    if ( modeTable[i] == 1 )
      continue;

    // initialize point list...
    pointCount = 0;

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    idxs = i * lN;
    for ( j = 0; j < lN; j++ )
      yk[j] = sdata[idxs+j];

    // Calculate the mean shift vector using the lattice
    // LatticeMSVector(Mh, yk); // modify to new
    /*****************************************************/
    // Initialize mean shift vector
    for ( j = 0; j < lN; j++ )
      Mh[j] = 0;
    wsuml = 0;
    // uniformLSearch(Mh, yk_ptr); // modify to new
    // find bucket of yk
    cBuck1 = ( int ) yk[0] + 1;
    cBuck2 = ( int ) yk[1] + 1;
    cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
    for ( j = 0; j < 27; j++ )
    {
      idxd = buckets[cBuck+bucNeigh[j]];
      // list parse, crt point is cHeadList
      while ( idxd >= 0 )
      {
        idxs = lN * idxd;
        // determine if inside search window
        el = sdata[idxs+0] - yk[0];
        diff = el * el;
        el = sdata[idxs+1] - yk[1];
        diff += el * el;

        if ( diff < 1.0 )
        {
          el = sdata[idxs+2] - yk[2];
          if ( yk[2] > hiLTr )
            diff = 4 * el * el;
          else
            diff = el * el;

          if ( N > 1 )
          {
            el = sdata[idxs+3] - yk[3];
            diff += el * el;
            el = sdata[idxs+4] - yk[4];
            diff += el * el;
          }

          if ( diff < 1.0 )
          {
            weight = 1 - weightMap[idxd];
            for ( k = 0; k < lN; k++ )
              Mh[k] += weight * sdata[idxs+k];
            wsuml += weight;

            //set basin of attraction mode table
            if ( diff < speedThreshold )
            {
              if ( modeTable[idxd] == 0 )
              {
                pointList[pointCount++] = idxd;
                modeTable[idxd] = 2;
              }
            }
          }
        }
        idxd = slist[idxd];
      }
    }
    if ( wsuml > 0 )
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = Mh[j] / wsuml - yk[j];
    }
    else
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
    }
    /*****************************************************/
    // Calculate its magnitude squared
    //mvAbs = 0;
    //for(j = 0; j < lN; j++)
    // mvAbs += Mh[j]*Mh[j];
    mvAbs = ( Mh[0] * Mh[0] + Mh[1] * Mh[1] ) * sigmaS * sigmaS;
    if ( N == 3 )
      mvAbs += ( Mh[2] * Mh[2] + Mh[3] * Mh[3] + Mh[4] * Mh[4] ) * sigmaR * sigmaR;
    else
      mvAbs += Mh[2] * Mh[2] * sigmaR * sigmaR;


    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // check to see if the current mode location is in the
      // basin of attraction...

      // calculate the location of yk on the lattice
      modeCandidateX = ( int ) ( sigmaS * yk[0] + 0.5 );
      modeCandidateY = ( int ) ( sigmaS * yk[1] + 0.5 );
      modeCandidate_i = modeCandidateY * width + modeCandidateX;

      // if mvAbs != 0 (yk did indeed move) then check
      // location basin_i in the mode table to see if
      // this data point either:

      // (1) has not been associated with a mode yet
      //     (modeTable[basin_i] = 0), so associate
      //     it with this one
      //
      // (2) it has been associated with a mode other
      //     than the one that this data point is converging
      //     to (modeTable[basin_i] = 1), so assign to
      //     this data point the same mode as that of basin_i

      if ( ( modeTable[modeCandidate_i] != 2 ) && ( modeCandidate_i != i ) )
      {
        // obtain the data point at basin_i to
        // see if it is within h*TC_DIST_FACTOR of
        // of yk
        diff = 0;
        idxs = lN * modeCandidate_i;
        for ( k = 2; k < lN; k++ )
        {
          el = sdata[idxs+k] - yk[k];
          diff += el * el;
        }

        // if the data point at basin_i is within
        // a distance of h*TC_DIST_FACTOR of yk
        // then depending on modeTable[basin_i] perform
        // either (1) or (2)
        if ( diff < speedThreshold )
        {
          // if the data point at basin_i has not
          // been associated to a mode then associate
          // it with the mode that this one will converge
          // to
          if ( modeTable[modeCandidate_i] == 0 )
          {
            // no mode associated yet so associate
            // it with this one...
            pointList[pointCount++]  = modeCandidate_i;
            modeTable[modeCandidate_i] = 2;

          } else
          {

            // the mode has already been associated with
            // another mode, thererfore associate this one
            // mode and the modes in the point list with
            // the mode associated with data[basin_i]...

            // store the mode info into yk using msRawData...
            for ( j = 0; j < N; j++ )
              yk[j+2] = msRawData[modeCandidate_i*N+j] / sigmaR;

            // update mode table for this data point
            // indicating that a mode has been associated
            // with it
            modeTable[i] = 1;

            // indicate that a mode has been associated
            // to this data point (data[i])
            mvAbs = -1;

            // stop mean shift calculation...
            break;
          }
        }
      }

      // Calculate the mean shift vector at the new
      // window location using lattice
      // Calculate the mean shift vector using the lattice
      // LatticeMSVector(Mh, yk); // modify to new
      /*****************************************************/
      // Initialize mean shift vector
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
      wsuml = 0;
      // uniformLSearch(Mh, yk_ptr); // modify to new
      // find bucket of yk
      cBuck1 = ( int ) yk[0] + 1;
      cBuck2 = ( int ) yk[1] + 1;
      cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
      cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      for ( j = 0; j < 27; j++ )
      {
        idxd = buckets[cBuck+bucNeigh[j]];
        // list parse, crt point is cHeadList
        while ( idxd >= 0 )
        {
          idxs = lN * idxd;
          // determine if inside search window
          el = sdata[idxs+0] - yk[0];
          diff = el * el;
          el = sdata[idxs+1] - yk[1];
          diff += el * el;

          if ( diff < 1.0 )
          {
            el = sdata[idxs+2] - yk[2];
            if ( yk[2] > hiLTr )
              diff = 4 * el * el;
            else
              diff = el * el;

            if ( N > 1 )
            {
              el = sdata[idxs+3] - yk[3];
              diff += el * el;
              el = sdata[idxs+4] - yk[4];
              diff += el * el;
            }

            if ( diff < 1.0 )
            {
              weight = 1 - weightMap[idxd];
              for ( k = 0; k < lN; k++ )
                Mh[k] += weight * sdata[idxs+k];
              wsuml += weight;

              //set basin of attraction mode table
              if ( diff < speedThreshold )
              {
                if ( modeTable[idxd] == 0 )
                {
                  pointList[pointCount++] = idxd;
                  modeTable[idxd] = 2;
                }
              }

            }
          }
          idxd = slist[idxd];
        }
      }
      if ( wsuml > 0 )
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = Mh[j] / wsuml - yk[j];
      }
      else
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = 0;
      }
      /*****************************************************/

      // Calculate its magnitude squared
      //mvAbs = 0;
      //for(j = 0; j < lN; j++)
      // mvAbs += Mh[j]*Mh[j];
      mvAbs = ( Mh[0] * Mh[0] + Mh[1] * Mh[1] ) * sigmaS * sigmaS;
      if ( N == 3 )
        mvAbs += ( Mh[2] * Mh[2] + Mh[3] * Mh[3] + Mh[4] * Mh[4] ) * sigmaR * sigmaR;
      else
        mvAbs += Mh[2] * Mh[2] * sigmaR * sigmaR;

      // Increment iteration count
      iterationCount++;

    }

    // if a mode was not associated with this data point
    // yet associate it with yk...
    if ( mvAbs >= 0 )
    {
      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // update mode table for this data point
      // indicating that a mode has been associated
      // with it
      modeTable[i] = 1;

    }

    for ( k = 0; k < N; k++ )
      yk[k+2] *= sigmaR;

    // associate the data point indexed by
    // the point list with the mode stored
    // by yk
    for ( j = 0; j < pointCount; j++ )
    {
      // obtain the point location from the
      // point list
      modeCandidate_i = pointList[j];

      // update the mode table for this point
      modeTable[modeCandidate_i] = 1;

      //store result into msRawData...
      for ( k = 0; k < N; k++ )
        msRawData[N*modeCandidate_i+k] = ( float ) ( yk[k+2] );
    }

    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;
  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif
  // de-allocate memory
  delete [] buckets;
  delete [] slist;
  delete [] sdata;

  delete [] yk;
  delete [] Mh;

  // done.
  return;

}

void msImageProcessor::NewNonOptimizedFilter ( float sigmaS, float sigmaR )
{

  // Declare Variables
  int   iterationCount, i, j, k;
  double mvAbs, diff, el;

  //make sure that a lattice height and width have
  //been defined...
  if ( !height )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "LFilter", ( char* ) "Lattice height and width are undefined." );
    return;
  }

  //re-assign bandwidths to sigmaS and sigmaR
  if ( ( ( h[0] = sigmaS ) <= 0 ) || ( ( h[1] = sigmaR ) <= 0 ) )
  {
    ErrorHandler ( ( char* ) "msImageProcessor", ( char* ) "Segment", ( char* ) "sigmaS and/or sigmaR is zero or negative." );
    return;
  }

  //define input data dimension with lattice
  int lN = N + 2;

  // Traverse each data point applying mean shift
  // to each data point

  // Allcocate memory for yk
  double *yk  = new double [lN];

  // Allocate memory for Mh
  double *Mh  = new double [lN];

  // let's use some temporary data
  double* sdata;
  sdata = new double[lN*L];

  // copy the scaled data
  int idxs, idxd;
  idxs = idxd = 0;
  if ( N == 3 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else if ( N == 1 )
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i / width ) / sigmaS;
      sdata[idxs++] = data[idxd++] / sigmaR;
    }
  } else
  {
    for ( i = 0; i < L; i++ )
    {
      sdata[idxs++] = ( i % width ) / sigmaS;
      sdata[idxs++] = ( i % width ) / sigmaS;
      for ( j = 0; j < N; j++ )
        sdata[idxs++] = data[idxd++] / sigmaR;
    }
  }
  // index the data in the 3d buckets (x, y, L)
  int* buckets;
  int* slist;
  slist = new int[L];
  int bucNeigh[27];

  double sMins; // just for L
  double sMaxs[3]; // for all
  sMaxs[0] = width / sigmaS;
  sMaxs[1] = height / sigmaS;
  sMins = sMaxs[2] = sdata[2];
  idxs = 2;
  double cval;
  for ( i = 0; i < L; i++ )
  {
    cval = sdata[idxs];
    if ( cval < sMins )
      sMins = cval;
    else if ( cval > sMaxs[2] )
      sMaxs[2] = cval;

    idxs += lN;
  }

  int nBuck1, nBuck2, nBuck3;
  int cBuck1, cBuck2, cBuck3, cBuck;
  nBuck1 = ( int ) ( sMaxs[0] + 3 );
  nBuck2 = ( int ) ( sMaxs[1] + 3 );
  nBuck3 = ( int ) ( sMaxs[2] - sMins + 3 );
  buckets = new int[nBuck1*nBuck2*nBuck3];
  for ( i = 0; i < ( nBuck1*nBuck2*nBuck3 ); i++ )
    buckets[i] = -1;

  idxs = 0;
  for ( i = 0; i < L; i++ )
  {
    // find bucket for current data and add it to the list
    cBuck1 = ( int ) sdata[idxs] + 1;
    cBuck2 = ( int ) sdata[idxs+1] + 1;
    cBuck3 = ( int ) ( sdata[idxs+2] - sMins ) + 1;
    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );

    slist[i] = buckets[cBuck];
    buckets[cBuck] = i;

    idxs += lN;
  }
  // init bucNeigh
  idxd = 0;
  for ( cBuck1 = -1; cBuck1 <= 1; cBuck1++ )
  {
    for ( cBuck2 = -1; cBuck2 <= 1; cBuck2++ )
    {
      for ( cBuck3 = -1; cBuck3 <= 1; cBuck3++ )
      {
        bucNeigh[idxd++] = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      }
    }
  }
  double wsuml, weight;
  double hiLTr = 80.0 / sigmaR;
  // done indexing/hashing

  // proceed ...
#ifdef PROMPT
  printf ( ( char* ) "done.\nApplying mean shift (Using Lattice)... " );
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\n 0%%" );
#endif
#endif

  for ( i = 0; i < L; i++ )
  {

    // Assign window center (window centers are
    // initialized by createLattice to be the point
    // data[i])
    idxs = i * lN;
    for ( j = 0; j < lN; j++ )
      yk[j] = sdata[idxs+j];

    // Calculate the mean shift vector using the lattice
    // LatticeMSVector(Mh, yk);
    /*****************************************************/
    // Initialize mean shift vector
    for ( j = 0; j < lN; j++ )
      Mh[j] = 0;
    wsuml = 0;
    // uniformLSearch(Mh, yk_ptr); // modify to new
    // find bucket of yk
    cBuck1 = ( int ) yk[0] + 1;
    cBuck2 = ( int ) yk[1] + 1;
    cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
    cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
    for ( j = 0; j < 27; j++ )
    {
      idxd = buckets[cBuck+bucNeigh[j]];
      // list parse, crt point is cHeadList
      while ( idxd >= 0 )
      {
        idxs = lN * idxd;
        // determine if inside search window
        el = sdata[idxs+0] - yk[0];
        diff = el * el;
        el = sdata[idxs+1] - yk[1];
        diff += el * el;

        if ( diff < 1.0 )
        {
          el = sdata[idxs+2] - yk[2];
          if ( yk[2] > hiLTr )
            diff = 4 * el * el;
          else
            diff = el * el;

          if ( N > 1 )
          {
            el = sdata[idxs+3] - yk[3];
            diff += el * el;
            el = sdata[idxs+4] - yk[4];
            diff += el * el;
          }

          if ( diff < 1.0 )
          {
            weight = 1 - weightMap[idxd];
            for ( k = 0; k < lN; k++ )
              Mh[k] += weight * sdata[idxs+k];
            wsuml += weight;
          }
        }
        idxd = slist[idxd];
      }
    }
    if ( wsuml > 0 )
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = Mh[j] / wsuml - yk[j];
    }
    else
    {
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
    }
    /*****************************************************/

    // Calculate its magnitude squared
    mvAbs = 0;
    for ( j = 0; j < lN; j++ )
      mvAbs += Mh[j] * Mh[j];

    // Keep shifting window center until the magnitude squared of the
    // mean shift vector calculated at the window center location is
    // under a specified threshold (Epsilon)

    // NOTE: iteration count is for speed up purposes only - it
    //       does not have any theoretical importance
    iterationCount = 1;
    while ( ( mvAbs >= EPSILON2 ) && ( iterationCount < LIMIT ) )
    {

      // Shift window location
      for ( j = 0; j < lN; j++ )
        yk[j] += Mh[j];

      // Calculate the mean shift vector at the new
      // window location using lattice
      // LatticeMSVector(Mh, yk);
      /*****************************************************/
      // Initialize mean shift vector
      for ( j = 0; j < lN; j++ )
        Mh[j] = 0;
      wsuml = 0;
      // uniformLSearch(Mh, yk_ptr); // modify to new
      // find bucket of yk
      cBuck1 = ( int ) yk[0] + 1;
      cBuck2 = ( int ) yk[1] + 1;
      cBuck3 = ( int ) ( yk[2] - sMins ) + 1;
      cBuck = cBuck1 + nBuck1 * ( cBuck2 + nBuck2 * cBuck3 );
      for ( j = 0; j < 27; j++ )
      {
        idxd = buckets[cBuck+bucNeigh[j]];
        // list parse, crt point is cHeadList
        while ( idxd >= 0 )
        {
          idxs = lN * idxd;
          // determine if inside search window
          el = sdata[idxs+0] - yk[0];
          diff = el * el;
          el = sdata[idxs+1] - yk[1];
          diff += el * el;

          if ( diff < 1.0 )
          {
            el = sdata[idxs+2] - yk[2];
            if ( yk[2] > hiLTr )
              diff = 4 * el * el;
            else
              diff = el * el;

            if ( N > 1 )
            {
              el = sdata[idxs+3] - yk[3];
              diff += el * el;
              el = sdata[idxs+4] - yk[4];
              diff += el * el;
            }

            if ( diff < 1.0 )
            {
              weight = 1 - weightMap[idxd];
              for ( k = 0; k < lN; k++ )
                Mh[k] += weight * sdata[idxs+k];
              wsuml += weight;
            }
          }
          idxd = slist[idxd];
        }
      }
      if ( wsuml > 0 )
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = Mh[j] / wsuml - yk[j];
      }
      else
      {
        for ( j = 0; j < lN; j++ )
          Mh[j] = 0;
      }
      /*****************************************************/

      // Calculate its magnitude squared
      //mvAbs = 0;
      //for(j = 0; j < lN; j++)
      // mvAbs += Mh[j]*Mh[j];
      mvAbs = ( Mh[0] * Mh[0] + Mh[1] * Mh[1] ) * sigmaS * sigmaS;
      if ( N == 3 )
        mvAbs += ( Mh[2] * Mh[2] + Mh[3] * Mh[3] + Mh[4] * Mh[4] ) * sigmaR * sigmaR;
      else
        mvAbs += Mh[2] * Mh[2] * sigmaR * sigmaR;

      // Increment interation count
      iterationCount++;
    }

    // Shift window location
    for ( j = 0; j < lN; j++ )
      yk[j] += Mh[j];

    //store result into msRawData...
    for ( j = 0; j < N; j++ )
      msRawData[N*i+j] = ( float ) ( yk[j+2] * sigmaR );

    // Prompt user on progress
#ifdef SHOW_PROGRESS
    percent_complete = ( float ) ( i / ( float ) ( L ) ) * 100;
    printf ( ( char* ) "\r%2d%%", ( int ) ( percent_complete + 0.5 ) );
#endif

    // Check to see if the algorithm has been halted
    if ( ( i % PROGRESS_RATE == 0 ) && ( ( ErrorStatus = msSys.Progress ( ( float ) ( i / ( float ) ( L ) ) * ( float ) ( 0.8 ) ) ) ) == EL_HALT )
      break;
  }

  // Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
  printf ( ( char* ) "\r" );
#endif
  printf ( ( char* ) "done." );
#endif

  // de-allocate memory
  delete [] buckets;
  delete [] slist;
  delete [] sdata;

  delete [] yk;
  delete [] Mh;

  // done.
  return;

}

void msImageProcessor::SetSpeedThreshold ( float speedUpThreshold )
{
  speedThreshold = speedUpThreshold;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ END OF CLASS DEFINITION @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/
/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/