imagesegmentation.cpp

图像分割算法

				//ImWrite(energies(i),s);
				filt = DOOG2D(side, sigma, offset, angle, ratio);
				energies(i+angles) += AbsI(Conv2(image[c], 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);

	}



	void block_write(Matrix<int>& matrix, int x_ind, int y_ind, int *keys, int energy);
	double my_atan2(double y, double x);

	Matrix<int> CreateFlowImage(Matrix<float>& xFlow, Matrix<float>& yFlow)
	{
		if(xFlow.Rows() != yFlow.Rows() || xFlow.Columns() != yFlow.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Matrix dimensions does not match to each other!");
		}

		int display_map[8][15] ={37,38,39,40,41,42,43,33,23,13,12,51,59,67,66,
								10,20,30,40,50,60,70,61,52,43,34,69,68,67,66,
								13,22,31,40,49,58,67,59,51,43,34,57,47,37,28,
								16,24,32,40,48,56,64,55,46,37,28,65,66,67,68,
								43,42,41,40,39,38,37,47,57,67,68,29,21,13,14,
								70,60,50,40,30,20,10,19,28,37,46,11,12,13,14,
								67,58,49,40,31,22,13,21,29,37,46,23,33,43,52,
								64,56,48,40,32,24,16,25,34,43,52,15,14,13,12};
								
		int quant[17] = {0, 45, 45, 90, 90, 135, 135, 180, 180, 225, 225, 270, 270, 315, 315, 0, 0};

		Matrix<int> angles(xFlow.Rows(), xFlow.Columns());
		Matrix<float> raw_angles(xFlow.Rows(), xFlow.Columns());
		Matrix<int> flow(9*xFlow.Rows(), 9*xFlow.Columns(), 255);

		Matrix<float> energies = SqrtI(xFlow*xFlow + yFlow*yFlow);
		float maximum = Maximum(energies(":"));
		float minimum = Minimum(energies(":"));
		Matrix<int> scaled_energies(xFlow.Rows(), xFlow.Columns(), 0);
		if(maximum-minimum > 1e-9)
		{
			scaled_energies = 255 - ToInt((energies-minimum)*255.0f/(maximum-minimum));
			//scaled_energies = 255 - ToInt(energies*255.0f/maximum);
		}

		int disp_keys[15];
		for(int i=0;i<xFlow.Columns();i++){
			for(int j=0;j<xFlow.Rows();j++){
				float angle = (float)my_atan2( yFlow(j,i), xFlow(j,i) ) ;
				raw_angles(j,i) = angle;
				int q_angle = (int)(angle/22.5);
				q_angle = (q_angle+17)%17;
				angles(j,i) =  quant[q_angle];
				int l = quant[q_angle] / 45;
				for(int s=0;s<15;s++){
					disp_keys[s] = display_map[l][s];
				}
				int q_energy = scaled_energies(j,i);
				block_write(flow, i, j, disp_keys, q_energy);
			}
		}

		return flow;
	}

	void block_write(Matrix<int>& matrix, int x_ind, int y_ind, int *keys, int energy){
		for(int t=0;t<15;t++){
			int x_loc = 9*x_ind + (keys[t] % 9);
			int y_loc = 9*y_ind + (keys[t] / 9);
			matrix(y_loc,x_loc) = energy;
		}
	}

	double my_atan2(double y, double x){
		if(x>=0 && y>=0){
			return atan2(y,x)*180.0/PI;
		}
		else if(x>=0 && y<=0){
			return ( 2*PI + atan2(y,x) )*180.0/PI;
		}
		else if(x<=0 && y>=0){
			return atan2(y,x)*180.0/PI;
		}
		else if(x<=0 && y<=0){
			return (2*PI + atan2(y,x) )*180.0/PI;
		}
		else
		{
			return -1;
		}
	}



	Matrix<int> HystThreshold(Matrix<float>& m, float thHigh, float thLow)
	{
		Matrix<int> th = (m > thHigh);
		Matrix<int> th2 = (m > thLow);

		int count;
		do{
			count = 0;
			for(int i=0;i<m.Rows();i++)
			{
				for(int j=0;j<m.Columns();j++)
				{
					if(th2.ElemNC(i,j) == 0 || th.ElemNC(i,j) == 1)
					{
						continue;
					}

					if(i>0 && th.ElemNC(i-1,j)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(i<m.Rows()-1 && th.ElemNC(i+1,j)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(j>0 && th.ElemNC(i,j-1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(j<m.Columns()-1 && th.ElemNC(i,j+1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}

					if(i>0 && j>0 && th.ElemNC(i-1,j-1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(i<m.Rows()-1 && j>0 && th.ElemNC(i+1,j-1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(i> 0 && j<m.Columns()-1 && th.ElemNC(i-1,j+1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}
					if(i<m.Rows()-1 && j<m.Columns()-1 && th.ElemNC(i+1,j+1)==1)
					{
						th.ElemNC(i,j) = 1;
						count++;
					}

				}
			}
			//cout << count << endl;
		}while(count > 0);
		
		return th;

	}







	Matrix<float>& PMAnisoDiff(Matrix<float>& image, float K, int iterations)
	{
		float lambda = 0.25;
		float K2Inv = 1/(K*K);
		Matrix<float> gradX(image.Rows(), image.Columns());
		Matrix<float> gradY(image.Rows(), image.Columns());
		gradX(0, image.Rows()-1, 0, 0) = 0;
		gradY(0, 0, 0, image.Columns()-1) = 0;

		for(int k=0; k<iterations; k++)
		{
			// Calculate gradient
			for(int i=0;i<image.Rows();i++)
			{
				for(int j=1;j<image.Columns();j++)
				{
					gradX.ElemNC(i,j) = image.ElemNC(i,j-1) - image.ElemNC(i,j);
				}
			}

			for(int i=1;i<image.Rows();i++)
			{
				for(int j=0;j<image.Columns();j++)
				{
					gradY.ElemNC(i,j) = image.ElemNC(i-1,j) - image.ElemNC(i,j);
				}
			}


			// Calculate flux
			Matrix<float> fluxX = gradX*Exp(-K2Inv*Abs(gradX));
			Matrix<float> fluxY = gradY*Exp(-K2Inv*Abs(gradY));


			// Update image
			for(int i=0;i<image.Rows()-1;i++)
			{
				for(int j=0;j<image.Columns()-1;j++)
				{
					image.ElemNC(i,j) += lambda*(fluxX.ElemNC(i,j) - fluxX.ElemNC(i,j+1) 
						+ fluxY.ElemNC(i,j) - fluxY.ElemNC(i+1,j) );
				}
			}


			for(int j=0;j<image.Rows()-1;j++)
			{
				image.ElemNC(j, image.Columns()-1) += lambda*(fluxY.ElemNC(j, image.Columns()-1) 
					- fluxY.ElemNC(j+1, image.Columns()-1) );
			}

			for(int i=0;i<image.Columns()-1;i++)
			{
				image.ElemNC(image.Rows()-1, i) += lambda*(fluxX.ElemNC(image.Rows()-1, i) 
					- fluxX.ElemNC(image.Rows()-1, i+1) );
			}
		}
		return image;
	}

	Matrix<double>& PMAnisoDiff(Matrix<double>& image, double K, int iterations)
	{
		double lambda = 0.25;
		double K2Inv = 1/(K*K);
		Matrix<double> gradX(image.Rows(), image.Columns());
		Matrix<double> gradY(image.Rows(), image.Columns());
		gradX(0, image.Rows()-1, 0, 0) = 0;
		gradY(0, 0, 0, image.Columns()-1) = 0;

		for(int k=0; k<iterations; k++)
		{
			// Calculate gradient
			for(int i=0;i<image.Rows();i++)
			{
				for(int j=1;j<image.Columns();j++)
				{
					gradX.ElemNC(i,j) = image.ElemNC(i,j-1) - image.ElemNC(i,j);
				}
			}

			for(int i=1;i<image.Rows();i++)
			{
				for(int j=0;j<image.Columns();j++)
				{
					gradY.ElemNC(i,j) = image.ElemNC(i-1,j) - image.ElemNC(i,j);
				}
			}


			// Calculate flux
			Matrix<double> fluxX = gradX*Exp(-K2Inv*Abs(gradX));
			Matrix<double> fluxY = gradY*Exp(-K2Inv*Abs(gradY));


			// Update image
			for(int i=0;i<image.Rows()-1;i++)
			{
				for(int j=0;j<image.Columns()-1;j++)
				{
					image.ElemNC(i,j) += lambda*(fluxX.ElemNC(i,j) - fluxX.ElemNC(i,j+1) 
						+ fluxY.ElemNC(i,j) - fluxY.ElemNC(i+1,j) );
				}
			}


			for(int j=0;j<image.Rows()-1;j++)
			{
				image.ElemNC(j, image.Columns()-1) += lambda*(fluxY.ElemNC(j, image.Columns()-1) 
					- fluxY.ElemNC(j+1, image.Columns()-1) );
			}

			for(int i=0;i<image.Columns()-1;i++)
			{
				image.ElemNC(image.Rows()-1, i) += lambda*(fluxX.ElemNC(image.Rows()-1, i) 
					- fluxX.ElemNC(image.Rows()-1, i+1) );
			}
		}
		return image;
	}

	//MatrixList<float>& PMAnisoDiff(MatrixList<float>& image, float K, int iterations)
	//{
	//}

	//MatrixList<double>& PMAnisoDiff(MatrixList<double>& image, double K, int iterations)
	//{
	//}





	class CompareGrads
	{
	public:
		CompareGrads(){}
		~CompareGrads(){}

		static Matrix<float> gradMatrix;
		static int Compare( Point px, Point py )  
		{
			if(gradMatrix(px.y,px.x) > gradMatrix(py.y,py.x))
			{
				return 1;
			}
			else if(gradMatrix(px.y,px.x) == gradMatrix(py.y,py.x))
			{
				return 0;
			}
			else
			{
				return -1;
			}
		}

	};


	Matrix<float> CompareGrads::gradMatrix;

	Matrix<int> GetEgdes(Matrix<float> &grads, Matrix<int> &thick, Matrix<int> &suppressed)
	{
		Matrix<int> edges(grads.Rows(), grads.Columns(), 0);

		Vector<int> LabelCount(grads.Rows()*grads.Columns(), 0);

		Vector<float> LabelAvg(grads.Rows()*grads.Columns(), 0.0f);
		
		Matrix<int> labels(grads.Rows(), grads.Columns(), 0);

		Matrix<int> dummy(grads.Rows(), grads.Columns(), 0);

		int label = 0;


		// Point struct needs to be defined.

		queue<Point> outsQ;
		queue<Point> tmpQ;


		for (int x=0; x<grads.Columns(); x++) 
		{
			for (int y=0; y<grads.Rows(); y++) 
			{
				if(dummy(y,x) == 0 && thick(y,x) == 0)
				{
					label++;
					tmpQ.push(Point(x,y));
					dummy(y,x) = 1;
					labels(y,x) = label;
					LabelCount(label)++;