Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

matrixbk.hh

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /************************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/matrixbk.hh,v $
00032  * $Revision: 3.11 $                                                    *
00033  * $Date: 2005/04/11 13:47:45 $                                         *
00034  * $Author: kampes $                                                    *
00035  *                                                                      *
00036  *  template (base) class for matrices.                                 *
00037  *                                                                      *
00038  *  ====== Defines (compiler -D option, see Makefile) ======            *
00039  *  __USE_FFTW_LIBRARY__ do use the fftw lib (compile separately)       *
00040  *                          prefered over veclib fft.                   *
00041  *  __USE_VECLIB_LIBRARY__ use VECLIB (FFT & matrix multiplication)     *
00042  *                      else internal routines will be used             *
00043  *  __USE_LAPACK_LIBRARY__ use LAPACK (cholesky routines)               *
00044  *                      else internal routines will be used             *
00045  *  __DEBUGMAT1         check index in matrix, dimensions, etc.         *
00046  *  __DEBUGMAT2         give info, can savely be un-defined             *
00047  *                                                                      *
00048  * Better make a vector baseclass (stl?) and a 2dmatrix class,          *
00049  *  which are friends of eachother,                                     *
00050  *  then derive an image class and add functions there.                 *
00051  * Note that data has to be continuously in memory, because             *
00052  *  VECLIB/FFTW assumes this (and me as well...) so a vector of vectors *
00053  *  cannot be used.                                                     *
00054  ************************************************************************/
00055 
00056 #ifndef MATRIXBK_H                              // guard
00057 #define MATRIXBK_H
00058 
00059 #include "constants.hh"                         // typedefs, window
00060 #include <fstream>                              // ofstream type
00061 
00062 
00063 //#ifdef __USE_FFTW_LIBRARY__
00064 //  #include <fftw3.h>                                  // fftw types
00065 //#endif
00066 
00067 // ______ message objects, global, set in main ______
00068 extern bk_messages matERROR;
00069 extern bk_messages matDEBUG;
00070 
00071 
00072 /************************************************************************
00073 friend functions are defined inside the class,
00074 since if they are declared here, and defined in matrixbk.cc
00075 then gcc gives errors like (for friend functions):
00076 In file included from processor.c:23:
00077 matrixbk.h:162: warning: friend declaration `void myswap(matrix<Type> &, matrix<Type> &)'
00078 matrixbk.h:162: warning:   declares a non-template function
00079 matrixbk.h:162: warning:   (if this is not what you intended, make sure
00080 matrixbk.h:162: warning:   the function template has already been declared,
00081 matrixbk.h:162: warning:   and add <> after the function name here)
00082 matrixbk.h:162: warning:   -Wno-non-template-friend disables this warning.
00083 
00084 Adding <> only makes the code less readable,
00085 cause other compilers don't want it.
00086 BK 17-Aug-2000
00087 *************************************************************************/
00088 
00089 
00090 
00091 // ====== Define template functions (no member no friend) ======
00092 // ______ (matrix class is declared way below) ______
00093 template <class Type> class matrix;
00094 
00095 #ifdef __USE_VECLIB_LIBRARY__
00096 // ______ See: matrixspecs.c ______
00097 // ______ Compiler uses specialization if declared here ______
00098 // ______ and not the template definition in matrixbk.cc ______
00099 matrix<real4>    operator *  (const matrix<real4>& A,   const matrix<real4>& B);
00100 matrix<real8>    operator *  (const matrix<real8>& A,   const matrix<real8>& B);
00101 matrix<complr4>  operator *  (const matrix<complr4>& A, const matrix<complr4>& B);
00102 matrix<complr8>  operator *  (const matrix<complr8>& A, const matrix<complr8>& B);
00103 matrix<real4>    matTxmat    (const matrix<real4> &A,   const matrix<real4> &B);
00104 matrix<real8>    matTxmat    (const matrix<real8> &A,   const matrix<real8> &B);
00105 matrix<complr4>  matTxmat    (const matrix<complr4> &A, const matrix<complr4> &B);
00106 matrix<complr8>  matTxmat    (const matrix<complr8> &A, const matrix<complr8> &B);
00107 matrix<real4>    matxmatT    (const matrix<real4> &A,   const matrix<real4> &B);
00108 matrix<real8>    matxmatT    (const matrix<real8> &A,   const matrix<real8> &B);
00109 matrix<complr4>  matxmatT    (const matrix<complr4> &A, const matrix<complr4> &B);
00110 matrix<complr8>  matxmatT    (const matrix<complr8> &A, const matrix<complr8> &B);
00111 #endif
00112 
00113 
00114 // ______ File: matrixspecs.c ______
00115 //#if defined (__DEBUGMAT2) || defined (__DEBUGMAT1)    // extra checking
00116 //void matDEBUG(char ch[ONE27]);
00117 //#endif
00118 //void matERROR(char ch[ONE27]);
00119 void matassert(
00120         const ofstream &str,
00121         const char* ofilename,
00122         const char* callingfilename = "?",
00123         int32 linenumber = 0);
00124 void matassert(
00125         const ifstream &str,
00126         const char* ifilename,
00127         const char* callingfilename = "?",
00128         int32 linenumber = 0);
00129 
00130 
00131 
00132 // ====== Without VECLIB these are also implemented (slower) ======
00133 // ______ NOW used is www.fftw.org routines ______
00134 // ______ matrixspecs.c: VECLIB [cz]1dfft or internal routine ______
00135 // ______ matrixspecs.c: FFTW  is used since 16-sep-2003 ______
00136   void            fft           (matrix<complr4> &A, int32 dimension);
00137   void            ifft          (matrix<complr4> &A, int32 dimension);
00138   void            fft2d         (matrix<complr4> &A);
00139   void            ifft2d        (matrix<complr4> &A);
00140   //void          fft2d         (matrix<real4> &A, matrix<real4> &B);
00141   //void          ifft2d        (matrix<real4> &A, matrix<real4> &B);
00142 
00143   matrix<complr4> oversample    (matrix<complr4> A,             // not by ref.
00144                                  uint frow, uint fcol);
00145   matrix<real4>   oversample    (const matrix<real4>   &A,
00146                                  uint frow, uint fcol);
00147 
00148 // ______ LAPACK_LIBRARY ______
00149 // ______ Without LAPACK, internal routines can be used (slower, inaccurate?) ______
00150   void            choles        (matrix<real4> &A);
00151   void            invertchol    (matrix<real4> &A);
00152   void            solvechol     (const matrix<real4> &A, matrix<real4> &B);
00153   void            choles        (matrix<real8> &A);
00154   void            invertchol    (matrix<real8> &A);
00155   void            solvechol     (const matrix<real8> &A, matrix<real8> &B);
00156 
00157 
00158   matrix<real4>   intensity     (const matrix<complr4> &A);
00159   matrix<real4>   magnitude     (const matrix<complr4> &A);
00160 
00161   matrix<real4>   real          (const matrix<complr4> &A);
00162   matrix<real4>   imag          (const matrix<complr4> &A);
00163 
00164   real4           norm2         (const matrix<complr4> &A);
00165   real4           norm2         (const matrix<real4> &A);
00166   matrix<complr4> norm          (const matrix<complr4> &A);
00167 
00168   matrix<real4>   abs           (const matrix<real4> &A);
00169   matrix<real8>   abs           (const matrix<real8> &A);
00170 
00171 // ______ Read file complex<short> into matrix complex<real4> ______
00172   void fileci2tomatcr4(
00173         matrix<complr4>         &Result,
00174         const char              *file, 
00175         uint                     filelines,
00176         window                   win,
00177         window                   winoffset);
00178 
00179 // ______ Casts ______
00180 // operator complr4 (matrix<real4>) seems better... (but how?)
00181 //  matrix<complr4> mat2cr4(
00182 //      const matrix<compli16>  &A);
00183   matrix<complr4> mat2cr4(
00184         const matrix<real4>     &A);
00185   matrix<complr4> mat2cr4(
00186         const matrix<real4>& A,
00187         const matrix<real4>& B);
00188   matrix<complr4> mat2cr4(
00189         const matrix<real8>& A,
00190         const matrix<real8>& B);
00191 
00192 // ______ phase ______
00193   matrix<real4>   angle(
00194         const matrix<complr4>   &A);    // phase        
00195   matrix<complr4> angle2cmplx(
00196         const matrix<real4>     &A);    // phasor, a=1
00197   void dotmultconjphase(
00198         matrix<complr4>         &complexinterferogram,//        by ref.
00199         const matrix<real4>     &refphase);     // phasor, a=1
00200 
00201 
00202 // ______ complex coherence ______
00203   matrix<complr4> coherence(
00204         const matrix<complr4>   &complex_interferogram,
00205         const matrix<complr4>   &norm_image1_and_2,
00206         uint                     estimatorwinsizeL,
00207         uint                     estimatorwinsizeP);
00208 
00209 // ______ complex coherence ______
00210   matrix<real4> coherence2(
00211         const matrix<complr4>   &complex_interferogram,
00212         const matrix<complr4>   &norm_image1_and_2,
00213         uint                     estimatorwinsizeL,
00214         uint                     estimatorwinsizeP);
00215 
00216 // should be in matrix class?
00217 // ______ sort rows of matrix on some column; uses qsort ______
00218 //    void mysort1(matrix<real4> &A);   // sort matrix based on col. number;
00219 // ______ (ascending) sort rows of matrix on first col, then second; uses qsort ______
00220   void mysort2(matrix<real4> &A);       // mysort(A);
00221   void mysort2(matrix<int32> &A);       // mysort(A);
00222 
00223 
00224 // ______ multiply strike x memory with values in B ______
00225 // ______ to multiply a column of A by vector B use:
00226 // ______ complr4 pntA=&A[0][c], strike=numpixels(A), length(B)=numrows.
00227   void dotmult (complr4 *startaddress, const matrix<real4> &B, int32 strike);
00228   //void dotmultrow (matrix &, row, matrix<real4> &B);
00229   //void dotmultcol (matrix &, col, const matrix<real4> &B);
00230   //matrix<complr4> dotmult(const matrix<complr4> &B, const matrix<real4> &A);
00231   //matrix<complr4> dotmult(const matrix<real4> &A,   const matrix<complr4> &B);
00232 
00233 // ______ Some trigonometric functions #%// BK 09-Nov-2000 ______
00234 // add lookup table for fast trig! (someday..)
00235 // Bert Kampes, 10-Apr-2005
00236   matrix<real4> cos (const matrix<real4> &A);
00237   matrix<real8> cos (const matrix<real8> &A);
00238   matrix<real4> sin (const matrix<real4> &A);
00239   matrix<real8> sin (const matrix<real8> &A);
00240 
00241 // ______ Some casts for complex (Re,Im) ______ HOW?
00242 //matrix<complr4> operator complr4  (const matrix<real4>& A);
00243 //matrix<complr4> operator complr4  (const matrix<real8>& A,   const matrix<real8>& B);
00244 //matrix<complr8> operator complr8  (const matrix<real8>& A,   const matrix<real8>& B);
00245 
00246 
00247 // ______ See: matrixbk.cc for definition (bottom) ______
00248 // ______ functions are no friends no members ______
00249 template <class Type>
00250   matrix<Type>  operator * (const matrix<Type>& A, const matrix<Type>& B);
00251 template <class Type>
00252   matrix<Type>  operator * (const matrix<Type>& A, Type scalar);
00253 template <class Type>
00254   matrix<Type>  operator *  (Type  scalar, const matrix<Type> &A);
00255 template <class Type>
00256   matrix<Type>  operator / (const matrix<Type>& A, Type B);
00257 template <class Type>
00258   matrix<Type>  operator / (const matrix<Type>& A, const matrix<Type>& B);
00259 template <class Type>
00260   matrix<Type>  operator - (const matrix<Type>& A, const matrix<Type>& B);
00261 template <class Type>
00262   matrix<Type>  operator - (const matrix<Type>& A, Type scalar);
00263 template <class Type>
00264   matrix<Type>  operator + (const matrix<Type>& A, const matrix<Type>& B);
00265 template <class Type>
00266   matrix<Type>  operator + (const matrix<Type>& A, Type B);
00267 template <class Type>
00268   Type          max(const matrix<Type> &A);
00269 template <class Type>
00270   Type          max(const matrix<Type> &A, uint& line, uint& pixel);
00271 template <class Type>
00272   Type          min(const matrix<Type> &A);
00273 template <class Type>
00274   Type          min(const matrix<Type> &A, uint& line, uint& pixel);
00275 template <class Type>
00276   matrix<Type>  matTxmat(const matrix<Type> &A, const matrix<Type> &B);
00277 template <class Type>
00278   matrix<Type>  matxmatT(const matrix<Type> &A, const matrix<Type> &B);
00279 template <class Type>
00280   void dumpasc(const char *file, const matrix<Type>& A);
00281 template <class Type>
00282   matrix<Type>  dotmult     (const matrix<Type> &A, const matrix<Type> &B);
00283 template <class Type>
00284   matrix<Type>  dotdiv      (const matrix<Type> &A, const matrix<Type> &B);
00285 template <class Type>
00286   matrix<Type>  sqr         (const matrix<Type> &A);
00287 template <class Type>
00288   matrix<Type>  conj        (const matrix<Type> &A);
00289 template <class Type>
00290   matrix<Type>  diagxmat    (const matrix<Type> &diag, const matrix<Type> &B);
00291 // ______ Specialisation R4*CR4 ______
00292 matrix<complr4> diagxmat    (const matrix<real4> &diag, const matrix<complr4> &B);
00293 template <class Type>
00294   matrix<Type>  multilook   (const matrix<Type> &A, uint factorL, uint factorP);
00295 template <class Type>
00296   matrix<real4> correlate   (const matrix<Type> &A, matrix<Type> Mask);
00297 template <class Type>
00298   matrix<Type>  operator -  (const matrix<Type>& A);
00299 template <class Type>
00300   real8         mean        (const matrix<Type> &A);
00301 template <class Type>
00302   matrix<Type>  sum         (const matrix<Type> &A, int32 dim);
00303 //template <class Type>
00304 //  matrix<Type>  operator ^ (const matrix<Type>& A);
00305 //template <class Type>
00306 //  matrix<Type>  operator ^ (Type scalar);
00307 
00308 
00309 
00310 
00311 
00312 
00313 // ====== ====== ====== ====== ====== ======
00314 // ====== Start of class declaration: ======
00315 // ====== ====== ====== ====== ====== ======
00316 template <class Type>                           // class template
00317 class matrix
00318   {
00319   private:
00320     // ______ Private data ______
00321     Type **data;                                // two dimensional array
00322     uint nrows;                                 // number of rows (lines)
00323     uint ncols;                                 // number of rows (pixels)
00324     uint nsize;                                 // lines*pixels
00325 
00326     // ______ Private functions ______
00327     void   allocate                                         (uint l, uint p);
00328     void   initialize                                       (uint l, uint p);
00329 
00330     #ifdef __DEBUGMAT1
00331     void   checkindex                                       (uint l, uint p) const;
00332     #endif
00333 
00334   public:
00335     // ______ Constructors ______
00336     matrix                      ();
00337     matrix                      (uint l, uint p);
00338     matrix                      (const matrix<Type>& A);
00339     matrix                      (window w, const matrix<Type>& A);
00340 
00341     // ______ Destructor ______
00342     ~matrix                     ();
00343 
00344     // ______ Data functions ______
00345     void          setdata       (Type w);
00346     void          setdata       (uint l, uint p, const matrix<Type>& A);
00347     void          setdata       (window winin, const matrix<Type> &A, window winA);
00348     void          setdata       (const matrix<Type> &A, window winA);
00349     matrix<Type>  getdata       (window w)                      const;
00350     matrix<Type>  getrow        (uint l)                        const;
00351     matrix<Type>  getcolumn     (uint p)                        const;
00352     // later... index matrix, mask, ... functies findlt findgt, etc.
00353     // later,,, void operator ()   (indexmatrix, mask, or whatever)
00354     // BK 23-Oct-2000
00355     //matrix<int32> find          (Type w)                              const;
00356     void          showdata      ()                              const;
00357     bool          isvector      ()                              const;
00358     uint          lines         ()                              const;
00359     uint          pixels        ()                              const;
00360     uint          size          ()                              const;
00361     void          resize        (uint l, uint p);
00362     void          clean         ();
00363 
00364     void          setrow        (uint l, const matrix<Type> &L);
00365     void          setrow        (uint l, Type scalar);
00366     void          setcolumn     (uint p, const matrix<Type> &C);
00367     void          setcolumn     (uint p, Type scalar);
00368     void          fliplr        ();
00369     void          flipud        ();
00370 
00371     Type*         operator []   (uint l)                        const;
00372     Type&         operator ()   (uint l, uint p)                const;
00373     matrix<Type>  operator ()   (window win)                    const;
00374     matrix<Type>  operator ()   (uint l0,uint lN,uint p0,uint pN) const;
00375 
00376     matrix<Type>& operator =    (const matrix<Type>& A);
00377     matrix<Type>& operator =    (const Type scalar);
00378     matrix<Type>& operator -=   (Type scalar);
00379     matrix<Type>& operator -=   (const matrix<Type>& A);
00380     matrix<Type>& operator +=   (Type scalar);
00381     matrix<Type>& operator +=   (const matrix<Type>& A);
00382     matrix<Type>& operator *=   (Type scalar);
00383     matrix<Type>& operator *=   (const matrix<Type> &A);
00384     matrix<Type>& operator /=   (Type scalar);
00385     matrix<Type>& operator /=   (const matrix<Type> &A);
00386     bool          operator ==   (Type scalar)                   const;
00387     bool          operator ==   (const matrix<Type> &A)         const;
00388     bool          operator !=   (Type scalar)                   const;
00389     bool          operator !=   (const matrix<Type> &A)         const;
00390     void          conj          ();
00391     void          mypow         (Type s);               // better use operator ^
00392     // ______ Functions to, e.g., multiply a CR4 by a R4 matrix ______
00393     // ______ Declare for all, but only define for a few in matrixspecs.c ______
00394     // ______ Hopefully the template Type,Type is not? overwritten, is same ______
00395     template <class TypeB>
00396       matrix<Type>& operator *=   (const matrix<TypeB> &A);
00397     template <class TypeB>
00398       matrix<Type>& operator /=   (const matrix<TypeB> &A);
00399     template <class TypeB>
00400       matrix<Type>& operator +=   (const matrix<TypeB> &A);
00401     template <class TypeB>
00402       matrix<Type>& operator -=   (const matrix<TypeB> &A);
00403 
00404 
00405 
00406 // ================
00407 //  Template friend functions have to be defined within class
00408 //  (where they are declared) to satisfy with g++ preference (why??)
00409 // ================
00410 
00411 
00412 /****************************************************************
00413  * file << A;                                                   *
00414  * Bert Kampes, 14-Jan-1999                                     *
00415  ****************************************************************/
00416 friend ostream&      operator << (ostream& file, const matrix<Type>& A)
00417   {
00418   const uint sizeofType = sizeof(Type);
00419   for (register int32 i=0;i<A.nrows;i++)
00420     for (register int32 j=0;j<A.ncols;j++)
00421       file.write((char*)&A.data[i][j],sizeofType);
00422   return file;
00423   } // END <<
00424 
00425 
00426 
00427 /****************************************************************
00428  * file >> A;                                                   *
00429  * Bert Kampes, 14-Jan-1999                                     *
00430  ****************************************************************/
00431 //friend istream&      operator >> (istream& file, const matrix<Type>& A)
00432 // BK 25-Sep-2000
00433 friend istream&      operator >> (istream& file, matrix<Type>& A)
00434   {
00435   const uint sizeofType = sizeof(Type);
00436   for (register int32 i=0;i<A.nrows;i++)
00437     for (register int32 j=0;j<A.ncols;j++)
00438       file.read((char*)&A.data[i][j],sizeofType);
00439   // too slow:  
00440   // file.read((char*)&A.data[0][0],A.nrows*A.ncols*sizeof(Type));
00441   return file;
00442   } // END >>
00443 
00444 
00445  
00446 /****************************************************************
00447  * myswap(A,B)                                                  *
00448  *    interchange matrices of same size                         *
00449  *    somehow g++ does not like name swap, so myswap            *
00450  *    Bert Kampes, 14-Apr-1999                                  *
00451  ****************************************************************/
00452 friend void          myswap      (matrix<Type> &A, matrix<Type> &B)
00453   {
00454   #ifdef __DEBUGMAT2
00455     matDEBUG.print("myswap.");
00456   #endif
00457   #ifdef __DEBUGMAT1
00458   if (A.nrows  != B.nrows ||
00459       A.ncols != B.ncols  )
00460     matERROR.print("swap: matrices must be same size (for now).");
00461   #endif
00462 
00463   Type **pntA = A.data;
00464   Type **pntB = B.data;
00465   A.data = pntB;
00466   B.data = pntA;
00467   } // END myswap
00468 
00469 
00470 
00471 /****************************************************************
00472  * A = sqrt(B)                                                  *
00473  * Bert Kampes, 16-Feb-1999                                     *
00474  ****************************************************************/
00475 friend matrix<Type>  sqrt        (const matrix<Type> &A)
00476   {
00477   #ifdef __DEBUGMAT2
00478   matDEBUG.print("sqrt");
00479   #endif
00480 
00481   matrix<Type> Result(A.lines(),A.pixels());
00482   // ______Ensure stride one memory______
00483   Type *pntA   = A.data[0];
00484   Type *pntR   = Result.data[0];
00485   for (register int32 i=0; i<Result.nsize; i++)
00486     (*pntR++) = (sqrt((*pntA++)));
00487   return Result;
00488   } // END sqrt
00489 
00490 
00491 
00492 /****************************************************************
00493  * readfile                                                     *
00494  *    Type must be same as on disk. Result is filled.           *
00495  *    Windows must be specified in same system.                 *
00496  *    Either use win[0:N] and offset[0:N] (where 0 is the first *
00497  *    line), or use win[1:N] and offset[1:N] (where 1 indicates *
00498  *    the first line.                                           *
00499  *    Since the information in the resultfiles is written       *
00500  *    starting at line 1 (not 0) a 'currentwindow' can be read  *
00501  *    by a call like:                                           *
00502  *    readfile(mat,name,lines,winstartingatline1,currentwindow);*
00503  *    Result is filled at R(0:L,0:P), L=win.lines.              *
00504  * Bert Kampes, 14-Jan-1999                                     *
00505  ****************************************************************/
00506 friend void readfile(matrix<Type> &Result, const char *file,
00507   uint filelines, window win, window winoffset)
00508   {
00509   #ifdef __DEBUGMAT2
00510   matDEBUG.print("readfile (window).");
00511   #endif
00512 
00513 // ______First account for possible offsets of file______
00514 // BK 18/1/00: winoffset starts at line 1?
00515   win.linelo = win.linelo - winoffset.linelo + 1;
00516   win.linehi = win.linehi - winoffset.linelo + 1;
00517   win.pixlo  = win.pixlo  - winoffset.pixlo  + 1;
00518   win.pixhi  = win.pixhi  - winoffset.pixlo  + 1;
00519   //win.pixhi  -= (winoffset.pixlo - 1);
00520 
00521   // ______ Check input ______
00522   // ifstream ifile(file, ios::in | ios::ate | ios::binary);
00523   // g++ seems to have a prob. with this... BK 130700
00524   // ifstream ifile(file, ios::in | ios::app | ios::binary);
00525   //ifstream ifile(file, ios::in | ios::binary | ios::nocreate);
00526   // new compiler v3.2 handles nocreate automatically fine 
00527 #ifdef __NO_IOS_BINARY__
00528   ifstream ifile(file, ios::in);
00529 #else
00530   ifstream ifile(file, ios::in | ios::binary);
00531 #endif
00532   matassert(ifile,file,__FILE__,__LINE__);
00533   ifile.seekg(0,ios::end);                              // pointer to end...
00534   const uint sizefile     = ifile.tellg();              // opened ate
00535   const uint pixelsxbytes = sizefile/filelines;
00536   const uint filepixels   = pixelsxbytes/sizeof(Type);  // Type on disk.
00537   const uint sizepixel    = pixelsxbytes/filepixels;
00538 
00539   #ifdef __DEBUGMAT2
00540   if (win.linelo<=0 || win.pixlo<=0)
00541     matERROR.print("minimum line(pixel) is 0, (should be 1?).");
00542   if (win.linelo>win.linehi || win.pixlo>win.pixhi)
00543     matERROR.print("minimum line(pixel) is larger than max.");
00544   if (win.linehi>filelines || win.pixhi>filepixels)
00545     matERROR.print("max. line (pixel) is larger then on file.");
00546   if (sizeof(Type)!=sizepixel)
00547     matERROR.print("Type on disk is different than type of matrix.");
00548   if (sizefile==0)
00549     matERROR.print("filesize==0...");
00550   #endif
00551 
00552   const uint lines  = win.lines();
00553   const uint pixels = win.pixels();
00554   const uint start  = ((win.linelo-1)*filepixels+win.pixlo-1)*sizepixel;
00555 
00556   #ifdef __DEBUGMAT1
00557   if (Result.lines() < lines || Result.pixels() < pixels)
00558     matERROR.print("readfile: matrix too small to contain window from file.");
00559   if (Result.lines() != lines || Result.pixels() != pixels)
00560     matDEBUG.print("debug info: readfile: matrix is not fully filled.");
00561   #endif
00562   for (register int32 lin=0; lin<lines; ++lin)
00563     {
00564     // read data at row: lin
00565     ifile.seekg(start+filepixels*lin*sizepixel,ios::beg);
00566     ifile.read((char*)&Result.data[lin][0],pixels*sizepixel);
00567     }
00568   ifile.close();
00569   } // END readfile
00570 
00571 
00572 
00573 /****************************************************************
00574  * writefile                                                    *
00575  * write matrix partially to file.                              *
00576  * win 0 = first line (matrix index)                            *
00577  * Bert Kampes, 14-Jan-1999                                     *
00578  ****************************************************************/
00579 friend void          writefile   (
00580   ofstream &file, const matrix<Type> &tobwritten, window win)
00581   {
00582   #ifdef __DEBUGMAT2
00583   matDEBUG.print("writefile.");
00584   #endif
00585   #ifdef __DEBUGMAT1
00586   if (win.linelo>win.linehi || win.pixlo>win.pixhi)
00587     matERROR.print("minimum line(pixel) is larger than max.");
00588   if (win.linehi>=tobwritten.lines() || win.pixhi>=tobwritten.pixels())
00589     matERROR.print("window not ok with matrix.");
00590   #endif
00591   matassert(file,"writefile: file unknown",__FILE__,__LINE__);
00592   const uint pixels = win.pixels();
00593   const uint size   = pixels*sizeof(Type);
00594   for (int32 line=win.linelo; line<=win.linehi; ++line)
00595     file.write((char*)&tobwritten.data[line][win.pixlo],size);
00596   } // END writefile
00597 
00598 
00599 
00600 /****************************************************************
00601  * fftshift(A)                                                  *
00602  *    fftshift of vector A is returned in A by reference        *
00603  *    shift DC term to middle. p=ceil(m/2); A=A[p:m-1 0:p-1];   *
00604  *    Bert Kampes, 24-Mar-2000                                  *
00605  ****************************************************************/
00606 friend void          fftshift    (matrix<Type> &A)
00607   {
00608   #ifdef __DEBUGMAT2
00609   matDEBUG.print("fftshift");
00610   #endif
00611   #ifdef __DEBUGMAT1
00612   if (!A.isvector())
00613     matERROR.print("fftshift: only vectors");
00614   #endif
00615   matrix<Type> Res(A.nrows,A.ncols);
00616   const int32 start  = int32(ceil(real8(A.nsize)/2));
00617   memcpy(Res.data[0],A.data[0]+start,sizeof(Type)*(A.nsize-start));
00618   memcpy(Res.data[0]+A.nsize-start,A.data[0],sizeof(Type)*start);
00619   myswap(A,Res);                // prevent copy
00620   } // END fftshift
00621 
00622 
00623 
00624 /****************************************************************
00625  * ifftshift(A)                                                 *
00626  *    ifftshift of vector A is returned in A by reference       *
00627  *    undo effect of fftshift. ?p=floor(m/2); A=A[p:m-1 0:p-1]; *
00628  *    Bert Kampes, 24-Mar-2000                                  *
00629  ****************************************************************/
00630 friend void          ifftshift   (matrix<Type> &A)
00631   {
00632   #ifdef __DEBUGMAT2
00633   matDEBUG.print("ifftshift");
00634   #endif
00635   #ifdef __DEBUGMAT1
00636   if (!A.isvector())
00637     matERROR.print("ifftshift: only vectors");
00638   #endif
00639   matrix<Type> Res(A.nrows,A.ncols);
00640   //const int32 length = A.size();
00641   const int32 start  = int32(floor(real8(A.nsize)/2));
00642   memcpy(Res.data[0],A.data[0]+start,sizeof(Type)*(A.nsize-start));
00643   memcpy(Res.data[0]+A.nsize-start,A.data[0],sizeof(Type)*start);
00644   myswap(A,Res);
00645   } // END ifftshift
00646 
00647 
00648 
00649 /****************************************************************
00650  * wshift(A,n)                                                  *
00651  *    circular shift of vector A by n pixels. positive n for    *
00652  *    right to left shift.                                      *
00653  *    implementation: WSHIFT(A,n) == WSHIFT(A,n-sizeA);         *
00654  *    A is changed itself!                                      *
00655  *    Bert Kampes, 02-Nov-2000                                  *
00656  ****************************************************************/
00657 friend void          wshift       (matrix<Type> &A, int32 n)
00658   {
00659   #ifdef __DEBUGMAT2
00660   matDEBUG.print("wshift");
00661   #endif
00662   #ifdef __DEBUGMAT1
00663   if (abs(n)>=A.nsize)
00664     matERROR.print("wshift: shift larger than matrix not implemented.");
00665   if (!A.isvector())
00666     matERROR.print("wshift: only vectors");
00667   #endif
00668   // positive only, use rem!  n = n%A.nsize;
00669   if (n==0) return;
00670   if (n<0)  n += A.nsize;
00671   matrix<Type> Res(A.nrows,A.ncols);
00672   // ______ n always >0 here ______
00673   memcpy(Res.data[0],A.data[0]+n,sizeof(Type)*(A.nsize-n));
00674   memcpy(Res.data[0]+A.nsize-n,A.data[0],sizeof(Type)*n);
00675   myswap(A,Res);                // prevent copy
00676   } // END wshift
00677 
00678 
00679   }; // END matrix class
00680 
00681 
00682 
00683 
00684 
00685 
00686 // ______ Compilation with g++ seems impossible any other way? ______
00687 #include "matrixbk.cc"
00688 
00689 #endif // MATRIXBK_H guard
00690 

Generated on Fri Apr 22 15:57:55 2005 for Doris by doxygen 1.3.6