combinatorial-leaf-fitting.h

Go to the documentation of this file.

/*
Copyright (C) 2012 Babette Dellen

Contour fitting algorithm 

*/

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

image<rgb> *combinatorial_leaf_fitting(int ** clusters, image<rgb> *model,image<rgb> *model_area, int num_clusters, double *& fit_score_list) 
{
  int width = model->width();
  int height = model->height();
  int ** segment_contoursX; 
  int ** segment_contoursY; 
  int ** segment_contoursLabel;
  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_contoursLabel = (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);
  segment_contoursLabel[u_label] = (int*)malloc(sizeof(int)*1);
  }
  
  
  int contour_size;
  
  int encountered_neighbor_labels [8]; 
  int number_encountered_labels [8]; 
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
      int n_count = 0
       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;
              encountered_neighbor_labels[n_count] = clusters[u_y][u_x];
              n_count = n_count +1;
              }
              
              }
              
              
              
              }
       }
       
       
       encountered_neighbor_labels = {0,0,0,0,0,0,0,0};
       number_encountered_labels = {0,0,0,0,0,0,0,0};
       
      
      
        if (edge_flag==1)
        {
          for (int j_n1=1;j_n1<n_count;j_n1++)
          {
            encountered_label_1 = encountered_neighbor_labels[j_n1];
              if (encountered_label_1>0)
                {
            
                  
                  for (int j_n2=1;j_n2<n_count;j_n2++)
                    {
                      encountered_label_2 = encountered_neighbor_labels[j_n2];
                      if (encountered_label_2==encountered_label_1)
                      {
                      number_encountered_labels[j_n1] = number_encountered_labels[j_n1] +1;
                      }
                      
                    }
            
                }
          }
          
           int max_encountered_label=0;
           int label_of_neighbor;
           for (int j_n3=1;j_n3<n_count;j_n3++)
                    {
                      if (number_encountered_labels[j_n3]>max_encountered_label)
                      {max_encountered_label = number_encountered_labels[j_n3];
                      label_of_neighbor = encountered_neighbor_labels[j_n3];
                      }
                      
                      
                    }
          
          
        //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_contoursLabel[current_label-1] = (int*) realloc(segment_contoursLabel[current_label-1],(contour_size+1)*sizeof(int));
        segment_contoursLabel[current_label-1][contour_size] = label_of_neighbor;
        segment_contour_size[current_label-1] = contour_size+1;
        
        
                
        
        
        
        
       
        }
        
      
      }
    
    
    
    }
  
  
   }
 
 
int is_boundary;

//int is_neighbor;

int ** NeighborGraph;
int j_count = num_clusters;
int delta = 1;
double depth;
double depth_neighbor;
double depth_dist;
double depth_thres = 5;
int label_neighbor;
int k_1;
int k_2;
  NeighborGraph=(int**)malloc(sizeof(int*)*j_count);
    
        for (int k=0; k<j_count; k++)
        {
                NeighborGraph[k] = (int*)malloc(sizeof(int)*j_count);
                
                
        }


for (k_1=0; k_1<j_count; k_1++)
        {
          for (k_2=0; k_2<j_count; k_2++)
        
          
          {
        NeighborGraph[k_1][k_2]=0;
        }
        
        }
          
          

int y;
int x;
int j_y;
int j_x;
for (y = 0; y < height; y++) {
    for (x = 0; x < width; x++) {
             
      is_boundary = 0;
             label = clusters[y][x];
            //if (comp_size_list[label]>size_thres){
             depth = (double)(imRef(fitted_depth, x, y).b);
             
         //   is_neighbor = 0;
           for (j_y = y-1; j_y<y+2; j_y++) {
              for (j_x = x-1; j_x<x+2; j_x++) {
            
              if ((j_y>=0)&&(j_x>=0)&&(j_x<width)&&(j_y<height))
              {
              label_neighbor = clusters[j_y][j_x];
              //printf ("%8d\n",j_x);
              if (label!=label_neighbor)
              { 
              is_boundary = 1;
              }
              
              }
                
             
              }
           }
      
            //}
      if (is_boundary==1)
      {
         for (j_y = y-delta; j_y< y+1+delta; j_y++) {
              for (j_x = x-delta; j_x < x+1+delta; j_x++) {
            
              if ((j_y>=0)&&(j_x>=0)&&(j_x<width)&&(j_y<height))
              {
                // check if neighbor
                
                label_neighbor = clusters[j_y][j_x];
                //if (comp_size_list[label_neighbor]>size_thres){
             //printf ("%8d\n",j_x);
                //if (label_neighbor>label){
                depth_neighbor = (double)(imRef(fitted_depth, j_x, j_y).b);
                depth_dist = abs(depth - depth_neighbor);
                if ((label!=label_neighbor)&&(depth_dist<depth_thres))
                {NeighborGraph[label][label_neighbor]=1;
                NeighborGraph[label_neighbor][label]=1;
                
                }
                
                //}
                
              //}
              }
              
              }
         }
        
        
      
      }
      
 
 
      
  }
}




vector<int> graph_label_1;
vector<int> graph_label_2;


for (k_1=0; k_1<j_count; k_1++) 
{
          for (k_2=0; k_2<j_count; k_2++)
          {
            if (k_1<k_2){
                
              if (NeighborGraph[k_1][k_2]==1)
                {
                 
                  graph_label_1.push_back (k_1);
                  graph_label_2.push_back (k_2);
                
                }
            }
        
        }
        
}





// sort graph links in ascending order      
      
//  vector<int> graph_index;
//  vector<int> graph_label_1_sort;
//  vector<int> graph_label_2_sort;
//  
//  BubbleSort(&graph_links,&graph_index);
//  ReArrange(&graph_label_1_sort,&graph_label_1,&graph_index);
//  ReArrange(&graph_label_2_sort,&graph_label_2,&graph_index);
      

//printf ("Links Sorted!");     

// ////// Do graph-based merging


int * labels_new;
labels_new = (int*)malloc(sizeof(int)*j_count);

int * comp_size_list_new;
comp_size_list_new = (int*)malloc(sizeof(int)*j_count);

double * comp_depth_list_new;
comp_depth_list_new = (double*)malloc(sizeof(double)*j_count);


for (int label_j=0;label_j<j_count-1;label_j++)
{labels_new[label_j]=label_j;
comp_size_list_new[label_j]=comp_size_list[label_j];//label_j;
comp_depth_list_new[label_j]=comp_depth_list[label_j];
}


int ** label_list;
label_list = (int**)malloc(sizeof(int*)*j_count);


int * number_of_labels_in_list;
number_of_labels_in_list = (int*)malloc(sizeof(int)*j_count);
for (int k_b=0; k_b<j_count; k_b++)
        {
                label_list[k_b] = (int*)malloc(sizeof(int));
                
        }

for (int label_j=0;label_j<j_count-1;label_j++)
{label_list[label_j][0]=label_j;
number_of_labels_in_list[label_j] = 1;
}

// int aa** = (int**) new int*[10];
// for (i=0..9)
//   aa[i]=(int*) new int[100];
// 
// for (i=0..9)
//   delete [] aa[i];
// delete [] aa;


int u_1;
int u_2;
int label_1;
int label_2;
int error_1;
int error_2;
int size_1;
int size_2;
int label_1_total_size;
int label_2_total_size;
double thres_merge = 15;
int number_merged_labels;
double min_error;
int number_of_links = graph_links.size();
int label_keep;
int label_replace;


for (int link_index=0;link_index<number_of_links;link_index++)
{
u_1 = graph_label_1.at(link_index);
u_2 = graph_label_2.at(link_index);

label_1=labels_new[u_1];
label_2=labels_new[u_2];



// merged label list
number_merged_labels = number_of_labels_in_list[label_1]+number_of_labels_in_list[label_2];
int * merged_labels = (int*)malloc(sizeof(int)*number_merged_labels);
 
for (int j_merge_1=0;j_merge_1<number_of_labels_in_list[label_1];j_merge_1++)
{
merged_labels[j_merge_1]= label_list[label_1][j_merge_1]; 
}

for (int j_merge_2=0;j_merge_2<number_of_labels_in_list[label_2];j_merge_2++)
{
merged_labels[j_merge_2+number_of_labels_in_list[label_1]]= label_list[label_2][j_merge_2]; 
}

number_of_labels_in_list[label_1]=number_merged_labels;
number_of_labels_in_list[label_2]=number_merged_labels;

if (label_1<label_2)
{
  label_replace = label_2;
  label_keep = label_1;
}
else
  {
  label_replace = label_1;
  label_keep = label_2;
}
  
// update label_new
for (int w_2=0;w_2<j_count-1;w_2++)
{
if (labels_new[w_2]==label_replace)
{
 labels_new[w_2]=label_keep;
}
} 

//update label_list
label_list[label_keep] = (int*) realloc (label_list[label_keep], number_merged_labels * sizeof(int));
for (int j_merge=0;j_merge<number_merged_labels;j_merge++)
{
label_list[label_keep][j_merge]=merged_labels[j_merge];
}



// free number_merged_labels
free (merged_labels);

}



  





// select best merges using an ordered list




// from accepted merges construct new segmentation 








  

  
// 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)){
//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;
    //final_min = min_fit_value [1];
   
    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 = 0;//(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)(0);
     imRef(output, j_x, j_y).g = (uchar)(0);
     imRef(output, j_x, j_y).b = (uchar)(0);
     
     }
       
  }
  }
  
   
//   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