mean_shift.cpp

mean-shift. pointer sample

	if (up)    gray = max(gray,get_persistence_of_pair(l,l_up));
	if (down)  gray = max(gray,get_persistence_of_pair(l,l_down));
	if (left)  gray = max(gray,get_persistence_of_pair(l,l_left));
	if (right) gray = max(gray,get_persistence_of_pair(l,l_right));

	gray = (gray < 0) ? 1.0 : pow(gray / max_persistence,0.5f);	
#endif

	double r,g,b;

	if (is_boundary){
	  gray = 0.1 + 0.9 * gray;
	  r = gray;
	  g = gray;
	  b = gray;
	}
	else{
	  
#ifdef COLOR_MODE
	  const vector_color_type& v = input(PIXEL) * shade_factor;
	  RGB_converter_to_Lab<real_type>::
	    from_cartesian_to_RGB(v[0],v[1],v[2],&r,&g,&b);	  
#endif

#ifdef GRAY_MODE
	  const real_type v = input(PIXEL) * shade_factor;
	  RGB_converter_to_Lab<real_type>::
	    from_cartesian_to_RGB(v,0,0,&r,&g,&b);
#endif
	
	}

	output(PIXEL) = rgb_color_type(r,g,b);
	  	
#ifdef MOVIE_MODE
      } // END OF for f
#endif
    } // END OF for y
  } // END OF for x
}





void color_cluster(const SR_domain_type&    SR_dom,
		   const output_label_type& label_data,
		   const size_type          size_threshold,
		   /*const*/ label_list_type&   relabel,
		   rgb_data_type* const     out){

#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::dense_type& SR_domain = *(SR_dom.dense);
  
#else
  
  const SR_domain_type& SR_domain = SR_dom;
  
#endif
  
  rgb_data_type& output = *out;
  
  output.assign(SIZES,rgb_color_type::black);

  size_type colored_pixels = 0;  

  
  FOR_ALL_PIXELS {
    
    label_type l = get_label(label_data(PIXEL),relabel);

#ifdef APPLY_SIZE_THRESHOLD
    if ((x > 0) && (l > 0) && (cluster_size[l] < size_threshold)){
      const label_type ll = get_label(label_data(PREVIOUS_PIXEL),relabel);
      relabel[l] = ll;
      l = ll;
    }
#endif
    
    if (l > 0){
      
#ifdef COLOR_CODED_CLUSTERS
      
      unsigned char r,g,b;

#ifdef OUTPUT_HIERARCHY
      
      b = l % 256; l /= 256;
      g = l % 256; l /= 256;
      r = l % 256;      
      
#else
      
      Color_tools::RGB_from_index(l,&r,&g,&b);

#endif // OUTPUT_HIERARCHY

      output(PIXEL) = rgb_color_type(r/255.0,g/255.0,b/255.0);
 
#else // COLOR_CODED_CLUSTERS
      
      double r,g,b;
      const SR_value_type v = SR_domain[cluster[l].summit_index];
      const real_type     d = DENSITY(v);
      const realN_type    p = point_from_index(v / d);
      
#ifdef GRAY_MODE
      
      RGB_converter_to_Lab<real_type>::
	from_cartesian_to_RGB(p[S_DIMENSION],0,0,&r,&g,&b);

      output(PIXEL) = rgb_color_type(r,g,b);
      
#endif // GRAY_MODE

#ifdef COLOR_MODE
      RGB_converter_to_Lab<real_type>::
	from_cartesian_to_RGB(p[S_DIMENSION + 0],
			      p[S_DIMENSION + 1],
			      p[S_DIMENSION + 2],&r,&g,&b);

      output(PIXEL) = rgb_color_type(r,g,b);
#endif // COLOR_MODE
      
#endif // COLOR_CODED_CLUSTERS
      
	  
      colored_pixels++;
    }	
  } 
}









void report_occupancy_ratio(const SR_domain_type& SR_domain){

#if (ARRAY_DENSITY == DENSE) 
  
  size_type count = 0;

  for(SR_domain_type::const_iterator i = SR_domain.begin();
      i != SR_domain.end();i++){

    count += (DENSITY(*i) != 0) ? 1 : 0;
  }

  cout<<"Occupancy ratio: "<<count<<"/"<<SR_domain.size()<<" = "
      <<(100*count/SR_domain.size())<<"%"<<endl;
  
#endif


  
#if (ARRAY_DENSITY == SPARSE)
  cout<<"Occupancy ratio: "<<SR_domain.size()<<"/"<<n_bins<<" = "
      <<(100*SR_domain.size()/n_bins)<<"% [sparse structure]"<<endl;  
#endif

  
  
#if (ARRAY_DENSITY == ADAPTIVE)

  if (SR_domain.dense != NULL){
    size_type count = 0;

    for(SR_domain_type::dense_type::const_iterator i = SR_domain.dense->begin();
	i != SR_domain.dense->end();i++){

      count += (DENSITY(*i) != 0) ? 1 : 0;
    }

    cout<<"Occupancy ratio: "<<count<<"/"<<SR_domain.dense->size()<<" = "
	<<(100*count/SR_domain.dense->size())<<"%"<<endl;
  }
  else{
    cout<<"Occupancy ratio: "<<SR_domain.sparse->size()<<"/"<<n_bins<<" = "
	<<(100*SR_domain.sparse->size()/n_bins)<<"% [sparse structure]"<<endl;  
  }
  
#endif
}






void write_rgb_data(const string         header,
		    const rgb_data_type& data){

#ifdef IMAGE_MODE
  PPM_tools::write_ppm((header + string(".ppm")).c_str(),data);
#endif

#ifdef MOVIE_MODE
  
  rgb_image_type output(width,height);
  
  for(size_type f = 0,f_end = depth;f < f_end;f++){
    
    for(size_type x = 0,x_end = width;x < x_end;x++){
      for(size_type y = 0,y_end = height;y < y_end;y++){

	output(x,y) = data(x,y,f);
      }
    }

    ostringstream name_out;
    name_out<<"movie/"<<header<<setw(4)<<setfill('0')<<f<<".png";
    PPM_tools::write_ppm(name_out.str().c_str(),output);    
  }
  
#endif
}




#if (ARRAY_DENSITY == ADAPTIVE)

void sparse_to_dense(SR_domain_type* const      SR_dom){

  SR_domain_type& SR_domain      = *SR_dom;
  
  SR_domain.dense = new SR_domain_type::dense_type(size);

  SR_domain_type::dense_type&  dense  = *(SR_domain.dense);
  SR_domain_type::sparse_type& sparse = *(SR_domain.sparse);

  static intN_type all_one(vector<key_coef_type>(space_dimension,1));
  static const intN_type corner1 = intN_type();
  static const intN_type corner2 = size - all_one;
  
  for(SR_domain_type::sparse_type::iterator
	s     = sparse.begin(),
	s_end = sparse.end();
      s != s_end;
      ++s){

    const intN_type& index = s->first;

    if (Geom_tools::inside_bounding_box(corner1,corner2,index)){
      dense[s->first] = s->second;
    }
  }

  delete SR_domain.sparse;
  SR_domain.sparse = NULL;
}

#endif




int main(int argc,char** argv){

#ifdef MOVIE_MODE
  const size_type time_proxy = 1;
#else
  const size_type time_proxy = 0;
#endif  


  if (argc != 5+time_proxy){
    cerr<<"\n"
	<<"usage: "<<argv[0]
	<<" file_name"
	<<" range_sigma"
#ifdef MOVIE_MODE
	<<" time_sigma"
#endif  
	<<" space_sigma"
	<<" persistence_threshold"
	<<"\n"
	<<endl;

    exit(1);
  }
  
  display_options();

  string input_name;
  real_type space_sigma,range_sigma,pers_thr;
  real_type sampling_factor = 1.0;
  real_type n_sub = 2.0;
  
#ifdef MOVIE_MODE
  real_type time_sigma;
#endif

  cluster.push_back(cluster_type());
  
  istringstream input_name_in(argv[1]);
  input_name_in >> input_name;

  istringstream range_sigma_in(argv[2]);
  range_sigma_in >> range_sigma;

#ifdef MOVIE_MODE
  istringstream time_sigma_in(argv[3]);
  time_sigma_in >> time_sigma;
#endif  
  
  istringstream space_sigma_in(argv[3+time_proxy]);
  space_sigma_in >> space_sigma;
  
  istringstream persistence_in(argv[4+time_proxy]);
  persistence_in >> pers_thr;

  
  
#ifdef IMAGE_MODE
  const size_type size_threshold =
    static_cast<size_type>(space_sigma * space_sigma / 2.0);
#endif

  
#ifdef MOVIE_MODE
  const size_type size_threshold =
    static_cast<size_type>(space_sigma * space_sigma * time_sigma / 2.0);
#endif

  
  
  if (verbose) cout<<"Load the input data... "<<endl;

  input_data_type input;

#ifdef IMAGE_MODE
  
  rgb_image_type* rgb_input = new rgb_image_type;
  PPM_tools::read_ppm(input_name.c_str(),rgb_input);

  width  = rgb_input->width();
  height = rgb_input->height();

  input.resize(SIZES);

#ifdef GRAY_MODE

  
  for(size_type x = 0;x < width;x++){
    for(size_type y = 0;y < height;y++){
      real_type r,g,b,L,A,B;
      (*rgb_input)(x,y).get_RGB(&r,&g,&b);
      RGB_converter_to_Lab<real_type>::from_RGB_to_cartesian(r,g,b,&L,&A,&B);
// 	cout<<r<<" "<<g<<" "<<b<<"\t"<<L<<" "<<A<<" "<<B<<endl;
      input(x,y) = L;
    }
  }
        
#endif
  
#ifdef COLOR_MODE

  for(size_type x = 0;x < width;x++){
    for(size_type y = 0;y < height;y++){
      for(size_type c = 0;c < 3;c++){

	real_type r,g,b,L,A,B;
	(*rgb_input)(x,y).get_RGB(&r,&g,&b);
	RGB_converter_to_Lab<real_type>::from_RGB_to_cartesian(r,g,b,&L,&A,&B);
	input(x,y)[0] = L;
	input(x,y)[1] = A;
	input(x,y)[2] = B;
      }
    }
  }
  
#endif
  
  delete rgb_input;

  

  if (verbose) cout<<"  File: "<<input_name<<" ("<<input.width()<<" x "
		   <<input.height()<<" = "<<input.size()<<")"<<endl;
  
#endif // END OF image mode


#ifdef MOVIE_MODE
  
  ifstream in(input_name.c_str());

  string path;
  getline(in,path);
  
  vector<string> file_name;
  string line;
  while(getline(in,line)){
    file_name.push_back(line);
  }

  depth = file_name.size();
  
  
  rgb_image_type* rgb_frame = new rgb_image_type;

  for(size_type f = 0;f < depth;f++){

    const string name = path + file_name[f];
    
    PPM_tools::read_ppm(name.c_str(),rgb_frame);

    if (f == 0){
      width  = rgb_frame->width();
      height = rgb_frame->height();

      input.resize(SIZES);
    }

    for(size_type x = 0;x < width;x++){
      for(size_type y = 0;y < height;y++){
	
	real_type r,g,b,L,A,B;
	(*rgb_frame)(x,y).get_RGB(&r,&g,&b);
	RGB_converter_to_Lab<real_type>::from_RGB_to_cartesian(r,g,b,&L,&A,&B);

#ifdef GRAY_MODE
	input(x,y,f) = L;
#endif

#ifdef COLOR_MODE
	input(x,y,f)[0] = L;
	input(x,y,f)[1] = A;
	input(x,y,f)[2] = B;
#endif
	
      } // END OF for y
    } // END OF for x

    if (verbose) cout<<"\r  "<<(f+1)<<" / "<<depth<<" frames loaded"<<flush;
  } // END OF for f

  if (verbose) cout<<"\n";

  delete rgb_frame;

  if (verbose) {
    cout<<"  File list: "<<input_name
	<<" ["<<input.width()<<" x "<<input.height()<<" x "
	<<input.depth()<<" = "<<input.size()<<" pixels";
    if (input.size() > 1e6){
      cout<<" ("<<((input.size() / 100000) / 10.0)<<" Mpixels)";
    }
    cout<<"]"<<endl;
  }
  
#endif // END OF movie mode
  
  if (verbose) cout<<"Done\n"<<endl;

  

  

  if (verbose) cout<<"Create the SxR domain... "<<endl;

  SR_domain_type      SR_domain;
  bool_SR_domain_type occupancy;
  feature_array_type  feature_image;
  
  
  init(input,
       space_sigma,
#ifdef MOVIE_MODE
       time_sigma,
#endif      
       range_sigma,
       sampling_factor,n_sub,
       &feature_image);
  
  n_bins = 1;

  if (verbose) {
    cout<<"  "<<init_chrono.report()<<endl;
    cout<<"  "<<"Domain size: ";
    for(size_type n = 0;n < space_dimension;n++){
      cout<<static_cast<int>(size[n]);
      n_bins *= size[n];
      if (n != space_dimension-1)
	cout<<" x ";
    }
    cout<<" = "<<n_bins<<" bins";
    if (n_bins > 1e6) {
      cout<<" ("<<((n_bins / 100000) / 10.0)<<" Mbins)";
    }
    cout<<endl;
  }
  
  
  fill_SR_domain(feature_image,&SR_domain,&occupancy);

  
  if (verbose) {
    cout<<"  "; report_occupancy_ratio(SR_domain);
    cout<<"  "<<creation_chrono.report()<<endl;
    
    cout<<"Done\n"<<endl;
  }
  


  


  if (verbose) cout<<"Filter the SxR domain... "<<endl;

  filter_chrono.start();
  
  filter_SR_domain(&SR_domain,sampling_factor);

  filter_chrono.stop();
  
  if (verbose) {
    cout<<"  ";  report_occupancy_ratio(SR_domain);
    cout<<"  "<<filter_chrono.report()<<endl;  
    cout<<"Done\n"<<endl;
  }



  


  label_SR_domain_type label;
  

  if (verbose) cout<<"Apply safe labels... "<<endl;

  index_list_type index_list;
  
  sort_chrono.start();

  build_sorted_bin_list(SR_domain,occupancy,&index_list);
  
#if (ARRAY_DENSITY == ADAPTIVE)

  sparse_to_dense(&SR_domain);

#endif
  
  sort_chrono.stop();
  if (verbose) cout<<"  "<<sort_chrono.report()<<endl;

  safe_labels_chrono.start();

  label_list_type safe_relabel;

  safe_label(SR_domain,index_list,
#ifdef LABEL_ON_THE_FLY
	     pers_thr,
#endif
	     &label,&safe_relabel);
  
  safe_labels_chrono.stop();
  
  if (verbose) {
    cout<<"  "<<safe_labels_chrono.report()<<endl;
    cout<<"Done\n"<<endl;
  }




  
  if (verbose) cout<<"Early classification...\n";
  
  output_label_type label_data;
  rgb_data_type     output;

  label_list_type no_relabel;
  init_relabel(&no_relabel);

  label_pixels(input,SR_domain,label,false,&label_data);
  
#ifdef LABEL_ON_THE_FLY

  color_cluster(SR_domain,label_data,0,safe_relabel,&output);

#else
  
  color_cluster(SR_domain,label_data,0,no_relabel,&output);

#endif
  
  write_rgb_data("partial",output);

  if (verbose) cout<<"Done\n"<<endl;

  
  



  

  if (verbose) cout<<"Label remaining pixels... "<<endl;

  pixel_chrono.start();

  label_pixels(input,SR_domain,label,true,&label_data);
  
  pixel_chrono.stop();
  
  if (verbose) {
    cout<<"  "<<pixel_chrono.report()<<endl;
    cout<<"Done\n"<<endl;
  }

#ifdef LABEL_ON_THE_FLY
  
  color_cluster(SR_domain,label_data,size_threshold,safe_relabel,&output);
  
#else
  
  label_list_type relabel;
  init_relabel(&relabel);

#ifdef OUTPUT_HIERARCHY
  color_cluster(SR_domain,label_data,size_threshold,relabel,&output);
#endif
  
  prune_persistence(SR_domain,pers_thr,&relabel);
  
#ifndef OUTPUT_HIERARCHY  
  color_cluster(SR_domain,label_data,size_threshold,relabel,&output);
#endif
  
#endif
   
  write_rgb_data("full",output);

  rgb_data_type boundaries;

#ifdef IMAGE_MODE


#ifdef LABEL_ON_THE_FLY

  draw_boundaries(input,SR_domain,label_data,safe_relabel,&boundaries);
  
#else

  draw_boundaries(input,SR_domain,label_data,relabel,&boundaries);

#endif

  write_rgb_data("boundaries",boundaries);

#endif

//   draw_boundaries(input,SR_domain,label_data,safe_relabel,&boundaries);
//   write_rgb_data("boundaries",boundaries);

  
  const real_type init_time     = init_chrono.time_in_seconds();
  const real_type creation_time = creation_chrono.time_in_seconds();
  const real_type filter_time   = filter_chrono.time_in_seconds();
  const real_type sort_time     = sort_chrono.time_in_seconds();
  const real_type cluster_time  = safe_labels_chrono.time_in_seconds();
  const real_type pixel_time    = pixel_chrono.time_in_seconds();
  real_type total_time    =
    init_time + creation_time + filter_time + sort_time + cluster_time + pixel_time;

#ifndef LABEL_ON_THE_FLY
  const real_type prune_time    = prune_chrono.time_in_seconds();
  total_time += prune_time;
#endif
  
#define TIME_OUTPUT(t) (t)<<"s\t("<<static_cast<int>(0.5+100*(t)/total_time)<<"%)"
  
  cout<<"\nRunning time summary\n--------------------\n";
  cout<<"  init:          "<<TIME_OUTPUT(init_time)<<endl;
  cout<<"  creation:      "<<TIME_OUTPUT(creation_time)<<endl;
  cout<<"  filter:        "<<TIME_OUTPUT(filter_time)<<endl;
  cout<<"  sort:          "<<TIME_OUTPUT(sort_time)<<endl;
  cout<<"  cluster label: "<<TIME_OUTPUT(cluster_time)<<endl;
#ifndef LABEL_ON_THE_FLY
  cout<<"  pruning:       "<<TIME_OUTPUT(prune_time)<<endl;
#endif
  cout<<"  pixel label:   "<<TIME_OUTPUT(pixel_time)<<endl;
  cout<<"  TOTAL TIME:    "<<total_time<<"s\n"<<endl;

#undef TIME_OUTPUT
  
}