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

matrixbk.cc

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.cc,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) ======                          *
00039  *  __USE_VECLIB_LIBRARY__      for using VECLIB library                *
00040  *  __USE_LAPACK_LIBRARY__      for using LAPACK library                *
00041  *  __USE_FFTW_LIBRARY__        for using FFTW library                  *
00042  *  __DEBUGMAT1         index checking, dimensions etc.                 *
00043  *  __DEBUGMAT2         info, can savely be un-defined                  *
00044  *                                                                      *
00045  * Better make a vector baseclass (stl?) and a 2dmatrix class,          *
00046  *  which are friends of eachother,                                     *
00047  *  then derive an image class and add functions there.                 *
00048  * Note that data has to be continuously in memory, because             *
00049  *  VECLIB/FFTW assumes this (and me as well...) so a vector of vectors *
00050  *  cannot be used.                                                     *
00051  ************************************************************************/
00052 
00053 #include "constants.hh"                 // typedefs, window
00054 
00055 #include <iostream>                     // cout etc.
00056 #include <fstream>                      // ofstream type
00057 #include <strstream>                    // memory stream
00058 #include <iomanip>                      // setw etc.
00059 #include <algorithm>                    // max
00060 #include <cstring>                      // memset according to g++ (?)
00061 #include <complex>
00062 #ifdef __DEBUGMAT1                      // use index checking, alloc
00063   #include <new>                        // bad_alloc
00064 #endif
00065 
00066 // ______ Keep track of total allocated memory ______
00067 // ______ This will work upto 2^31 bytes (2GB) ______
00068 // ______ there is a problem here, but how to declare a global? ______
00069 // ______ if called from different files, not correct bookkeeping ______
00070 // ______ therefor we define it above main in processor.cc ______
00071 #ifdef __DEBUGMAT2                              // use index checking, alloc
00072   extern uint totalallocated;                   // [B]
00073 #endif
00074 
00075 // ______ message objects, global, set in main ______
00076 extern bk_messages matERROR;
00077 extern bk_messages matDEBUG;
00078 
00079 
00080 
00081 
00082 // ====== Start of class definition ======
00083 // ====== Private functions ======
00084 /****************************************************************
00085  * allocate                                                     *
00086  *    allocate 2d array in major row order, continuously in     *
00087  *    memory (required by VECLIB/FFTW). consider rewriting to vector    *
00088  *    of vectors of stl, but may not be cont. in mem.           *
00089  *    g++ w/o exception handling, ...                           *
00090  *    Bert Kampes, 01-Feb-1999                                  *
00091  ****************************************************************/
00092 template <class Type>
00093 void matrix<Type>::allocate(uint numlines, uint numpixels)      // allocator
00094   {
00095 #ifdef __DEBUGMAT1
00096   if (numlines==0 || numpixels==0)
00097     {
00098     matERROR << "Allocation impossible: size (l,p): "
00099          << numlines << ", " << numpixels;
00100     matERROR.print();
00101     }
00102 #endif
00103 
00104   nrows = numlines;
00105   ncols = numpixels;
00106   nsize = numlines*numpixels;
00107 
00108 // Bert Kampes, 07-Apr-2005: try/catch should work by now...
00109 //#if defined (__DEBUGMAT1) && !defined (__GNUC__) // try not ok with g++ ??
00110   try
00111     {
00112 //#endif
00113   data = new Type*[numlines];                   // get memory
00114 //#if defined (__DEBUGMAT1) && !defined (__GNUC__) // try not ok with g++ ??
00115     }
00116   catch (bad_alloc)
00117     {
00118     matERROR << "code 502: allocation failed: 1size: "
00119          << numlines << ", " << numpixels;
00120     matERROR.print();
00121     }
00122   try
00123     {
00124 //#endif
00125   data[0] = new Type[nsize];                    // get memory
00126 //#if defined (__DEBUGMAT1) && !defined (__GNUC__) // try not ok with g++ ??
00127     }
00128   catch(bad_alloc)
00129     {
00130     matERROR << "code 502: allocation failed: 2size: "
00131          << numlines << ", " << numpixels;
00132     matERROR.print();
00133     }
00134 //#endif
00135   for (register int32 i=1; i<int32(numlines); i++)
00136     data[i] = data[i-1]+numpixels;              // start at 0,0
00137 
00138 #ifdef __DEBUGMAT2
00139   // --- do this in main program? ---
00140   //uint totalallocated=0;              // [B]
00141   uint allocated = sizeof(Type)*nsize;         // [B]
00142   totalallocated += allocated;                  // [B]
00143   matDEBUG << "allocated   matrix("
00144        << numlines << "," << numpixels << ") at: " << &data[0][0] << " ("
00145        << setw(10) << allocated << "B, total: "
00146        << setw(10) << totalallocated << "B)";
00147   matDEBUG.print();
00148 #endif
00149   } // END allocate
00150 
00151 
00152 
00153 /****************************************************************
00154  * initialize = allocate + set 0                                *
00155  *    Bert Kampes, 01-Feb-1999                                  *
00156  ****************************************************************/
00157 template <class Type>
00158 void matrix<Type>::initialize(uint numlines, uint numpixels)
00159   {
00160   #ifdef __DEBUGMAT2
00161   matDEBUG.print("initialization matrix.");
00162   #endif
00163   allocate(numlines,numpixels);
00164   clean();
00165   } // END initialize
00166 
00167 
00168 
00169 #ifdef __DEBUGMAT1
00170 /****************************************************************
00171  * checkindex (l,p)                                             *
00172  *    Bert Kampes, 01-Feb-1999                                  *
00173  ****************************************************************/
00174 template <class Type>
00175 void matrix<Type>::checkindex(uint line, uint pixel) const
00176   {
00177   if (line > nrows-1 || pixel > ncols-1)
00178     {
00179     matERROR << "Wrong index (l,p)=(" << line << "," << pixel
00180          << "); Matrix(" << nrows << "," << ncols
00181          << ") at " << &data[0][0];
00182     matERROR.print();
00183     }
00184   } // END checkindex
00185 #endif
00186 
00187 
00188 
00189 // ====== Public functions ======
00190 // ====== Constructors ======
00191 /****************************************************************  \\
00192  * matrix<real8> A;                                             *  \\
00193  * Bert Kampes, 11-Dec-1998                                     *  \\
00194  ****************************************************************/
00195 template <class Type>
00196 matrix<Type>::matrix()                          // constructor (0 arg)
00197   {
00198   nrows = 0;
00199   ncols = 0;
00200   nsize = 0;
00201   data  = 0;
00202   } // END constructor
00203 
00204 
00205 
00206 /****************************************************************
00207  * matrix<real8> A(3,3);                                        *
00208  * Bert Kampes, 11-Dec-1998                                     *
00209  ****************************************************************/
00210 template <class Type>
00211 matrix<Type>::matrix(uint lines, uint pixels)
00212   {
00213   initialize(lines,pixels);                             // set to 0
00214   } // END constructor
00215 
00216 
00217 
00218 /****************************************************************
00219  * matrix<real8> A=B;                                           *
00220  * copy constructor; avoids default for dynamical memory:       * 
00221  * bitwise copy.                                                *
00222  * Bert Kampes, 11-Dec-1998                                     *
00223  ****************************************************************/
00224 template <class Type>
00225 matrix<Type>::matrix(const matrix<Type>& A)
00226   {
00227   #ifdef __DEBUGMAT2
00228   matDEBUG.print("copy constructor.");
00229   #endif
00230   if (A.nsize)
00231     {
00232     allocate(A.nrows,A.ncols);
00233     memcpy(data[0],A.data[0],nsize*sizeof(Type));
00234     }
00235   } // END constructor
00236 
00237 
00238 
00239 /****************************************************************
00240  * matrix<int32> A(win,B)                                       *
00241  * Constructor (part)                                           *
00242  *    Bert Kampes, 01-Feb-1999                                  *
00243  *    Bert Kampes, 31-Mar-1999 memcpy ipv. for_j                *
00244  * window starts at 0 (matrix index)                            *
00245  #%// BK 25-Oct-2000                                            *
00246  ****************************************************************/
00247 template <class Type>
00248 matrix<Type>::matrix (window win, const matrix<Type>& A)
00249   {
00250   #ifdef __DEBUGMAT2
00251   matDEBUG.print("constructor as part.");
00252   #endif
00253   // ______ Check arguments ______
00254   #ifdef __DEBUGMAT1
00255   if (win.linehi<win.linelo)
00256     matERROR.print("constructor (4uint,matrix): win.linehi<linelo ?");
00257   if (win.pixhi<win.pixlo)
00258     matERROR.print("constructor (4uint,matrix): win.pixhi<pixlo ?");
00259   A.checkindex(win.linehi,win.pixhi);
00260   #endif
00261   // ______ Allocate new matrix and fill ______
00262   const uint numlin  = win.lines();
00263   const uint numpix  = win.pixels();
00264   const uint sizelin = numpix*sizeof(Type);
00265   allocate(numlin,numpix);
00266   for(register int32 i=0; i<numlin; i++)
00267     memcpy(data[i],A[win.linelo+i]+win.pixlo,sizelin);
00268   } // END constructor
00269 
00270 
00271 
00272 // ======Destructor======
00273 /****************************************************************
00274  * Destructor                                                   *
00275  *    Bert Kampes, 01-Feb-1999                                  *
00276  ****************************************************************/
00277 template <class Type>
00278 matrix<Type>::~matrix()
00279   {
00280   if (data==0) return;
00281   #ifdef __DEBUGMAT2
00282   uint deallocated  = sizeof(Type)*nsize;       // [B]
00283   totalallocated   -= deallocated;              // [B]
00284   matDEBUG << "deallocated matrix("
00285        << nrows << "," << ncols << ") at: " << data[0] << " ("
00286        << setw(10) << deallocated    << "B; total: "
00287        << setw(10) << totalallocated << "B)";
00288   matDEBUG.print();
00289   #endif
00290 
00291   delete [] data[0];                    // deallocate
00292   delete [] data;                       // deallocate
00293   data=0;                               // set to null pointer
00294   } // END destructor
00295 
00296 
00297 
00298 // ======Data functions======
00299 /****************************************************************
00300  * A.setdata(w)                                                 *
00301  *    Bert Kampes, 01-Feb-1999                                  *
00302  ****************************************************************/
00303 template <class Type>
00304 void matrix<Type>::setdata(Type w)      // sets 2 w
00305   {
00306   #ifdef __DEBUGMAT2
00307   matDEBUG.print("A.setdata(w)");
00308   #endif
00309   //  memset(data[0],w,nsize*sizeof(Type));
00310   Type *pnt = data[0];
00311   for (register int32 i=0;i<nsize;i++)
00312     (*pnt++) = Type(w);
00313 
00314   // ______ not faster, and more difficult to read ______
00315   //  for (register Type *pnt=&data[0][0]; 
00316   //     pnt<=&data[nrows-1][ncols-1];
00317   //     ++*pnt++=Type(w)); // do nothing
00318   } // END setdata
00319 
00320 
00321 
00322 /****************************************************************
00323  * setdata(i,j,A)                                               *
00324  *    put matrix A on l1,p1                                     *
00325  *    Bert Kampes, 01-Feb-1999                                  *
00326  ****************************************************************/
00327 template <class Type>
00328 void matrix<Type>::setdata(uint l1, uint p1, const matrix<Type>& A)
00329   {
00330   #ifdef __DEBUGMAT1
00331   checkindex(l1+A.nrows-1,p1+A.ncols-1);        // far most corner 
00332   #endif
00333   const uint sizelin = A.ncols*sizeof(Type);
00334   for (register int32 i=0;i<A.nrows;i++)
00335     memcpy(data[i+l1]+p1,A[i],sizelin);
00336   } // END setdata
00337 
00338 
00339 
00340 /****************************************************************
00341  * B.setdata(win, A, winA):                                     *
00342  *  set win of B to winA of A                                   *
00343  * if win==0 defaults to totalB, winA==0 defaults to totalA     *
00344  * first line matrix =0 (?)
00345  *    Bert Kampes, 01-Feb-1999                                  *
00346  ****************************************************************/
00347 template <class Type>
00348 void matrix<Type>::setdata(window winin, const matrix<Type> &A, window winA)
00349   {
00350   #ifdef __DEBUGMAT2
00351   matDEBUG.print("setdata (win,A,win)");
00352   #endif
00353 
00354   // ______Check default request______
00355   if (winin.linehi == 0 && winin.pixhi == 0)
00356     {winin.linehi = nrows  -1;
00357      winin.pixhi  = ncols -1;}
00358   if  (winA.linehi == 0 &&  winA.pixhi == 0)
00359     {winA.linehi  = A.lines()  -1;
00360      winA.pixhi   = A.pixels() -1;}
00361 
00362   #ifdef __DEBUGMAT1
00363   if (((winin.linehi - winin.linelo) != (winA.linehi - winA.linelo)) ||
00364        ((winin.pixhi  - winin.pixlo) != (winA.pixhi  - winA.pixlo))    )
00365     matERROR.print("code 901: wrong input.");
00366   if (winin.linehi<winin.linelo || winin.pixhi<winin.pixlo)
00367     matERROR.print("code 901: wrong input.1");
00368   if ((winin.linehi > nrows  -1) ||
00369       (winin.pixhi  > ncols -1)   )
00370     matERROR.print("code 901: wrong input.2");
00371   if ((winA.linehi > A.lines()  -1) ||
00372       (winA.pixhi  > A.pixels() -1)   )
00373     matERROR.print("code 901: wrong input.3");
00374   #endif
00375 
00376   // ______ Fill data ______
00377   const uint sizelin = (winA.pixhi - winA.pixlo + 1)*sizeof(Type);
00378   for(register int32 i=winin.linelo; i<=winin.linehi;i++)
00379     memcpy(data[i]+winin.pixlo,A[i-winin.linelo+winA.linelo]+winA.pixlo,sizelin);
00380   } // END setdata
00381 
00382 
00383 
00384 /****************************************************************
00385  * B.setdata(A, winA):                                          *
00386  *    set total of B to winA of A                               *
00387  *    Bert Kampes, 17-Mar-1999                                  *
00388  ****************************************************************/
00389 template <class Type>
00390 void matrix<Type>::setdata(const matrix<Type> &A, window winA)
00391   {
00392   #ifdef __DEBUGMAT2
00393   matDEBUG.print("setdata (A,win)");
00394   #endif
00395   #ifdef __DEBUGMAT1
00396   if (((nrows  -1) != (winA.linehi - winA.linelo)) ||
00397        ((ncols -1) != (winA.pixhi  - winA.pixlo))    )
00398     matERROR.print("code 901: wrong input.");
00399   if ((winA.linehi > A.lines()  -1) ||
00400       (winA.pixhi  > A.pixels() -1)   )
00401     matERROR.print("code 901: wrong input.3");
00402   #endif
00403   
00404   // ______ Fill data ______
00405   const uint sizelin = ncols*sizeof(Type);
00406   for(register int32 i=0; i<nrows;i++)
00407     memcpy(data[i],A[i+winA.linelo]+winA.pixlo,sizelin);
00408   } // END setdata
00409 
00410 
00411 
00412 
00413 /****************************************************************
00414  * A=B.getdata(win)                                             *
00415  *    Bert Kampes, 01-Feb-1999                                  *
00416  ****************************************************************/
00417 template <class Type>
00418 matrix<Type> matrix<Type>::getdata(window win) const
00419   {
00420   #ifdef __DEBUGMAT2
00421   matDEBUG.print("getdata.");
00422   #endif
00423   #ifdef __DEBUGMAT1
00424   if (win.linehi<win.linelo || win.pixhi<win.pixlo)
00425     matERROR.print(
00426     "code 501: matrix::getdata (win): arguments are wrong, l1<l2,p1<p2");
00427   checkindex(win.linehi,win.pixhi);
00428   #endif
00429   //  uint numlin = win.lines();
00430   //  uint numpix = win.pixels();
00431   uint numlin = win.linehi - win.linelo + 1;
00432   uint numpix = win.pixhi -  win.pixlo  + 1;
00433   matrix<Type> Result(numlin,numpix);                   // =data(;
00434   for(register int32 i=0; i<numlin; i++)
00435     memcpy(Result[i],data[i+win.linelo]+win.pixlo,numpix*sizeof(Type));
00436   return Result;
00437   } // END getdata
00438 
00439 
00440 /****************************************************************
00441  * A=B.getrow(row)                                              *
00442  * rows: 0 1 2 ..                                               *       
00443  *    Bert Kampes, 01-Feb-1999                                  *
00444  ****************************************************************/
00445 template <class Type>
00446 matrix<Type> matrix<Type>::getrow(uint line) const
00447   {
00448   #ifdef __DEBUGMAT1
00449   checkindex(line,0);
00450   #endif
00451   matrix<Type> Result(1,ncols);
00452   memcpy(Result[0],data[line],ncols*sizeof(Type));
00453   return Result;
00454   } // END getrow
00455 
00456 
00457 /****************************************************************
00458  * A=B.getcolumn(col)                                           *
00459  * cols: 0 1 2 ..                                               *       
00460  *    Bert Kampes, 01-Feb-1999                                  *
00461  ****************************************************************/
00462 template <class Type>
00463 matrix<Type> matrix<Type>::getcolumn(uint pixel) const
00464   {
00465   #ifdef __DEBUGMAT1
00466     checkindex(0,pixel);
00467   #endif
00468   matrix<Type> Result(nrows,1);
00469 
00470 // ok, but can faster?
00471 //    for (register int32 i=0; i<nrows; ++i)
00472 //      Result(i,0) = data[i][pixel];
00473 
00474 // NOT ok??
00475 //    for (register Type *pntA  = data[0]+pixel;
00476 //                        pntA  < data[nrows];
00477 //                        pntA += ncols)
00478 //      (*pntR++) = *pntA;
00479 
00480   Type *pntA = data[0]+pixel;
00481   Type *pntR = Result[0];
00482   for (register int32 i=0; i<nrows; ++i)
00483     {
00484     (*pntR++) = *pntA;
00485     pntA += ncols;
00486     }
00487 
00488   return Result;
00489   } // END getcolumn
00490 
00491 
00492 /****************************************************************
00493  * B.showdata()                                                 *
00494  *    Bert Kampes, 01-Feb-1999                                  *
00495  ****************************************************************/
00496 template <class Type>
00497 void matrix<Type>::showdata() const                     // show all data in matrix
00498   {
00499   #ifdef __DEBUGMAT1
00500     matDEBUG.print("showdata.");
00501   #endif
00502   #ifdef __DEBUGMAT2
00503   if (nrows>100 || ncols>15)
00504     {
00505     matDEBUG << "matrix ("
00506          << nrows << "," << ncols
00507          << "); only showing data (0:99,0:14).";
00508     matDEBUG.print();
00509     }
00510   #endif
00511   const uint L = (nrows<=100)  ? nrows  : 100;
00512   const uint P = (ncols<=15) ? ncols : 15;
00513   for (register int32 i=0; i<L; i++) { 
00514     for (register int32 j=0; j<P; j++) { 
00515       cout << setw(5) << data[i][j] << " "; }
00516     cout << endl; } 
00517   } // END showdata
00518 
00519 
00520 
00521 /****************************************************************
00522  * bool v = A.isvector()                                        *
00523  *    Bert Kampes, 01-Feb-1999                                  *
00524  ****************************************************************/
00525 template <class Type>
00526 bool matrix<Type>::isvector() const
00527   {
00528   return (nrows==1 || ncols==1) ? true : false;
00529   } // END isvector
00530 
00531 
00532 
00533 /****************************************************************
00534  * uint l = A.lines()                                           *
00535  *    Bert Kampes, 01-Feb-1999                                  *
00536  ****************************************************************/
00537 template <class Type>
00538 uint matrix<Type>::lines() const                // return number of lines
00539   {
00540   return nrows;
00541   } // END lines
00542 
00543 
00544 /****************************************************************
00545  * uint p = A.pixels()                                          *
00546  *    Bert Kampes, 01-Feb-1999                                  *
00547  ****************************************************************/
00548 template <class Type>
00549 uint matrix<Type>::pixels() const               // return number of pixels
00550   {
00551   return ncols;
00552   } // END pixels
00553 
00554 
00555 /****************************************************************
00556  * uint s = A.size()                                            *
00557  *    Bert Kampes, 01-Feb-1999                                  *
00558  ****************************************************************/
00559 template <class Type>
00560 uint matrix<Type>::size() const                         // return nsize
00561   {
00562   return nsize;
00563   } // END size
00564 
00565 
00566 /****************************************************************
00567  * A.resize(l,p)                                                *
00568  *    Bert Kampes, 01-Feb-1999                                  *
00569  ****************************************************************/
00570 template <class Type>
00571 void matrix<Type>::resize(uint l1, uint p1)
00572   {
00573   #ifdef __DEBUGMAT2
00574   matDEBUG.print("resize.");
00575   #endif
00576   if (l1 == nrows && p1 == ncols) return;
00577   else if (data!=0)                             // check for allocated memory
00578     {
00579     #ifdef __DEBUGMAT2
00580     uint deallocated  = sizeof(Type)*nsize;       // [B]
00581     totalallocated   -= deallocated;            // [B]
00582     matDEBUG << "deallocated matrix("
00583              << nrows << "," << ncols << ") at: " << data[0] << " ("
00584              << setw(10) << deallocated << "B; total: "
00585              << setw(10) << totalallocated << "B)";
00586     matDEBUG.print();
00587     #endif
00588     delete [] data[0];
00589     delete [] data;
00590     data=0;
00591     //#ifdef __USE_FFTW_LIBRARY__
00592     //fftw_fwd_plan = NULL;                  // plan only used for complex float!
00593     //fftw_inv_plan = NULL;                  // plan only used for complex float!
00594     //#endif
00595     }
00596   initialize(l1,p1);                            // set to 0
00597   } // END resize
00598 
00599 
00600 /****************************************************************
00601  * A.clean()                                                    *
00602  *    Bert Kampes, 01-Feb-1999                                  *
00603  ****************************************************************/
00604 template <class Type>
00605 void matrix<Type>::clean()              // sets 2 zero
00606   {
00607   memset(data[0],0,nsize*sizeof(Type));
00608   } // END clean
00609 
00610 
00611 
00612 
00613 /****************************************************************
00614  * B.setrow(row, data)                                          *
00615  * rows: 0 1 2 .., should fit exactly                           *       
00616  * orientation of vector is disregarded.                        *
00617  *    Bert Kampes, 12-Oct-1999                                  *
00618  ****************************************************************/
00619 template <class Type>
00620 void matrix<Type>::setrow(uint line, const matrix<Type> &LINE)
00621   {
00622   #ifdef __DEBUGMAT1
00623   checkindex(line,0);
00624   if (!(LINE.nrows == 1 || LINE.ncols == 1))
00625     matERROR.print("setrow: only vector input.");
00626   if (LINE.nsize != ncols)
00627     matERROR.print("setrow: sizeofvector should be same as matrix.");
00628   #endif
00629   memcpy(data[line],LINE[0],ncols*sizeof(Type));
00630   } // END setrow
00631 
00632 
00633 
00634 /****************************************************************
00635  * B.setrow(row, scalar)                                        *
00636  * rows: 0 1 2 .., should fit exactly                           *       
00637  * orientation of vector is disregarded.                        *
00638  *    Bert Kampes, 12-Oct-1999                                  *
00639  ****************************************************************/
00640 template <class Type>
00641 void matrix<Type>::setrow(uint line, Type scalar)
00642   {
00643   #ifdef __DEBUGMAT1
00644   checkindex(line,0);
00645   #endif
00646   for (register Type *pntB =&data[line][0];
00647                       pntB<=&data[line][ncols-1];
00648                       pntB++)
00649     *pntB = scalar;
00650   } // END setrow
00651 
00652 
00653 
00654 /****************************************************************
00655  * B.setcolumn(col, COL)                                        *
00656  * cols: 0 1 2 ..                                               *       
00657  * orientation of vector is disregarded.                        *
00658  *    Bert Kampes, 12-Oct-1999                                  *
00659  ****************************************************************/
00660 template <class Type>
00661 void matrix<Type>::setcolumn(uint pixel, const matrix<Type> &COLUMN)
00662   {
00663   #ifdef __DEBUGMAT1
00664   checkindex(0,pixel);
00665   if (!(COLUMN.nrows == 1 || COLUMN.ncols == 1))
00666     matERROR.print("setcolumn: only vector input.");
00667   if (COLUMN.nsize != nrows)
00668     matERROR.print("setcolumn: sizeofvector should be same as matrix.");
00669   #endif
00670   Type *pntCOL = COLUMN[0];
00671   for (register Type *pntB =&data[0][pixel];
00672                       pntB<=&data[nrows-1][pixel];
00673                       pntB+=ncols)
00674     *pntB = *pntCOL++;
00675   } // END setcolumn
00676 
00677 
00678 /****************************************************************
00679  * B.setcolumn(col, scalar)                                     *
00680  * cols: 0 1 2 ..                                               *       
00681  #%// BK 25-Sep-2000                                            *
00682  ****************************************************************/
00683 template <class Type>
00684 void matrix<Type>::setcolumn(uint pixel, Type scalar)
00685   {
00686   #ifdef __DEBUGMAT1
00687   checkindex(0,pixel);
00688   #endif
00689   for (register Type *pntB =&data[0][pixel];
00690                       pntB<=&data[nrows-1][pixel];
00691                       pntB+=ncols)
00692     *pntB = scalar;
00693   } // END setcolumn
00694 
00695 
00696 
00697 /****************************************************************
00698  * B.fliplr()                                                   *
00699  * Mirror in center vertical (flip left right).                 *
00700  *    Bert Kampes, 23-Mar-2000                                  *
00701  ****************************************************************/
00702 template <class Type>
00703 void matrix<Type>::fliplr()
00704   {
00705   #ifdef __DEBUGMAT2
00706   matDEBUG.print("fliplr.");
00707   #endif
00708   if (nrows==1)
00709     {
00710     Type *pnt1 = data[0];                                       // first one
00711     Type *pnt2 = data[0]+ncols-1;                       // last one
00712     Type tmp   = *pnt1;
00713     for (register int32 i=0; i<int32(ncols/2); ++i)     // floor
00714       {
00715       *pnt1++  = *pnt2;
00716       *pnt2--  = tmp;
00717       tmp      = *pnt1;
00718       //Type tmp = data[i];
00719       //data[i] = data[ncols-i-1];
00720       //data[ncols-i-1] = tmp;
00721       }
00722     }
00723   else
00724     {
00725     for (register int32 i=0; i<int32(ncols/2); ++i)     // floor
00726       {
00727       matrix<Type> tmp1 = getcolumn(i);
00728       matrix<Type> tmp2 = getcolumn(ncols-i-1);
00729       setcolumn(i,tmp2);
00730       setcolumn(ncols-i-1,tmp1);
00731       }
00732     }
00733   } // END fliplr
00734 
00735 
00736 
00737 /****************************************************************
00738  * B.flipud()                                                   *
00739  * Mirror in center vertical (flip left right).                 *
00740  *    Actually move data around, not pointers, to be sure       *
00741  *    veclib works ok, data is cont. in memory.                 *
00742  *    Bert Kampes, 23-Mar-2000                                  *
00743  ****************************************************************/
00744 template <class Type>
00745 void matrix<Type>::flipud()
00746   {
00747   #ifdef __DEBUGMAT2
00748   matDEBUG.print("flipud.");
00749   #endif
00750   if (ncols==1)
00751     {
00752     Type *pnt1 = data[0];                                       // first one
00753     Type *pnt2 = data[0]+nrows-1;                       // last one
00754     Type tmp   = *pnt1;
00755     for (register int32 i=0; i<int32(nrows/2); ++i)     // floor
00756       {
00757       *pnt1++  = *pnt2;
00758       *pnt2--  = tmp;
00759       tmp      = *pnt1;
00760       }
00761     }
00762   else
00763     {
00764     for (register int32 i=0; i<int32(ncols/2); ++i)     // floor
00765       {
00766       matrix<Type> tmp1 = getrow(i);
00767       matrix<Type> tmp2 = getrow(nrows-i-1);
00768       setrow(i,tmp2);
00769       setrow(nrows-i-1,tmp1);
00770       }
00771     }
00772   } // END flipud
00773 
00774 
00775 
00776 /****************************************************************
00777  * OVERLOADED OPS                                               *
00778  ****************************************************************/
00779 
00780 
00781 
00782 /****************************************************************
00783  * a = A[5] ?;                                                  *
00784  * Bert Kampes, 14-Jan-1999                                     *
00785  ****************************************************************/
00786 template <class Type>
00787 Type* matrix<Type>::operator [] (uint line) const
00788   {
00789   #ifdef __DEBUGMAT1
00790   checkindex(line,0);
00791   #endif
00792   return data[line];
00793   } // END []
00794 
00795 
00796 /****************************************************************
00797  * a = A(i,j);                                                  *
00798  * Bert Kampes, 14-Jan-1999                                     *
00799  ****************************************************************/
00800 template <class Type>
00801 Type& matrix<Type>::operator () (uint line, uint pixel) const
00802   {
00803   #ifdef __DEBUGMAT1
00804   checkindex(line,pixel);
00805   #endif
00806   return data[line][pixel];
00807   } // END ()
00808 
00809 
00810 
00811 /****************************************************************
00812  * matrix<T> B = A(window);                                     *
00813  * may not be to efficient cause twice allocation?              *
00814  * Bert Kampes, 31-Mar-2000                                     *
00815  ****************************************************************/
00816 template <class Type>
00817 matrix<Type> matrix<Type>::operator () (window win) const
00818   {
00819   matrix<Type> Res(win,*this);
00820   return Res;
00821   //return (win,*this);
00822   } // END (win)
00823 
00824 
00825 
00826 /****************************************************************
00827  * matrix<T> B = A(uint,uint,uint,uint);                        *
00828  * Bert Kampes, 31-Mar-2000                                     *
00829  ****************************************************************/
00830 template <class Type>
00831 matrix<Type> matrix<Type>::operator () (uint l0,uint lN,uint p0,uint pN) const
00832   {
00833   const window win(l0,lN,p0,pN);
00834   matrix<Type> Res(win,*this);
00835   return Res;
00836   } // END (4uint)
00837 
00838 
00839 
00840 /****************************************************************
00841  *  =                                                           *
00842  * Bert Kampes, 14-Jan-1999                                     *
00843  ****************************************************************/
00844 template <class Type>
00845 matrix<Type>& matrix<Type>::operator = (const matrix<Type>& A)
00846   {
00847   #ifdef __DEBUGMAT2
00848   matDEBUG.print("operator =");
00849   #endif
00850   if (this != &A)                               // prevent copy to itself
00851     {
00852     if (A.nsize)                                // if allocated
00853       {
00854       if (data != 0)                            // if allocated
00855         {
00856         #ifdef __DEBUGMAT2
00857         uint deallocated  = sizeof(Type)*nsize; // [B]
00858         totalallocated   -= deallocated;        // [B]
00859         matDEBUG << "deallocated matrix("
00860              << nrows << "," << ncols << ") at: " << data[0] << " ("
00861              << setw(10) << deallocated << "B; total: "
00862              << setw(10) << totalallocated << "B)";
00863         matDEBUG.print();
00864         #endif
00865 
00866         delete [] data[0];
00867         delete [] data;
00868         data = 0;
00869         }
00870       allocate(A.nrows,A.ncols);
00871       memcpy(data[0],A.data[0],nsize*sizeof(Type));
00872       }
00873     }
00874   return *this;
00875   } // END =
00876 
00877 
00878 /****************************************************************
00879  *  =                                                           *
00880  #%// BK 09-Nov-2000                                            *
00881  ****************************************************************/
00882 template <class Type>
00883 matrix<Type>& matrix<Type>::operator = (const Type scalar)
00884   {
00885   #ifdef __DEBUGMAT2
00886   matDEBUG.print("operator = (scalar)");
00887   #endif
00888   setdata(scalar);
00889   return *this;
00890   } // END = (scalar)
00891 
00892 
00893 
00894 
00895 /****************************************************************
00896  * C *= 5.0;                                                    *
00897  * Bert Kampes, 14-Jan-1999                                     *
00898  ****************************************************************/
00899 template <class Type>
00900 matrix<Type>& matrix<Type>::operator *= (Type scalar)
00901   {
00902   #ifdef __DEBUGMAT2
00903   matDEBUG.print("*=");
00904   #endif
00905   #ifdef __DEBUGMAT1
00906   if (!nsize)
00907     matERROR.print("matrixbaseclass: *= with empty matrix.");
00908   #endif
00909 
00910   Type *pntmat = data[0];
00911   for (register int32 i=0; i<nsize; ++i)
00912     (*pntmat++) *= scalar;      
00913   return *this;
00914   } // END *= scalar
00915 
00916 
00917 
00918 
00919 /****************************************************************
00920  * C *= A;      pointwise multiplication, a,c same size         *
00921  * Bert Kampes, 06-Oct-1999                                     *
00922  ****************************************************************/
00923 template <class Type>
00924 matrix<Type>& matrix<Type>::operator *= (const matrix<Type> &A)
00925   {
00926   #ifdef __DEBUGMAT2
00927   matDEBUG.print("*= pointwise");
00928   #endif
00929   #ifdef __DEBUGMAT1
00930   if (nrows != A.lines() || ncols != A.pixels())
00931     matERROR.print("matrixbaseclass: *= matrices must be same size.");
00932   #endif
00933 
00934   Type *pntmat = data[0];
00935   Type *pntA   = A[0];
00936   for (register int32 i=0; i<nsize; i++)
00937     (*pntmat++) *= (*pntA++);   
00938   return *this;
00939   } // END *= matrices
00940 
00941 
00942 
00943 
00944 /****************************************************************
00945  * C /= 5.0;                                                    *
00946  * Bert Kampes, 12-Oct-1999                                     *
00947  ****************************************************************/
00948 template <class Type>
00949 matrix<Type>& matrix<Type>::operator /= (Type scalar)
00950   {
00951   #ifdef __DEBUGMAT2
00952   matDEBUG.print("/= scalar");
00953   #endif
00954   #ifdef __DEBUGMAT1
00955   if (!nsize)
00956     matERROR.print("matrixbaseclass: /= with empty matrix.");
00957   #endif
00958     Type *pntmat = data[0];
00959     for (register int32 i=0; i<nsize; i++)
00960       (*pntmat++) /= scalar;    
00961     return *this;
00962   } // END /= scalar
00963 
00964 
00965 
00966 
00967 /****************************************************************
00968  * C /= A;      pointwise division, a,c same size               *
00969  * Bert Kampes, 06-Oct-1999                                     *
00970  ****************************************************************/
00971 template <class Type>
00972 matrix<Type>& matrix<Type>::operator /= (const matrix<Type> &A)
00973   {
00974   #ifdef __DEBUGMAT2
00975   matDEBUG.print("/=");
00976   #endif
00977   #ifdef __DEBUGMAT1
00978   if (nrows != A.lines() || ncols != A.pixels())
00979     matERROR.print("matrixbaseclass: /= matrices must be same size.");
00980   #endif
00981 
00982   Type *pntmat = data[0];
00983   Type *pntA   = A[0];
00984   for (register int32 i=0; i<nsize; i++)
00985     (*pntmat++) /= (*pntA++);   
00986   return *this;
00987   } // END /= matrices
00988 
00989 
00990 
00991 /****************************************************************
00992  * C -= A;                                                      *
00993  * Bert Kampes, 14-Jan-1999                                     *
00994  ****************************************************************/
00995 template <class Type>
00996 matrix<Type>& matrix<Type>::operator -= (const matrix<Type>& A)
00997   {
00998   #ifdef __DEBUGMAT2
00999   matDEBUG.print("-= mat");
01000   #endif
01001   #ifdef __DEBUGMAT1
01002   if (nrows != A.nrows || ncols != A.ncols)
01003     matERROR.print("error dimensions.");
01004   if (!A.nsize)
01005     matERROR.print("matrixbaseclass: -= with empty matrices.");
01006   #endif
01007 
01008   Type *pntmat = data[0];
01009   Type *pntA   = A.data[0];
01010   for (register int32 i=0; i<nsize; i++)
01011     (*pntmat++) -= (*pntA++);   
01012   return *this;
01013   } // END -=
01014 
01015 
01016 
01017 /****************************************************************
01018  * C -= 5.0;                                                    *
01019  * Bert Kampes, 26-Jan-1999                                     *
01020  ****************************************************************/
01021 template <class Type>
01022 matrix<Type>& matrix<Type>::operator -= (Type scalar)
01023   {
01024   #ifdef __DEBUGMAT2
01025   matDEBUG.print("-= scalar");
01026   #endif
01027   #ifdef __DEBUGMAT1
01028   if (!nsize)
01029     matERROR.print("matrixbaseclass: -= with empty matrix.");
01030   #endif
01031 
01032   Type *pntmat = data[0];
01033   for (register int32 i=0; i<nsize; i++)
01034     (*pntmat++) -= scalar;      
01035   return *this;
01036   } // END -= scalar
01037 
01038 
01039 
01040 
01041 /****************************************************************
01042  * C += A;                                                      *
01043  * Bert Kampes, 14-Jan-1999                                     *
01044  #%// Bert Kampes, 10-Apr-2005  (why const?)
01045  ****************************************************************/
01046 template <class Type>
01047 matrix<Type>& matrix<Type>::operator += (const matrix<Type>& A)
01048   {
01049   #ifdef __DEBUGMAT2
01050   matDEBUG.print("+=");
01051   #endif
01052   #ifdef __DEBUGMAT1
01053   if (nrows != A.nrows || ncols != A.ncols)
01054     matERROR.print("error dimensions.");
01055   if (!A.nsize)
01056     matERROR.print("matrixbaseclass: += with empty matrices.");
01057   #endif
01058 
01059   Type *pntmat = data[0];
01060   Type *pntA   = A.data[0];
01061   for (register int32 i=0; i<nsize; i++)
01062     (*pntmat++) += (*pntA++);   
01063   return *this;
01064   } // END +=
01065 
01066 
01067 
01068 
01069 /****************************************************************
01070  * C += 5.0;                                                    *
01071  * Bert Kampes, 26-Jan-1999                                     *
01072  ****************************************************************/
01073 template <class Type>
01074 matrix<Type>& matrix<Type>::operator += (Type scalar)
01075   {
01076   #ifdef __DEBUGMAT2
01077   matDEBUG.print("+= scalar");
01078   #endif
01079   #ifdef __DEBUGMAT1
01080   if (!nsize)
01081     matERROR.print("matrixbaseclass: += with empty matrix.");
01082   #endif
01083 
01084   Type *pntmat = data[0];
01085   for (register int32 i=0; i<nsize; i++)
01086     (*pntmat++) += scalar;      
01087   return *this;
01088   } // END += scalar
01089 
01090 
01091 
01092 /****************************************************************
01093  * A.conj();complex conjugated                                  *
01094  *  can be more efficient, but complex<float>::conj is unknown? *
01095  * Bert Kampes, 18-Oct-1999                                     *
01096  ****************************************************************/
01097 template <class Type>
01098 void matrix<Type>::conj()
01099   {
01100   #ifdef __DEBUGMAT2
01101   matDEBUG.print("conj"); 
01102   #endif
01103   //    // not ok..
01104   //    for (register Type *pnt=&data[0][0];
01105   //         pnt <= &data[nrows-1][ncols-1];
01106   //         *pnt = (conj(*pnt)++)); // do nothing
01107   // this works ok, but is not very fast...
01108   for (register int32 i=0;i<nrows;i++)
01109     for (register int32 j=0;j<ncols;j++)
01110       data[i][j] = Type((data[i][j]).real(),-(data[i][j]).imag());
01111   } // END conj()
01112 
01113 
01114 
01115 /****************************************************************
01116  * if (A==scalar) all elements equal scalar                     *
01117  * Bert Kampes, 08-Oct-1999                                     *
01118  ****************************************************************/
01119 template <class Type>
01120 bool matrix<Type>::operator == (Type scalar) const
01121   {
01122   #ifdef __DEBUGMAT2
01123   matDEBUG.print("== (scalar)");
01124   #endif
01125   if (nsize == 0)
01126     return false;
01127   bool same = true;
01128   Type *pnt = data[0];
01129   for (register int32 i=0; i<nsize; ++i)
01130     {
01131     if ((*pnt) != scalar)
01132       {
01133       same = false;
01134       break;
01135       }
01136     }
01137   return same;
01138   } // END ==
01139 
01140 
01141 
01142 /****************************************************************
01143  * if (A==B)                                                    *
01144  * Bert Kampes, 08-Oct-1999                                     *
01145  ****************************************************************/
01146 template <class Type>
01147 bool matrix<Type>::operator == (const matrix<Type> &A) const
01148   {
01149 #ifdef __DEBUGMAT2
01150   matDEBUG.print("==");
01151 #endif
01152   if ((A.lines() != nrows) || (A.pixels() != ncols))
01153     return false;
01154   bool same  = true;
01155   Type *pnt  = data[0];
01156   Type *pntA = A[0];
01157   for (register int32 i=0; i<nsize; ++i)
01158     {
01159     if ((*pnt) != (*pntA))
01160       {
01161       same = false;
01162       break;
01163       }
01164     }
01165   return same;
01166   } // END ==
01167 
01168 
01169 
01170 /****************************************************************
01171  * if (A!=scalar)                                               *
01172  * Bert Kampes, 08-Oct-1999                                     *
01173  ****************************************************************/
01174 template <class Type>
01175 bool matrix<Type>::operator != (Type scalar) const
01176   {
01177 #ifdef __DEBUGMAT2
01178   matDEBUG.print("!= (scalar)");
01179 #endif
01180   if (*this == scalar)
01181     return false;
01182   else
01183     return true;
01184   } // END !=
01185 
01186 
01187 
01188 /****************************************************************
01189  * if (A!=B)                                                    *
01190  * Bert Kampes, 08-Oct-1999                                     *
01191  ****************************************************************/
01192 template <class Type>
01193 bool matrix<Type>::operator != (const matrix<Type> &A) const
01194   {
01195 #ifdef __DEBUGMAT2
01196   matDEBUG.print("!=");
01197 #endif
01198   if (*this == A)
01199     return false;
01200   else
01201     return true;
01202   } // END !=
01203 
01204 
01205 
01206 
01207 // ++++++++++++++++++++++++++
01208 // template functions, stupid place, but else not found?
01209 // specializations in matrixspecs.c are NOT? used before these?
01210 // BK 16-Aug-2000
01211 // ++++++++++++++++++++++++++
01212 /****************************************************************
01213  * C = A * B;                                                   *
01214  * Bert Kampes, 14-Jan-1999                                     *
01215  ****************************************************************/
01216 template <class Type>
01217 matrix<Type> operator * (const matrix<Type>& A, const matrix<Type>& B)
01218   {
01219   #ifdef __DEBUGMAT2
01220   matDEBUG.print("matrices * (no veclib)");
01221   #endif
01222   #ifdef __DEBUGMAT1
01223   if (A.pixels() != B.lines())
01224     matERROR.print("matrix::operator *: multiplication not possible");
01225   if (!A.size())
01226     matERROR.print("matrixbaseclass: operator * with empty matrices.");
01227   #endif
01228   // ______ Straightforward, no veclib _____
01229   matrix<Type>  Result(A.lines(),B.pixels());
01230   register Type sum = Type(0.);
01231   for (register int32 i=0; i<Result.lines(); i++) 
01232     {
01233     for (register int32 j=0; j<Result.pixels(); j++) 
01234       {
01235       for (register int32 k=0; k<A.pixels(); k++) 
01236         {
01237         //sum += A.data[i][k] * B.data[k][j]; 
01238         sum += A(i,k) * B(k,j); 
01239         }
01240       //Result.data[i][j] = sum; 
01241       Result(i,j) = sum; 
01242       sum = Type(0.);                                   // complex requires this
01243       }
01244     }
01245   return Result; 
01246   } // END *
01247 
01248 
01249 /****************************************************************
01250  * C = 5.0 * B;                                                 *
01251  * Bert Kampes, 14-Jan-1999                                     *
01252  ****************************************************************/
01253 template <class Type>
01254 matrix<Type> operator * (const matrix<Type>& A, Type scalar)
01255   {
01256   #ifdef __DEBUGMAT2
01257   matDEBUG.print("scalar *");
01258   #endif
01259   // ______Perform multiplication______
01260   matrix<Type>  Result=A;
01261   return Result *= scalar;              // checks in *=
01262   } // END * scalar
01263 
01264 
01265 
01266 /****************************************************************
01267  * C = B * 5.0;                                                 *
01268  * Bert Kampes, 14-Jan-1999                                     *
01269  ****************************************************************/
01270 template <class Type>
01271 matrix<Type> operator *  (Type  scalar, const matrix<Type> &A)
01272   {
01273   #ifdef __DEBUGMAT2
01274   matDEBUG.print("* scalar.");
01275   #endif
01276   return A*scalar;                      // checks in *=
01277   } // END scalar *
01278 
01279 
01280 
01281 /****************************************************************
01282  * C = A / 5;                                                   *
01283  * Bert Kampes, 04-Apr-2000                                     *
01284  ****************************************************************/
01285 template <class Type>
01286 matrix<Type> operator / (const matrix<Type>& A, Type scalar)
01287   {
01288   #ifdef __DEBUGMAT2
01289   matDEBUG.print("/");
01290   #endif
01291   matrix<Type> Result = A;
01292   return Result      /= scalar;         // checks are performed here.
01293   } // END /
01294 
01295 
01296 
01297 /****************************************************************
01298  * C = A / B;                                                   *
01299  * Bert Kampes, 04-Apr-2000                                     *
01300  ****************************************************************/
01301 template <class Type>
01302 matrix<Type> operator / (const matrix<Type>& A, const matrix<Type>& B)
01303   {
01304   #ifdef __DEBUGMAT2
01305   matDEBUG.print("/");
01306   #endif
01307   matrix<Type> Result = A;
01308   return Result      /= B;      // checks are performed here.
01309   } // END /
01310 
01311 
01312 
01313 /****************************************************************
01314  * C = A - B;                                                   *
01315  * Bert Kampes, 14-Jan-1999                                     *
01316  ****************************************************************/
01317 template <class Type>
01318 matrix<Type> operator - (const matrix<Type>& A, const matrix<Type>& B)
01319   {
01320   #ifdef __DEBUGMAT2
01321   matDEBUG.print("-");
01322   #endif
01323   matrix<Type> Result = A;
01324   return Result      -= B;              // checks are performed here.
01325   } // END - (binary)
01326 
01327 
01328 
01329 /****************************************************************
01330  * C = A - 5;                                                   *
01331  * Bert Kampes, 04-Apr-2000                                     *
01332  ****************************************************************/
01333 template <class Type>
01334 matrix<Type> operator - (const matrix<Type>& A, Type scalar)
01335   {
01336   #ifdef __DEBUGMAT2
01337   matDEBUG.print("-");
01338   #endif
01339   matrix<Type> Result = A;
01340   return Result      -= scalar;         // checks are performed here.
01341   } // END - (binary)
01342 
01343 
01344 
01345 /****************************************************************
01346  * C = A + B;                                                   *
01347  * Bert Kampes, 14-Jan-1999                                     *
01348  ****************************************************************/
01349 template <class Type>
01350 matrix<Type> operator + (const matrix<Type>& A, const matrix<Type>& B)
01351   {
01352   #ifdef __DEBUGMAT2
01353   matDEBUG.print("+");
01354   #endif
01355   matrix<Type> Result = B;
01356   return Result      += A;              // checks are in +=
01357   } // END + (binary)
01358 
01359 
01360 
01361 /****************************************************************
01362  * C = A + 5;                                                   *
01363  * Bert Kampes, 04-Apr-2000                                     *
01364  ****************************************************************/
01365 template <class Type>
01366 matrix<Type> operator + (const matrix<Type>& A, Type scalar)
01367   {
01368   #ifdef __DEBUGMAT2
01369   matDEBUG.print("+");
01370   #endif
01371   matrix<Type> Result = A;
01372   return Result      += scalar;          // checks are performed here.
01373   } // END + (binary)
01374 
01375 
01376 
01377 /****************************************************************
01378  * a = max(A)                                                   *
01379  * Bert Kampes, 14-Jan-1999                                     *
01380  ****************************************************************/
01381 template <class Type>
01382 Type max(const matrix<Type> &A)
01383   {
01384   #ifdef __DEBUGMAT2
01385   matDEBUG.print("max");
01386   #endif
01387   Type m=A(0,0);
01388   for (register int32 i=0; i<A.lines(); ++i)
01389     for (register int32 j=0; j<A.pixels(); ++j)
01390       if (A(i,j)>m) m=A(i,j);
01391   return m;     
01392   } // END max
01393 
01394 
01395 
01396 /****************************************************************
01397  * a = max(A,linemax,pixelmax)                                  *
01398  * Bert Kampes, 14-Jan-1999                                     *
01399  ****************************************************************/
01400 template <class Type>
01401 Type max(const matrix<Type> &A, uint& line, uint& pixel)
01402   {
01403   #ifdef __DEBUGMAT2
01404   matDEBUG.print("max");
01405   #endif
01406   Type m=A(0,0);
01407   for (register int32 i=0; i<A.lines(); ++i)
01408     for (register int32 j=0; j<A.pixels(); ++j)
01409       if (A(i,j)>=m)
01410         {
01411         m     = A(i,j);
01412         line  = i;
01413         pixel = j;
01414         }
01415   return m;     
01416   } // END max
01417 
01418 
01419 
01420 /****************************************************************
01421  * a = min(A)                                                   *
01422  * Bert Kampes, 14-Jan-1999                                     *
01423  ****************************************************************/
01424 template <class Type>
01425 Type min(const matrix<Type> &A)
01426   {
01427   #ifdef __DEBUGMAT2
01428   matDEBUG.print("min");
01429   #endif
01430   Type m=A(0,0);
01431   for (register int32 i=0; i<A.lines(); ++i)
01432     for (register int32 j=0; j<A.pixels(); ++j)
01433       if (A(i,j)<m) m=A(i,j);
01434   return m;     
01435   } // END min
01436 
01437 
01438 
01439 /****************************************************************
01440  * a = min(A,linemax,pixelmax)                                  *
01441  * Bert Kampes, 14-Jan-1999                                     *
01442  ****************************************************************/
01443 template <class Type>
01444 Type min(const matrix<Type> &A, uint& line, uint& pixel)
01445   {
01446   #ifdef __DEBUGMAT2
01447   matDEBUG.print("min");
01448   #endif
01449   Type m=A(0,0);
01450   for (register int32 i=0; i<A.lines(); i++)
01451     for (register int32 j=0; j<A.pixels(); j++)
01452       if (A(i,j)<=m) 
01453         {
01454         m     = A(i,j);
01455         line  = i;
01456         pixel = j;
01457         }
01458   return m;     
01459   } // END min
01460 
01461 
01462 
01463 /****************************************************************
01464  * C=matTxmat(A,B) C=trans(A)*B; specialized for veclib         *
01465  *    Bert Kampes, 22-Feb-1999                                  *
01466  ****************************************************************/
01467 template <class Type>
01468 matrix<Type> matTxmat(const matrix<Type> &A, const matrix<Type> &B)
01469   {
01470   #ifdef __DEBUGMAT2
01471     matDEBUG.print("matTxmat: no veclib");
01472   #endif
01473   #ifdef __DEBUGMAT1
01474   if (A.lines() != B.lines())
01475     matERROR.print("matTxmat: size A,B: input is A,B; computed is trans(A)*B.");
01476   #endif
01477 
01478   matrix<Type> Result(A.pixels(),B.pixels());
01479   register Type sum = Type(0.);
01480   for (register int32 i=0; i<Result.lines(); i++)
01481     {
01482     for (register int32 j=0; j<Result.pixels(); j++)
01483       {
01484       for (register int32 k=0; k<A.lines(); k++)
01485         {
01486         sum += A(k,i) * B(k,j);
01487         //sum += A.data[k][i] * B.data[k][j];
01488         }
01489       //Result.data[i][j] = sum;
01490       Result(i,j) = sum;
01491       sum = Type(0.);
01492       }
01493     }
01494   return Result;
01495   } // END matTxmat
01496 
01497 
01498 
01499 /****************************************************************
01500  * C=matxmatT(A,B) C=A*trans(B); specialized for veclib         *
01501  *    Bert Kampes, 22-Feb-1999                                  *
01502  ****************************************************************/
01503 template <class Type>
01504 matrix<Type> matxmatT(const matrix<Type> &A, const matrix<Type> &B)
01505   {
01506   #ifdef __DEBUGMAT2
01507     matDEBUG.print("matxmatT: no veclib");
01508   #endif
01509   #ifdef __DEBUGMAT1
01510   if (A.pixels() != B.pixels())
01511     matERROR.print("matxmatT: size A,B: input is A,B; computed is A*trans(B).");
01512   #endif
01513 
01514   Type sum = Type(0.);
01515   matrix<Type> Result(A.lines(),B.lines());
01516   for (register int32 i=0; i<Result.lines(); i++)
01517     {
01518     for (register int32 j=0; j<Result.pixels(); j++)
01519       {
01520       for (register int32 k=0; k<A.pixels(); k++)
01521         {
01522         //sum += A.data[i][k] * B.data[j][k];
01523         sum += (A(i,k) * B(j,k));
01524         }
01525       //Result.data[i][j] = sum;
01526       Result(i,j) = sum;
01527       sum = Type(0.);
01528       }
01529     }
01530   return Result;
01531   } // END matxmatT
01532 
01533 
01534 
01535 /****************************************************************
01536  * dumpasc(file,A);                                             *
01537  * Bert Kampes, 14-Jan-1999                                     *
01538  ****************************************************************/
01539 template <class Type>
01540 void dumpasc(const char *file, const matrix<Type>& A)
01541   {
01542   #ifdef __DEBUGMAT2
01543   matDEBUG.print("dumpasc to file");
01544   #endif
01545   ofstream fo(file,ios::out | ios::trunc);
01546   matassert(fo,file,__FILE__,__LINE__);
01547   fo.precision(3);
01548   fo.width(11);
01549   // this is ok, but test if not: fo.setf(ios::left);
01550   fo.setf(ios::fixed);
01551 
01552   for (register int32 i=0; i<A.lines(); ++i)
01553     {
01554     for (register int32 j=0; j<A.pixels(); ++j)
01555       {
01556       //fo << A.data[i][j] << " ";
01557       fo << A(i,j) << " ";
01558       }
01559     fo << endl;
01560     }
01561   fo.close();
01562   } // END dumpasc
01563 
01564 
01565 
01566 
01567 /****************************************************************
01568  * C = dotmult(A,B) = A .* B                                    *
01569  * Bert Kampes, 26-Jan-1999                                     *
01570  ****************************************************************/
01571 template <class Type>
01572 matrix<Type>  dotmult     (const matrix<Type> &A, const matrix<Type> &B)
01573   {
01574   #ifdef __DEBUGMAT2
01575   matDEBUG.print("dotmult");
01576   #endif
01577   //#ifdef __DEBUGMAT1
01578   //if (A.lines() != B.lines() || A.pixels() != B.pixels())
01579     //matERROR.print("code 901: wrong input.");
01580   //#endif
01581 
01582   matrix<Type> Result = A;
01583   Result             *= B;                      // checks are here
01584   return Result;
01585 
01586   //Type *pntR   = Result.data[0];
01587   //Type *pntB   = B.data[0];
01588   //for (register int32 i=0; i<A.nsize; i++)
01589     //(*pntR++) *= (*pntB++);
01590   //return Result;
01591   } // END dotmult
01592 
01593 
01594 
01595 /****************************************************************
01596  * C = dotdiv(A,B) = A/B                                        *
01597  * Bert Kampes, 26-Jan-1999                                     *
01598  ****************************************************************/
01599 template <class Type>
01600 matrix<Type>  dotdiv      (const matrix<Type> &A, const matrix<Type> &B)
01601   {
01602   #ifdef __DEBUGMAT2
01603   matDEBUG.print("dotdiv");
01604   #endif
01605   //#ifdef __DEBUGMAT1
01606   //if (A.lines() != B.lines() || A.pixels() != B.pixels())
01607   //  matERROR.print("code 901: wrong input.");
01608   //#endif
01609 
01610   matrix<Type> Result = A;
01611   Result             /= B;                      // checks are here
01612   return Result;
01613 
01614   //Type *pntR   = Result.data[0];
01615   //Type *pntB   = B.data[0];
01616   //for (register int32 i=0; i<A.nsize; i++)
01617   //  (*pntR++) /= (*pntB++);
01618   //return Result;
01619   } // END dotdiv
01620 
01621 
01622 
01623 /****************************************************************
01624  * A = sqr(B)                                                   *
01625  * Bert Kampes, 16-Feb-1999                                     *
01626  ****************************************************************/
01627 template <class Type>
01628 matrix<Type>  sqr         (const matrix<Type> &A)
01629   {
01630   #ifdef __DEBUGMAT2
01631   matDEBUG.print("sqr");
01632   #endif
01633   matrix<Type> Result = dotmult(A,A);
01634 
01635 //    matrix<Type> Result(A.lines(),A.pixels());
01636 //    Type *pntR   = Result.data[0];
01637 //   
01638 //    // ______Ensure stride one memory______
01639 //    Type *pntA   = A.data[0];
01640 //    for (register int32 i=0; i<Result.nsize; i++)
01641 //      (*pntR++) = (sqr((*pntA++)));
01642 
01643   return Result;
01644   } // END sqr
01645 
01646 
01647 /****************************************************************
01648  *    conj                                                      *
01649  *    Bert Kampes, 02-Mar-1999                                  *
01650  ****************************************************************/
01651 template <class Type>
01652 matrix<Type> conj       (const matrix<Type> &A)
01653   {
01654   #ifdef __DEBUGMAT2
01655   matDEBUG.print("conj: repair stupid conj not known with gcc ?? 2.95.2");
01656   #endif
01657   // ______Ensure stride one memory______
01658   matrix<Type> Result(A.lines(),A.pixels());
01659   Type *pntA = A[0];
01660   Type *pntR = Result[0];
01661   for (register int32 i=0; i<int32(A.size()); i++)
01662     //#ifdef __GplusplusCOMPILER__         // conj not found?? with g++ 2.95.2?
01663     #ifdef __GNUC__         // conj not found?? with g++ 2.95.2?
01664       //(*pntR++) = Type::conj(*pntA++);
01665       //(*pntR++) = Type(*pntA).real(),-(*pntA++).imag());
01666       (*pntR++) = Type(real(*pntA),-imag(*pntA++));
01667     #else
01668       (*pntR++) = conj(*pntA++);
01669     #endif
01670     // not ok: (*pntR++) = Type(*pntA.real(),-*pntA.imag()++);
01671   return Result;
01672   } // END conj
01673 
01674 
01675 
01676 /****************************************************************
01677  * C=diagxmat(vec,B) C=diag(vec) * B;                           *
01678  *    Bert Kampes, 22-Feb-1999                                  *
01679  ****************************************************************/
01680 template <class Type>
01681 matrix<Type>  diagxmat    (const matrix<Type> &diag, const matrix<Type> &B)
01682   {
01683   #ifdef __DEBUGMAT2
01684   matDEBUG.print("diagxmat");
01685   #endif
01686   #ifdef __DEBUGMAT1
01687   if (min(diag.lines(),diag.pixels()) != 1)
01688     matERROR.print("diagxmat: sizes A,B: diag is vector.");
01689   if (diag.size() != B.lines())
01690     matERROR.print("diagxmat: sizes A,B: input is vector, matrix.");
01691   #endif
01692  
01693   matrix<Type> Result=B;
01694   if (diag.lines() != 1)        // standing
01695     {
01696     for (register int32 i=0; i<int32(Result.lines()); i++)
01697       for (register int32 j=0; j<int32(Result.pixels()); j++)
01698         Result(i,j) *= diag(i,0);
01699         //Result.data[i][j] *= diag.data[i][0];
01700     }
01701   else
01702     {
01703     for (register int32 i=0; i<int32(Result.lines()); i++)
01704       for (register int32 j=0; j<int32(Result.pixels()); j++)
01705         Result(i,j) *= diag(0,i);
01706         //Result.data[i][j] *= diag.data[0][i];
01707     }
01708   return Result;
01709   } // END diagxmat
01710 
01711 
01712 
01713 /****************************************************************
01714  * multilook A with factors l,p                                 *
01715  * lines(A) have to be multiple of factorL.                     *
01716  * size of output is (Al/fl,Ap/fp)                              *
01717  *  multilooked is averaged by factorLP.                        *
01718  * Bert Kampes, 19-Apr-1999                                     *
01719  ****************************************************************/
01720 template <class Type>
01721 matrix<Type>  multilook   (const matrix<Type> &A, uint factorL, uint factorP)
01722   {
01723   #ifdef __DEBUGMAT2
01724   matDEBUG.print("multilook.");
01725   #endif
01726   #ifdef __DEBUGMAT1
01727   if (A.lines()%factorL)
01728     matERROR.print("lines A must be multiple of factorL.");
01729   #endif
01730 
01731   if (factorL==1 && factorP==1)
01732     {
01733     matrix<Type> R=A;
01734     return R;
01735     }
01736  
01737   Type sum;
01738   Type factorLP = Type(factorL * factorP);
01739   matrix<Type> Result(A.lines()/factorL,A.pixels()/factorP);
01740   for (register int32 i=0; i<int32(Result.lines()); i++)
01741     {
01742     for (register int32 j=0; j<int32(Result.pixels()); j++)
01743       {
01744       sum = Type(0.);
01745       for (register int32 k=i*factorL; k<(i+1)*factorL; k++)
01746         {
01747         for (register int32 l=j*factorP; l<(j+1)*factorP; l++)
01748           {
01749           sum += A(k,l);
01750           }
01751         }
01752       Result(i,j) = sum / factorLP;
01753       }
01754     }
01755   return Result;
01756   } // END multilook
01757 
01758 
01759 /****************************************************************
01760  * Correlate A with maskB, return C(sizeA)                      *
01761  *  egde is set to zero, not pad with zeros                     *
01762  * Bert Kampes, 14-Jan-1999                                     *
01763  ****************************************************************/
01764 template <class Type>
01765 matrix<real4> correlate   (const matrix<Type> &A, matrix<Type> Mask)
01766   {
01767   #ifdef __DEBUGMAT2
01768   matDEBUG.print("Correlate.");
01769   matDEBUG.print("not yet correct for complex?: returns real4");
01770   if (Mask.lines()<2 || Mask.pixels()<2)
01771     matERROR.print("very small mask.");
01772   #endif
01773   #ifdef __DEBUGMAT1
01774   if (A.lines()<Mask.lines() || A.pixels()<Mask.pixels())
01775     matERROR.print("matrix input smaller than mask.");
01776   #endif
01777 
01778   real8 varM = 0.;                                // variance of Mask
01779   Mask      -= mean(Mask);
01780   //Type *pntMsk = Mask.data[0];
01781   Type *pntMsk = Mask[0];
01782   for (register int32 ii=0; ii<int32(Mask.size()); ii++)
01783     varM += sqr(*pntMsk++);                                     // 1/N later
01784 
01785 // ______Compute correlation at these points______
01786   uint beginl = (Mask.lines()-1)  / 2 ;                   // floor
01787   uint beginp = (Mask.pixels()-1) / 2 ;                   // floor
01788   matrix<real4> Result(A.lines(),A.pixels());       // init to 0
01789 
01790 // ______First window of A, updated at end of loop______
01791   window winA   (0, Mask.lines()-1, 0, Mask.pixels()-1);
01792   window windef (0,0,0,0);// defaults to total Am
01793 
01794 // ______Correlate part of Result______
01795   matrix<Type> Am(Mask.lines(),Mask.pixels());
01796   for (register int32 i=beginl; i<A.lines()-beginl; i++)
01797     {
01798     for (register int32 j=beginp; j<A.pixels()-beginp; j++)
01799       {
01800       Am.setdata(windef,A,winA);                // Am no allocs.
01801       Am -= mean(Am);                           // center around mean
01802       real8 covAM  = 0.;                        // covariance A,Mask
01803       real8 varA   = 0.;                        // variance of A(part)
01804       Type  *pntM  = Mask[0];
01805       Type  *pntAm = Am[0];
01806       for (register int32 l=0; l<Mask.size(); l++)
01807         {
01808         covAM += ((*pntM++) * (*pntAm));        // wait for move pnt
01809         varA  += sqr(*pntAm++);                 // pnt ++
01810         }
01811       Result(i,j) = covAM / sqrt(varM*varA);
01812       winA.pixlo++;
01813       winA.pixhi++;
01814       }
01815     winA.linelo++;
01816     winA.linehi++;
01817     winA.pixlo = 0;
01818     winA.pixhi = winA.pixlo + Mask.pixels() - 1;
01819     }
01820   return Result;
01821   } // END correlate
01822 
01823 
01824 
01825 /****************************************************************
01826  * C = -A;                                                      *
01827  * Bert Kampes, 14-Jan-1999                                     *
01828  ****************************************************************/
01829 template <class Type>
01830 matrix<Type>  operator - (const matrix<Type>& A)
01831   {
01832   #ifdef __DEBUGMAT2
01833   matDEBUG.print("operator -");
01834   #endif
01835   #ifdef __DEBUGMAT1
01836   if (!A.size())
01837     matERROR.print("matrixbaseclass: unary minus with empty matrix");
01838   #endif
01839   matrix<Type> Result(A.lines(),A.pixels());
01840   Type *pntA = A[0];
01841   Type *pntR = Result[0];
01842   for (register int32 i=0; i<Result.size(); ++i)
01843     (*pntR++) = -(*pntA++);
01844   return Result;
01845   } // END - (unary)
01846 
01847 
01848 
01849 /****************************************************************
01850  * a = mean(A)                                                  *
01851  * Bert Kampes, 14-Jan-1999                                     *
01852  ****************************************************************/
01853 template <class Type>
01854 real8         mean        (const matrix<Type> &A)
01855   {
01856   #ifdef __DEBUGMAT2
01857   matDEBUG.print("mean");
01858   #endif
01859   real8 sum=0.;
01860   // ______Ensure stride one memory______
01861   Type *pntA   = A[0];
01862   for (register int32 i=0; i<A.size(); ++i)
01863     sum += (*pntA++);
01864   return sum/real8(A.size());
01865   } // END mean
01866 
01867 
01868 
01869 /****************************************************************
01870  * matrix<Type> S = sum(A,dim)                                  *
01871  * return sum over dim of A. dim=1 by default.                  *
01872  * always returns a matrix, not very handy...                   *
01873  * A=[1 2 3]  then: sum(A,1) = [5 7 9]; sum(A,2)=[6]            *
01874  *   [4 5 6]                                     [15]           *
01875  * Bert Kampes, 28-Mar-2000                                     *
01876  ****************************************************************/
01877 template <class Type>
01878 matrix<Type>  sum         (const matrix<Type> &A, int32 dim)
01879   {
01880   #ifdef __DEBUGMAT2
01881   matDEBUG.print("sum");
01882   #endif
01883   register int32 i,j;
01884   Type sum = Type(0);
01885   matrix<Type> Res;                             // may be 1x1 ...
01886   if (A.isvector())
01887     {
01888     Res.resize(1,1);
01889     Type *pntA   = A[0];
01890     for (i=0; i<A.size(); i++)
01891       sum += (*pntA++);
01892     Res(0,0) = sum;
01893     }
01894   else // no vector
01895     {
01896     switch (dim)
01897       {
01898       // ______ sum over rows ______
01899       case 1:
01900         {
01901         Res.resize(1,A.pixels());
01902         for (i=0; i<A.pixels(); ++i)
01903           {
01904           for (j=0; j<A.lines(); ++j)
01905             {
01906             //sum += A.data[j][i];
01907             sum += A(j,i);
01908             }
01909           //Res.data[0][i] = sum;
01910           Res(0,i) = sum;
01911           sum = Type(0);
01912           }
01913         }
01914         break;
01915       // ______ sum over columns, may be done by pointers for speed ______
01916       case 2:
01917         {
01918         Res.resize(A.lines(),1);
01919         for (i=0; i<A.lines(); ++i)
01920           {
01921           for (j=0; j<A.pixels(); ++j)
01922             {
01923             //sum += A.data[i][j];
01924             sum += A(i,j);
01925             }
01926           //Res.data[i][0] = sum;
01927           Res(i,0) = sum;
01928           sum = Type(0);
01929           }
01930         }
01931         break;
01932       default: matERROR.print("sum: dim!={1,2}");
01933       }
01934     } // ifelse vector
01935   return Res;
01936   } // END sum
01937 
01938 
01939 
01940 /****************************************************************
01941  * A.mypow(scalar)                                              *
01942  * NOmatrix<Type> S = pow(A,scalar)                             *
01943  * return pointwize power.                                      *
01944  * always returns a matrix, not very handy...                   *
01945  #%// BK 26-Oct-2000
01946  ****************************************************************/
01947 template <class Type>
01948 void matrix<Type>::mypow(Type s)
01949   {
01950   for (register int32 i=0; i<nrows; ++i)
01951     for (register int32 j=0; j<ncols; ++j)
01952       data[i][j] = pow(data[i][j],s);
01953   } // END sum
01954 
01955 
01956 
01957 

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