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
00054
00055
00056 #ifndef MATRIXBK_H // guard
00057 #define MATRIXBK_H
00058
00059 #include "constants.hh"
00060 #include <fstream>
00061
00062
00063
00064
00065
00066
00067
00068 extern bk_messages matERROR;
00069 extern bk_messages matDEBUG;
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 template <class Type> class matrix;
00094
00095 #ifdef __USE_VECLIB_LIBRARY__
00096
00097
00098
00099 matrix<real4> operator * (const matrix<real4>& A, const matrix<real4>& B);
00100 matrix<real8> operator * (const matrix<real8>& A, const matrix<real8>& B);
00101 matrix<complr4> operator * (const matrix<complr4>& A, const matrix<complr4>& B);
00102 matrix<complr8> operator * (const matrix<complr8>& A, const matrix<complr8>& B);
00103 matrix<real4> matTxmat (const matrix<real4> &A, const matrix<real4> &B);
00104 matrix<real8> matTxmat (const matrix<real8> &A, const matrix<real8> &B);
00105 matrix<complr4> matTxmat (const matrix<complr4> &A, const matrix<complr4> &B);
00106 matrix<complr8> matTxmat (const matrix<complr8> &A, const matrix<complr8> &B);
00107 matrix<real4> matxmatT (const matrix<real4> &A, const matrix<real4> &B);
00108 matrix<real8> matxmatT (const matrix<real8> &A, const matrix<real8> &B);
00109 matrix<complr4> matxmatT (const matrix<complr4> &A, const matrix<complr4> &B);
00110 matrix<complr8> matxmatT (const matrix<complr8> &A, const matrix<complr8> &B);
00111 #endif
00112
00113
00114
00115
00116
00117
00118
00119 void matassert(
00120 const ofstream &str,
00121 const char* ofilename,
00122 const char* callingfilename = "?",
00123 int32 linenumber = 0);
00124 void matassert(
00125 const ifstream &str,
00126 const char* ifilename,
00127 const char* callingfilename = "?",
00128 int32 linenumber = 0);
00129
00130
00131
00132
00133
00134
00135
00136 void fft (matrix<complr4> &A, int32 dimension);
00137 void ifft (matrix<complr4> &A, int32 dimension);
00138 void fft2d (matrix<complr4> &A);
00139 void ifft2d (matrix<complr4> &A);
00140
00141
00142
00143 matrix<complr4> oversample (matrix<complr4> A,
00144 uint frow, uint fcol);
00145 matrix<real4> oversample (const matrix<real4> &A,
00146 uint frow, uint fcol);
00147
00148
00149
00150 void choles (matrix<real4> &A);
00151 void invertchol (matrix<real4> &A);
00152 void solvechol (const matrix<real4> &A, matrix<real4> &B);
00153 void choles (matrix<real8> &A);
00154 void invertchol (matrix<real8> &A);
00155 void solvechol (const matrix<real8> &A, matrix<real8> &B);
00156
00157
00158 matrix<real4> intensity (const matrix<complr4> &A);
00159 matrix<real4> magnitude (const matrix<complr4> &A);
00160
00161 matrix<real4> real (const matrix<complr4> &A);
00162 matrix<real4> imag (const matrix<complr4> &A);
00163
00164 real4 norm2 (const matrix<complr4> &A);
00165 real4 norm2 (const matrix<real4> &A);
00166 matrix<complr4> norm (const matrix<complr4> &A);
00167
00168 matrix<real4> abs (const matrix<real4> &A);
00169 matrix<real8> abs (const matrix<real8> &A);
00170
00171
00172 void fileci2tomatcr4(
00173 matrix<complr4> &Result,
00174 const char *file,
00175 uint filelines,
00176 window win,
00177 window winoffset);
00178
00179
00180
00181
00182
00183 matrix<complr4> mat2cr4(
00184 const matrix<real4> &A);
00185 matrix<complr4> mat2cr4(
00186 const matrix<real4>& A,
00187 const matrix<real4>& B);
00188 matrix<complr4> mat2cr4(
00189 const matrix<real8>& A,
00190 const matrix<real8>& B);
00191
00192
00193 matrix<real4> angle(
00194 const matrix<complr4> &A);
00195 matrix<complr4> angle2cmplx(
00196 const matrix<real4> &A);
00197 void dotmultconjphase(
00198 matrix<complr4> &complexinterferogram,
00199 const matrix<real4> &refphase);
00200
00201
00202
00203 matrix<complr4> coherence(
00204 const matrix<complr4> &complex_interferogram,
00205 const matrix<complr4> &norm_image1_and_2,
00206 uint estimatorwinsizeL,
00207 uint estimatorwinsizeP);
00208
00209
00210 matrix<real4> coherence2(
00211 const matrix<complr4> &complex_interferogram,
00212 const matrix<complr4> &norm_image1_and_2,
00213 uint estimatorwinsizeL,
00214 uint estimatorwinsizeP);
00215
00216
00217
00218
00219
00220 void mysort2(matrix<real4> &A);
00221 void mysort2(matrix<int32> &A);
00222
00223
00224
00225
00226
00227 void dotmult (complr4 *startaddress, const matrix<real4> &B, int32 strike);
00228
00229
00230
00231
00232
00233
00234
00235
00236 matrix<real4> cos (const matrix<real4> &A);
00237 matrix<real8> cos (const matrix<real8> &A);
00238 matrix<real4> sin (const matrix<real4> &A);
00239 matrix<real8> sin (const matrix<real8> &A);
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 template <class Type>
00250 matrix<Type> operator * (const matrix<Type>& A, const matrix<Type>& B);
00251 template <class Type>
00252 matrix<Type> operator * (const matrix<Type>& A, Type scalar);
00253 template <class Type>
00254 matrix<Type> operator * (Type scalar, const matrix<Type> &A);
00255 template <class Type>
00256 matrix<Type> operator / (const matrix<Type>& A, Type B);
00257 template <class Type>
00258 matrix<Type> operator / (const matrix<Type>& A, const matrix<Type>& B);
00259 template <class Type>
00260 matrix<Type> operator - (const matrix<Type>& A, const matrix<Type>& B);
00261 template <class Type>
00262 matrix<Type> operator - (const matrix<Type>& A, Type scalar);
00263 template <class Type>
00264 matrix<Type> operator + (const matrix<Type>& A, const matrix<Type>& B);
00265 template <class Type>
00266 matrix<Type> operator + (const matrix<Type>& A, Type B);
00267 template <class Type>
00268 Type max(const matrix<Type> &A);
00269 template <class Type>
00270 Type max(const matrix<Type> &A, uint& line, uint& pixel);
00271 template <class Type>
00272 Type min(const matrix<Type> &A);
00273 template <class Type>
00274 Type min(const matrix<Type> &A, uint& line, uint& pixel);
00275 template <class Type>
00276 matrix<Type> matTxmat(const matrix<Type> &A, const matrix<Type> &B);
00277 template <class Type>
00278 matrix<Type> matxmatT(const matrix<Type> &A, const matrix<Type> &B);
00279 template <class Type>
00280 void dumpasc(const char *file, const matrix<Type>& A);
00281 template <class Type>
00282 matrix<Type> dotmult (const matrix<Type> &A, const matrix<Type> &B);
00283 template <class Type>
00284 matrix<Type> dotdiv (const matrix<Type> &A, const matrix<Type> &B);
00285 template <class Type>
00286 matrix<Type> sqr (const matrix<Type> &A);
00287 template <class Type>
00288 matrix<Type> conj (const matrix<Type> &A);
00289 template <class Type>
00290 matrix<Type> diagxmat (const matrix<Type> &diag, const matrix<Type> &B);
00291
00292 matrix<complr4> diagxmat (const matrix<real4> &diag, const matrix<complr4> &B);
00293 template <class Type>
00294 matrix<Type> multilook (const matrix<Type> &A, uint factorL, uint factorP);
00295 template <class Type>
00296 matrix<real4> correlate (const matrix<Type> &A, matrix<Type> Mask);
00297 template <class Type>
00298 matrix<Type> operator - (const matrix<Type>& A);
00299 template <class Type>
00300 real8 mean (const matrix<Type> &A);
00301 template <class Type>
00302 matrix<Type> sum (const matrix<Type> &A, int32 dim);
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 template <class Type>
00317 class matrix
00318 {
00319 private:
00320
00321 Type **data;
00322 uint nrows;
00323 uint ncols;
00324 uint nsize;
00325
00326
00327 void allocate (uint l, uint p);
00328 void initialize (uint l, uint p);
00329
00330 #ifdef __DEBUGMAT1
00331 void checkindex (uint l, uint p) const;
00332 #endif
00333
00334 public:
00335
00336 matrix ();
00337 matrix (uint l, uint p);
00338 matrix (const matrix<Type>& A);
00339 matrix (window w, const matrix<Type>& A);
00340
00341
00342 ~matrix ();
00343
00344
00345 void setdata (Type w);
00346 void setdata (uint l, uint p, const matrix<Type>& A);
00347 void setdata (window winin, const matrix<Type> &A, window winA);
00348 void setdata (const matrix<Type> &A, window winA);
00349 matrix<Type> getdata (window w) const;
00350 matrix<Type> getrow (uint l) const;
00351 matrix<Type> getcolumn (uint p) const;
00352
00353
00354
00355
00356 void showdata () const;
00357 bool isvector () const;
00358 uint lines () const;
00359 uint pixels () const;
00360 uint size () const;
00361 void resize (uint l, uint p);
00362 void clean ();
00363
00364 void setrow (uint l, const matrix<Type> &L);
00365 void setrow (uint l, Type scalar);
00366 void setcolumn (uint p, const matrix<Type> &C);
00367 void setcolumn (uint p, Type scalar);
00368 void fliplr ();
00369 void flipud ();
00370
00371 Type* operator [] (uint l) const;
00372 Type& operator () (uint l, uint p) const;
00373 matrix<Type> operator () (window win) const;
00374 matrix<Type> operator () (uint l0,uint lN,uint p0,uint pN) const;
00375
00376 matrix<Type>& operator = (const matrix<Type>& A);
00377 matrix<Type>& operator = (const Type scalar);
00378 matrix<Type>& operator -= (Type scalar);
00379 matrix<Type>& operator -= (const matrix<Type>& A);
00380 matrix<Type>& operator += (Type scalar);
00381 matrix<Type>& operator += (const matrix<Type>& A);
00382 matrix<Type>& operator *= (Type scalar);
00383 matrix<Type>& operator *= (const matrix<Type> &A);
00384 matrix<Type>& operator /= (Type scalar);
00385 matrix<Type>& operator /= (const matrix<Type> &A);
00386 bool operator == (Type scalar) const;
00387 bool operator == (const matrix<Type> &A) const;
00388 bool operator != (Type scalar) const;
00389 bool operator != (const matrix<Type> &A) const;
00390 void conj ();
00391 void mypow (Type s);
00392
00393
00394
00395 template <class TypeB>
00396 matrix<Type>& operator *= (const matrix<TypeB> &A);
00397 template <class TypeB>
00398 matrix<Type>& operator /= (const matrix<TypeB> &A);
00399 template <class TypeB>
00400 matrix<Type>& operator += (const matrix<TypeB> &A);
00401 template <class TypeB>
00402 matrix<Type>& operator -= (const matrix<TypeB> &A);
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416 friend ostream& operator << (ostream& file, const matrix<Type>& A)
00417 {
00418 const uint sizeofType = sizeof(Type);
00419 for (register int32 i=0;i<A.nrows;i++)
00420 for (register int32 j=0;j<A.ncols;j++)
00421 file.write((char*)&A.data[i][j],sizeofType);
00422 return file;
00423 }
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 friend istream& operator >> (istream& file, matrix<Type>& A)
00434 {
00435 const uint sizeofType = sizeof(Type);
00436 for (register int32 i=0;i<A.nrows;i++)
00437 for (register int32 j=0;j<A.ncols;j++)
00438 file.read((char*)&A.data[i][j],sizeofType);
00439
00440
00441 return file;
00442 }
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452 friend void myswap (matrix<Type> &A, matrix<Type> &B)
00453 {
00454 #ifdef __DEBUGMAT2
00455 matDEBUG.print("myswap.");
00456 #endif
00457 #ifdef __DEBUGMAT1
00458 if (A.nrows != B.nrows ||
00459 A.ncols != B.ncols )
00460 matERROR.print("swap: matrices must be same size (for now).");
00461 #endif
00462
00463 Type **pntA = A.data;
00464 Type **pntB = B.data;
00465 A.data = pntB;
00466 B.data = pntA;
00467 }
00468
00469
00470
00471
00472
00473
00474
00475 friend matrix<Type> sqrt (const matrix<Type> &A)
00476 {
00477 #ifdef __DEBUGMAT2
00478 matDEBUG.print("sqrt");
00479 #endif
00480
00481 matrix<Type> Result(A.lines(),A.pixels());
00482
00483 Type *pntA = A.data[0];
00484 Type *pntR = Result.data[0];
00485 for (register int32 i=0; i<Result.nsize; i++)
00486 (*pntR++) = (sqrt((*pntA++)));
00487 return Result;
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506 friend void readfile(matrix<Type> &Result, const char *file,
00507 uint filelines, window win, window winoffset)
00508 {
00509 #ifdef __DEBUGMAT2
00510 matDEBUG.print("readfile (window).");
00511 #endif
00512
00513
00514
00515 win.linelo = win.linelo - winoffset.linelo + 1;
00516 win.linehi = win.linehi - winoffset.linelo + 1;
00517 win.pixlo = win.pixlo - winoffset.pixlo + 1;
00518 win.pixhi = win.pixhi - winoffset.pixlo + 1;
00519
00520
00521
00522
00523
00524
00525
00526
00527 #ifdef __NO_IOS_BINARY__
00528 ifstream ifile(file, ios::in);
00529 #else
00530 ifstream ifile(file, ios::in | ios::binary);
00531 #endif
00532 matassert(ifile,file,__FILE__,__LINE__);
00533 ifile.seekg(0,ios::end);
00534 const uint sizefile = ifile.tellg();
00535 const uint pixelsxbytes = sizefile/filelines;
00536 const uint filepixels = pixelsxbytes/sizeof(Type);
00537 const uint sizepixel = pixelsxbytes/filepixels;
00538
00539 #ifdef __DEBUGMAT2
00540 if (win.linelo<=0 || win.pixlo<=0)
00541 matERROR.print("minimum line(pixel) is 0, (should be 1?).");
00542 if (win.linelo>win.linehi || win.pixlo>win.pixhi)
00543 matERROR.print("minimum line(pixel) is larger than max.");
00544 if (win.linehi>filelines || win.pixhi>filepixels)
00545 matERROR.print("max. line (pixel) is larger then on file.");
00546 if (sizeof(Type)!=sizepixel)
00547 matERROR.print("Type on disk is different than type of matrix.");
00548 if (sizefile==0)
00549 matERROR.print("filesize==0...");
00550 #endif
00551
00552 const uint lines = win.lines();
00553 const uint pixels = win.pixels();
00554 const uint start = ((win.linelo-1)*filepixels+win.pixlo-1)*sizepixel;
00555
00556 #ifdef __DEBUGMAT1
00557 if (Result.lines() < lines || Result.pixels() < pixels)
00558 matERROR.print("readfile: matrix too small to contain window from file.");
00559 if (Result.lines() != lines || Result.pixels() != pixels)
00560 matDEBUG.print("debug info: readfile: matrix is not fully filled.");
00561 #endif
00562 for (register int32 lin=0; lin<lines; ++lin)
00563 {
00564
00565 ifile.seekg(start+filepixels*lin*sizepixel,ios::beg);
00566 ifile.read((char*)&Result.data[lin][0],pixels*sizepixel);
00567 }
00568 ifile.close();
00569 }
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 friend void writefile (
00580 ofstream &file, const matrix<Type> &tobwritten, window win)
00581 {
00582 #ifdef __DEBUGMAT2
00583 matDEBUG.print("writefile.");
00584 #endif
00585 #ifdef __DEBUGMAT1
00586 if (win.linelo>win.linehi || win.pixlo>win.pixhi)
00587 matERROR.print("minimum line(pixel) is larger than max.");
00588 if (win.linehi>=tobwritten.lines() || win.pixhi>=tobwritten.pixels())
00589 matERROR.print("window not ok with matrix.");
00590 #endif
00591 matassert(file,"writefile: file unknown",__FILE__,__LINE__);
00592 const uint pixels = win.pixels();
00593 const uint size = pixels*sizeof(Type);
00594 for (int32 line=win.linelo; line<=win.linehi; ++line)
00595 file.write((char*)&tobwritten.data[line][win.pixlo],size);
00596 }
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 friend void fftshift (matrix<Type> &A)
00607 {
00608 #ifdef __DEBUGMAT2
00609 matDEBUG.print("fftshift");
00610 #endif
00611 #ifdef __DEBUGMAT1
00612 if (!A.isvector())
00613 matERROR.print("fftshift: only vectors");
00614 #endif
00615 matrix<Type> Res(A.nrows,A.ncols);
00616 const int32 start = int32(ceil(real8(A.nsize)/2));
00617 memcpy(Res.data[0],A.data[0]+start,sizeof(Type)*(A.nsize-start));
00618 memcpy(Res.data[0]+A.nsize-start,A.data[0],sizeof(Type)*start);
00619 myswap(A,Res);
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 friend void ifftshift (matrix<Type> &A)
00631 {
00632 #ifdef __DEBUGMAT2
00633 matDEBUG.print("ifftshift");
00634 #endif
00635 #ifdef __DEBUGMAT1
00636 if (!A.isvector())
00637 matERROR.print("ifftshift: only vectors");
00638 #endif
00639 matrix<Type> Res(A.nrows,A.ncols);
00640
00641 const int32 start = int32(floor(real8(A.nsize)/2));
00642 memcpy(Res.data[0],A.data[0]+start,sizeof(Type)*(A.nsize-start));
00643 memcpy(Res.data[0]+A.nsize-start,A.data[0],sizeof(Type)*start);
00644 myswap(A,Res);
00645 }
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657 friend void wshift (matrix<Type> &A, int32 n)
00658 {
00659 #ifdef __DEBUGMAT2
00660 matDEBUG.print("wshift");
00661 #endif
00662 #ifdef __DEBUGMAT1
00663 if (abs(n)>=A.nsize)
00664 matERROR.print("wshift: shift larger than matrix not implemented.");
00665 if (!A.isvector())
00666 matERROR.print("wshift: only vectors");
00667 #endif
00668
00669 if (n==0) return;
00670 if (n<0) n += A.nsize;
00671 matrix<Type> Res(A.nrows,A.ncols);
00672
00673 memcpy(Res.data[0],A.data[0]+n,sizeof(Type)*(A.nsize-n));
00674 memcpy(Res.data[0]+A.nsize-n,A.data[0],sizeof(Type)*n);
00675 myswap(A,Res);
00676 }
00677
00678
00679 };
00680
00681
00682
00683
00684
00685
00686
00687 #include "matrixbk.cc"
00688
00689 #endif // MATRIXBK_H guard
00690