mean_shift.cpp

mean-shift. pointer sample

    }

    mask = buffer;
  }
}




void filter_SR_domain(SR_domain_type* const SR_dom,
		      const real_type       /*sampling_factor*/){

  
#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::sparse_type& SR_domain = *(SR_dom->sparse);

  typedef SR_domain_type::sparse_type buffer_type;
  
#else
  
  SR_domain_type& SR_domain = *SR_dom;

  typedef SR_domain_type buffer_type;
  
#endif

  
#if ARRAY_DENSITY == DENSE
  
  const real_type delta = 0.5 / n_substeps;
  
  SR_domain_type buffer(size);

  intN_type start_index(vector<key_coef_type>(space_dimension,1));

  const intN_type end_index = size - start_index;

  for(size_type n = 0;n < SR_domain_type::dimension;n++){

    intN_type zero_index,offset_index;
    offset_index[n] = 1;

    const SR_domain_type::difference_type offset
      = &(SR_domain[offset_index]) - &(SR_domain[zero_index]);

    for(size_type iter = 0;iter < n_iterations;iter++){

      memcpy(&(*buffer.begin()),
	     &(*SR_domain.begin()),
	     SR_domain.size() * sizeof(SR_domain_type::value_type));
	
      intN_type index = start_index;
      do{
	
	const SR_value_type& b = buffer[index];

	if (DENSITY(b) > 0){
	  const SR_value_type v  = b * delta;
	  
	  SR_value_type* const ptr = &(SR_domain[index]);
	  *(ptr-offset) += v;
	  *(ptr+offset) += v;
	  *ptr          -= v * 2.0;
	}
      }

      while(SR_domain.advance(&index,start_index,end_index));

    } // END OF for iter
  } // END OF for n

#endif



  
#if (ARRAY_DENSITY == SPARSE) || (ARRAY_DENSITY == ADAPTIVE)
  
  buffer_type buffer;
  
  buffer.resize(size);
  vector<real_type> mask;
  create_filter_mask(n_iterations,&mask);
  
  for(size_type n=0;n<buffer_type::dimension;n++){

    SR_domain.swap(buffer); SR_domain.clear();

    for(buffer_type::iterator
	  i     = buffer.begin(),
	  i_end = buffer.end();
	i != i_end;
	i++){
      
      intN_type index        = i->first;
      const SR_value_type& v = i->second;

      index[n] -= n_iterations;
      
      for(size_type
	    o     = 0,
	    o_end = 2 * n_iterations + 1;
	  o <= o_end;o++,index[n]++){

	SR_domain[index] += v * mask[o];
      } // END OF for o      
    } // END OF for i   
  } // END OF for n

#endif
  
}






void build_sorted_bin_list(const SR_domain_type&      SR_dom,
#ifdef PROCESS_EMPTY_BINS
			   const bool_SR_domain_type&,
#else
			   const bool_SR_domain_type& occ,
#endif
			   
			   index_list_type* const     bin_list){

  index_list_type& index_list = *bin_list;

  valued_index_list_type valued_index_list;

#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::sparse_type& SR_domain      = *(SR_dom.sparse);
  bool_SR_domain_type::sparse_type& occupancy = *(occ.sparse);

  typedef SR_domain_type::sparse_type domain_type;
  
#else
  
  const SR_domain_type&      SR_domain = SR_dom;
  const bool_SR_domain_type& occupancy = occ;

  typedef SR_domain_type domain_type;

#endif


  
#if (ARRAY_DENSITY == DENSE)
  
  static intN_type all_one(vector<key_coef_type>(space_dimension,1));

  const intN_type start_index = all_one;
  const intN_type end_index   = size - all_one;
  
  intN_type index = start_index;
  do{

//     const real_type v = SR_domain[index][value_dimension-1];
    const real_type v = DENSITY(SR_domain[index]);

#ifdef PROCESS_EMPTY_BINS
    if (v > 0)
#else
      if ((v > 0) && occupancy[index])
#endif
      {
	const valued_index_type valued_index(v,index);
	valued_index_list.push_back(valued_index);
      }
  }
  while(SR_domain.advance(&index,start_index,end_index));
  
#endif


  
  
#if (ARRAY_DENSITY == SPARSE) || (ARRAY_DENSITY == ADAPTIVE)

  for(domain_type::iterator
	i     = const_cast<domain_type&>(SR_domain).begin(),
	i_end = const_cast<domain_type&>(SR_domain).end();
      i != i_end;i++){

//     const real_type v = i->second[value_dimension-1];
    const real_type v = DENSITY(i->second);
    
#ifdef PROCESS_EMPTY_BINS
    if (v > 0)
#else
      if ((v > 0) && occupancy[i->first])
#endif
      {
	const valued_index_type valued_index(v,i->first);
	valued_index_list.push_back(valued_index);
      }
    
  }
  
#endif

  sort(valued_index_list.begin(),valued_index_list.end());

  index_list.clear();
  index_list.reserve(valued_index_list.size());

  for(valued_index_list_type::reverse_iterator
	i     = valued_index_list.rbegin(),
	i_end = valued_index_list.rend();
      i != i_end;
      i++){

    index_list.push_back(i->data);
  }
}












label_type classify_point(
  
#if (ARRAY_DENSITY == ADAPTIVE)
  const SR_domain_type::dense_type&       SR_domain,
  const label_SR_domain_type::dense_type& label,
#else
  const SR_domain_type&       SR_domain,
  const label_SR_domain_type& label,
#endif
  const realN_type&           point){

  label_type l = NO_LABEL;
  realN_type p = point;
  realN_type ms,prev_ms;
  
  size_type n = 0;

#ifdef DENSITY_ONLY
  const size_type max_iterations = 20;
#else
  const size_type max_iterations = 20;
#endif
  
  do{

#ifdef DENSITY_ONLY

    Math_tools::Nlinear_interpolation_gradient(SR_domain,p,&ms);
    ms.normalize();
    ms *= 0.2;

    
#else
    
    const SR_value_type mean =
      Math_tools::Nlinear_interpolation(SR_domain,p);

    const real_type density = DENSITY(mean);
    
    for(size_type n=0;n<space_dimension;n++){
      ms[n] = mean[n] / density - p[n];
    }

    
#endif
    
    const real_type dot = ms * prev_ms;
    
    if (dot < 0){
      ms -= prev_ms * dot;
    }
    
    const real_type sq_norm = ms.square_norm();
    const real_type thr = 0.2;
    if (sq_norm < thr*thr){
      ms /= 1.0 / thr * sqrt(sq_norm);
    }
    
    p += ms;
    prev_ms = ms.unit();

    l = label[index_of_real_index(p)];
  }
  while((l <= 0) && (++n < max_iterations));
  return l;
}






void build_neighbor_list(index_list_type* const neighbor){

  index_list_type& n = *neighbor;

#ifdef USE_LARGE_NEIGHBORHOOD
  
  const size_type size = Math_tools::power<space_dimension>(3) - 1;
  n.resize(size);
  
  intN_type all_one(vector<key_coef_type>(space_dimension,1));
  intN_type all_three(vector<key_coef_type>(space_dimension,3));

  Array_ND<space_dimension,bool> proxy(all_three);

  size_type i = 0;
  intN_type index;
  do{
    if (index != all_one){
      n[i] = index - all_one;
      i++;
    }
  }
  while(proxy.advance(&index));
  
#else
  
  const size_type size = space_dimension * 2;
  n.resize(size);

  for(size_type d = 0;d < space_dimension;d++){

    intN_type index;
    index[d] = -1;
    n[2*d]   = index;
    
    index[d] = 1;
    n[2*d+1] = index;
  }

#endif
}




void safe_label(const SR_domain_type&       SR_dom,
		const index_list_type&      index_list,
#ifdef LABEL_ON_THE_FLY
		const real_type             persistence_threshold,
#endif
		label_SR_domain_type* const lab,
		label_list_type* const      relab){

#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::dense_type& SR_domain   = *(SR_dom.dense);

  lab->sparse = NULL;
  lab->dense = new label_SR_domain_type::dense_type;
  label_SR_domain_type::dense_type& label = *(lab->dense);
  
#else
  
  const SR_domain_type& SR_domain = SR_dom;
  label_SR_domain_type& label     = *lab;
  
#endif
  
  static intN_type all_one(vector<key_coef_type>(space_dimension,1));
  static intN_type all_two(vector<key_coef_type>(space_dimension,2));
  
  label_list_type& relabel = *relab;
  relabel.clear();
  relabel.push_back(0);

  size_type labelled_bins = 0;
  size_type easy_bins = 0;
  size_type saved_bins = 0;

  label.resize(size);
  
  label_type max_ink = 0;

  const size_type n_neighbors = Math_tools::power<space_dimension>(3)-1;

  index_list_type neighbor_offset(n_neighbors);
  build_neighbor_list(&neighbor_offset);

  for(index_list_type::const_iterator
	i     = index_list.begin(),
	i_end = index_list.end();
      i != i_end;
      i++){

    const intN_type& index = *i;
    
    valued_label_list_type neighbor;
    label_set_type touching_cluster;
      
    for(size_type o=0 ; o<n_neighbors ; ++o){
    
      const intN_type sub_index = index + neighbor_offset[o];
    
      const label_type l = label[sub_index];

      if (l > 0){
	const real_type d   = DENSITY(SR_domain[sub_index]);
	neighbor.push_back(valued_label_type(d,l));
	touching_cluster.insert(l);
      }
    }

    label_type ink = NO_LABEL;

    switch (touching_cluster.size()){
      
    case 0: // new summit found, create a new cluster
      max_ink++;
      ink = max_ink;
      cluster.push_back(cluster_type(ink,index,
				     DENSITY(SR_domain[index])));
      relabel.push_back(max_ink);
      break;
      
    case 1: // non-ambiguous labeling, just copy the label
      ink = neighbor[0].data;
      break;

    default: // ambiguous labeling, two or more clusters nearby
    {
      const label_type default_label
	= max_element(neighbor.begin(),neighbor.end())->data;
      ink = - 2 * default_label;

      label_set_type::iterator l1 = touching_cluster.begin();
      label_set_type::iterator post_l1 = l1;
      post_l1++;
	
      label_set_type::iterator l1_end = touching_cluster.end();
      for(; post_l1 != l1_end ; l1 = post_l1 , ++post_l1){

	cluster_type& c1 = cluster[*l1];
	  
	label_set_type::iterator l2 = l1;
	l2++;
	for(label_set_type::iterator l2_end = touching_cluster.end();
	    l2 != l2_end;l2++){

	  cluster_type& c2 = cluster[*l2];

#ifdef LABEL_ON_THE_FLY
	  if ((get_label(c1.label,relabel) == get_label(c2.label,relabel))){
	    ink = - 2 * c1.label - 1;
	  }
	  else 
#endif
	  if (!cluster_type::are_neighbors(c1,c2,relabel)){
	    
	    saddle_type s(c1,c2,SR_domain[index]);

#ifdef LABEL_ON_THE_FLY
	    update_saddle(SR_domain,relabel,&s);
	    
	    if (persistence(SR_domain,s) < persistence_threshold){
	      relabel[s.minor] = s.major;
	      ink = - 2 * s.major - 1; // indicate that it does not need to be refined
	      if (c2 < c1) c1.neighbor.insert(c1.neighbor.end(),
					      c2.neighbor.begin(),
					      c2.neighbor.end());
	      else c2.neighbor.insert(c2.neighbor.end(),
				      c1.neighbor.begin(),
				      c1.neighbor.end());
	    }
	    else{
	      saddle.push_back(s);	      
	      c1.neighbor.push_back(c2.label);
	      c2.neighbor.push_back(c1.label);
	    }
#else			    
	    saddle.push_back(s);	      
	    c1.neighbor.push_back(c2.label);
	    c2.neighbor.push_back(c1.label);
#endif
	  } // END OF if !are_neighbors
	} // END OF for l2
      } // END OF for l1
            
    } // END OF default case 
    } // END OF switch touching_cluster.size()
    
    if (ink > 0) labelled_bins++;
    
    label[index] = ink;
  }
  
  if (verbose) cout<<"  Created clusters: "<<max_ink<<'\n'
		   <<"  Labelled bins: "<<labelled_bins<<" on "<<index_list.size()
		   <<" ("<<(100*labelled_bins/index_list.size())<<"%)\n"
		   <<"  Easy bins:     "<<easy_bins
		   <<" ("<<(100*easy_bins/index_list.size())<<"%)\n"
		   <<"  Saved bins:    "<<saved_bins
		   <<" ("<<(100*saved_bins/index_list.size())<<"%)"
		   <<endl;
}






void label_pixels(const input_data_type&      input,
		  const SR_domain_type&       SR_dom,
		  const label_SR_domain_type& lab,
		  const bool                  fill_holes,
		  output_label_type* const    out){

#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::dense_type& SR_domain   = *(SR_dom.dense);
  label_SR_domain_type::dense_type& label = *(lab.dense);
  
#else
  
  const SR_domain_type& SR_domain   = SR_dom;
  const label_SR_domain_type& label = lab;
  
#endif

  output_label_type& label_image = *out;
  label_image.assign(SIZES,NO_LABEL);

  intN_type index;
  size_type labelled_pixels = 0;

  cluster_size.assign(cluster.size(),0);
    
  if (fill_holes){
    
    FOR_ALL_PIXELS {	  

      input_feature_type input_feature;
      input_feature_from_pixel(PIXEL,input,&input_feature);

      GET_FEATURE(input_feature);

      const realN_type r = real_index_of_feature(feature);
      const intN_type i  = index_of_real_index(r);
      
      label_type l = label[i];

      if (l <= 0){
	
	label_type point_label;
	if (l%2 == 0){
	  point_label = classify_point(SR_domain,label,r);
	}
	else{
	  point_label = -l/2;
	}

	l = (point_label > 0) ? point_label :
	  ((point_label < 0 ) ? -point_label/2 : -l/2);
      }
      
      if (l > 0){
	labelled_pixels++;
	cluster_size[l]++;
	label_image(PIXEL) = l;
      }

    } 
  }
  else{

    FOR_ALL_PIXELS {	  
      input_feature_type input_feature;
      input_feature_from_pixel(PIXEL,input,&input_feature);

      GET_FEATURE(input_feature);

      const realN_type r = real_index_of_feature(feature);
      const intN_type i  = index_of_real_index(r);	  
      label_type l       = label[i];      

      if ((l<0) && (l%2 != 0)){
	l = -l/2;
      }
      
      if (l>0){
	labelled_pixels++;
	cluster_size[l]++;
	label_image(PIXEL) = l;
      }     
    } 
  }

  if (verbose) cout<<"  Labelled pixels: "<<labelled_pixels
		   <<" on "<<input.size()<<" ("
		   <<(static_cast<long int>(labelled_pixels)*100/input.size())<<"%)"
		   <<endl;  
}





void draw_boundaries(const input_data_type&   input,
		     const SR_domain_type&    SR_dom,
		     const output_label_type& label_data,
		     const label_list_type&   relabel,
		     rgb_data_type* const     out){

  const real_type shade_factor = 0.0;
  
#if (ARRAY_DENSITY == ADAPTIVE)
  
  SR_domain_type::dense_type& SR_domain = *(SR_dom.dense);
  
#else
  
  const SR_domain_type& SR_domain = SR_dom;
  
#endif
  
  real_type max_persistence = 0;
  
  for(saddle_list_type::iterator
	s     = saddle.begin(),
	s_end = saddle.end();
      s != s_end;s++){

    update_saddle(SR_domain,relabel,&(*s));
    max_persistence = max(persistence(SR_domain,*s),max_persistence);
  }

  
  rgb_data_type& output = *out;
  output.resize(SIZES);

  real_type gray = 1.0;
  
  for(size_type x = 0,x_end = width;x < x_end;x++){

    const size_type xm = (x > 0) ? (x - 1) : x; 
    const size_type xp = (x < width-1) ? (x + 1) : x; 
    
    for(size_type y = 0,y_end = height;y < y_end;y++){

      const size_type ym = (y > 0) ? (y - 1) : y; 
      const size_type yp = (y < height-1) ? (y + 1) : y; 
      
#ifdef MOVIE_MODE
      for(size_type f = 0,f_end = depth;f < f_end;f++){

// 	const size_type fm = (f > 0) ? (f - 1) : f; 
// 	const size_type fp = (f < depth-1) ? (f + 1) : f; 

	const label_type l       = get_label(label_data(x,y, f),relabel);
	const label_type l_up    = get_label(label_data(x,yp,f),relabel);
	const label_type l_down  = get_label(label_data(x,ym,f),relabel);
	const label_type l_left  = get_label(label_data(xp,y,f),relabel);
	const label_type l_right = get_label(label_data(xm,y,f),relabel);

	const bool up    = (l != l_up)    && (yp != height-1) && (ym != 0);
	const bool down  = (l != l_down)  && (yp != height-1) && (ym != 0);
	const bool left  = (l != l_left)  && (xp != width-1) && (xm != 0);
	const bool right = (l != l_right) && (xp != width-1) && (xm != 0);
	
	const bool is_boundary = (up || down || left || right);
	
	gray = -1.0;

	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

#ifdef IMAGE_MODE	
	const label_type l       = get_label(label_data(x,y),relabel);
	const label_type l_up    = get_label(label_data(x,yp),relabel);
	const label_type l_down  = get_label(label_data(x,ym),relabel);
	const label_type l_left  = get_label(label_data(xp,y),relabel);
	const label_type l_right = get_label(label_data(xm,y),relabel);

	const bool up    = (l != l_up)    && (yp != height-1) && (ym != 0);
	const bool down  = (l != l_down)  && (yp != height-1) && (ym != 0);
	const bool left  = (l != l_left)  && (xp != width-1) && (xm != 0);
	const bool right = (l != l_right) && (xp != width-1) && (xm != 0);
	
	const bool is_boundary = (up || down || left || right);

	gray = -1.0;