imagesegmentation.cpp

图像分割算法

					LabelAvg(label) += grads(y,x);

					while(tmpQ.size() > 0)
					{
						Point p = tmpQ.front();
						tmpQ.pop();

						for (int i=-1; i<=1; i++) 
						{
							for (int j=-1; j<=1; j++) 
							{
								if(i == 0 || j == 0) 
								{
									if(p.x+i>=0 && p.x+i<grads.Columns() && p.y+j>=0 && p.y+j<grads.Rows()) 
									{
										if(dummy(p.y+j,p.x+i) == 0 && thick(p.y+j,p.x+i) == 0)
										{
											tmpQ.push(Point(p.x+i,p.y+j));
											dummy(p.y+j,p.x+i) = 1;
											labels(p.y+j,p.x+i) = label;
											LabelCount(label)++;
											LabelAvg(label) += grads(p.y+j,p.x+i);
										}
									}
								}
							}
						}
					}

				}
				else if(thick(y,x) == 1)
				{
					outsQ.push(Point(x,y));
				}
			}
		}


		for(int i=1;i<=label;i++)
		{
			LabelAvg(i) /= LabelCount(i);
			//					Console.WriteLine(i + ": " + LabelCount[i] + " : " + LabelAvg[i]);
		}
		




		// Utilize suppressed edges and use sorting for gradients...
		// Region growing here...

		CompareGrads::gradMatrix = grads;


		Matrix<int> tmp(grads.Rows(), grads.Columns(), 0);
		//Hashtable Marked = new Hashtable();
		bool greedy = false;
		int counter = 0;
		while(outsQ.size() > 0)
		{
			counter++;
			//Marked.Clear();
			int len = outsQ.size();
			//cout << "Length: " << len << endl;

			// copy queue to a list
			vector<Point> psA;
			for(int i=0; i<len; i++)
			{
				psA.push_back(outsQ.front());
				outsQ.pop();
			}
			std::sort( psA.begin( ), psA.end( ), CompareGrads::Compare );

			int marked = 0;
			for(int i=0; i<psA.size(); i++)
			{
				Point p = psA[i];
				if(labels(p.y,p.x) == 0)
				{
					vector<Point> A;


					for (int i=-1; i<=1; i++) 
					{
						for (int j=-1; j<=1; j++) 
						{
							if(i == 0 || j == 0) 
							{
								if(p.x+i>=0 && p.x+i<grads.Columns() && p.y+j>=0 && p.y+j<grads.Rows()) 
								{
									if(labels(p.y+j,p.x+i) != 0)
									{
										A.push_back(Point(p.x+i,p.y+j));
									}
								}
							}
						}
					}

					if(A.size() == 0)
					{
						outsQ.push(p);
					}
					else if(A.size() == 1)
					{
						Point p2 = A[0];
						if(suppressed(p.y,p.x) == 0 || greedy)
						{
							labels(p.y,p.x) = labels(p2.y,p2.x);
							marked++;
						}
						else
						{
							outsQ.push(p);
						}
					}
					else
					{
						Point pp = A[0];
						bool equal = true;
						for(int j=0;j<A.size();j++)
						{
							Point pp2 = A[j];
							if(labels(pp.y,pp.x) != labels(pp2.y,pp2.x))
							{
								equal = false;
							}
						}

						if(equal)
						{
							if(suppressed(p.y,p.x) == 0 || greedy)
							{
								labels(p.y,p.x) = labels(pp.y,pp.x);
								marked++;
							}
							else
							{
								outsQ.push(p);
							}
						}
					}
				}
				else
				{
					cout << "Something is wrong: A labeled point is added to the Queue!!!" << endl;
				}
			}


			if(marked == 0 && greedy)
			{
				break;
			}

			if(marked == 0)
			{
				greedy = true;
			}
		}
		






		for (int x=0; x<grads.Columns(); x++) 
		{
			for (int y=0; y<grads.Rows(); y++) 
			{
				if(labels(y,x) == 0)
				{
					edges(y,x) = 1;
				}
			}
		}

		cout << "Extracting boundaries finished!" << endl;

		return edges;


	
	
	}


	float Angle(float x1, float y1, float x2,  float y2)
	{
		float scalar = x1*x2+y1*y2;
		float m1 = sqrt(x1*x1+y1*y1);
		float m2 = sqrt(x2*x2+y2*y2);
		float ac = scalar/(m1*m2);
		ac = ac>1?1:ac;
		ac = ac<-1?-1:ac;
		return acos(ac)*180/PI;
	}


//The following code for RGB to Lab conversion is adapted from code 
//written by Yossi Rubner and Mark Ruzon. Following comments belong to them.
//Baris Sumengen 09/09/2005
//
//	Convert between RGB and CIE-Lab color spaces
//	Uses ITU-R recommendation BT.709 with D65 as reference white.
//	Yossi Rubner 2/24/98
//
//	Mark Ruzon 3/18/99 -- Added thresholds to truncate parts of the space
//	where changing a value has no visible effect on the color.


	void RGB2Lab(double R, double G, double B, double &L, double &a, double &b)
	{
		const double BLACK = 20;
		const double YELLOW = 70;
		double X, Y, Z, fX, fY, fZ;

		X = 0.412453*R + 0.357580*G + 0.180423*B;
		Y = 0.212671*R + 0.715160*G + 0.072169*B;
		Z = 0.019334*R + 0.119193*G + 0.950227*B;

		X /= (255 * 0.950456);
		Y /=  255;
		Z /= (255 * 1.088754);

		if(Y > 0.008856)
		{
			fY = pow(Y, 1.0/3.0);
			L = 116.0*fY - 16.0;
		}
		else
		{
			fY = 7.787*Y + 16.0/116.0;
			L = 903.3*Y;
		}

		if(X > 0.008856)
		{
			fX = pow(X, 1.0/3.0);
		}
		else
		{
			fX = 7.787*X + 16.0/116.0;
		}

		if(Z > 0.008856)
		{
			fZ = pow(Z, 1.0/3.0);
		}
		else
		{
			fZ = 7.787*Z + 16.0/116.0;
		}

		a = 500.0*(fX - fY);
		b = 200.0*(fY - fZ);

		if(L < BLACK) 
		{
			a *= exp((L - BLACK) / (BLACK / 4));
			b *= exp((L - BLACK) / (BLACK / 4));
			L = BLACK;
		}
		if(b > YELLOW)
		{
			b = YELLOW;
		}
	}


	void Lab2RGB(double L, double a, double b, double &R, double &G, double &B)
	{
		const double BLACK = 20;
		const double YELLOW = 70;
		double X, Y, Z, fX, fY, fZ;
		double RR, GG, BB;

		fY = pow((L + 16.0) / 116.0, 3.0);
		if(fY < 0.008856)
		{
			fY = L / 903.3;
		}
		Y = fY;

		if(fY > 0.008856)
		{
			fY = pow(fY, 1.0/3.0);
		}
		else
		{
			fY = 7.787 * fY + 16.0/116.0;
		}

		fX = a / 500.0 + fY;      
		if(fX > 0.206893)
		{
			X = pow(fX, 3.0);
		}
		else
		{
			X = (fX - 16.0/116.0) / 7.787;
		}

		fZ = fY - b /200.0;      
		if(fZ > 0.206893)
		{
			Z = pow(fZ, 3.0);
		}
		else
		{
			Z = (fZ - 16.0/116.0) / 7.787;
		}

		X *= (0.950456 * 255);
		Y *=             255;
		Z *= (1.088754 * 255);

		RR =  3.240479*X - 1.537150*Y - 0.498535*Z;
		GG = -0.969256*X + 1.875992*Y + 0.041556*Z;
		BB =  0.055648*X - 0.204043*Y + 1.057311*Z;

		R = (double)(RR < 0 ? 0 : RR > 255 ? 255 : RR);
		G = (double)(GG < 0 ? 0 : GG > 255 ? 255 : GG);
		B = (double)(BB < 0 ? 0 : BB > 255 ? 255 : BB);

	}


	MatrixList<float> RGB2Lab(MatrixList<float> input)
	{
		if(input.Planes() != 3)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Length of MatrixList should be 3 for an RGB image!");
		}
		MatrixList<float> tmp(3, input.Rows(),input.Columns());
		for(int i=0;i<input.Rows();i++)
		{
			for(int j=0;j<input.Columns();j++)
			{
				double L,a,b;
				L = a = b = 0;
				RGB2Lab(input[0].ElemNC(i,j),input[1].ElemNC(i,j),input[2].ElemNC(i,j), L, a, b);
				tmp[0].ElemNC(i,j) = (float)L;
				tmp[1].ElemNC(i,j) = (float)a;
				tmp[2].ElemNC(i,j) = (float)b;
			}
		}
		return tmp;

	}


	MatrixList<float> Lab2RGB(MatrixList<float> input)
	{
		if(input.Planes() != 3)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Length of MatrixList should be 3 for an Lab image!");
		}
		MatrixList<float> tmp(3, input.Rows(),input.Columns());
		for(int i=0;i<input.Rows();i++)
		{
			for(int j=0;j<input.Columns();j++)
			{
				double R,G,B;
				R = G = B = 0;
				Lab2RGB(input[0].ElemNC(i,j),input[1].ElemNC(i,j),input[2].ElemNC(i,j), R, G, B);
				tmp[0].ElemNC(i,j) = (float)R;
				tmp[1].ElemNC(i,j) = (float)G;
				tmp[2].ElemNC(i,j) = (float)B;
			}
		}
		return tmp;

	}


	extern "C" float *poisson(float *tmp, int X_DIM, int Y_DIM);

	Matrix<int> SegmentEF(Matrix<float> &im, bool normalized, float initScale, float scaleJump, 
		float endScale, float angleLimit, float ratioLimit, float smoothWeighting, float stopError, int accuracy)
	{
		float diag = sqrt(float(im.Rows()*im.Rows() + im.Columns()*im.Columns()));
		float atomic = diag/1000.0;
		cout << "Unit scale: " << atomic << " pixels" << endl << endl;
		float scale = initScale*atomic;

		cout << "Flow field at scale: " << scale << endl;
		MatrixList<float> flow = EdgeflowVectorField(im, 8, scale, scale, 1.0f, normalized);
		endScale = endScale*atomic;
		scaleJump = scaleJump*atomic;
		scale += scaleJump;
		while( scale <=  endScale)
		{
			Matrix<float> efMag = Sqrt(flow(0)*flow(0)+flow(1)*flow(1));
			float magLimit = Maximum(efMag(":"))/ratioLimit;

			cout << "Flow field at scale: " << scale << " pixels" << endl;
			MatrixList<float> tempflow = EdgeflowVectorField(im, 4, scale, scale, 1.0f, normalized);
			Matrix<float> tempMag = Sqrt(tempflow(0)*tempflow(0)+tempflow(1)*tempflow(1));

			for (int i=0; i<flow.Rows(); i++) {
				for (int j=0; j<flow.Columns(); j++) {
					if(efMag(i,j) < magLimit)
					{
						flow(0)(i,j) = tempflow(0)(i,j);
						flow(1)(i,j) = tempflow(1)(i,j);
					}
					else if( Angle(flow(0)(i,j),flow(1)(i,j),tempflow(0)(i,j),tempflow(1)(i,j)) < angleLimit )
					{
						flow(0)(i,j) += tempflow(0)(i,j);
						flow(1)(i,j) += tempflow(1)(i,j);
					}
				}
			}
			scale += scaleJump;
		}


		Matrix<float> BB = Divergence(flow(0), flow(1));
		float *tmp = poisson((-BB).Data(), BB.Columns(), BB.Rows());
		Matrix<float> CC(BB.Rows(), BB.Columns());
		for (int j=0; j<CC.Rows(); j++) {
			for (int i=0; i<CC.Columns(); i++) {
				CC(j,i) = tmp[i+j*CC.Columns()];
			}
		}
		ImWrite(CC, "edgefunction.png");
		float min = Minimum(CC(":"));
		float max = Maximum(CC(":"));
		CC -= min;
		CC *= 2.0f/(max-min);
		flow(0) *= 40.0f/(max-min);
		flow(1) *= 40.0f/(max-min);

		Matrix<float> b(im.Rows()+6, im.Columns()+6,0);
		b(3,im.Rows()+2,3,im.Columns()+2) = smoothWeighting*CC;

		Matrix<float> u(im.Rows()+6, im.Columns()+6,0);
		u(3,im.Rows()+2,3,im.Columns()+2) = flow(0);

		Matrix<float> v(im.Rows()+6, im.Columns()+6,0);
		v(3,im.Rows()+2,3,im.Columns()+2) = flow(1);

		Matrix<float> phi(im.Rows()+6, im.Columns()+6);
		phi(3,im.Rows()+2,3,im.Columns()+2) = im;

		Matrix<float> delta(im.Rows()+6, im.Columns()+6);
		Matrix<float> delta2(im.Rows()+6, im.Columns()+6);
		Matrix<float> old = im.Clone();

		PerfTimer pt200;
		pt200.Tic();
		float oldChange = 1e10;
		int k=0;
		while(1)
		{
			ExtendConst2D(phi,3);
			Evolve2DKappa(phi, 1, 1, 1, 1, b, delta);
			switch( accuracy )
			{
				case 1:
					Evolve2DVectorENO1(phi, 1, 1, u, v, delta2);
					break;
				case 2:
					Evolve2DVectorENO2(phi, 1, 1, u, v, delta2);
					break;
				case 3:
					Evolve2DVectorENO3(phi, 1, 1, u, v, delta2);
					break;
				case 5:
					Evolve2DVectorWENO(phi, 1, 1, u, v, delta2);
					break;
				default:
					cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
					Utility::RunTimeError("accuracy parameter can only be 1, 2, 3, or 5!");
			}

			float dt = Getdt2DVectorKappa(0.95f, 1, 1, 1, 1, u, v, b);
			AddPhi2D(phi,delta,dt);
			SubtractPhi2D(phi,delta2,dt);
			if(k%200 == 199 ){
				cout << "EF: " << k+1 << endl;
				pt200.Toc();
				char str[100];
				sprintf(str, "phi_%03d.bmp", k);
				for(int i=0; i<phi.Numel(); i++)
				{
					if(phi(i) < 0)
					{
						phi(i) = 0;
					}
					if(phi(i) > 255)
					{
						phi(i) = 255;
					}
				}
				Matrix<float> diffused = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
				//	ImWrite(diffused, str);
				float change = Sum(Vector<float>(Abs(old-diffused)))/im.Numel();
				cout << "Error: "  << change << endl << endl;
				if(change < stopError || k > 10000)
				{
					break;
				}
				if(change > oldChange)
				{

					Utility::Warning("!!!!! Instaability detected during diffusion... Exiting diffusion stage!");
					break;
				}
				oldChange = change;
				old = diffused.Clone();
				pt200.Tic();
			}
			k++;
		}

		phi = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
		for(int i=0; i<phi.Numel(); i++)
		{
			if(phi(i) < 0)
			{
				phi(i) = 0;
			}
			if(phi(i) > 255)
			{
				phi(i) = 255;
			}
		}

		ImWrite(phi, "diffused.png");

		Matrix<float> gx = FilterFDGaussian2D(phi, atomic, 0);
		Matrix<float> gy = FilterFDGaussian2D(phi, atomic, 90);
		Matrix<float> gradIm = Sqrt(SQR(gx) + SQR(gy));
		Matrix<float> gradSupp = NonMaximaSuppress(gradIm, gx, gy);
		Matrix<int> edges = GetEgdes(gradIm, gradIm >