levelsetmethods.cpp

图像分割算法

				double phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j));

				double phi_y_m = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy);
				double phi_y_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi;
			}
		}
	}

	void Evolve2DNormalENO2(Matrix<double>& phi, double dx, double dy, Matrix<double>& Vn, Matrix<double>& delta, Matrix<double>& H1Abs, Matrix<double>& H2Abs)
	{
		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 phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				double phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j));

				double phi_y_m = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				double phi_y_p = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi;
			}
		}
	}

	void Evolve2DNormalENO3(Matrix<double>& phi, double dx, double dy, Matrix<double>& Vn, Matrix<double>& delta, Matrix<double>& H1Abs, Matrix<double>& H2Abs)
	{
		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 phi_x_m = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				double phi_x_p = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j));

				double phi_y_m = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				double phi_y_p = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi;
			}
		}
	}

	void Evolve2DNormalWENO(Matrix<double>& phi, double dx, double dy, Matrix<double>& Vn, Matrix<double>& delta, Matrix<double>& H1Abs, Matrix<double>& H2Abs)
	{
		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 phi_x_m = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				double phi_x_p = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j));

				double phi_y_m = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				double phi_y_p = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi;
			}
		}
	}






















	void Evolve2DVectorENO1(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				}
				else
				{
					phi_x = 0;
				}

				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}


	void Evolve2DVectorENO2(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else
				{
					phi_x = 0;
				}

				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}


	void Evolve2DVectorENO3(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				}
				else
				{
					phi_x = 0;
				}


				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}

	void Evolve2DVectorWENO(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				}
				else
				{
					phi_x = 0;
				}


				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}






	void Evolve2DVectorENO1(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				}
				else
				{
					phi_x = 0;
				}

				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}


	void Evolve2DVectorENO2(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else
				{
					phi_x = 0;
				}

				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}


	void Evolve2DVectorENO3(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				}
				else
				{
					phi_x = 0;
				}


				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}

	void Evolve2DVectorWENO(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 phi_x, phi_y;
				if(u.ElemNC(i,j)>0)
				{
					phi_x = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				}
				else if(u.ElemNC(i,j)<0)
				{
					phi_x = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx);
				}
				else
				{
					phi_x = 0;
				}


				if(v.ElemNC(i,j)>0)
				{
					phi_y = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy);
				}
				else if(v.ElemNC(i,j)<0)
				{
					phi_y = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy);
				}
				else
				{
					phi_y = 0;
				}

				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y;
			}
		}
	}