imagesegmentation.cpp

图像分割算法

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.


#include "./ImageSegmentation.h"
#include <queue>
#include <vector>





namespace ImageProcessing
{


// Gaussian Filter
	Matrix<float> Gaussian2D(int side, float sigma_x, float angle, float ratio)
	{
		Matrix<float> temp(2*side+1, 2*side+1); 

		float sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = 1.0/(2.0*PI*sigma_x*sigma_y)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = (float)(sum_G/T);
			}
		}
		return temp;
	}



	Matrix<double> Gaussian2D(int side, double sigma_x, double angle, double ratio)
	{
		Matrix<double> temp(2*side+1, 2*side+1); 

		double sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = 1.0/(2.0*PI*sigma_x*sigma_y)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = sum_G/T;
			}
		}
		return temp;
	}





// First Directional Derivative of Gaussian Filter


	Matrix<float> FDGaussian2D(int side, float sigma_x, float angle, float ratio)
	{
		Matrix<float> temp(2*side+1, 2*side+1);

		float sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = -x_new/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = (float)(sum_G/T);
			}
		}
		return temp;
	}



	Matrix<double> FDGaussian2D(int side, double sigma_x, double angle, double ratio)
	{
		Matrix<double> temp(2*side+1, 2*side+1);

		double sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = -x_new/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = sum_G/T;
			}
		}
		return temp;
	}


// Second Directional Derivative of Gaussian Filter


	Matrix<float> SDGaussian2D(int side, float sigma_x, float angle, float ratio)
	{
		Matrix<float> temp(2*side+1, 2*side+1);

		float sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);
		}
		else
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = (x_new*x_new-sigma_x*sigma_x)/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = (float)(sum_G/T);
			}
		}
		return temp;
	}



	Matrix<double> SDGaussian2D(int side, double sigma_x, double angle, double ratio)
	{
		Matrix<double> temp(2*side+1, 2*side+1);

		double sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = (x_new*x_new-sigma_x*sigma_x)/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = sum_G/T;
			}
		}
		return temp;
	}





// Laplacian of Gaussian Filter
	Matrix<float> LOG(int side, float sigma_x, float angle, float ratio)
	{
		Matrix<float> temp(2*side+1, 2*side+1); 

		float sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = (x_new*x_new+y_new*y_new-sigma_x*sigma_x-sigma_y*sigma_y)
							/(2.0*PI*sigma_x*sigma_y*(sigma_x*sigma_x*sigma_x*sigma_x+sigma_y*sigma_y*sigma_y*sigma_y))
							*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = (float)(sum_G/T);
			}
		}
		return temp;
	}



	Matrix<double> LOG(int side, double sigma_x, double angle, double ratio)
	{
		Matrix<double> temp(2*side+1, 2*side+1); 

		double sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}

		double T = (double) (2*sample+1)*(2*sample+1);

		for (int i=0;i<2*side+1;i++) 
		{
			for (int j=0;j<2*side+1;j++) 
			{
				double sum_G = 0.0;
				double y = (double) (j-side)-d*sample-d;
				for (int sx=0;sx<2*sample+1;sx++) 
				{
					y += d;
					double x = (double) (i-side)-d*sample-d;
					for (int sy=0;sy<2*sample+1;sy++) 
					{
						x += d;

						double x_new = x*cos_angle + y*sin_angle;
						double y_new = -x*sin_angle + y*cos_angle;
						double g = (x_new*x_new+y_new*y_new-sigma_x*sigma_x-sigma_y*sigma_y)
							/(2.0*PI*sigma_x*sigma_y*(sigma_x*sigma_x*sigma_x*sigma_x+sigma_y*sigma_y*sigma_y*sigma_y))
							*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
						sum_G += g;
					}
				}

				temp.ElemNC(j,i) = sum_G/T;
			}
		}
		return temp;
	}





// Difference of offset Gaussians filter
	Matrix<float> DOOG2D(int side, float sigma_x, float offset, float angle, float ratio)
	{
		Matrix<float> temp(2*side+1, 2*side+1);

		float sigma_y = sigma_x*ratio;

		double sin_angle = sin(angle*PI/180);
		double cos_angle = cos(angle*PI/180);

		int sample;
		double d;
		if (sigma_x < 2.0 || sigma_y < 2.0) 
		{
			// use denser sample for better filter quality 
			sample = 5;
			d = 1.0/(2.0*sample+1.0);      
		}
		else 
		{
			sample = 0;
			d = 1.0;
		}