imagesegmentation.cpp

图像分割算法

			}
		}


		//Case 8
		for(int i=0;i<theta.Columns()-1;i++)
		{
			for(int j=1;j<theta.Rows();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;
					}
				}
			}
		}

		return mask;

	}


	Matrix<int> NonMaximaMask(Matrix<float>& edges, Matrix<float>& vectorX, Matrix<float>& 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;
					}
				}
			}
		}


		//Case 8
		for(int i=0;i<theta.Columns()-1;i++)
		{
			for(int j=1;j<theta.Rows();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;
					}
				}
			}
		}

		return mask;

	}






	double Direction(double y, double x)
	{
		if(x>=0 && y>=0)
		{
			return atan2(y,x);
		}
		else if(x<0 && y == 0)
		{
			return atan2(y,x) - PI;
		}
		else if(x<0 && y>0)
		{
			return atan2(y,x);
		}
		else if(x<=0 && y<0)
		{
			return PI + atan2(y,x) ;
		}
		else if(x>0 && y<0)
		{
			return PI + atan2(y,x) ;
		}
		else 
		{
			cout << "Something went wrong with Direction(x,y) function.. Returning -1..\n" << endl;
			return -1;
		}
	}

	float Direction(float y, float x)
	{
		if(x>=0 && y>=0)
		{
			return atan2(y,x);
		}
		else if(x<0 && y == 0)
		{
			return atan2(y,x) - PI;
		}
		else if(x<0 && y>0)
		{
			return atan2(y,x);
		}
		else if(x<=0 && y<0)
		{
			return PI + atan2(y,x) ;
		}
		else if(x>0 && y<0)
		{
			return PI + atan2(y,x) ;
		}
		else 
		{
			cout << "Something went wrong with Direction(x,y) function.. Returning -1..\n" << endl;
			return -1;
		}
	}



// Edgeflow

// both grayscale and multi-valued
	// Outout is two dimensional (x and y components of the vector field)
	MatrixList<float> EdgeflowVectorField(Matrix<float>& image, int angles, float sigma, float offset, float ratio, bool normalized)
	{

		Vector<float> radAngles(2*angles);
		Vector<float> cosAngles(2*angles);
		Vector<float> sinAngles(2*angles);
		for(int i=0; i<angles; i++)
		{
			radAngles[i] = (float)(i*PI/angles);
			cosAngles[i] = cos(radAngles[i]);
			sinAngles[i] = sin(radAngles[i]);

			radAngles[i+angles] = (float)((i+angles)*PI/angles);
			cosAngles[i+angles] = cos(radAngles[i+angles]);
			sinAngles[i+angles] = sin(radAngles[i+angles]);

		}

		// Find energies at directions theta and theta+pi
		Vector<Matrix<float> > energies(2*angles);
		float sigma_y = ratio*sigma;
		int side = (int)((3*sigma+offset>=3*sigma_y?3*sigma+offset:3*sigma_y)+0.5);
		for(int i=0; i<angles; i++)
		{
			float angle = (float)(i*(180.0/angles));
			
			//char s[200];
			Matrix<float> filt = DOOG2D(side, sigma, offset, 180+angle, ratio);
			energies(i) = AbsI(Conv2(image, filt));
			//sprintf(s, "%02d.0.bmp", i);
			//ImWrite(energies(i),s);
			filt = DOOG2D(side, sigma, offset, angle, ratio);
			energies(i+angles) = AbsI(Conv2(image, filt));
			//sprintf(s, "%02d.1.bmp", i);
			//ImWrite(energies(i+angles),s);

			//Matrix<float> filt = DOOG2DCentered(side, sigma, offset, angle, ratio);
			//char s[200];
			//Matrix<float> filtOut = Conv2(image, filt, Full);
			//
			//int offsetX = (int)(offset*cos(radAngles[i])/2+0.5);
			//offsetX = offsetX==0?1:offsetX;
			//int offsetY = (int)(offset*sin(radAngles[i])/2+0.5);
			//offsetY = offsetY==0?1:offsetY;
			//// cout << (fmod(2.0+cos((double)radAngles[i]),2.0)-1) << " : " << (fmod(2.0+sin((double)radAngles[i]),2.0)-1) << endl;

			//energies(i) = Abs(filtOut.Slice(side+offsetY, image.Rows()+side+offsetY-1, side+offsetX, image.Columns()+side+offsetX-1));
			//sprintf(s, "%02d.0.bmp", i);
			//ImWrite(energies(i), s);
			//
			//energies(i+angles) = Abs(filtOut.Slice(side-offsetY, image.Rows()+side-offsetY-1, side-offsetX, image.Columns()+side-offsetX-1));
			//sprintf(s, "%02d.1.bmp", i);
			//ImWrite(energies(i+angles), s);
		}
		//pf.Toc();

		// Sum energies over half circles
		//pf.Tic();
		Vector<Matrix<float> > sumEnergies(2*angles);
		for(int i=0; i<2*angles; i++)
		{
			sumEnergies(i) = Matrix<float>(image.Rows(), image.Columns(), 0.0f);
			int start = i-angles/2;
			int end = start+angles;
			for(int j=start; j<end; j++)
			{
				int ind = (j+2*angles)%(2*angles);
				sumEnergies(i) += energies(ind);
			}
		}
		//pf.Toc();

		// normalize summed energies to make them like probabilities.
		//pf.Tic();
		Vector<Matrix<float> > probabilities(2*angles);
		for(int i=0; i<angles; i++)
		{
			probabilities(i) = Matrix<float>(image.Rows(), image.Columns());
			probabilities(i+angles) = Matrix<float>(image.Rows(), image.Columns());
			Matrix<float> total = sumEnergies(i) + sumEnergies(i+angles);
			for(int r=0;r<image.Rows();r++)
			{
				for(int c=0;c<image.Columns();c++)
				{
					if(total(r,c) > 1e-9)
					{
						probabilities(i).Elem(r,c) = sumEnergies(i).Elem(r,c)/total(r,c);
						probabilities(i+angles).Elem(r,c) = sumEnergies(i+angles).Elem(r,c)/total(r,c);
					}
					else
					{
						probabilities(i).Elem(r,c) = 0.5;
						probabilities(i+angles).Elem(r,c) = 0.5;
					}
				}
			}
		}
		//pf.Toc();

		//pf.Tic();
		Matrix<float> xFlow(image.Rows(), image.Columns());
		Matrix<float> yFlow(image.Rows(), image.Columns());
		Matrix<float> maxEnergy(image.Rows(), image.Columns());
		for(int r=0;r<image.Rows();r++)
		{
			for(int c=0;c<image.Columns();c++)
			{

				float dirX = 0;
				float dirY = 0;
				//int count = 0;
				for(int k=0;k<2*angles;k++)
				{
					float v = probabilities(k).Elem(r,c);
					float vl = probabilities((k+2*angles-1)%(2*angles)).Elem(r,c);
					float vr = probabilities((k+2*angles+1)%(2*angles)).Elem(r,c);
					if(v>=vl && v>=vr)
					{
						float orientation, strength;
						if(vl+vr != 2*v)
						{
							orientation = 0.5f*(vl - vr)/(vl + vr - 2*v);
							strength = v + 0.5f*( (vr - vl)*orientation 
								+ (vr + vl - 2*v)*orientation*orientation );
							orientation = (float)((k+orientation)*PI/angles);
						}
						else
						{
							strength = v;
							orientation = (float)(k*PI/angles);
						}
						dirX += cos(orientation)*strength;
						dirY += sin(orientation)*strength;
						//maxEnergy(r,c) += strength;
						//count++;
					}
				}

				//maxEnergy(r,c) /= (float)count;

				float dir = sqrt(dirX*dirX+dirY*dirY);
				dirX /= dir;
				dirY /= dir;


				float value;
				int ang;
				if(probabilities(0).Elem(r,c)>probabilities(angles).Elem(r,c))
				{
					value = probabilities(0).Elem(r,c);
					ang = 0;
				}
				else
				{
					value = probabilities(angles).Elem(r,c);
					ang = angles;
				}


				float maximum = value;
				int maxIndex = ang;
				float minimum = value;
				int minIndex = ang;
				int maxOrientation = 0;
				int minOrientation = 0;

				for(int k=1;k<angles;k++)
				{
					if(probabilities(k).Elem(r,c)>probabilities(k+angles).Elem(r,c))
					{
						value = probabilities(k).Elem(r,c);
						ang = k;
					}
					else
					{
						value = probabilities(k+angles).Elem(r,c);
						ang = k+angles;
					}

					if(value > maximum)
					{
						maximum = value;
						maxIndex = ang;
						maxOrientation = k;
					}
					if(value < minimum)
					{
						minimum = value;
						minIndex = ang;
						minOrientation = k;
					}
				}
				if(normalized)
				{
					xFlow(r,c) = dirX*probabilities(maxIndex)(r,c);
					yFlow(r,c) = dirY*probabilities(maxIndex)(r,c);
				}
				else
				{
					xFlow(r,c) = dirX*sumEnergies(maxIndex)(r,c);
					yFlow(r,c) = dirY*sumEnergies(maxIndex)(r,c);
				}
				maxEnergy(r,c) = sumEnergies(maxIndex)(r,c);
			}
		}
		//if(normalized)
		//{
		//	xFlow *= ToFloat(maxEnergy > (float)angles);
		//	yFlow *= ToFloat(maxEnergy > (float)angles);
		//}

		//ImWrite(xFlow, "xflow.bmp");
		//ImWrite(yFlow, "yflow.bmp");
		//ImWrite(maxEnergy, "maxEnergy.bmp");
		return MatrixList<float>(xFlow, yFlow);

	}

	//MatrixList<double> EdgeflowVectorField(Matrix<double> image, double sigma, double offset, double ratio, int angles)
	//{
	//}

	MatrixList<float> EdgeflowVectorField(MatrixList<float>& image, int angles, float sigma, float offset, float ratio, bool normalized)
	{

		Vector<float> radAngles(2*angles);
		Vector<float> cosAngles(2*angles);
		Vector<float> sinAngles(2*angles);
		for(int i=0; i<angles; i++)
		{
			radAngles[i] = (float)(i*PI/angles);
			cosAngles[i] = cos(radAngles[i]);
			sinAngles[i] = sin(radAngles[i]);

			radAngles[i+angles] = (float)((i+angles)*PI/angles);
			cosAngles[i+angles] = cos(radAngles[i+angles]);
			sinAngles[i+angles] = sin(radAngles[i+angles]);

		}

		// Find energies at directions theta and theta+pi
		Vector<Matrix<float> > energies(2*angles);
		for(int i=0;i<energies.Length();i++)
		{
			energies(i) = Matrix<float>(image.Rows(), image.Columns(), 0);
		}
		float sigma_y = ratio*sigma;
		int side = (int)((3*sigma+offset>=3*sigma_y?3*sigma+offset:3*sigma_y)+0.5);
		for(int i=0; i<angles; i++)
		{
			float angle = (float)(i*(180.0/angles));
			for(int c=0; c<image.Planes();c++)
			{
				//char s[200];
				Matrix<float> filt = DOOG2D(side, sigma, offset, 180+angle, ratio);
				energies(i) += AbsI(Conv2(image[c], filt));
				//sprintf(s, "%02d.0.bmp", i);