levelsetmethods.cpp

图像分割算法

	// for (+) derivative v1 = phi_i, for (-) derivative v2 = phi_i
	inline float DerENO1Minus(float v1, float v2, float dx)
	{
		return (v2-v1)/dx;
	}
	inline float DerENO1Plus(float v1, float v2, float dx)
	{
		return (v2-v1)/dx;
	}

	// v3 = phi_i
	inline float DerENO2Minus(float v1, float v2, float v3, float v4, float dx)
	{
		float D1_1 = (v2-v1);
		float D1_2 = (v3-v2);
		float D1_3 = (v4-v3);

		float D2_1 = 0.5f*(D1_2-D1_1);
		float D2_2 = 0.5f*(D1_3-D1_2);

		float absD2_1 = fabs(D2_1);
		float absD2_2 = fabs(D2_2);
		float Q2p;
		if(absD2_1 <= absD2_2)
		{
			Q2p = D2_1;
		}
		else
		{
			Q2p = D2_2;
		}

		return (D1_2 + Q2p)/dx;
	}

	// v2 = phi_i
	inline float DerENO2Plus(float v1, float v2, float v3, float v4, float dx)
	{
		float D1_1 = (v2-v1);
		float D1_2 = (v3-v2);
		float D1_3 = (v4-v3);

		float D2_1 = 0.5f*(D1_2-D1_1);
		float D2_2 = 0.5f*(D1_3-D1_2);

		float absD2_1 = fabs(D2_1);
		float absD2_2 = fabs(D2_2);
		float Q2p;
		if(absD2_1 <= absD2_2)
		{
			Q2p = -D2_1;
		}
		else
		{
			Q2p = -D2_2;
		}

		return (D1_2 + Q2p)/dx;
	}

	// v4 = phi_i
	inline float DerENO3Minus(float v1, float v2, float v3, float v4, float v5, float v6, float dx)
	{
		float D1_1 = (v2-v1);
		float D1_2 = (v3-v2);
		float D1_3 = (v4-v3);
		float D1_4 = (v5-v4);
		float D1_5 = (v6-v5);

		float D2_1 = 0.5f*(D1_2-D1_1);
		float D2_2 = 0.5f*(D1_3-D1_2);
		float D2_3 = 0.5f*(D1_4-D1_3);
		float D2_4 = 0.5f*(D1_5-D1_4);

		float D3_1 = (D2_2-D2_1)/3;
		float D3_2 = (D2_3-D2_2)/3;
		float D3_3 = (D2_4-D2_3)/3;

		float absD2_2 = fabs(D2_2);
		float absD2_3 = fabs(D2_3);
		float Q2p;
		float Q3p;
		if(absD2_2 <= absD2_3)
		{
			Q2p = D2_2;
			float cstar;
			float absD3_1 = fabs(D3_1);
			float absD3_2 = fabs(D3_2);
			if(absD3_1 <= absD3_2)
			{
				cstar = D3_1;
			}
			else
			{
				cstar = D3_2;
			}
			Q3p = cstar*2;
		}
		else
		{
			Q2p = D2_3;
			float cstar;
			float absD3_2 = fabs(D3_2);
			float absD3_3 = fabs(D3_3);
			if(absD3_2 <= absD3_3)
			{
				cstar = D3_2;
			}
			else
			{
				cstar = D3_3;
			}
			Q3p = -cstar;
		}

		return (D1_3 + Q2p + Q3p)/dx;
	}

	// v3 = phi_i
	inline float DerENO3Plus(float v1, float v2, float v3, float v4, float v5, float v6, float dx)
	{
		float D1_1 = (v2-v1);
		float D1_2 = (v3-v2);
		float D1_3 = (v4-v3);
		float D1_4 = (v5-v4);
		float D1_5 = (v6-v5);

		float D2_1 = 0.5f*(D1_2-D1_1);
		float D2_2 = 0.5f*(D1_3-D1_2);
		float D2_3 = 0.5f*(D1_4-D1_3);
		float D2_4 = 0.5f*(D1_5-D1_4);

		float D3_1 = (D2_2-D2_1)/3;
		float D3_2 = (D2_3-D2_2)/3;
		float D3_3 = (D2_4-D2_3)/3;

		float absD2_2 = fabs(D2_2);
		float absD2_3 = fabs(D2_3);
		float Q2p;
		float Q3p;
		if(absD2_2 <= absD2_3)
		{
			Q2p = -D2_2;
			float cstar;
			float absD3_1 = fabs(D3_1);
			float absD3_2 = fabs(D3_2);
			if(absD3_1 <= absD3_2)
			{
				cstar = D3_1;
			}
			else
			{
				cstar = D3_2;
			}
			Q3p = -cstar;
		}
		else
		{
			Q2p = -D2_3;
			float cstar;
			float absD3_2 = fabs(D3_2);
			float absD3_3 = fabs(D3_3);
			if(absD3_2 <= absD3_3)
			{
				cstar = D3_2;
			}
			else
			{
				cstar = D3_3;
			}
			Q3p = cstar*2;
		}

		return (D1_3 + Q2p + Q3p)/dx;
	}


	// vv4 = phi_i
	inline float DerWENOMinus(float vv1, float vv2, float vv3, float vv4, float vv5, float vv6, float dx)
	{
		float v1 = (vv2-vv1)/dx;
		float v2 = (vv3-vv2)/dx;
		float v3 = (vv4-vv3)/dx;
		float v4 = (vv5-vv4)/dx;
		float v5 = (vv6-vv5)/dx;

		float data_x_1 = (2*v1 - 7*v2 + 11*v3)/6;
		float data_x_2 = (-v2 + 5*v3 + 2*v4)/6;
		float data_x_3 = (2*v3 + 5*v4 - v5)/6;

		float dummy1 = v1-2*v2+v3;
		float dummy2 = v1-4*v2+3*v3;
		float dummy3 = v2-2*v3+v4;
		float dummy4 = v2-v4;
		float dummy5 = v3-2*v4+v5;
		float dummy6 = 3*v3-4*v4+v5;
		float dummy_c = 13/12;
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2;
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4;
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6;

		float m = v1*v1;
		float temp = v2*v2;
		m = m>=temp?m:temp;
		temp = v3*v3;
		m = m>=temp?m:temp;
		temp = v4*v4;
		m = m>=temp?m:temp;
		temp = v5*v5;
		m = m>=temp?m:temp;
		double epsil = m*1e-6 + 1e-199;

		double alpha1 = 0.1/((S1+epsil)*(S1+epsil));
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil));
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil));
		return (float)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3));
	}

	// v3 = phi_i
	inline float DerWENOPlus(float vv1, float vv2, float vv3, float vv4, float vv5, float vv6, float dx)
	{

		float v5 = (vv2-vv1)/dx;
		float v4 = (vv3-vv2)/dx;
		float v3 = (vv4-vv3)/dx;
		float v2 = (vv5-vv4)/dx;
		float v1 = (vv6-vv5)/dx;

		float data_x_1 = (2*v1 - 7*v2 + 11*v3)/6;
		float data_x_2 = (-v2 + 5*v3 + 2*v4)/6;
		float data_x_3 = (2*v3 + 5*v4 - v5)/6;

		float dummy1 = v1-2*v2+v3;
		float dummy2 = v1-4*v2+3*v3;
		float dummy3 = v2-2*v3+v4;
		float dummy4 = v2-v4;
		float dummy5 = v3-2*v4+v5;
		float dummy6 = 3*v3-4*v4+v5;
		float dummy_c = 13/12;
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2;
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4;
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6;

		float m = v1*v1;
		float temp = v2*v2;
		m = m>=temp?m:temp;
		temp = v3*v3;
		m = m>=temp?m:temp;
		temp = v4*v4;
		m = m>=temp?m:temp;
		temp = v5*v5;
		m = m>=temp?m:temp;
		double epsil = m*1e-6 + 1e-199;

		double alpha1 = 0.1/((S1+epsil)*(S1+epsil));
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil));
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil));
		return (float)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3));
	}






	// for (+) derivative v1 = phi_i, for (-) derivative v2 = phi_i
	inline double DerENO1Minus(double v1, double v2, double dx)
	{
		return (v2-v1)/dx;
	}
	inline double DerENO1Plus(double v1, double v2, double dx)
	{
		return (v2-v1)/dx;
	}

	// v3 = phi_i
	inline double DerENO2Minus(double v1, double v2, double v3, double v4, double dx)
	{
		double D1_1 = (v2-v1);
		double D1_2 = (v3-v2);
		double D1_3 = (v4-v3);

		double D2_1 = 0.5f*(D1_2-D1_1);
		double D2_2 = 0.5f*(D1_3-D1_2);

		double absD2_1 = fabs(D2_1);
		double absD2_2 = fabs(D2_2);
		double Q2p;
		if(absD2_1 <= absD2_2)
		{
			Q2p = D2_1;
		}
		else
		{
			Q2p = D2_2;
		}

		return (D1_2 + Q2p)/dx;
	}

	// v2 = phi_i
	inline double DerENO2Plus(double v1, double v2, double v3, double v4, double dx)
	{
		double D1_1 = (v2-v1);
		double D1_2 = (v3-v2);
		double D1_3 = (v4-v3);

		double D2_1 = 0.5f*(D1_2-D1_1);
		double D2_2 = 0.5f*(D1_3-D1_2);

		double absD2_1 = fabs(D2_1);
		double absD2_2 = fabs(D2_2);
		double Q2p;
		if(absD2_1 <= absD2_2)
		{
			Q2p = -D2_1;
		}
		else
		{
			Q2p = -D2_2;
		}

		return (D1_2 + Q2p)/dx;
	}

	// v4 = phi_i
	inline double DerENO3Minus(double v1, double v2, double v3, double v4, double v5, double v6, double dx)
	{
		double D1_1 = (v2-v1);
		double D1_2 = (v3-v2);
		double D1_3 = (v4-v3);
		double D1_4 = (v5-v4);
		double D1_5 = (v6-v5);

		double D2_1 = 0.5f*(D1_2-D1_1);
		double D2_2 = 0.5f*(D1_3-D1_2);
		double D2_3 = 0.5f*(D1_4-D1_3);
		double D2_4 = 0.5f*(D1_5-D1_4);

		double D3_1 = (D2_2-D2_1)/3;
		double D3_2 = (D2_3-D2_2)/3;
		double D3_3 = (D2_4-D2_3)/3;

		double absD2_2 = fabs(D2_2);
		double absD2_3 = fabs(D2_3);
		double Q2p;
		double Q3p;
		if(absD2_2 <= absD2_3)
		{
			Q2p = D2_2;
			double cstar;
			double absD3_1 = fabs(D3_1);
			double absD3_2 = fabs(D3_2);
			if(absD3_1 <= absD3_2)
			{
				cstar = D3_1;
			}
			else
			{
				cstar = D3_2;
			}
			Q3p = cstar*2;
		}
		else
		{
			Q2p = D2_3;
			double cstar;
			double absD3_2 = fabs(D3_2);
			double absD3_3 = fabs(D3_3);
			if(absD3_2 <= absD3_3)
			{
				cstar = D3_2;
			}
			else
			{
				cstar = D3_3;
			}
			Q3p = -cstar;
		}

		return (D1_3 + Q2p + Q3p)/dx;
	}

	// v3 = phi_i
	inline double DerENO3Plus(double v1, double v2, double v3, double v4, double v5, double v6, double dx)
	{
		double D1_1 = (v2-v1);
		double D1_2 = (v3-v2);
		double D1_3 = (v4-v3);
		double D1_4 = (v5-v4);
		double D1_5 = (v6-v5);

		double D2_1 = 0.5f*(D1_2-D1_1);
		double D2_2 = 0.5f*(D1_3-D1_2);
		double D2_3 = 0.5f*(D1_4-D1_3);
		double D2_4 = 0.5f*(D1_5-D1_4);

		double D3_1 = (D2_2-D2_1)/3;
		double D3_2 = (D2_3-D2_2)/3;
		double D3_3 = (D2_4-D2_3)/3;

		double absD2_2 = fabs(D2_2);
		double absD2_3 = fabs(D2_3);
		double Q2p;
		double Q3p;
		if(absD2_2 <= absD2_3)
		{
			Q2p = -D2_2;
			double cstar;
			double absD3_1 = fabs(D3_1);
			double absD3_2 = fabs(D3_2);
			if(absD3_1 <= absD3_2)
			{
				cstar = D3_1;
			}
			else
			{
				cstar = D3_2;
			}
			Q3p = -cstar;
		}
		else
		{
			Q2p = -D2_3;
			double cstar;
			double absD3_2 = fabs(D3_2);
			double absD3_3 = fabs(D3_3);
			if(absD3_2 <= absD3_3)
			{
				cstar = D3_2;
			}
			else
			{
				cstar = D3_3;
			}
			Q3p = cstar*2;
		}

		return (D1_3 +