www.pudn.com > Image_segment.rar > PDE.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 "./PDE.h" 
 
 
 
 
 
namespace PDE 
{ 
 
	Vector tridiagonal_solve(Vector& a, Vector& b, Vector& c, Vector& w) 
	{ 
		// Finds the tridiagonal solution for Mu = w, M is NxN 
		// a,b,c are the diagonals from left to right all length N 
		// with a(0) and c(end) are irrelevant; 
		 
		int n = b.Numel(); 
		Vector x(n); 
		Vector y(n); 
		Vector u(n); 
		 
		x[n-1] = -a[n-1]/b[n-1]; 
		y[n-1] = w[n-1]/b[n-1]; 
		 
		for(int i=n-2;i>=1;i--) 
		{ 
			x[i] = -a[i]/(b[i] + c[i]*x[i+1]); 
			y[i] = (w[i] - c[i]*y[i+1])/(b[i] + c[i]*x[i+1]); 
		} 
		 
		x[0] = 0; 
		y[0] = (w[0] - c[0]*y[1])/(b[0] + c[0]*x[1]); 
		 
		u[0] = y[0]; 
		for(int i=1; i<=n-1; i++) 
		{ 
			u[i] = x[i]*u[i-1]+y[i]; 
		} 
		 
		return u; 
	} 
 
 
 
	Vector Poisson1D(Vector& v, float dx, BoundaryCondition boundary) 
	{ 
 
		//tridiagonals 
		Vector a; 
		Vector b; 
		Vector c; 
 
		Vector ou;  
		if(boundary == Neumann) 
		{ 
			a = Vector(v.Numel(), 1); 
			b = Vector(v.Numel(), -2); 
			c = Vector(v.Numel(), 1); 
 
			b(0) = -1; 
			b(b.Numel()-1) = -1; 
			ou = tridiagonal_solve(a, b, c, v*(dx*dx)); 
		} 
		else 
		{ 
			a = Vector(v.Numel()-2, 1); 
			b = Vector(v.Numel()-2, -2); 
			c = Vector(v.Numel()-2, 1); 
			ou = Vector(v.Numel(), 0); 
			ou(1, v.Numel()-2) = tridiagonal_solve(a, b, c, v.Slice(1, v.Numel()-2)*(dx*dx)); 
		} 
 
		//Vector a(v.Numel()-2, 1); 
		//Vector b(v.Numel()-2, -2); 
		//Vector c(v.Numel()-2, 1); 
 
		return ou; 
	} 
 
	Vector Poisson1DFFT(Vector& v, float dx, BoundaryCondition boundary) 
	{ 
		Vector V; 
		if(boundary == Neumann) 
		{ 
			V = FFTCos(v); 
		} 
		else 
		{ 
			V = FFTSin(v); 
		} 
 
		V(0) = 0; 
		for(int i=1; i poi; 
		if(boundary == Neumann) 
		{ 
			poi = IFFTCos(V); 
		} 
		else 
		{ 
			poi = IFFTSin(V); 
		} 
		return poi; 
	} 
 
	Matrix Poisson2DFFT(Matrix& v, float dx, BoundaryCondition boundary) 
	{ 
		Matrix V; 
		if(boundary == Neumann) 
		{ 
			V = FFT2Cos(v); 
		} 
		else 
		{ 
			//V = FFTSin(v); 
		} 
 
		for(int i=0; i poi; 
		if(boundary == Neumann) 
		{ 
			poi = IFFT2Cos(V); 
		} 
		else 
		{ 
			//poi = IFFTSin(V); 
		} 
		return poi; 
	} 
 
	//Matrix Poisson2D(Matrix& v, float dx, BoundaryCondition boundary) 
	//{ 
 
	//} 
 
	Vector FFTSin(Vector& m) 
	{ 
		int end = m.Numel()-1; 
		Vector mm(2*end); 
		mm(0) = m(0); 
		mm(end) = m(0); 
		for(int i=1; i MM = FFT(mm); 
		 
		Vector M(end+1); 
		M(0) = 0; 
		M(end) = 0; 
		for(int i=1; i IFFTSin(Vector& M) 
	{ 
		int end = M.Numel()-1; 
		Vector MM(2*end); 
		MM(0) = 0; 
		MM(end) = 0; 
		for(int i=1; i mm = IFFT(MM); 
		 
		Vector m(end+1); 
		m(0) = 0; 
		m(end) = 0; 
		for(int i=1; i FFTCos(Vector& m) 
	{ 
		int end = m.Numel()-1; 
		Vector mm(2*end); 
		mm(0) = m(0); 
		mm(end) = m(end); 
		for(int i=1; i MM = FFT(mm); 
		 
		Vector M(end+1); 
		M(0) = real(MM(0)); 
		M(end) = real(MM(end)); 
		for(int i=1; i IFFTCos(Vector& M) 
	{ 
		int end = M.Numel()-1; 
		Vector MM(2*end); 
		MM(0) = M(0); 
		MM(end) = M(end); 
		for(int i=1; i mm = IFFT(MM); 
		 
		Vector m(end+1); 
		m(0) = real(mm(0)); 
		m(end) = real(mm(end)); 
		for(int i=1; i FFT2Cos(Matrix& m) 
	{ 
		int endR = m.Rows()-1; 
		int endC = m.Columns()-1; 
		Matrix mm = Matrix::Cat(1, (m , FlipLRI(m.Slice(0,endR,1,endC-1))), (FlipUDI(m.Slice(1,endR-1,0,endC)), FlipUDI(FlipLRI(m.Slice(1,endR-1,1,endC-1)))) ); 
 
		Matrix MM = FFT2(ToComplexFloat(mm)); 
		 
		Matrix M(m.Rows(), m.Columns()); 
		for(int i=1; i IFFT2Cos(Matrix& M) 
	{ 
		int endR = M.Rows()-1; 
		int endC = M.Columns()-1; 
 
		Matrix MM(M.Rows(),M.Columns()); 
		for(int i=1; i::Cat(1, (MM , FlipLRI(MM.Slice(0,endR,1,endC-1))), (FlipUDI(MM.Slice(1,endR-1,0,endC)), FlipUDI(FlipLRI(MM.Slice(1,endR-1,1,endC-1)))) ); 
 
		Matrix mm = Real(IFFT2(ToComplexFloat(MM))); 
		return mm.Slice(0,endR,0,endC); 
	} 
 
 
 
};