00001 /* LIBGUL - Geometry Utility Library 00002 * Copyright (C) 1998-1999 Norbert Irmer 00003 * 00004 * This library is free software; you can redistribute it and/or 00005 * modify it under the terms of the GNU Library General Public 00006 * License as published by the Free Software Foundation; either 00007 * version 2 of the License, or (at your option) any later version. 00008 * 00009 * This library is distributed in the hope that it will be useful, 00010 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 * Library General Public License for more details. 00013 * 00014 * You should have received a copy of the GNU Library General Public 00015 * License along with this library; if not, write to the 00016 * Free Software Foundation, Inc., 59 Temple Place - Suite 330, 00017 * Boston, MA 02111-1307, USA. 00018 */ 00019 00020 #ifndef GUMA_LINEQ_H 00021 #define GUMA_LINEQ_H 00022 00023 namespace guma { 00024 00025 using gul::Ptr; 00026 using gul::rtr; 00027 using gul::Square; 00028 using gul::point; 00029 00030 /*-----------------------------------------------------------------------*//** 00031 Decompose a matrix into lower and upper triangular matrix. the function 00032 throws an exception when an singular matrix is encountered (for 00033 convenience, so that not only the calling program can react to this, 00034 but also functions higher in the call stack) */ 00035 /*--------------------------------------------------------------------------*/ 00036 template< class T > 00037 GULAPI bool LUDecomposition( int nRowsCols, Ptr< Ptr<T> >& A, 00038 int *perm_sign, Ptr<int>& index ); 00039 00040 /*---------------------------------------------------------------*//** 00041 Solve system of linear equations with help of LU-Decomposition 00042 (given in A). */ 00043 /*------------------------------------------------------------------*/ 00044 template< class T > 00045 GULAPI void LUBackSubstitution( int nRowsCols, Ptr<int>& index, 00046 Ptr< Ptr<T> >& A, Ptr<T>& b ); 00047 00048 /*-------------------------------------------------------------------------*//** 00049 Multiply a band matrix and a vector: b = A*x, m1 = number of subdiagonal 00050 elements, m2 = number of superdiagonal elements, n = number of rows. 00051 It demonstrates the indexing scheme for band matrices */ 00052 /*----------------------------------------------------------------------------*/ 00053 template< class T > 00054 GULAPI void BandMVProduct( int n, int m1, int m2, 00055 Ptr< Ptr<T> > a, Ptr<T> x, Ptr<T> b); 00056 00057 /*-----------------------------------------------------------------------*//** 00058 Decompose a band matrix into lower and upper triangular matrix: A = L*U. 00059 m1 = number of subdiagonal elements, m2 = number of superdiagonal elements, 00060 n = number of rows and columns. L is returned in L[0..n-1][0..m1-1], 00061 U in the input matrix A; the rowwise permutation from pivoting is returned 00062 in 'index', its sign in 'perm_sign'. */ 00063 /*--------------------------------------------------------------------------*/ 00064 template< class T > 00065 GULAPI void BandDecomposition( 00066 int n, int m1, int m2, Ptr< Ptr<T> >& A, 00067 Ptr< Ptr<T> >& L, int *perm_sign, Ptr<int>& index ); 00068 00069 /*-----------------------------------------------------------------------*//** 00070 Solve system of linear equations with help of LU-Decomposition of 00071 a band matrix (given by L and A). */ 00072 /*--------------------------------------------------------------------------*/ 00073 template< class T > 00074 GULAPI void BandBackSubstitution( 00075 int n, int m1, int m2, const Ptr< Ptr<T> >& A, 00076 const Ptr< Ptr<T> >& L, const Ptr<int>& index, Ptr<T>& b ); 00077 00078 00079 /*----------------------------------------------------------------*//** 00080 Solve systems of linear equations with Gauss Jordan Elimination. 00081 00082 (A)nn * (b')nm = (b)nm 00083 (A)nn * (A')nn = (1)nn 00084 00085 After that, matrix 'A' contains the inverse matrix of 'A', 00086 and the columns of 'b' contain the solution vectors of the 00087 different right sides in 'b'. (see "Numerical Recipes in C") */ 00088 /*-------------------------------------------------------------------*/ 00089 template< class T > 00090 GULAPI bool GaussJordan( int nRowsCols, int nRightSides, 00091 Ptr< Ptr<T> >& A, Ptr< Ptr<T> > &b ); 00092 00093 00094 /*---------------------------------------------------------------------*//** 00095 Calculate (a^2 + b^2)^{1/2}, with precautions to avoid destructive 00096 over- or underflow. */ 00097 /*------------------------------------------------------------------------*/ 00098 template< class T > 00099 inline T Pythagoras( const T& a, const T& b ) 00100 { 00101 T absa,absb; 00102 00103 absa = rtr<T>::fabs(a); 00104 absb = rtr<T>::fabs(b); 00105 if( absa > absb ) 00106 return absa * rtr<T>::sqrt((T)1+Square(absb/absa)); 00107 00108 return absb == (T)0 ? (T)0 : absb * rtr<T>::sqrt((T)1+Square(absa/absb)); 00109 } 00110 00111 /*------------------------------------------------------------------------*//** 00112 Singular Value Decomposition. 00113 ============================= 00114 00115 A = U * W * V^{T} with 00116 00117 A = m x n Matrix (m >= n) 00118 U = m x n Matrix 00119 W = n x n Matrix 00120 V = n x n Matrix 00121 00122 U,V orthogonal ( i.e. U^{T} * U = (1), V^{T} * V = (1) ) 00123 W diagonal matrix, containing the singular values. */ 00124 /*---------------------------------------------------------------------------*/ 00125 template< class T > 00126 GULAPI void SVDecomposition( int m, int n, Ptr< Ptr<T> >& a, 00127 Ptr<T>& w, Ptr< Ptr<T> >& v ); 00128 00129 /*-------------------------------------------------------------------*//** 00130 Solve system of linear equations, when a Singular Value Decomposition 00131 is given */ 00132 /*----------------------------------------------------------------------*/ 00133 template< class T > 00134 GULAPI void SVBackSubstitution( int m, int n, 00135 Ptr< Ptr<T> >& u, Ptr<T>& w, 00136 Ptr< Ptr<T> >& v, Ptr<T>& b, 00137 Ptr<T>& x ); 00138 template< class T > 00139 GULAPI void SVBackSubstitution( int m, int n, 00140 Ptr< Ptr<T> >& u, Ptr<T>& w, 00141 Ptr< Ptr<T> >& v, Ptr<T>& b ); 00142 00143 /*---------------------------------------------------------------------*//** 00144 Sets very small singular values to zero. This avoids disastrous 00145 rounding errors, when the coefficient matrix of a system of 00146 linear equations is ill-conditioned. 00147 If some singular values were zero, or had to be set to zero, the function 00148 returns false, else it returns true */ 00149 /*------------------------------------------------------------------------*/ 00150 template< class T > 00151 GULAPI bool SVZero( int n, Ptr<T>& w ); 00152 00153 00154 /*-----------------------------------------------------------------------*//** 00155 Intersect two lines in 3-dimensions, using singular decomposition. 00156 (remark: this is a overkill, if you alredy know that the two lines 00157 intersect (but on the other hand this is the first oppurtinity 00158 to do something senseful with the above SVD algorithm:))))) 00159 00160 retcode = 0 => intersection point found 00161 retcode = 1 => lines are parallel 00162 retcode = 2 => lines are not parallel, but don't intersect either 00163 (what's the english word for "windschief" ?) 00164 In the last two cases the solution minimizes |x|^2 or |Ax + b|, 00165 and the midst of the two points on the two lines is returned. */ 00166 /*---------------------------------------------------------------------------*/ 00167 template< class T > 00168 int SVDIntersectLines( point<T> P1, point<T> T1, 00169 point<T> P2, point<T> T2, 00170 T *alpha1, T *alpha2, point<T> *P ); 00171 00172 } 00173 00174 #endif 00175 00176 00177 00178 00179 00180 00181 00182 00183 00184 00185 00186 00187 00188 00189 00190 00191 00192 00193 00194 00195 00196 00197 00198