gcoptimization.cpp

Graph Cut algorithm implementation. Includes MATLAB compiled codes.

	for ( i = size - 1; i >= 0; i-- )
	{
		pix = m_lookupPixVar[i];
		y = pix/m_width;
		x = pix - y*m_width;

		m_lookupPixVar[pix] = i;

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;

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

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;

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

		if ( y > 0 )
		{
			nPix = pix - m_width;
	
			if ( m_labeling[nPix] == alpha_label )
			   e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label),
			   	                 m_smoothcost(m_labeling[pix],alpha_label));
		}	
			
	}
	
}

/**********************************************************************************************/
/* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */
/* specified by a function */
void GCoptimization::set_up_expansion_energy_G_ARRAY_VW_pix(int size, LabelType alpha_label,Energy *e,
													  Energy::Var *variables, PixelType *pixels, int num )
{
	int i,nPix,pix,x,y;
	EnergyTermType weight;


	for ( i = 0; i < num; i++ )
	{
		pix = pixels[i];
		if ( m_labeling[pix]!= alpha_label )
		{
			y = pix/m_width;
			x = pix - y*m_width;

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

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

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

/**********************************************************************************************/
/* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */
/* specified by a function */
void GCoptimization::set_up_expansion_energy_G_ARRAY_pix(int size, LabelType alpha_label,Energy *e,
													 Energy::Var *variables, PixelType *pixels,
													 int num)
{
	int i,nPix,pix,x,y;


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


		if ( m_labeling[pix]!= alpha_label )
		{
			if ( x < m_width - 1 )
			{
				nPix = pix + 1;

				if ( m_labeling[nPix] != alpha_label && m_lookupPixVar[pix] != -1)
					e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
					          m_smoothcost(alpha_label,alpha_label),
							  m_smoothcost(alpha_label,m_labeling[nPix]),
							  m_smoothcost(m_labeling[pix],alpha_label),
							  m_smoothcost(m_labeling[pix],m_labeling[nPix]));
				else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),
				                 m_smoothcost(m_labeling[pix],m_labeling[nPix]));
			}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;

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

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


/**********************************************************************************************/
/* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */
/* specified by a function */
void GCoptimization::set_up_expansion_energy_G_FnPix(int size, LabelType alpha_label,Energy *e,
													 Energy::Var *variables )
{
	int i,nPix,pix,x,y;


	for ( i = size - 1; i >= 0; i-- )
	{
		pix = m_lookupPixVar[i];
		y = pix/m_width;
		x = pix - y*m_width;

		m_lookupPixVar[pix] = i;

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;

			if ( m_labeling[nPix] != alpha_label )
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
					          m_smoothFnPix(pix,nPix,alpha_label,alpha_label),
							  m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
							  m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label),
							  m_smoothFnPix(pix,nPix,m_labeling[pix],m_labeling[nPix]));
			else   e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),
				                 m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;

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

		if ( y > 0 )
		{
			nPix = pix - m_width;
	
			if ( m_labeling[nPix] == alpha_label )
			   e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,alpha_label),
			   	                 m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));
		}	
			
	}
	
}

/**********************************************************************************************/
/* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */
/* specified by a function */
void GCoptimization::set_up_expansion_energy_G_FnCoord(int size, LabelType alpha_label,Energy *e,
													 Energy::Var *variables )
{
	int i,nPix,pix,x,y;


	for ( i = size - 1; i >= 0; i-- )
	{
		pix = m_lookupPixVar[i];
		y = pix/m_width;
		x = pix - y*m_width;

		m_lookupPixVar[pix] = i;

		if ( x < m_width - 1 )
		{
			nPix = pix + 1;

			if ( m_labeling[nPix] != alpha_label )
				e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],
					          m_horz_cost(x,y,alpha_label,alpha_label),
							  m_horz_cost(x,y,alpha_label,m_labeling[nPix]),
							  m_horz_cost(x,y,m_labeling[pix],alpha_label),
							  m_horz_cost(x,y,m_labeling[pix],m_labeling[nPix]));
			else   e ->add_term1(variables[i],m_horz_cost(x,y,alpha_label,alpha_label),
				                 m_horz_cost(x,y,m_labeling[pix],alpha_label));
		}	

		if ( y < m_height - 1 )
		{
			nPix = pix + m_width;

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

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



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

GCoptimization::EnergyType GCoptimization::oneExpansionIteration()
{
	int next;
   
	terminateOnError( m_dataType == NONE,"You have to set up the data cost before running optimization");
	terminateOnError( m_smoothType == NONE,"You have to set up the smoothness cost before running optimization");

	if (m_random_label_order) scramble_label_table();
	

	for (next = 0;  next < m_num_labels;  next++ )
		perform_alpha_expansion(m_labelTable[next]);
	
	
	return(compute_energy());
}

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

void GCoptimization::setNeighbors(PixelType pixel1, int pixel2, EnergyTermType weight)
{

	assert(pixel1 < m_num_pixels && pixel1 >= 0 && pixel2 < m_num_pixels && pixel2 >= 0);
	assert(m_grid_graph == 0);

	Neighbor *temp1 = (Neighbor *) new Neighbor;
	Neighbor *temp2 = (Neighbor *) new Neighbor;

	temp1->weight  = weight;
	temp1->to_node = pixel2;

	temp2->weight  = weight;
	temp2->to_node = pixel1;

	m_neighbors[pixel1].addFront(temp1);
	m_neighbors[pixel2].addFront(temp2);
	
}

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

void GCoptimization::setNeighbors(PixelType pixel1, int pixel2)
{

	assert(pixel1 < m_num_pixels && pixel1 >= 0 && pixel2 < m_num_pixels && pixel2 >= 0);
	assert(m_grid_graph == 0);
	

	Neighbor *temp1 = (Neighbor *) new Neighbor;
	Neighbor *temp2 = (Neighbor *) new Neighbor;

	temp1->weight  = (EnergyTermType) 1;
	temp1->to_node = pixel2;

	temp2->weight  = (EnergyTermType) 1;
	temp2->to_node = pixel1;

	m_neighbors[pixel1].addFront(temp1);
	m_neighbors[pixel2].addFront(temp2);
	
}

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

GCoptimization::~GCoptimization()
{

    // <!-- bagon
    // mexWarnMsgTxt("Calling destructor"); /* bagon added */
    if ( m_dataType == ARRAY )
        delete [] m_datacost;
    if ( m_smoothType == ARRAY )
        delete [] m_smoothcost;
    if ( m_varying_weights == 1 ) {
        delete [] m_vertWeights;
        delete [] m_horizWeights;
    }
    // bagon -->
    
	if ( m_deleteLabeling ) 
		delete [] m_labeling;


	if ( m_dataInput == SET_INDIVIDUALLY )
		delete [] m_datacost;

	if ( m_smoothInput == SET_INDIVIDUALLY )
			delete [] m_smoothcost;
 		
	if ( ! m_grid_graph )
		delete [] m_neighbors;			


	delete [] m_labelTable;
	delete [] m_lookupPixVar;
			
}


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

GCoptimization::EnergyType GCoptimization::swap(int max_num_iterations)
{
	return(start_swap(max_num_iterations)); 
}

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


GCoptimization::EnergyType GCoptimization::swap()
{
	return(start_swap(MAX_INTT));
}

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


GCoptimization::EnergyType GCoptimization::start_swap(int max_num_iterations )
{
	
	int curr_cycle = 1;
	EnergyType new_energy,old_energy;
	

	new_energy = compute_energy();

	old_energy = (new_energy+1)*10; // BAGON changed init value to exceed current energy by factor of 10 (thanks to A. Khan)    

	while ( old_energy > new_energy  && curr_cycle <= max_num_iterations)
	{
		old_energy = new_energy;
		new_energy = oneSwapIteration();
		
		curr_cycle++;	
	}

	return(new_energy);
}

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

GCoptimization::EnergyType GCoptimization::oneSwapIteration()
{
	int next,next1;
   
	if (m_random_label_order) scramble_label_table();
		

	for (next = 0;  next < m_num_labels;  next++ )
		for (next1 = m_num_labels - 1;  next1 >= 0;  next1-- )
			if ( m_labelTable[next] < m_labelTable[next1] )
				perform_alpha_beta_swap(m_labelTable[next],m_labelTable[next1]);

	return(compute_energy());
}

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

GCoptimization::EnergyType GCoptimization::alpha_beta_swap(LabelType alpha_label, LabelType beta_label)
{