imagesegmentation.cpp

图像分割算法

		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/(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)
							- 1/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new-offset)*(x_new-offset)/(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> DOOG2D(int side, double sigma_x, double offset, 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/(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)
							- 1/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new-offset)*(x_new-offset)/(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 (centered at the middle of two Gaussians)
	Matrix<float> DOOG2DCentered(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);
		double half_offset = offset/2.0;

		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/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new+half_offset)*(x_new+half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
							- 1/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new-half_offset)*(x_new-half_offset)/(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> DOOG2DCentered(int side, double sigma_x, double offset, 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);
		double half_offset = offset/2.0;

		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/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new+half_offset)*(x_new+half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
							- 1/(2.0*PI*sigma_x*sigma_y)
							*exp(-((x_new-half_offset)*(x_new-half_offset)/(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;	
	}








// Filter image with these filters

	Matrix<float> FilterGaussian2D(Matrix<float>& image, float sigma_x, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<float> filt = Gaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<double> FilterGaussian2D(Matrix<double>& image, double sigma_x, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<double> filt = Gaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}


	// First derivative of Gaussian
	Matrix<float> FilterFDGaussian2D(Matrix<float>& image, float sigma_x, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<float> filt = FDGaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<double> FilterFDGaussian2D(Matrix<double>& image, double sigma_x, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<double> filt = FDGaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}


	// Second derivative of Gaussian
	Matrix<float> FilterSDGaussian2D(Matrix<float>& image, float sigma_x, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<float> filt = SDGaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<double> FilterSDGaussian2D(Matrix<double>& image, double sigma_x, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<double> filt = SDGaussian2D(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}


	// Laplacian of Gaussian
	Matrix<float> FilterLOG(Matrix<float>& image, float sigma_x, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<float> filt = LOG(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}
	Matrix<double> FilterLOG(Matrix<double>& image, double sigma_x, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
		int side = (int)(3*sigma_max+0.5);
		Matrix<double> filt = LOG(side, sigma_x, angle, ratio);
		return Conv2(image, filt);
	}

	
	// Difference of offset Gaussians
	Matrix<float> FilterDOOG2D(Matrix<float>& image, float sigma_x, float offset, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		int side = (int)((3*sigma_x+offset>=3*sigma_y?3*sigma_x+offset:3*sigma_y)+0.5);
		Matrix<float> filt = DOOG2D(side, sigma_x, offset, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<double> FilterDOOG2D(Matrix<double>& image, double sigma_x, double offset, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		int side = (int)((3*sigma_x+offset>=3*sigma_y?3*sigma_x+offset:3*sigma_y)+0.5);
		Matrix<double> filt = DOOG2D(side, sigma_x, offset, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<float> FilterDOOG2DCentered(Matrix<float>& image, float sigma_x, float offset, float angle, float ratio)
	{
		float sigma_y = ratio*sigma_x;
		int side = (int)((3*sigma_x+offset/2>=3*sigma_y?3*sigma_x+offset/2:3*sigma_y)+0.5);
		Matrix<float> filt = DOOG2DCentered(side, sigma_x, offset, angle, ratio);
		return Conv2(image, filt);
	}

	Matrix<double> FilterDOOG2DCentered(Matrix<double>& image, double sigma_x, double offset, double angle, double ratio)
	{
		double sigma_y = ratio*sigma_x;
		int side = (int)((3*sigma_x+offset/2>=3*sigma_y?3*sigma_x+offset/2:3*sigma_y)+0.5);
		Matrix<double> filt = DOOG2DCentered(side, sigma_x, offset, angle, ratio);
		return Conv2(image, filt);
	}


// Various Edge Detection schemes
	Matrix<double> NonMaximaSuppress(Matrix<double>& edgesMain, Matrix<double>& theta)
	{
		Matrix<double> vectorX = Cos(theta);
		Matrix<double> vectorY = Sin(theta);
		return NonMaximaSuppress(edgesMain, vectorX, vectorY);
	}

	Matrix<float> NonMaximaSuppress(Matrix<float>& edgesMain, Matrix<float>& theta)
	{
		Matrix<float> vectorX = Cos(theta);
		Matrix<float> vectorY = Sin(theta);
		return NonMaximaSuppress(edgesMain, vectorX, vectorY);
	}


	Matrix<double> NonMaximaSuppress(Matrix<double>& edgesMain, Matrix<double>& vectorX, Matrix<double>& vectorY)
	{
		Matrix<double> edges = edgesMain.Clone();
		Matrix<int> mask = NonMaximaMask(edges - Minimum(edges(":")), vectorX, vectorY);

		for(int i=0;i<edges.Columns();i++)
		{
			for(int j=0;j<edges.Rows();j++)
			{
				if(mask(j,i) == 0)
				{
					edges(j,i) = 0;
				}
			}
		}

		return edges;
	}

	Matrix<float> NonMaximaSuppress(Matrix<float>& edgesMain, Matrix<float>& vectorX, Matrix<float>& vectorY)
	{
		Matrix<float> edges = edgesMain.Clone();
		Matrix<int> mask = NonMaximaMask(edges - Minimum(edges(":")), vectorX, vectorY);

		for(int i=0;i<edges.Columns();i++)
		{
			for(int j=0;j<edges.Rows();j++)
			{
				if(mask(j,i) == 0)
				{
					edges(j,i) = 0;
				}
			}
		}

		return edges;
	}


	Matrix<int> NonMaximaMask(Matrix<double>& edges, Matrix<double>& vectorX, Matrix<double>& vectorY)
	{
		Matrix<int> mask(edges.Rows(), edges.Columns(), 1);
		Matrix<double> theta(edges.Rows(), edges.Columns(), 0);
		
		Matrix<int> mask15(edges.Rows(), edges.Columns(), 0);
		Matrix<int> mask26(edges.Rows(), edges.Columns(), 0);
		Matrix<int> mask37(edges.Rows(), edges.Columns(), 0);
		Matrix<int> mask48(edges.Rows(), edges.Columns(), 0);

		for(int i=0;i<theta.Columns();i++)
		{
			for(int j=0;j<theta.Rows();j++)
			{
				theta(j,i) = Direction(vectorY(j,i), vectorX(j,i));

				if( theta(j,i) >= 0 && theta(j,i) < PI/4 )
				{
					mask15(j,i) = 1;
				}

				else if( theta(j,i) >= PI/4 && theta(j,i) < PI/2 )
				{
					mask26(j,i) = 1;
				}

				else if( theta(j,i) >= PI/2 && theta(j,i) < PI*3/4 )
				{
					mask37(j,i) = 1;
				}

				else if( theta(j,i) >= PI*3/4 && theta(j,i) < PI )
				{
					mask48(j,i) = 1;
				}

			}

		}



		//Case 1
		for(int i=0;i<theta.Columns()-1;i++)
		{
			for(int j=0;j<theta.Rows()-1;j++)
			{
				if(mask15(j,i) == 1)
				{
					double d = tan(theta(j,i));
					double interpolated = edges(j,i+1)*(1-d) + edges(j+1,i+1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}

		//Case 5
		for(int i=1;i<theta.Columns();i++)
		{
			for(int j=1;j<theta.Rows();j++)
			{
				if(mask15(j,i) == 1)
				{
					double d = tan(theta(j,i));
					double interpolated = edges(j,i-1)*(1-d) + edges(j-1,i-1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}


		//Case 2
		for(int i=0;i<theta.Columns()-1;i++)
		{
			for(int j=0;j<theta.Rows()-1;j++)
			{
				if(mask26(j,i) == 1)
				{
					double d = tan(PI/2 - theta(j,i));
					double interpolated = edges(j+1,i)*(1-d) + edges(j+1,i+1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}

		//Case 6
		for(int i=1;i<theta.Columns();i++)
		{
			for(int j=1;j<theta.Rows();j++)
			{
				if(mask26(j,i) == 1)
				{
					double d = tan(PI/2 - theta(j,i));
					double interpolated = edges(j-1,i)*(1-d) + edges(j-1,i-1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}




		//Case 3
		for(int i=1;i<theta.Columns();i++)
		{
			for(int j=0;j<theta.Rows()-1;j++)
			{
				if(mask37(j,i) == 1)
				{
					double d = tan(theta(j,i) - PI/2);
					double interpolated = edges(j+1,i)*(1-d) + edges(j+1,i-1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}


		//Case 7
		for(int i=0;i<theta.Columns()-1;i++)
		{
			for(int j=1;j<theta.Rows();j++)
			{
				if(mask37(j,i) == 1)
				{
					double d = tan(theta(j,i) - PI/2);
					double interpolated = edges(j-1,i)*(1-d) + edges(j-1,i+1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}
			}
		}


		//Case 4
		for(int i=1;i<theta.Columns();i++)
		{
			for(int j=0;j<theta.Rows()-1;j++)
			{
				if(mask48(j,i) == 1)
				{
					double d = tan(PI - theta(j,i));
					double interpolated = edges(j,i-1)*(1-d) + edges(j+1,i-1)*d;
					if(edges(j,i)<interpolated)
					{
						mask(j,i) = 0;
					}
				}