gcoptimization.cpp

Graph Cut algorithm implementation. Includes MATLAB compiled codes.

	terminateOnError( alpha_label < 0 || alpha_label >= m_num_labels || beta_label < 0 || beta_label >= m_num_labels,
		"Illegal Label to Expand On");
	perform_alpha_beta_swap(alpha_label,beta_label);
	return(compute_energy());
}
/**************************************************************************************/

void GCoptimization::add_t_links_ARRAY_swap(Energy *e,Energy::Var *variables,int size,
											LabelType alpha_label, LabelType beta_label,
											PixelType *pixels)
{
	for ( int i = 0; i < size; i++ )
		e -> add_term1(variables[i], m_datacost(pixels[i],alpha_label),
									 m_datacost(pixels[i],beta_label));

}
	
/**************************************************************************************/

void GCoptimization::add_t_links_FnPix_swap(Energy *e,Energy::Var *variables,int size,
											LabelType alpha_label, LabelType beta_label,
											PixelType *pixels)
{
	for ( int i = 0; i < size; i++ )
		e -> add_term1(variables[i], m_dataFnPix(pixels[i],alpha_label),
									 m_dataFnPix(pixels[i],beta_label));

}
/**************************************************************************************/

void GCoptimization::add_t_links_FnCoord_swap(Energy *e,Energy::Var *variables,int size,
											  LabelType alpha_label, LabelType beta_label,
											  PixelType *pixels)
{
	int x,y;

	for ( int i = 0; i < size; i++ )
	{

		y = pixels[i]/m_width;
		x = pixels[i] - y*m_width;

		e -> add_term1(variables[i], m_dataFnCoord(x,y,alpha_label),m_dataFnCoord(x,y,beta_label));
	}
}


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

void  GCoptimization::perform_alpha_beta_swap(LabelType alpha_label, LabelType beta_label)
{
	PixelType i,size = 0;
	Energy *e = new Energy(mexErrMsgTxt);
	PixelType *pixels = new PixelType[m_num_pixels];
	

	for ( i = 0; i < m_num_pixels; i++ )
	{
		if ( m_labeling[i] == alpha_label || m_labeling[i] == beta_label)
		{
			pixels[size]    = i;
			m_lookupPixVar[i] = size;
			size++;
		}
	}

	if ( size == 0 )
	{
		delete e;
		delete [] pixels;
		return;
	}


	Energy::Var *variables = (Energy::Var *) new Energy::Var[size];


	for ( i = 0; i < size; i++ )
		variables[i] = e ->add_variable();

	if ( m_dataType == ARRAY ) add_t_links_ARRAY_swap(e,variables,size,alpha_label,beta_label,pixels);
	else  if  ( m_dataType == FUNCTION_PIX ) add_t_links_FnPix_swap(e,variables,size,alpha_label,beta_label,pixels);
	else  add_t_links_FnCoord_swap(e,variables,size,alpha_label,beta_label,pixels);



	if ( m_grid_graph )
	{
		if ( m_smoothType == ARRAY )
		{
			if (m_varying_weights) set_up_swap_energy_G_ARRAY_VW(size,alpha_label,beta_label,pixels,e,variables);
			else set_up_swap_energy_G_ARRAY(size,alpha_label,beta_label,pixels,e,variables);
		}
		else if ( m_smoothType == FUNCTION_PIX ) set_up_swap_energy_G_FnPix(size,alpha_label,beta_label,pixels,e,variables);
		else  set_up_swap_energy_G_FnCoord(size,alpha_label,beta_label,pixels,e,variables);
			
	}
	else
	{
		if ( m_smoothType == ARRAY ) set_up_swap_energy_NG_ARRAY(size,alpha_label,beta_label,pixels,e,variables);
		else if ( m_smoothType == FUNCTION_PIX ) set_up_swap_energy_NG_FnPix(size,alpha_label,beta_label,pixels,e,variables);
	}
		

	e -> minimize();

	for ( i = 0; i < size; i++ )
		if ( e->get_var(variables[i]) == 0 )
			m_labeling[pixels[i]] = alpha_label;
		else m_labeling[pixels[i]] = beta_label;


	delete [] variables;
	delete [] pixels;
	delete e;

}

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

void GCoptimization::set_up_swap_energy_NG_ARRAY(int size,LabelType alpha_label,LabelType beta_label,
												 PixelType *pixels,Energy* e, Energy::Var *variables)
{
	PixelType nPix,pix,i;
	EnergyTermType weight;
	Neighbor *tmp;
	


	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		if ( !m_neighbors[pix].isEmpty() )
		{
			m_neighbors[pix].setCursorFront();
			
			while ( m_neighbors[pix].hasNext() )
			{
				tmp = (Neighbor *) (m_neighbors[pix].next());
				nPix   = tmp->to_node;
				weight = tmp->weight;
				
				if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				{
					if ( pix < nPix )
						e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
							          m_smoothcost(alpha_label,alpha_label)*weight,
									  m_smoothcost(alpha_label,beta_label)*weight,
									  m_smoothcost(beta_label,alpha_label)*weight,
									  m_smoothcost(beta_label,beta_label)*weight);
				}
				else
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,
						                       m_smoothcost(beta_label,m_labeling[nPix])*weight);
			}
		}
	}
}

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

void GCoptimization::set_up_swap_energy_NG_FnPix(int size,LabelType alpha_label,LabelType beta_label,
												 PixelType *pixels,Energy* e, Energy::Var *variables)
{
	PixelType nPix,pix,i;
	Neighbor *tmp;
	

	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		if ( !m_neighbors[pix].isEmpty() )
		{
			m_neighbors[pix].setCursorFront();
			
			while ( m_neighbors[pix].hasNext() )
			{
				tmp = (Neighbor *) (m_neighbors[pix].next());
				nPix   = tmp->to_node;
				
				if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				{
					if ( pix < nPix )
						e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
							          m_smoothFnPix(pix,nPix,alpha_label,alpha_label),
									  m_smoothFnPix(pix,nPix,alpha_label,beta_label),
									  m_smoothFnPix(pix,nPix,beta_label,alpha_label),
									  m_smoothFnPix(pix,nPix,beta_label,beta_label) );
				}
				else
					e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
						                       m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));
			}
		}
	}
}

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

void GCoptimization::set_up_swap_energy_G_FnPix(int size,LabelType alpha_label,LabelType beta_label,
												PixelType *pixels,Energy* e, Energy::Var *variables)
{
	PixelType nPix,pix,i,x,y;

	
	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		y = pix/m_width;
		x = pix - y*m_width;

		if ( x > 0 )
		{
			nPix = pix - 1;
	
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothFnPix(pix,nPix,alpha_label,alpha_label),
							  m_smoothFnPix(pix,nPix,alpha_label,beta_label),
							  m_smoothFnPix(pix,nPix,beta_label,alpha_label),
							  m_smoothFnPix(pix,nPix,beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
 					                       m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));
	
		}	
		if ( y > 0 )
		{
			nPix = pix - m_width;
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothFnPix(pix,nPix,alpha_label,alpha_label),
							  m_smoothFnPix(pix,nPix,alpha_label,beta_label),
							  m_smoothFnPix(pix,nPix,beta_label,alpha_label),
							  m_smoothFnPix(pix,nPix,beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
 					                       m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));
		}	

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
					e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
 					                           m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
				e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
				                           m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));

		}
	}
}
/**************************************************************************************/

void GCoptimization::set_up_swap_energy_G_FnCoord(int size,LabelType alpha_label,LabelType beta_label,PixelType *pixels,
     											 Energy* e, Energy::Var *variables)
{  
	PixelType nPix,pix,i,x,y;

	
	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		y = pix/m_width;
		x = pix - y*m_width;

		if ( x > 0 )
		{
			nPix = pix - 1;
	
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_horz_cost(x-1,y,alpha_label,alpha_label),
							  m_horz_cost(x-1,y,alpha_label,beta_label),
							  m_horz_cost(x-1,y,beta_label,alpha_label),
							  m_horz_cost(x-1,y,beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_horz_cost(x-1,y,alpha_label,m_labeling[nPix]),
 					                           m_horz_cost(x-1,y,beta_label,m_labeling[nPix]));
	
		}	
		if ( y > 0 )
		{
			nPix = pix - m_width;
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_vert_cost(x,y-1,alpha_label,alpha_label),
							  m_vert_cost(x,y-1,alpha_label,beta_label),
							  m_vert_cost(x,y-1,beta_label,alpha_label),
							  m_vert_cost(x,y-1,beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_vert_cost(x,y-1,alpha_label,m_labeling[nPix]),
 					                       m_vert_cost(x,y-1,beta_label,m_labeling[nPix]));
		}	

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
					e ->add_term1(variables[i],m_horz_cost(x,y,alpha_label,m_labeling[nPix]),
 					                           m_horz_cost(x,y,beta_label,m_labeling[nPix]));
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
				e ->add_term1(variables[i],m_vert_cost(x,y,alpha_label,m_labeling[nPix]),
				                           m_vert_cost(x,y,beta_label,m_labeling[nPix]));

		}
	}
}

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

void GCoptimization::set_up_swap_energy_G_ARRAY_VW(int size,LabelType alpha_label,LabelType beta_label,
												   PixelType *pixels,Energy* e, Energy::Var *variables)
{
	PixelType nPix,pix,i,x,y;
	EnergyTermType weight;	



	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		y = pix/m_width;
		x = pix - y*m_width;

		if ( x > 0 )
		{
			nPix = pix - 1;
			weight = m_horizWeights[nPix];
	
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothcost(alpha_label,alpha_label)*weight,
							  m_smoothcost(alpha_label,beta_label)*weight,
							  m_smoothcost(beta_label,alpha_label)*weight,
							  m_smoothcost(beta_label,beta_label)*weight );
	
				else
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,
 					                       m_smoothcost(beta_label,m_labeling[nPix])*weight);
	
		}	
		if ( y > 0 )
		{
			nPix = pix - m_width;
			weight = m_vertWeights[nPix];

			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothcost(alpha_label,alpha_label)*weight,
							  m_smoothcost(alpha_label,beta_label)*weight,
							  m_smoothcost(beta_label,alpha_label)*weight,
							  m_smoothcost(beta_label,beta_label)*weight );
	
				else
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,
 					                       m_smoothcost(beta_label,m_labeling[nPix])*weight);
		}	

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_horizWeights[pix],
 					                           m_smoothcost(beta_label,m_labeling[nPix])*m_horizWeights[pix]);
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
				e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_vertWeights[pix],
				                           m_smoothcost(beta_label,m_labeling[nPix])*m_vertWeights[pix]);

		}
	}
}

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

void GCoptimization::set_up_swap_energy_G_ARRAY(int size,LabelType alpha_label,LabelType beta_label,
											   PixelType *pixels,Energy* e, Energy::Var *variables)

{
	PixelType nPix,pix,i,x,y;
	


	for ( i = 0; i < size; i++ )
	{
		pix = pixels[i];
		y = pix/m_width;
		x = pix - y*m_width;

		if ( x > 0 )
		{
			nPix = pix - 1;
	
	
			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothcost(alpha_label,alpha_label),
							  m_smoothcost(alpha_label,beta_label),
							  m_smoothcost(beta_label,alpha_label),
							  m_smoothcost(beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),
 					                       m_smoothcost(beta_label,m_labeling[nPix]));
	
		}	
		if ( y > 0 )
		{
			nPix = pix - m_width;

			if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
				              m_smoothcost(alpha_label,alpha_label),
							  m_smoothcost(alpha_label,beta_label),
							  m_smoothcost(beta_label,alpha_label),
							  m_smoothcost(beta_label,beta_label) );
	
				else
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),
 					                       m_smoothcost(beta_label,m_labeling[nPix]));
		}	

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )
					e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),
 					                           m_smoothcost(beta_label,m_labeling[nPix]));
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;
	
			if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label))
				e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),
				                           m_smoothcost(beta_label,m_labeling[nPix]));

		}
	}
}


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

void GCoptimization::setLabelOrder(bool RANDOM_LABEL_ORDER)
{
	m_random_label_order = RANDOM_LABEL_ORDER;
}

/****************************************************************************/
/* This procedure checks if an error has occured, terminates program if yes */

void GCoptimization::terminateOnError(bool error_condition,const char *message)

{ 
   if  (error_condition) 
   {
       mexErrMsgIdAndTxt("GraphCut:internal_error","\nGCoptimization error: %s\n", message);
    }
}