levelsetmethods.cpp

图像分割算法

			}
		}
		return alpha/(maxs+(maxs>1e-19?0:1));
	}

	double Getdt2DNormalVectorKappa(double alpha, double dx, double dy, double dx2, double dy2, Matrix<double>& H1Abs, Matrix<double>& H2Abs, Matrix<double>& b)
	{
		int endR = b.Rows()-1;
		int endC = b.Columns()-1;
		double maxs = 0.0;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				double temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy + 
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2;
				maxs = maxs>=temp?maxs:temp;
			}
		}
		return alpha/(maxs+(maxs>1e-19?0:1));
	}













	void Evolve2DKappa(Matrix<float>& phi, float dx, float dy, float dx2, float dy2, Matrix<float>& b, Matrix<float>& delta)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				float kappaAbsPhi = KappaAbsPhi2D(phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+1,j+1), 
					phi.ElemNC(i,j+1), phi.ElemNC(i-1,j+1), phi.ElemNC(i-1,j), phi.ElemNC(i-1,j-1), 
					phi.ElemNC(i,j-1), phi.ElemNC(i+1,j-1), dx, dy, dx2, dy2);
				delta.ElemNC(i,j) = b.ElemNC(i,j)*kappaAbsPhi;
			}
		}
	}




	inline float KappaAbsPhi2D(float phi_i_j, float phi_ip1_j, float phi_ip1_jp1, 
		float phi_i_jp1, float phi_im1_jp1, float phi_im1_j, float phi_im1_jm1, 
		float phi_i_jm1, float phi_ip1_jm1, float dx, float dy, float dx2, float dy2)
	{
		float phi_x = (phi_i_jp1 - phi_i_jm1)/(2*dx);
		float phi_y = (phi_ip1_j - phi_im1_j)/(2*dy);
		float phi_xx = (phi_i_jp1 - 2*phi_i_j + phi_i_jm1)/dx2;
		float phi_yy = (phi_ip1_j - 2*phi_i_j + phi_im1_j)/dy2;
		float phi_xy = (phi_ip1_jp1 + phi_im1_jm1 - phi_ip1_jm1 - phi_im1_jp1)/(4*dx*dy);
		float absGradPhiSq = (phi_x*phi_x + phi_y*phi_y);
		float kappaAbsPhi = (phi_xx*phi_y*phi_y - 2*phi_y*phi_x*phi_xy 
			+ phi_yy*phi_x*phi_x)/(absGradPhiSq + (float)(absGradPhiSq>1e-19?0:1));
		return kappaAbsPhi;

	}




	void Evolve2DKappa(Matrix<double>& phi, double dx, double dy, double dx2, double dy2, Matrix<double>& b, Matrix<double>& delta)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				double kappaAbsPhi = KappaAbsPhi2D(phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+1,j+1), 
					phi.ElemNC(i,j+1), phi.ElemNC(i-1,j+1), phi.ElemNC(i-1,j), phi.ElemNC(i-1,j-1), 
					phi.ElemNC(i,j-1), phi.ElemNC(i+1,j-1), dx, dy, dx2, dy2);
				delta.ElemNC(i,j) = b.ElemNC(i,j)*kappaAbsPhi;
			}
		}
	}




	inline double KappaAbsPhi2D(double phi_i_j, double phi_ip1_j, double phi_ip1_jp1, 
		double phi_i_jp1, double phi_im1_jp1, double phi_im1_j, double phi_im1_jm1, 
		double phi_i_jm1, double phi_ip1_jm1, double dx, double dy, double dx2, double dy2)
	{
		double phi_x = (phi_i_jp1 - phi_i_jm1)/(2*dx);
		double phi_y = (phi_ip1_j - phi_im1_j)/(2*dy);
		double phi_xx = (phi_i_jp1 - 2*phi_i_j + phi_i_jm1)/dx2;
		double phi_yy = (phi_ip1_j - 2*phi_i_j + phi_im1_j)/dy2;
		double phi_xy = (phi_ip1_jp1 + phi_im1_jm1 - phi_ip1_jm1 - phi_im1_jp1)/(4*dx*dy);
		double absGradPhiSq = (phi_x*phi_x + phi_y*phi_y);
		double kappaAbsPhi = (phi_xx*phi_y*phi_y - 2*phi_y*phi_x*phi_xy 
			+ phi_yy*phi_x*phi_x)/(absGradPhiSq + (double)(absGradPhiSq>1e-19?0:1));
		return kappaAbsPhi;

	}
















	void ReinitSignedDistanceENO1(Matrix<float>& phi, float dx, float dy, float alpha, int iterations)
	{
		Matrix<float> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<float> delta(phi.Rows(), phi.Columns());
		Matrix<float> H1Abs(phi.Rows(), phi.Columns());
		Matrix<float> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,1);
			Evolve2DNormalENO1(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			float dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}
	
	void ReinitSignedDistanceENO2(Matrix<float>& phi, float dx, float dy, float alpha, int iterations)
	{
		Matrix<float> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<float> delta(phi.Rows(), phi.Columns());
		Matrix<float> H1Abs(phi.Rows(), phi.Columns());
		Matrix<float> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,2);
			Evolve2DNormalENO2(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			float dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}

	void ReinitSignedDistanceENO3(Matrix<float>& phi, float dx, float dy, float alpha, int iterations)
	{
		Matrix<float> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<float> delta(phi.Rows(), phi.Columns());
		Matrix<float> H1Abs(phi.Rows(), phi.Columns());
		Matrix<float> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,3);
			Evolve2DNormalENO3(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			float dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}

	void ReinitSignedDistanceWENO(Matrix<float>& phi, float dx, float dy, float alpha, int iterations)
	{
		Matrix<float> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<float> delta(phi.Rows(), phi.Columns());
		Matrix<float> H1Abs(phi.Rows(), phi.Columns());
		Matrix<float> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,3);
			Evolve2DNormalWENO(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			float dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}





	void ReinitSignedDistanceENO1(Matrix<double>& phi, double dx, double dy, double alpha, int iterations)
	{
		Matrix<double> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<double> delta(phi.Rows(), phi.Columns());
		Matrix<double> H1Abs(phi.Rows(), phi.Columns());
		Matrix<double> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,1);
			Evolve2DNormalENO1(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			double dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}
	
	void ReinitSignedDistanceENO2(Matrix<double>& phi, double dx, double dy, double alpha, int iterations)
	{
		Matrix<double> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<double> delta(phi.Rows(), phi.Columns());
		Matrix<double> H1Abs(phi.Rows(), phi.Columns());
		Matrix<double> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,2);
			Evolve2DNormalENO2(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			double dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}

	void ReinitSignedDistanceENO3(Matrix<double>& phi, double dx, double dy, double alpha, int iterations)
	{
		Matrix<double> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<double> delta(phi.Rows(), phi.Columns());
		Matrix<double> H1Abs(phi.Rows(), phi.Columns());
		Matrix<double> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,3);
			Evolve2DNormalENO3(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			double dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}

	void ReinitSignedDistanceWENO(Matrix<double>& phi, double dx, double dy, double alpha, int iterations)
	{
		Matrix<double> S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx);
		Matrix<double> delta(phi.Rows(), phi.Columns());
		Matrix<double> H1Abs(phi.Rows(), phi.Columns());
		Matrix<double> H2Abs(phi.Rows(), phi.Columns());
		for(int k=0;k<iterations;k++)
		{
			Extend2D(phi,3);
			Evolve2DNormalWENO(phi, dx, dy, S_phi_0, delta, H1Abs, H2Abs);
			double dt = Getdt2DNormal(alpha, dx, dy, H1Abs, H2Abs);
			AddPhi2D(phi, S_phi_0-delta, dt);
		}
	}
















	void SubtractPhi2D(Matrix<float>& phi, Matrix<float>& delta, float dt)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				phi.ElemNC(i,j) -= dt*delta.ElemNC(i,j);
			}
		}
	}
	void AddPhi2D(Matrix<float>& phi, Matrix<float>& delta, float dt)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				phi.ElemNC(i,j) += dt*delta.ElemNC(i,j);
			}
		}
	}


	void SubtractPhi2D(Matrix<double>& phi, Matrix<double>& delta, double dt)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				phi.ElemNC(i,j) -= dt*delta.ElemNC(i,j);
			}
		}
	}
	void AddPhi2D(Matrix<double>& phi, Matrix<double>& delta, double dt)
	{
		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				phi.ElemNC(i,j) += dt*delta.ElemNC(i,j);
			}
		}
	}














	inline float SelectDerNormal(float phi_x_m, float phi_x_p, float Vn)
	{
		float phi_x;
		float VnPhi_m = Vn*phi_x_m;
		float VnPhi_p = Vn*phi_x_p;
		if(VnPhi_m <= 0 && VnPhi_p <= 0)
		{
			phi_x = phi_x_p;
		}
		else if(VnPhi_m >= 0 && VnPhi_p >= 0)
		{
			phi_x = phi_x_m;
		}
		else if(VnPhi_m <= 0 && VnPhi_p >= 0)
		{
			phi_x = 0;
		}
		else if(VnPhi_m >= 0 && VnPhi_p <= 0)
		{
			if(fabs(VnPhi_p) >= fabs(VnPhi_m))
			{
				phi_x = phi_x_p;
			}
			else
			{
				phi_x = phi_x_m;
			}
		}
		return phi_x;
	}


	inline float SelectDerNormalVectorSD(float phi_x_m, float phi_x_p, float u, float Vn)
	{
		float phi_x;

		float H1_m = u + Vn*phi_x_m;
		float H1_p = u + Vn*phi_x_p;

		if(H1_m <= 0 && H1_p <= 0)
		{
			phi_x = phi_x_p;
		}
		else if(H1_m >= 0 && H1_p >= 0)
		{
			phi_x = phi_x_m;
		}
		else if(H1_m <= 0 && H1_p >= 0)
		{
			phi_x = -u/Vn;
		}
		else if(H1_m >= 0 && H1_p <= 0)
		{
			if(fabs(H1_p) >= fabs(H1_m))
			{
				phi_x = phi_x_p;
			}
			else
			{
				phi_x = phi_x_m;
			}
		}
		return phi_x;
	}

	inline double SelectDerNormal(double phi_x_m, double phi_x_p, double Vn)
	{
		double phi_x;
		double VnPhi_m = Vn*phi_x_m;
		double VnPhi_p = Vn*phi_x_p;
		if(VnPhi_m <= 0 && VnPhi_p <= 0)
		{
			phi_x = phi_x_p;
		}
		else if(VnPhi_m >= 0 && VnPhi_p >= 0)
		{
			phi_x = phi_x_m;
		}
		else if(VnPhi_m <= 0 && VnPhi_p >= 0)
		{
			phi_x = 0;
		}
		else if(VnPhi_m >= 0 && VnPhi_p <= 0)
		{
			if(fabs(VnPhi_p) >= fabs(VnPhi_m))
			{
				phi_x = phi_x_p;
			}
			else
			{
				phi_x = phi_x_m;
			}
		}
		return phi_x;
	}


	inline double SelectDerNormalVectorSD(double phi_x_m, double phi_x_p, double u, double Vn)
	{
		double phi_x;

		double H1_m = u + Vn*phi_x_m;
		double H1_p = u + Vn*phi_x_p;

		if(H1_m <= 0 && H1_p <= 0)
		{
			phi_x = phi_x_p;
		}
		else if(H1_m >= 0 && H1_p >= 0)
		{
			phi_x = phi_x_m;
		}
		else if(H1_m <= 0 && H1_p >= 0)
		{
			phi_x = -u/Vn;
		}
		else if(H1_m >= 0 && H1_p <= 0)
		{
			if(fabs(H1_p) >= fabs(H1_m))
			{
				phi_x = phi_x_p;
			}
			else
			{
				phi_x = phi_x_m;
			}
		}
		return phi_x;
	}