levelsetmethods.cpp

图像分割算法

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#include "./LevelSetMethods.h"


namespace LevelSetMethods
{


	Matrix<int> ExtractCurve(Matrix<float> phi)
	{
		Matrix<int> curve(phi.Rows(), phi.Columns(), 0);

		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		for(int i=4;i<=endR-4;i++)
		{
			for(int j=4;j<=endC-4;j++)
			{
				if(phi.ElemNC(i+1,j)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i+1,j) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i-1,j)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i-1,j) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i,j+1)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i,j+1) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i,j-1)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i,j-1) = 1;
					curve.ElemNC(i,j) = 1;
				}
			}
		}

		return curve;
	}


	Matrix<int> ExtractCurve(Matrix<double> phi)
	{
		Matrix<int> curve(phi.Rows(), phi.Columns(), 0);

		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		for(int i=4;i<=endR-4;i++)
		{
			for(int j=4;j<=endC-4;j++)
			{
				if(phi.ElemNC(i+1,j)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i+1,j) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i-1,j)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i-1,j) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i,j+1)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i,j+1) = 1;
					curve.ElemNC(i,j) = 1;
				}
				if(phi.ElemNC(i,j-1)*phi.ElemNC(i,j) <= 0)
				{
					curve.ElemNC(i,j-1) = 1;
					curve.ElemNC(i,j) = 1;
				}
			}
		}

		return curve;
	}




	Matrix<float>& ExtendConst2D(Matrix<float>& phi, int size)
	{
		if(size < 1 || size >3)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Extension size should be between 1 and 3!");
		}

		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		// Extrapolate first layer
		for(int i=3;i<=endR-3;i++)
		{
			phi.ElemNC(i,2) = phi.ElemNC(i,3);
			phi.ElemNC(i,endC-2) = phi.ElemNC(i,endC-3);
		}
		for(int j=3;j<=endC-3;j++)
		{
			phi.ElemNC(2,j) = phi.ElemNC(3,j);
			phi.ElemNC(endR-2,j) = phi.ElemNC(endR-3,j);
		}
		// May be useful for kappa 
		phi.ElemNC(2,2) = phi.ElemNC(3,3);
		phi.ElemNC(2,endC-2) = phi.ElemNC(3,endC-3);
		phi.ElemNC(endR-2,2) = phi.ElemNC(endR-3,3);
		phi.ElemNC(endR-2,endC-2) = phi.ElemNC(endR-3,endC-3);

		// Extend second layer
		if(size>=2)
		{
			for(int i=2;i<=endR-2;i++)
			{
				phi.ElemNC(i,1) = phi.ElemNC(i,2);
				phi.ElemNC(i,endC-1) = phi.ElemNC(i,endC-2);
			}
			for(int j=2;j<=endC-2;j++)
			{
				phi.ElemNC(1,j) = phi.ElemNC(2,j);
				phi.ElemNC(endR-1,j) = phi.ElemNC(endR-2,j);
			}
		}


		// Extend third layer
		if(size>=3)
		{
			for(int i=1;i<=endR-1;i++)
			{
				phi.ElemNC(i,0) = phi.ElemNC(i,1);
				phi.ElemNC(i,endC-0) = phi.ElemNC(i,endC-1);
			}
			for(int j=1;j<=endC-1;j++)
			{
				phi.ElemNC(0,j) = phi.ElemNC(1,j);
				phi.ElemNC(endR-0,j) = phi.ElemNC(endR-1,j);
			}
		}


		return phi;
	}

	Matrix<float>& Extend2D(Matrix<float>& phi, int size)
	{
		if(size < 1 || size >3)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Extension size should be between 1 and 3!");
		}

		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		// Extrapolate first layer
		for(int i=3;i<=endR-3;i++)
		{
			phi.ElemNC(i,2) = 2*phi.ElemNC(i,3) - phi.ElemNC(i,4);
			phi.ElemNC(i,endC-2) = 2*phi.ElemNC(i,endC-3) - phi.ElemNC(i,endC-4);
		}
		for(int j=3;j<=endC-3;j++)
		{
			phi.ElemNC(2,j) = 2*phi.ElemNC(3,j) - phi.ElemNC(4,j);
			phi.ElemNC(endR-2,j) = 2*phi.ElemNC(endR-3,j) - phi.ElemNC(endR-4,j);
		}
		// May be useful for kappa 
		phi.ElemNC(2,2) = 2*phi.ElemNC(3,3) - phi.ElemNC(4,4);
		phi.ElemNC(2,endC-2) = 2*phi.ElemNC(3,endC-3) - phi.ElemNC(4,endC-4);
		phi.ElemNC(endR-2,2) = 2*phi.ElemNC(endR-3,3) - phi.ElemNC(endR-4,4);
		phi.ElemNC(endR-2,endC-2) = 2*phi.ElemNC(endR-3,endC-3) - phi.ElemNC(endR-4,endC-4);


		// Extrapolate second layer
		if(size>=2)
		{
			for(int i=2;i<=endR-2;i++)
			{
				phi.ElemNC(i,1) = 2*phi.ElemNC(i,2) - phi.ElemNC(i,3);
				phi.ElemNC(i,endC-1) = 2*phi.ElemNC(i,endC-2) - phi.ElemNC(i,endC-3);
			}
			for(int j=2;j<=endC-2;j++)
			{
				phi.ElemNC(1,j) = 2*phi.ElemNC(2,j) - phi.ElemNC(3,j);
				phi.ElemNC(endR-1,j) = 2*phi.ElemNC(endR-2,j) - phi.ElemNC(endR-3,j);
			}
		}


		// Extrapolate third layer
		if(size>=3)
		{
			for(int i=1;i<=endR-1;i++)
			{
				phi.ElemNC(i,0) = 2*phi.ElemNC(i,1) - phi.ElemNC(i,2);
				phi.ElemNC(i,endC-0) = 2*phi.ElemNC(i,endC-1) - phi.ElemNC(i,endC-2);
			}
			for(int j=1;j<=endC-1;j++)
			{
				phi.ElemNC(0,j) = 2*phi.ElemNC(1,j) - phi.ElemNC(2,j);
				phi.ElemNC(endR-0,j) = 2*phi.ElemNC(endR-1,j) - phi.ElemNC(endR-2,j);
			}
		}

		return phi;
	}


	Matrix<double>& Extend2D(Matrix<double>& phi, int size)
	{
		if(size < 1 || size >3)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Extension size should be between 1 and 3!");
		}

		int endR = phi.Rows()-1;
		int endC = phi.Columns()-1;

		// Extrapolate first layer
		for(int i=3;i<=endR-3;i++)
		{
			phi.ElemNC(i,2) = 2*phi.ElemNC(i,3) - phi.ElemNC(i,4);
			phi.ElemNC(i,endC-2) = 2*phi.ElemNC(i,endC-3) - phi.ElemNC(i,endC-4);
		}
		for(int j=3;j<=endC-3;j++)
		{
			phi.ElemNC(2,j) = 2*phi.ElemNC(3,j) - phi.ElemNC(4,j);
			phi.ElemNC(endR-2,j) = 2*phi.ElemNC(endR-3,j) - phi.ElemNC(endR-4,j);
		}
		// May be useful for kappa 
		phi.ElemNC(2,2) = 2*phi.ElemNC(3,3) - phi.ElemNC(4,4);
		phi.ElemNC(2,endC-2) = 2*phi.ElemNC(3,endC-3) - phi.ElemNC(4,endC-4);
		phi.ElemNC(endR-2,2) = 2*phi.ElemNC(endR-3,3) - phi.ElemNC(endR-4,4);
		phi.ElemNC(endR-2,endC-2) = 2*phi.ElemNC(endR-3,endC-3) - phi.ElemNC(endR-4,endC-4);


		// Extrapolate second layer
		if(size>=2)
		{
			for(int i=2;i<=endR-2;i++)
			{
				phi.ElemNC(i,1) = 2*phi.ElemNC(i,2) - phi.ElemNC(i,3);
				phi.ElemNC(i,endC-1) = 2*phi.ElemNC(i,endC-2) - phi.ElemNC(i,endC-3);
			}
			for(int j=2;j<=endC-2;j++)
			{
				phi.ElemNC(1,j) = 2*phi.ElemNC(2,j) - phi.ElemNC(3,j);
				phi.ElemNC(endR-1,j) = 2*phi.ElemNC(endR-2,j) - phi.ElemNC(endR-3,j);
			}
		}


		// Extrapolate third layer
		if(size>=3)
		{
			for(int i=1;i<=endR-1;i++)
			{
				phi.ElemNC(i,0) = 2*phi.ElemNC(i,1) - phi.ElemNC(i,2);
				phi.ElemNC(i,endC-0) = 2*phi.ElemNC(i,endC-1) - phi.ElemNC(i,endC-2);
			}
			for(int j=1;j<=endC-1;j++)
			{
				phi.ElemNC(0,j) = 2*phi.ElemNC(1,j) - phi.ElemNC(2,j);
				phi.ElemNC(endR-0,j) = 2*phi.ElemNC(endR-1,j) - phi.ElemNC(endR-2,j);
			}
		}

		return phi;
	}










	void Evolve2DNormalENO1(Matrix<float>& phi, float dx, float dy, 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 = SelectDerNormal(phi_x_m, phi_x_p, 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 = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = (float)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<float>& phi, float dx, float dy, 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 = SelectDerNormal(phi_x_m, phi_x_p, 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 = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = (float)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<float>& phi, float dx, float dy, 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 = SelectDerNormal(phi_x_m, phi_x_p, 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 = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = (float)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<float>& phi, float dx, float dy, 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 = SelectDerNormal(phi_x_m, phi_x_p, 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 = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j));

				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y);
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) );
				H2Abs.ElemNC(i,j) = (float)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 Evolve2DNormalENO1(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 = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx);