leaf-fitting.h

Go to the documentation of this file.

/*
Copyright (C) 2012 Babette Dellen

Contour fitting algorithm 

*/

#ifndef LEAF_FITTING
#define LEAF_FITTING
#include <iostream>
//#include "cv.h"
#include <cstdlib>
#include <image.h>
#include <misc.h>
#include <filter.h>
#include <vector>
#define PI 3.14159265
using std::vector;


double ContourFittingError(int * contour_pointX, int * contour_pointY,int number_contour_points,image<rgb> *model, double * current_trans)
{
  double dist;
  //dist=0;
  int j;
  int width = model->width();
  int height = model->height();
  int x;
  int y;
  double meanX=0;
  double meanY=0;
  double match_point_sum = number_contour_points;
  double rotX;
  double rotY;
  image<float> *dummy = new image<float>(width, height);
  // smooth each color channel  
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      imRef(dummy, x, y) = 0;
      
    }
  }
  
  
  //transform contour according to parameters
  for (j=0;j<number_contour_points;j++)
  {//shift contour 
  meanX = meanX + contour_pointX[j];
  meanY = meanY + contour_pointY[j];
  
  }
  
  meanX = meanX/((double)(number_contour_points));
  meanY = meanY/((double)(number_contour_points));
  
  
  for (j=0;j<number_contour_points;j++)
  {
    // rotating + scaling + shifting
    rotX = current_trans[0]*((contour_pointX[j]-meanX)*cos(current_trans[3]) - (contour_pointY[j]-meanY)*sin(current_trans[3])) + meanX + current_trans[1];
    rotY = current_trans[0]*((contour_pointX[j]-meanX)*sin(current_trans[3]) + (contour_pointY[j]-meanY)*cos(current_trans[3])) + meanY + current_trans[2];
  
    x = round(rotX);
    y = round(rotY);
    //imRef(dummy, x, y)
    
    if ((x>=0)&&(x<width)&&(y>=0)&&(y<height))
    {
//       printf("error %d value\n", x);
//       printf("error %d value\n", y);
//       printf("error %8f value\n",imRef(model,x,y).b);
      if (imRef(dummy,x,y)==0)
      {
      imRef(dummy,x,y)=1;
      match_point_sum = match_point_sum - (double)(imRef(model,x,y).b)/((double)(170));
      //printf("match %f value\n", (double)(imRef(model,x,y).b)/((double)(170)));
      }
      
    }
    
    
  }
  
  dist = match_point_sum/((double)(number_contour_points));
  delete dummy;
    
return dist;
}



//****************************************************************************80

void nelmin_contour (int n, double * start, double * xmin, 
  double * ynewlo, double reqmin, double * step, int konvge, int kcount, 
  int *icount, int *numres, int *ifault,int * contour_pointX, int * contour_pointY,int number_contour_points,image<rgb> *model)
  
  
//****************************************************************************80
//
//  Purpose:
//
//    NELMIN minimizes a function using the Nelder-Mead algorithm.
//
//  Discussion:
//
//    This routine seeks the minimum value of a user-specified function.
//
//    Simplex function minimisation procedure due to Nelder+Mead(1965),
//    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
//    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
//    25, 97) and Hill(1978, 27, 380-2)
//
//    The function to be minimized must be defined by a function of
//    the form
//
//      function fn ( x, f )
//      double fn
//      double x(*)
//
//    and the name of this subroutine must be declared EXTERNAL in the
//    calling routine and passed as the argument FN.
//
//    This routine does not include a termination test using the
//    fitting of a quadratic surface.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    27 February 2008
//
//  Author:
//
//    Original FORTRAN77 version by R ONeill.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    John Nelder, Roger Mead,
//    A simplex method for function minimization,
//    Computer Journal,
//    Volume 7, 1965, pages 308-313.
//
//    R ONeill,
//    Algorithm AS 47:
//    Function Minimization Using a Simplex Procedure,
//    Applied Statistics,
//    Volume 20, Number 3, 1971, pages 338-345.
//
//  Parameters:
//
//    Input, double FN ( double x[] ), the name of the routine which evaluates
//    the function to be minimized.
//
//    Input, int N, the number of variables.
//
//    Input/output, double START[N].  On input, a starting point
//    for the iteration.  On output, this data may have been overwritten.
//
//    Output, double XMIN[N], the coordinates of the point which
//    is estimated to minimize the function.
//
//    Output, double YNEWLO, the minimum value of the function.
//
//    Input, double REQMIN, the terminating limit for the variance
//    of function values.
//
//    Input, double STEP[N], determines the size and shape of the
//    initial simplex.  The relative magnitudes of its elements should reflect
//    the units of the variables.
//
//    Input, int KONVGE, the convergence check is carried out 
//    every KONVGE iterations.
//
//    Input, int KCOUNT, the maximum number of function 
//    evaluations.
//
//    Output, int *ICOUNT, the number of function evaluations 
//    used.
//
//    Output, int *NUMRES, the number of restarts.
//
//    Output, int *IFAULT, error indicator.
//    0, no errors detected.
//    1, REQMIN, N, or KONVGE has an illegal value.
//    2, iteration terminated because KCOUNT was exceeded without convergence.
//
{
  double ccoeff = 0.5;
  double del;
  double dn;
  double dnn;
  double ecoeff = 2.0;
  double eps = 0.001;
  int i;
  int ihi;
  int ilo;
  int j;
  int jcount;
  int l;
  int nn;
  double *p;
  double *p2star;
  double *pbar;
  double *pstar;
  double rcoeff = 1.0;
  double rq;
  double x;
  double *y;
  double y2star;
  double ylo;
  double ystar;
  double z;
  double * start_init;
//
//  Check the input parameters.
//
  if ( reqmin <= 0.0 )
  {
    *ifault = 1;
    return;
  }

  if ( n < 1 )
  {
    *ifault = 1;
    return;
  }

  if ( konvge < 1 )
  {
    *ifault = 1;
    return;
  }

  p = new double[n*(n+1)];
  pstar = new double[n];
  p2star = new double[n];
  pbar = new double[n];
  y = new double[n+1];
  //start_init = new double[n];
  
//   for ( i = 0; i < n; i++ )
//     { 
//       start_init[i] = start[i];
//     }
//   

  *icount = 0;
  *numres = 0;

  jcount = konvge; 
  dn = ( double ) ( n );
  nn = n + 1;
  dnn = ( double ) ( nn );
  del = 1.0;
  rq = reqmin * dn;
//
//  Initial or restarted loop.
//
  for ( ; ; )
  {
    for ( i = 0; i < n; i++ )
    { 
      p[i+n*n] = start[i];
    }
    y[n] = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,start);
    //FittingError(x_val,y_val,z_val,start,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( start );
    *icount = *icount + 1;

    for ( j = 0; j < n; j++ )
    {
      x = start[j];
      start[j] = start[j] + step[j] * del;
      for ( i = 0; i < n; i++ )
      {
//      if (i==j)
//      {start[i]=start_init[j] + step[j] * del;
//      p[i+j*n] =start_init[j] + step[j] * del;
//      }
//      else
//      {start[i]=start_init[i];
//      p[i+j*n]=start_init[i];
//      }
        
          
        p[i+j*n] = start[i];
      }
      y[j] = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,start);
      //FittingError(x_val,y_val,z_val,start,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( start );
      *icount = *icount + 1;
      start[j] = x;
    }
//                    
//  The simplex construction is complete.
//                    
//  Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
//  the vertex of the simplex to be replaced.
//                
    ylo = y[0];
    ilo = 0;

    for ( i = 1; i < nn; i++ )
    {
      if ( y[i] < ylo )
      {
        ylo = y[i];
        ilo = i;
      }
    }
//
//  Inner loop.
//
    for ( ; ; )
    {
      if ( kcount <= *icount )
      {
        break;
      }
      *ynewlo = y[0];
      ihi = 0;

      for ( i = 1; i < nn; i++ )
      {
        if ( *ynewlo < y[i] )
        {
          *ynewlo = y[i];
          ihi = i;
        }
      }
//
//  Calculate PBAR, the centroid of the simplex vertices
//  excepting the vertex with Y value YNEWLO.
//
      for ( i = 0; i < n; i++ )
      {
        z = 0.0;
        for ( j = 0; j < nn; j++ )
        { 
          z = z + p[i+j*n];
        }
        z = z - p[i+ihi*n];  
        pbar[i] = z / dn;
      }
//
//  Reflection through the centroid.
//
      for ( i = 0; i < n; i++ )
      {
        pstar[i] = pbar[i] + rcoeff * ( pbar[i] - p[i+ihi*n] );
      }
      ystar = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,pstar);
      //FittingError(x_val,y_val,z_val,p2star,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( pstar );
      *icount = *icount + 1;
//
//  Successful reflection, so extension.
//
      if ( ystar < ylo )
      {
        for ( i = 0; i < n; i++ )
        {
          p2star[i] = pbar[i] + ecoeff * ( pstar[i] - pbar[i] );
        }
        y2star = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,p2star);
        //FittingError(x_val,y_val,z_val,p2star,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( p2star );
        *icount = *icount + 1;
//
//  Check extension.
//
        if ( ystar < y2star )
        {
          for ( i = 0; i < n; i++ )
          {
            p[i+ihi*n] = pstar[i];
          }
          y[ihi] = ystar;
        }
//
//  Retain extension or contraction.
//
        else
        {
          for ( i = 0; i < n; i++ )
          {
            p[i+ihi*n] = p2star[i];
          }
          y[ihi] = y2star;
        }
      }
//
//  No extension.
//
      else
      {
        l = 0;
        for ( i = 0; i < nn; i++ )
        {
          if ( ystar < y[i] )
          {
            l = l + 1;
          }
        }

        if ( 1 < l )
        {
          for ( i = 0; i < n; i++ )
          {
            p[i+ihi*n] = pstar[i];
          }
          y[ihi] = ystar;
        }
//
//  Contraction on the Y(IHI) side of the centroid.
//
        else if ( l == 0 )
        {
          for ( i = 0; i < n; i++ )
          {
            p2star[i] = pbar[i] + ccoeff * ( p[i+ihi*n] - pbar[i] );
          }
          y2star = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,p2star);
          //FittingError(x_val,y_val,z_val,p2star,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( p2star );
          *icount = *icount + 1;
//
//  Contract the whole simplex.
//
          if ( y[ihi] < y2star )
          {
            for ( j = 0; j < nn; j++ )
            {
              for ( i = 0; i < n; i++ )
              {
                p[i+j*n] = ( p[i+j*n] + p[i+ilo*n] ) * 0.5;
                xmin[i] = p[i+j*n];
              }
              y[j] = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,xmin);
              //FittingError(x_val,y_val,z_val,xmin,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( xmin );
              *icount = *icount + 1;
            }
            ylo = y[0];
            ilo = 0;

            for ( i = 1; i < nn; i++ )
            {
              if ( y[i] < ylo )
              {
                ylo = y[i];
                ilo = i;
              }
            }
            continue;
          }
//
//  Retain contraction.
//
          else
          {
            for ( i = 0; i < n; i++ )
            {
              p[i+ihi*n] = p2star[i];
            }
            y[ihi] = y2star;
          }
        }
//
//  Contraction on the reflection side of the centroid.
//
        else if ( l == 1 )
        {
          for ( i = 0; i < n; i++ )
          {
            p2star[i] = pbar[i] + ccoeff * ( pstar[i] - pbar[i] );
          }
          y2star = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,p2star);
          //FittingError(x_val,y_val,z_val,p2star,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( p2star );
          *icount = *icount + 1;
//
//  Retain reflection?
//
          if ( y2star <= ystar )
          {
            for ( i = 0; i < n; i++ )
            {
              p[i+ihi*n] = p2star[i];
            }
            y[ihi] = y2star;
          }
          else
          {
            for ( i = 0; i < n; i++ )
            {
              p[i+ihi*n] = pstar[i];
            }
            y[ihi] = ystar;
          }
        }
      }
//
//  Check if YLO improved.
//
      if ( y[ihi] < ylo )
      {
        ylo = y[ihi];
        ilo = ihi;
      }
      jcount = jcount - 1;

      if ( 0 < jcount )
      {
        continue;
      }
//
//  Check to see if minimum reached.
//
      if ( *icount <= kcount )
      {
        jcount = konvge;

        z = 0.0;
        for ( i = 0; i < nn; i++ )
        {
          z = z + y[i];
        }
        x = z / dnn;

        z = 0.0;
        for ( i = 0; i < nn; i++ )
        {
          z = z + pow ( y[i] - x, 2 );
        }

        if ( z <= rq )
        {
          break;
        }
      }
    }
//
//  Factorial tests to check that YNEWLO is a local minimum.
//
    for ( i = 0; i < n; i++ )
    {
      xmin[i] = p[i+ilo*n];
    }
    *ynewlo = y[ilo];

    if ( kcount < *icount )
    {
      *ifault = 2;
      break;
    }

    *ifault = 0;

    for ( i = 0; i < n; i++ )
    {
      del = step[i] * eps;
      xmin[i] = xmin[i] + del;
      z = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,xmin);
      //FittingError(x_val,y_val,z_val,xmin,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( xmin );
      *icount = *icount + 1;
      if ( z < *ynewlo )
      {
        *ifault = 2;
        break;
      }
      xmin[i] = xmin[i] - del - del;
      z = ContourFittingError(contour_pointX, contour_pointY,number_contour_points,model,xmin);
      //FittingError(x_val,y_val,z_val,xmin,number_of_elements_in_segment,current_label_list, number_current_labels);//num_p,current_label);//fn ( xmin );
      *icount = *icount + 1;
      if ( z < *ynewlo )
      {
        *ifault = 2;
        break;
      }
      xmin[i] = xmin[i] + del;
    }

    if ( *ifault == 0 )
    {
      break;
    }
//
//  Restart the procedure.
//
    for ( i = 0; i < n; i++ )
    {
      start[i] = xmin[i];
    }
    del = eps;
    *numres = *numres + 1;
  }
  delete [] p;
  delete [] pstar;
  delete [] p2star;
  delete [] pbar;
  delete [] y;

  return;
}



/*
 * Segment an image
 *
 * Returns a color image representing the segmentation.
 *
 * im: image to segment.
 * sigma: to smooth the image.
 * c: constant for treshold function.
 * min_size: minimum component size (enforced by post-processing stage).
 * num_ccs: number of connected components in the segmentation.
 */
// image<float> *segment_image(image<rgb> *im, image<rgb> *pc, float sigma, float c, int min_size,
//                        int *num_ccs) {
image<rgb> *leaf_fitting(int ** clusters, image<rgb> *model,image<rgb> *model_area,int num_clusters, double *& fit_score_list, double **& transformation_params) 
{
  int width = model->width();
  int height = model->height();
  int ** segment_contoursX; 
  int ** segment_contoursY; 
  //double ** transformation_params;
  //double ** found_transformation_params;
  int current_label;
  int * segment_contour_size;
  int edge_flag;
  int j_x;
  int j_y;
  int u_x;
  int u_y;
  int length_model_bound = 0;
  
  double center_modelX=0;
  double center_modelY=0;
  double model_size = 0;
  for (j_x=0;j_x<width;j_x++)
  {
  for (j_y=0;j_y<height;j_y++)
  {
    edge_flag = 0;
    //printf("%8d\n",current_label);
    if (imRef(model_area,j_x,j_y).b>10)
    {
       center_modelX = center_modelX + j_x;
       center_modelY = center_modelY + j_y;
       model_size = model_size +1;
       
      //test if edge point
      
       for (u_y = j_y-1; u_y< j_y+1+1; u_y++) 
       {
              for (u_x = j_x-1; u_x < j_x+1+1; u_x++) 
              {
             
              if ((u_y>=0)&&(u_x>=0)&&(u_x<width)&&(u_y<height))
              {
              if (imRef(model_area,u_x,u_y).b<10)
              {
              edge_flag = 1;
              }
              
              }
              
              
              
              }
       }
      
        
      }
    
        if (edge_flag==1)
        {
          length_model_bound = length_model_bound + 1;
       
        }
        
    
    
    
    }
  
  
   }
  
  center_modelX = center_modelX/model_size;
  center_modelY = center_modelY/model_size;
  
  segment_contoursX = (int**)malloc(sizeof(int*)*(num_clusters));
  segment_contoursY = (int**)malloc(sizeof(int*)*(num_clusters));
  segment_contour_size = (int*)malloc(sizeof(int)*num_clusters);
  double * segment_centerX = (double*)malloc(sizeof(double)*num_clusters);
  double * segment_centerY = (double*)malloc(sizeof(double)*num_clusters);
  double * segment_size = (double*)malloc(sizeof(double)*num_clusters);
  transformation_params = (double**)malloc(sizeof(double*)*(num_clusters));
  
  int u_label;
  for (u_label=0;u_label<num_clusters;u_label++)
  {
  segment_contour_size[u_label]=0;
  segment_size[u_label]=0;
  transformation_params[u_label] = (double*)malloc(sizeof(double)*4);
  segment_contoursX[u_label] = (int*)malloc(sizeof(int)*1);
  segment_contoursY[u_label] = (int*)malloc(sizeof(int)*1);
  }
  
  
  int contour_size;
for (j_x=0;j_x<width;j_x++)
  {
  for (j_y=0;j_y<height;j_y++)
  {
    current_label = clusters[j_y][j_x];
    //printf("%8d\n",current_label);
    if (current_label>0)
    {
      segment_centerX[current_label-1] =  segment_centerX[current_label-1] + j_x;
      segment_centerY[current_label-1] =  segment_centerY[current_label-1] + j_y;
      segment_size[current_label-1] =  segment_size[current_label-1] + 1;
       edge_flag = 0;
      //test if edge point
      
       for (u_y = j_y-1; u_y< j_y+1+1; u_y++) 
       {
              for (u_x = j_x-1; u_x < j_x+1+1; u_x++) 
              {
             
              if ((u_y>=0)&&(u_x>=0)&&(u_x<width)&&(u_y<height))
              {
              if (clusters[u_y][u_x]!=current_label)
              {
              edge_flag = 1;
              }
              
              }
              
              
              
              }
       }
      
      
        if (edge_flag==1)
        {//add to boundary list
        contour_size = segment_contour_size[current_label-1];
        //printf("%8d\n",current_label);
        segment_contoursX[current_label-1] = (int*) realloc(segment_contoursX[current_label-1],(contour_size+1)*sizeof(int));
        segment_contoursX[current_label-1][contour_size] = j_x;
        segment_contoursY[current_label-1] = (int*) realloc(segment_contoursY[current_label-1],(contour_size+1)*sizeof(int));
        segment_contoursY[current_label-1][contour_size] = j_y;
        segment_contour_size[current_label-1] = contour_size+1;
       
        }
        
      
      }
    
    
    
    }
  
  
   }
  
  
// fit extracted contours to model contour

 
double initial_parameters_trans [4] = {1,0,0,0};//mean_depth};scale,trans_x,trans_y,angle
double parameters_at_minimum [4] = {100,100,100,100};
double min_fit_value [1] = {100}; 
double required_min = 0.1;
double initial_simplex_step [4] = {1,10,10,0.2};
int konvergence_check = 100;
int max_number_evaluations = 4000;
int used_number_evaluations [1] = {0};
int number_restarts [1] = {0};
int error_detected [1] = {0};
int parameter_trans = 4;  
double final_params [4]; 
double final_min;
double dist;
double x_shift;
double y_shift;

fit_score_list = (double*)malloc(sizeof(double)*num_clusters);
 for (u_label=0;u_label<num_clusters-1;u_label++)
  {
    //if (((model_size/segment_size[u_label])<10)&&((model_size/segment_size[u_label])>0.1)){
    if (((model_size/segment_size[u_label])<2)&&((model_size/segment_size[u_label])>0.5)){
//dist = ContourFittingError(segment_contoursX[u_label], segment_contoursY[u_label],segment_contour_size[u_label],model, initial_parameters_trans);
      final_min = 1000;
      final_params[0]=0;
      final_params[1]=0;
      final_params[2]=0;
      final_params[3]=0;
      for (int u_rot=0;u_rot<18;u_rot++)
      {
        //int u_rot = 0;
        x_shift = center_modelX - segment_centerX[u_label]/segment_size[u_label];
        y_shift = center_modelY - segment_centerY[u_label]/segment_size[u_label];
        initial_parameters_trans[0]=model_size/segment_size[u_label];
    
      //printf("scale_param %8f value\n",model_size);
      //printf("scale_param %8f value\n",initial_parameters_trans[0]);

        initial_parameters_trans[1]= x_shift;
        initial_parameters_trans[2]= y_shift;
        initial_parameters_trans[3]=PI*(u_rot*20)/((double)(180));
    
      //printf("int angle %f value\n",initial_parameters_trans[3]);
      //dist = ContourFittingError(segment_contoursX[u_label], segment_contoursY[u_label],segment_contour_size[u_label],model, initial_parameters_trans);
      //printf("error %8f value\n", dist);
    
        nelmin_contour(parameter_trans,initial_parameters_trans, parameters_at_minimum, min_fit_value, required_min, initial_simplex_step, konvergence_check, max_number_evaluations,used_number_evaluations,number_restarts, error_detected,segment_contoursX[u_label],segment_contoursY[u_label],segment_contour_size[u_label],model);
      //printf("error %8f value\n", min_fit_value[0]);
              if (min_fit_value[0]<final_min)
              {
                final_min = min_fit_value[0];
                final_params[0]=parameters_at_minimum[0];
                final_params[1]=parameters_at_minimum[1];
                final_params[2]=parameters_at_minimum[2];
                final_params[3]=parameters_at_minimum[3];
              }
      }
  
    //printf("final error %f value\n", final_min);
    //printf("size %f value\n", segment_size[u_label]);

      fit_score_list[u_label] = final_min;
      transformation_params[u_label][0]=final_params[0];
      transformation_params[u_label][1]=final_params[1];
      transformation_params[u_label][2]=final_params[2];
      transformation_params[u_label][3]=final_params[3];
    //printf("scale %f value\n",final_params[0]);
    //printf("angle %f value\n",final_params[3]);
    //printf("shift x %f value\n",final_params[1]);
    //printf("shift y %f value\n",final_params[2]);
    }
    else
    {
      fit_score_list[u_label] = 1;}
    }
    
    image<rgb> *output = new image<rgb>(width, height); 
  
    for (j_x=0;j_x<width;j_x++)
    {
      for (j_y=0;j_y<height;j_y++)
      {
        current_label = clusters[j_y][j_x];
        if (current_label>0)
        {
          imRef(output, j_x, j_y).r = (uchar)((1-fit_score_list[current_label-1])*250);
          imRef(output, j_x, j_y).g = (uchar)(current_label);//(uchar)((1-fit_score_list[current_label-1]))*250;
          imRef(output, j_x, j_y).b = 0;//(uchar)((1-fit_score_list[current_label-1]))*250;
        }
        else 
        {
          imRef(output, j_x, j_y).r = (uchar)(255);
          imRef(output, j_x, j_y).g = (uchar)(255);
          imRef(output, j_x, j_y).b = (uchar)(255);
        }
      }
    }
  
   
//   image<rgb> *output = new image<rgb>(width, height);
// 
//   // pick random colors for each component
//   rgb *colors = new rgb[width*height];
//   
//   c.r = (uchar)random();
//   c.g = (uchar)random();
//   c.b = (uchar)random();
// 
//   return c;
// }
//   for (int y = 0; y < height; y++) {
//     for (int x = 0; x < width; x++) {
//       int comp = ClusterArray[y][x]-1;
// //       int comp = u->find(y * width + x);
// //       printf ("%8d\n",comp);
//       int comp_fin = labels_new[comp];
// //       printf ("%8d\n",comp_fin);
//       imRef(output, x, y) = clusters[y][x];//colors[comp_fin];
//     }
//   }  
/*
  delete [] colors;  
  delete u;*/

  //return fitted_depth;
  return(output);
}

#endif