levelsetmethods.cpp

图像分割算法

	void Evolve2DNormalVectorENO1SignedDistance(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, Matrix<float>& Vn, Matrix<float>& delta, Matrix<float>& H1Abs, Matrix<float>& 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++)
			{
				float phi_x_m = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx);
				float phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j));

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

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

	void Evolve2DNormalVectorENO2SignedDistance(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, Matrix<float>& Vn, Matrix<float>& delta, Matrix<float>& H1Abs, Matrix<float>& 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++)
			{
				float phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				float phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx);
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j));

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

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


	void Evolve2DNormalVectorENO3SignedDistance(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, Matrix<float>& Vn, Matrix<float>& delta, Matrix<float>& H1Abs, Matrix<float>& 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++)
			{
				float 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);
				float 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);
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j));

				float 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);
				float 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);
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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


	void Evolve2DNormalVectorWENOSignedDistance(Matrix<float>& phi, float dx, float dy, Matrix<float>& u, Matrix<float>& v, Matrix<float>& Vn, Matrix<float>& delta, Matrix<float>& H1Abs, Matrix<float>& 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++)
			{
				float 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);
				float 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);
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j));

				float 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);
				float 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);
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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





	void Evolve2DNormalVectorENO1SignedDistance(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx);
				double phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx);
				double phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), 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 = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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

	void Evolve2DNormalVectorENO2SignedDistance(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), 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 = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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


	void Evolve2DNormalVectorENO3SignedDistance(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), 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 = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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


	void Evolve2DNormalVectorWENOSignedDistance(Matrix<double>& phi, double dx, double dy, Matrix<double>& u, Matrix<double>& v, 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 = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), 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 = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j));

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
















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

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


	float Getdt2DVector(float alpha, float dx, float dy, Matrix<float>& u, Matrix<float>& v)
	{
		int endR = u.Rows()-1;
		int endC = u.Columns()-1;
		float maxs = 0.0;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				float temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy;
				maxs = maxs>=temp?maxs:temp;
			}
		}
		return (float)(alpha/(maxs+(maxs>1e-19?0:1)));
	}



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

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

	float Getdt2DNormalKappa(float alpha, float dx, float dy, float dx2, float dy2, Matrix<float>& H1Abs, Matrix<float>& H2Abs, Matrix<float>& b)
	{
		int endR = b.Rows()-1;
		int endC = b.Columns()-1;
		float maxs = 0.0;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				float 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 (float)(alpha/(maxs+(maxs>1e-19?0:1)));
	}

	float Getdt2DNormalVectorKappa(float alpha, float dx, float dy, float dx2, float dy2, Matrix<float>& H1Abs, Matrix<float>& H2Abs, Matrix<float>& b)
	{
		int endR = b.Rows()-1;
		int endC = b.Columns()-1;
		float maxs = 0.0;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				float 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 (float)(alpha/(maxs+(maxs>1e-19?0:1)));
	}








	double Getdt2DNormal(double alpha, double dx, double dy, Matrix<double>& H1Abs, Matrix<double>& H2Abs)
	{
		int endR = H1Abs.Rows()-1;
		int endC = H1Abs.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;
				maxs = maxs>=temp?maxs:temp;
			}
		}
		return alpha/(maxs+(maxs>1e-19?0:1));
	}

	double Getdt2DNormalVector(double alpha, double dx, double dy, Matrix<double>& H1Abs, Matrix<double>& H2Abs)
	{
		int endR = H1Abs.Rows()-1;
		int endC = H1Abs.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;
				maxs = maxs>=temp?maxs:temp;
			}
		}
		return alpha/(maxs+(maxs>1e-19?0:1));
	}


	double Getdt2DVector(double alpha, double dx, double dy, Matrix<double>& u, Matrix<double>& v)
	{
		int endR = u.Rows()-1;
		int endC = u.Columns()-1;
		double maxs = 0.0;
		for(int j=3;j<=endC-3;j++)
		{
			for(int i=3;i<=endR-3;i++)
			{
				double temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy;
				maxs = maxs>=temp?maxs:temp;
			}
		}
		return alpha/(maxs+(maxs>1e-19?0:1));
	}



	double Getdt2DKappa(double alpha, double dx2, double dy2, 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 = (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));
	}

	double Getdt2DVectorKappa(double alpha, double dx, double dy, double dx2, double dy2, Matrix<double>& u, Matrix<double>& v, 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 = fabs(u.ElemNC(i,j))/dx + fabs(v.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));
	}

	double Getdt2DNormalKappa(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;