00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "constants.hh"
00054
00055 #include <iostream>
00056 #include <fstream>
00057 #include <strstream>
00058 #include <iomanip>
00059 #include <algorithm>
00060 #include <cstring>
00061 #include <complex>
00062 #ifdef __DEBUGMAT1 // use index checking, alloc
00063 #include <new>
00064 #endif
00065
00066
00067
00068
00069
00070
00071 #ifdef __DEBUGMAT2 // use index checking, alloc
00072 extern uint totalallocated;
00073 #endif
00074
00075
00076 extern bk_messages matERROR;
00077 extern bk_messages matDEBUG;
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 template <class Type>
00093 void matrix<Type>::allocate(uint numlines, uint numpixels)
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
00109
00110 try
00111 {
00112
00113 data = new Type*[numlines];
00114
00115 }
00116 catch (bad_alloc)
00117 {
00118 matERROR << "code 502: allocation failed: 1size: "
00119 << numlines << ", " << numpixels;
00120 matERROR.print();
00121 }
00122 try
00123 {
00124
00125 data[0] = new Type[nsize];
00126
00127 }
00128 catch(bad_alloc)
00129 {
00130 matERROR << "code 502: allocation failed: 2size: "
00131 << numlines << ", " << numpixels;
00132 matERROR.print();
00133 }
00134
00135 for (register int32 i=1; i<int32(numlines); i++)
00136 data[i] = data[i-1]+numpixels;
00137
00138 #ifdef __DEBUGMAT2
00139
00140
00141 uint allocated = sizeof(Type)*nsize;
00142 totalallocated += allocated;
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 }
00150
00151
00152
00153
00154
00155
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 }
00166
00167
00168
00169 #ifdef __DEBUGMAT1
00170
00171
00172
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 }
00185 #endif
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 template <class Type>
00196 matrix<Type>::matrix()
00197 {
00198 nrows = 0;
00199 ncols = 0;
00200 nsize = 0;
00201 data = 0;
00202 }
00203
00204
00205
00206
00207
00208
00209
00210 template <class Type>
00211 matrix<Type>::matrix(uint lines, uint pixels)
00212 {
00213 initialize(lines,pixels);
00214 }
00215
00216
00217
00218
00219
00220
00221
00222
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 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
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
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
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 }
00269
00270
00271
00272
00273
00274
00275
00276
00277 template <class Type>
00278 matrix<Type>::~matrix()
00279 {
00280 if (data==0) return;
00281 #ifdef __DEBUGMAT2
00282 uint deallocated = sizeof(Type)*nsize;
00283 totalallocated -= deallocated;
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];
00292 delete [] data;
00293 data=0;
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303 template <class Type>
00304 void matrix<Type>::setdata(Type w)
00305 {
00306 #ifdef __DEBUGMAT2
00307 matDEBUG.print("A.setdata(w)");
00308 #endif
00309
00310 Type *pnt = data[0];
00311 for (register int32 i=0;i<nsize;i++)
00312 (*pnt++) = Type(w);
00313
00314
00315
00316
00317
00318 }
00319
00320
00321
00322
00323
00324
00325
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);
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 }
00337
00338
00339
00340
00341
00342
00343
00344
00345
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
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
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 }
00381
00382
00383
00384
00385
00386
00387
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
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 }
00409
00410
00411
00412
00413
00414
00415
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
00430
00431 uint numlin = win.linehi - win.linelo + 1;
00432 uint numpix = win.pixhi - win.pixlo + 1;
00433 matrix<Type> Result(numlin,numpix);
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 }
00438
00439
00440
00441
00442
00443
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 }
00455
00456
00457
00458
00459
00460
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
00471
00472
00473
00474
00475
00476
00477
00478
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 }
00490
00491
00492
00493
00494
00495
00496 template <class Type>
00497 void matrix<Type>::showdata() const
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 }
00518
00519
00520
00521
00522
00523
00524
00525 template <class Type>
00526 bool matrix<Type>::isvector() const
00527 {
00528 return (nrows==1 || ncols==1) ? true : false;
00529 }
00530
00531
00532
00533
00534
00535
00536
00537 template <class Type>
00538 uint matrix<Type>::lines() const
00539 {
00540 return nrows;
00541 }
00542
00543
00544
00545
00546
00547
00548 template <class Type>
00549 uint matrix<Type>::pixels() const
00550 {
00551 return ncols;
00552 }
00553
00554
00555
00556
00557
00558
00559 template <class Type>
00560 uint matrix<Type>::size() const
00561 {
00562 return nsize;
00563 }
00564
00565
00566
00567
00568
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)
00578 {
00579 #ifdef __DEBUGMAT2
00580 uint deallocated = sizeof(Type)*nsize;
00581 totalallocated -= deallocated;
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
00592
00593
00594
00595 }
00596 initialize(l1,p1);
00597 }
00598
00599
00600
00601
00602
00603
00604 template <class Type>
00605 void matrix<Type>::clean()
00606 {
00607 memset(data[0],0,nsize*sizeof(Type));
00608 }
00609
00610
00611
00612
00613
00614
00615
00616
00617
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 }
00631
00632
00633
00634
00635
00636
00637
00638
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 }
00651
00652
00653
00654
00655
00656
00657
00658
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 }
00676
00677
00678
00679
00680
00681
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 }
00694
00695
00696
00697
00698
00699
00700
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];
00711 Type *pnt2 = data[0]+ncols-1;
00712 Type tmp = *pnt1;
00713 for (register int32 i=0; i<int32(ncols/2); ++i)
00714 {
00715 *pnt1++ = *pnt2;
00716 *pnt2-- = tmp;
00717 tmp = *pnt1;
00718
00719
00720
00721 }
00722 }
00723 else
00724 {
00725 for (register int32 i=0; i<int32(ncols/2); ++i)
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 }
00734
00735
00736
00737
00738
00739
00740
00741
00742
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];
00753 Type *pnt2 = data[0]+nrows-1;
00754 Type tmp = *pnt1;
00755 for (register int32 i=0; i<int32(nrows/2); ++i)
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)
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 }
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
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 }
00794
00795
00796
00797
00798
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 }
00808
00809
00810
00811
00812
00813
00814
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
00822 }
00823
00824
00825
00826
00827
00828
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 }
00837
00838
00839
00840
00841
00842
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)
00851 {
00852 if (A.nsize)
00853 {
00854 if (data != 0)
00855 {
00856 #ifdef __DEBUGMAT2
00857 uint deallocated = sizeof(Type)*nsize;
00858 totalallocated -= deallocated;
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 }
00876
00877
00878
00879
00880
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 }
00891
00892
00893
00894
00895
00896
00897
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 }
00915
00916
00917
00918
00919
00920
00921
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 }
00940
00941
00942
00943
00944
00945
00946
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 }
00963
00964
00965
00966
00967
00968
00969
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 }
00988
00989
00990
00991
00992
00993
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 }
01014
01015
01016
01017
01018
01019
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 }
01037
01038
01039
01040
01041
01042
01043
01044
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 }
01065
01066
01067
01068
01069
01070
01071
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 }
01089
01090
01091
01092
01093
01094
01095
01096
01097 template <class Type>
01098 void matrix<Type>::conj()
01099 {
01100 #ifdef __DEBUGMAT2
01101 matDEBUG.print("conj");
01102 #endif
01103
01104
01105
01106
01107
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 }
01112
01113
01114
01115
01116
01117
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 }
01139
01140
01141
01142
01143
01144
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 }
01167
01168
01169
01170
01171
01172
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 }
01185
01186
01187
01188
01189
01190
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 }
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
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
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
01238 sum += A(i,k) * B(k,j);
01239 }
01240
01241 Result(i,j) = sum;
01242 sum = Type(0.);
01243 }
01244 }
01245 return Result;
01246 }
01247
01248
01249
01250
01251
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
01260 matrix<Type> Result=A;
01261 return Result *= scalar;
01262 }
01263
01264
01265
01266
01267
01268
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;
01277 }
01278
01279
01280
01281
01282
01283
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;
01293 }
01294
01295
01296
01297
01298
01299
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;
01309 }
01310
01311
01312
01313
01314
01315
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;
01325 }
01326
01327
01328
01329
01330
01331
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;
01341 }
01342
01343
01344
01345
01346
01347
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;
01357 }
01358
01359
01360
01361
01362
01363
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;
01373 }
01374
01375
01376
01377
01378
01379
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 }
01393
01394
01395
01396
01397
01398
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 }
01417
01418
01419
01420
01421
01422
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 }
01436
01437
01438
01439
01440
01441
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 }
01460
01461
01462
01463
01464
01465
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
01488 }
01489
01490 Result(i,j) = sum;
01491 sum = Type(0.);
01492 }
01493 }
01494 return Result;
01495 }
01496
01497
01498
01499
01500
01501
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
01523 sum += (A(i,k) * B(j,k));
01524 }
01525
01526 Result(i,j) = sum;
01527 sum = Type(0.);
01528 }
01529 }
01530 return Result;
01531 }
01532
01533
01534
01535
01536
01537
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
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
01557 fo << A(i,j) << " ";
01558 }
01559 fo << endl;
01560 }
01561 fo.close();
01562 }
01563
01564
01565
01566
01567
01568
01569
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
01578
01579
01580
01581
01582 matrix<Type> Result = A;
01583 Result *= B;
01584 return Result;
01585
01586
01587
01588
01589
01590
01591 }
01592
01593
01594
01595
01596
01597
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
01606
01607
01608
01609
01610 matrix<Type> Result = A;
01611 Result /= B;
01612 return Result;
01613
01614
01615
01616
01617
01618
01619 }
01620
01621
01622
01623
01624
01625
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
01636
01637
01638
01639
01640
01641
01642
01643 return Result;
01644 }
01645
01646
01647
01648
01649
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
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
01663 #ifdef __GNUC__ // conj not found?? with g++ 2.95.2?
01664
01665
01666 (*pntR++) = Type(real(*pntA),-imag(*pntA++));
01667 #else
01668 (*pntR++) = conj(*pntA++);
01669 #endif
01670
01671 return Result;
01672 }
01673
01674
01675
01676
01677
01678
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)
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
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
01707 }
01708 return Result;
01709 }
01710
01711
01712
01713
01714
01715
01716
01717
01718
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 }
01757
01758
01759
01760
01761
01762
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.;
01779 Mask -= mean(Mask);
01780
01781 Type *pntMsk = Mask[0];
01782 for (register int32 ii=0; ii<int32(Mask.size()); ii++)
01783 varM += sqr(*pntMsk++);
01784
01785
01786 uint beginl = (Mask.lines()-1) / 2 ;
01787 uint beginp = (Mask.pixels()-1) / 2 ;
01788 matrix<real4> Result(A.lines(),A.pixels());
01789
01790
01791 window winA (0, Mask.lines()-1, 0, Mask.pixels()-1);
01792 window windef (0,0,0,0);
01793
01794
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);
01801 Am -= mean(Am);
01802 real8 covAM = 0.;
01803 real8 varA = 0.;
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));
01809 varA += sqr(*pntAm++);
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 }
01822
01823
01824
01825
01826
01827
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 }
01846
01847
01848
01849
01850
01851
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
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 }
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
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;
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
01895 {
01896 switch (dim)
01897 {
01898
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
01907 sum += A(j,i);
01908 }
01909
01910 Res(0,i) = sum;
01911 sum = Type(0);
01912 }
01913 }
01914 break;
01915
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
01924 sum += A(i,j);
01925 }
01926
01927 Res(i,0) = sum;
01928 sum = Type(0);
01929 }
01930 }
01931 break;
01932 default: matERROR.print("sum: dim!={1,2}");
01933 }
01934 }
01935 return Res;
01936 }
01937
01938
01939
01940
01941
01942
01943
01944
01945
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 }
01954
01955
01956
01957