diff --git a/src/mc/mcmain1d/mr1d_gmca.cc b/src/mc/mcmain1d/mr1d_gmca.cc index b65cd6b..412fc57 100644 --- a/src/mc/mcmain1d/mr1d_gmca.cc +++ b/src/mc/mcmain1d/mr1d_gmca.cc @@ -27,6 +27,7 @@ #include "Array.h" #include "NR.h" #include "IM_Obj.h" +#include "IM_IO.h" #include "GMCA.h" #include "MR1D1D.h" #include "IM_Noise.h" @@ -34,8 +35,8 @@ #include "MR1D_Filter.h" /****************************************************************************/ - -class GMCA_1D: public MR1D1D, public GMCA + +class GMCA_1D: public MR1D1D, public GMCA { public: Bool UseRMSMap; @@ -46,7 +47,7 @@ class GMCA_1D: public MR1D1D, public GMCA void inpainting_run(fltarray &TabCannels, fltarray & InpData); void transrecons_sources(fltarray &TabVect,fltarray &Recdata); void transform_sources(fltarray &Data,fltarray &TabVect); - void recons_sources(fltarray &DataIn, fltarray &EstSources); // Applying the matrix on the data + void recons_sources(dblarray &DataIn, dblarray &EstSources); // Applying the matrix on the data void HT_Sources(fltarray &TabSource, float &KThrd) ; ~GMCA_1D() {} ; @@ -170,12 +171,12 @@ void GMCA_1D::transrecons_sources(fltarray &TabVect,fltarray &Recdata) /****************************************************************************/ -void GMCA_1D::recons_sources(fltarray &DataIn, fltarray &EstSources) // Applying the matrix on the data +void GMCA_1D::recons_sources(dblarray &DataIn, dblarray &EstSources) // Applying the matrix on the data { int Nx = DataIn.nx(); int Ny = DataIn.ny(); int i,k,l; - fltarray RefData,RefSources; + dblarray RefData,RefSources; int Deb = 0; // cout << "NEW recons_sources " << endl; @@ -588,7 +589,7 @@ int test_main(int argc, char *argv[]) int main(int argc, char *argv[]) { - fltarray Dat; + dblarray Dat; /* Get command line arguments, open input file(s) if necessary */ fitsstruct Header; char Cmd[512]; @@ -624,7 +625,7 @@ int main(int argc, char *argv[]) if (Verbose == True) cout << "\n Reading the data"<< endl; - fits_read_fltarr(Name_Cube_In, Dat, &Header); + fits_read_dblarr(Name_Cube_In, Dat); int Nx = Dat.nx(); int Ny = Dat.ny(); @@ -676,8 +677,8 @@ int main(int argc, char *argv[]) // WT.write(Name_Out); // Compute the 2D1D transform - fltarray TabVect; - WT.transform_to_vectarray(Dat, TabVect); + dblarray TabVect; + WT.transform_to_vectdblarray(Dat, TabVect); // fits_write_fltarr ("xx_tabvect.fits", TabVect); // Initalize the class for GMCA @@ -699,7 +700,7 @@ int main(int argc, char *argv[]) WT.GlobThrd = GThrd; WT.SVConst = UsePCA; WT.MatNbrScale1D = Nbr_Plan; - fltarray QSVec; + dblarray QSVec; if (UsePCA == True) { @@ -709,13 +710,13 @@ int main(int argc, char *argv[]) if (UseMask == True) { - fits_read_fltarr (Name_Mask, WT.Mask); + fits_read_dblarr (Name_Mask, WT.Mask); } if (UseKnownColomn == True) { - fits_read_fltarr (Name_KnowColumn, WT.MatKnownColumn); + fits_read_dblarr (Name_KnowColumn, WT.MatKnownColumn); if (WT.MatKnownColumn.naxis() == 1) WT.NbrKnownColumn = 1; else WT.NbrKnownColumn = WT.MatKnownColumn.axis(2); } @@ -726,7 +727,7 @@ int main(int argc, char *argv[]) if (EstimNbSources == False) // THE NUMBER OF SOURCES IS FIXED { - fltarray TabSource; + dblarray TabSource; int NbrCoef = TabVect.nx(); TabSource.alloc(NbrCoef,NbrSources); WT.GMCA::Verbose = Verbose; @@ -749,7 +750,7 @@ int main(int argc, char *argv[]) NbrSources++; if (Verbose == True) cout << "Running GMCA ... Number of Estimated Sources : " << NbrSources << endl; WT.NbrSources = NbrSources; - fltarray TabSource; + dblarray TabSource; int NbrCoef = TabVect.nx(); TabSource.alloc(NbrCoef,NbrSources); WT.GMCA::Verbose = Verbose; @@ -782,7 +783,7 @@ int main(int argc, char *argv[]) // Reconstruction : if (Verbose == True) cout << "Reconstruction ... "<< endl; - fltarray EstSources; + dblarray EstSources; // cout << "GO REC" << endl; WT.recons_sources(Dat,EstSources); @@ -795,16 +796,16 @@ int main(int argc, char *argv[]) // fits_write_fltarr ("xx_InvMixingMat.fits", WT.InvMixingMat); // Header.origin = Cmd; - fits_write_fltarr(Name_Out, EstSources); - if (WriteMixing == True) fits_write_fltarr (Name_Out_2, WT.RecMixingMat); - if (WriteChannels == True) + fits_write_dblarr(Name_Out, EstSources); + if (WriteMixing == True) fits_write_dblarr (Name_Out_2, WT.RecMixingMat); + if (WriteChannels == True) { - fltarray EstChannels, TranspMixingMat; + dblarray EstChannels, TranspMixingMat; MatOper MAT; // See file $Tools/MatrixOper.cc and .h MAT.transpose(WT.MixingMat,TranspMixingMat); WT.apply_mat(EstSources, TranspMixingMat, EstChannels); // WT.apply_mat(EstSources, WT.MixingMat, EstChannels); - fits_write_fltarr (Name_Out_3, EstChannels); + fits_write_dblarr (Name_Out_3, EstChannels); } exit(0); } diff --git a/src/mc/mcmain2d/mr_gmca.cc b/src/mc/mcmain2d/mr_gmca.cc index c5baeb0..96d950a 100644 --- a/src/mc/mcmain2d/mr_gmca.cc +++ b/src/mc/mcmain2d/mr_gmca.cc @@ -229,20 +229,18 @@ static void usage(char *argv[]) fprintf(OUTMAN, " Apply a l_1 constraint also on the mixing matrix. Default is no. \n"); fprintf(OUTMAN, " [-d]\n"); fprintf(OUTMAN, " Estimate the number of sources. Default is no. \n"); - fprintf(OUTMAN, " [-d]\n"); - fprintf(OUTMAN, " Estimate the number of sources. Default is no. \n"); fprintf(OUTMAN, " [-m]\n"); fprintf(OUTMAN, " Mad-based stopping criterion when the number of sources is estimated. Default is 5 - default criterion is l2-based. \n"); fprintf(OUTMAN, " [-L]\n"); fprintf(OUTMAN, " L2-based stopping criterion when the number of sources is estimated. Default is 40 (in dB). \n"); - fprintf(OUTMAN, " [-D] \n"); - fprintf(OUTMAN, " Spectra with disjoint supports for thresholds higher than 7 Mad. \n"); // Should be an option +// fprintf(OUTMAN, " [-D] \n"); +// fprintf(OUTMAN, " Spectra with disjoint supports for thresholds higher than 7 Mad. \n"); // Should be an option fprintf(OUTMAN, " [-K Last K-Mad]\n"); fprintf(OUTMAN, " Last value of K for K-Mad Thresholding. \n"); fprintf(OUTMAN, " [-G Global Thresholding]\n"); - fprintf(OUTMAN, " [-O]\n"); - fprintf(OUTMAN, " Orthogonalization of the spectra\n"); - verbose_usage(); +// fprintf(OUTMAN, " [-O]\n"); +// fprintf(OUTMAN, " Orthogonalization of the spectra\n"); + verbose_usage(); vm_usage(); manline(); exit(-1); @@ -435,7 +433,7 @@ static void transinit(int argc, char *argv[]) int main(int argc, char *argv[]) { - fltarray Dat; + dblarray Dat; /* Get command line arguments, open input file(s) if necessary */ fitsstruct Header; char Cmd[512]; @@ -466,38 +464,44 @@ int main(int argc, char *argv[]) if (Verbose == True) cout << "\n Reading the data"<< endl; - io_3d_read_data(Name_Cube_In, Dat, &Header); - + // io_3d_read_data(Name_Cube_In, Dat, &Header); + fits_read_dblarr (Name_Cube_In, Dat); + int Nx = Dat.nx(); int Ny = Dat.ny(); int Nz = Dat.nz(); - if (Verbose == True) cout << "Nx = " << Dat.nx() << " Ny = " << Dat.ny() << " Nz = " << Dat.nz() << endl; - + if (Verbose == True) + cout << "Nx = " << Dat.nx() << " Ny = " << Dat.ny() << " Nz = " << Dat.nz() << endl; + // Dat.info("READ data"); + + if (Normalize == True) { double Mean = Dat.mean(); double Sigma = Dat.sigma(); + // cout << "Data mean = " << Mean << " Sigma = " << Sigma << endl; + // printf("Sigmad=f = %f", Sigma); for (int i=0;i ErrRate) L0n++; return L0n; @@ -982,11 +1010,11 @@ float GMCA::calc_L0norm(fltarray &TabSource) /****************************************************************************/ -float GMCA::l2_error(fltarray & TabData, fltarray & TabSource) -{ - fltarray TabChannel; - fltarray TranspMixingMat; - fltarray Mat= MixingMat; +double GMCA::l2_error(dblarray & TabData, dblarray & TabSource) +{ + dblarray TabChannel; + dblarray TranspMixingMat; + dblarray Mat= MixingMat; MAT.transpose(Mat,TranspMixingMat); MAT.mat_mult(TranspMixingMat, TabSource, TabChannel); TabChannel -= TabData; @@ -995,15 +1023,15 @@ float GMCA::l2_error(fltarray & TabData, fltarray & TabSource) /****************************************************************************/ -float GMCA::CountResiThrd(fltarray & TabData, fltarray & TabSource, float & Thrd) -{ - fltarray TabChannel; - fltarray TranspMixingMat; - fltarray Mat= MixingMat; +double GMCA::CountResiThrd(dblarray & TabData, dblarray & TabSource, float & Thrd) +{ + dblarray TabChannel; + dblarray TranspMixingMat; + dblarray Mat= MixingMat; MAT.transpose(Mat,TranspMixingMat); MAT.mat_mult(TranspMixingMat, TabSource, TabChannel); TabChannel -= TabData; - float count=0; + double count=0; int i,j; for (i=0; i Thrd) count++; return count; @@ -1011,17 +1039,17 @@ float GMCA::CountResiThrd(fltarray & TabData, fltarray & TabSource, float & Thrd /****************************************************************************/ -float GMCA::MyKurt(fltarray & TabData) -{ +double GMCA::MyKurt(dblarray & TabData) +{ int Nx = TabData.nx(); int Ny = TabData.ny(); - float KurtV,MeanVal,SigVal; + double KurtV,MeanVal,SigVal; int i,j,p=0; - fltarray V(Ny); + dblarray V(Ny); int NCol = Nx* Ny; - fltarray W(NCol); - - for (i=0; i 0) // Initialize with a known sub-matrix { if (L1_Matrix_Constraint == True) @@ -1133,7 +1162,7 @@ void GMCA::init_mixmat(fltarray &TabCannels) // Devrait être fait dans le const if (NbrColumnInit > 1) { // Transform the known spectra - fltarray TempM(NbrColumnInit,NbrCannels); + dblarray TempM(NbrColumnInit,NbrCannels); MAT.transpose(MatColumnInit,TempM); mat_wt1d_trans(TempM); MAT.transpose(TempM,MatColumnInit); @@ -1141,7 +1170,7 @@ void GMCA::init_mixmat(fltarray &TabCannels) // Devrait être fait dans le const if (NbrColumnInit == 1) { - for (k=0;k 1) { // Transform the known spectra - fltarray TempM(NbrKnownColumn,NbrCannels); + dblarray TempM(NbrKnownColumn,NbrCannels); MAT.transpose(MatKnownColumn,TempM); mat_wt1d_trans(TempM); MAT.transpose(TempM,MatKnownColumn); @@ -1192,7 +1221,8 @@ void GMCA::init_mixmat(fltarray &TabCannels) // Devrait être fait dans le const } - + // cout << " IN init_mixmat " << NbrKnownColumn << endl; + { for (i=0;i < NbrKnownColumn;i++) { @@ -1200,14 +1230,15 @@ void GMCA::init_mixmat(fltarray &TabCannels) // Devrait être fait dans le const if (NbrKnownColumn == 1) for (k=0;k < NbrCannels;k++) Temp(k,i)=MatKnownColumn(k); } } - + // cout << " IN init_mixmat " << NbrKnownColumn << endl; + MixingMat = Temp; RecMixingMat = Temp; } /****************************************************************************/ -void GMCA::mat_wt1d_trans(fltarray &Mat) // Il vaudrait mieux mettre un fltarray en entree car ce sont les colonnes de TabData que nous allons transformer +void GMCA::mat_wt1d_trans(dblarray &Mat) // Il vaudrait mieux mettre un dblarray en entree car ce sont les colonnes de TabData que nous allons transformer { int i,k,s; fltarray Vect(NbrCannels); @@ -1226,7 +1257,7 @@ void GMCA::mat_wt1d_trans(fltarray &Mat) // Il vaudrait mieux mettre un fltarra { for (k=0; k < MR_Mat.size_scale_np (s); k++) { - Mat(i,ind+k) = MR_Mat(s,k); + Mat(i,ind+k) = (double) MR_Mat(s,k); } @@ -1237,7 +1268,7 @@ void GMCA::mat_wt1d_trans(fltarray &Mat) // Il vaudrait mieux mettre un fltarra /****************************************************************************/ -void GMCA::mat_wt1d_recons(fltarray &Mat) +void GMCA::mat_wt1d_recons(dblarray &Mat) { int i,k; fltarray Vect(NbrCannels); @@ -1246,21 +1277,21 @@ void GMCA::mat_wt1d_recons(fltarray &Mat) { int ind=0; for (int s=0; s < MR_Mat.nbr_band (); s++) - for (k=0; k < MR_Mat.size_scale_np (s); k++) MR_Mat(s,k) = Mat(ind++,i); + for (k=0; k < MR_Mat.size_scale_np (s); k++) MR_Mat(s,k) = (float) Mat(ind++,i); MR_Mat.recons(Vect); - for (k=0;k X^T MAT.mat_mult(Mat,transpM,cM); // X -> X^T*X @@ -1274,12 +1305,12 @@ void GMCA::leftpi_mat(fltarray &Mat,fltarray &LPIMat) // compute the inverse ma /****************************************************************************/ -void GMCA::rightpi_mat(fltarray &Mat, fltarray &RPIMat) // compute the inverse matrix (should be great to define a similar function to compute pseudo-inverses (left and right)) +void GMCA::rightpi_mat(dblarray &Mat, dblarray &RPIMat) // compute the inverse matrix (should be great to define a similar function to compute pseudo-inverses (left and right)) { - fltarray transpM; - fltarray cM; - fltarray IcM; - fltarray tempA; + dblarray transpM; + dblarray cM; + dblarray IcM; + dblarray tempA; MAT.transpose(Mat,transpM); MAT.mat_mult(Mat,transpM,cM); @@ -1292,11 +1323,11 @@ void GMCA::rightpi_mat(fltarray &Mat, fltarray &RPIMat) // compute the inverse /******************** CalcMadMat ******************************/ -float GMCA::CalcMadMat(fltarray &TabSource) +double GMCA::CalcMadMat(dblarray &TabSource) { int i,k; int Nxd = TabSource.nx(),Nyd = TabSource.ny(); - fltarray V(Nxd*Nyd); + dblarray V(Nxd*Nyd); int Val = 0; for (i=0;i < Nxd;i++) @@ -1312,13 +1343,13 @@ float GMCA::CalcMadMat(fltarray &TabSource) /******************** Hard Threshold the sources******************************/ -void GMCA::HT_Sources(fltarray &TabSource, float &KThrd) // Applying the matrix on the data +void GMCA::HT_Sources(dblarray &TabSource, float &KThrd) // Applying the matrix on the data { int NbrCoef = TabSource.nx(); int Nyd = TabSource.ny(); int i,k; float Mad; - fltarray V(NbrCoef); + dblarray V(NbrCoef); float Thrd; if (GlobThrd == True) @@ -1350,13 +1381,13 @@ void GMCA::HT_Sources(fltarray &TabSource, float &KThrd) // Applying the matrix /******************** Hard Threshold the spectra******************************/ -void GMCA::HT_MixMat(fltarray &Mat, float &KThrd) // Applying the matrix on the data +void GMCA::HT_MixMat(dblarray &Mat, float &KThrd) // Applying the matrix on the data { int NLin_Mat = Mat.axis(1); int NCol_Mat = Mat.axis(2); int i,k; float Mad; - fltarray V; + dblarray V; V.alloc(1,NCol_Mat-KStart); @@ -1377,12 +1408,12 @@ void GMCA::HT_MixMat(fltarray &Mat, float &KThrd) // Applying the matrix on the /****************************************************************************/ -void GMCA::apply_mat(fltarray & TabData, fltarray &Mat, fltarray & Result) // Applying the matrix on the data +void GMCA::apply_mat(dblarray & TabData, dblarray &Mat, dblarray & Result) // Applying the matrix on the data { int NbrCoef = TabData.nx(); int Nyd = TabData.ny(); int i,k,NCol_Mat = Mat.axis(2); - fltarray V,R; + dblarray V,R; V.alloc(1,Nyd); R.alloc(1,NCol_Mat); @@ -1401,13 +1432,13 @@ void GMCA::apply_mat(fltarray & TabData, fltarray &Mat, fltarray & Result) // A /******************** Hard Threshold the sources******************************/ -void GMCA::recons_sources(fltarray &DataIn, fltarray &EstSources) // Applying the matrix on the data +void GMCA::recons_sources(dblarray &DataIn, dblarray &EstSources) // Applying the matrix on the data { int Nx = DataIn.nx(); int Ny = DataIn.ny(); int Nz = DataIn.nz(); int i,k,l; - fltarray RefData,RefSources; + dblarray RefData,RefSources; int Deb = 0; // Reform the data @@ -1442,7 +1473,7 @@ void GMCA::recons_sources(fltarray &DataIn, fltarray &EstSources) // Applying t /**********************************************************************************/ - void GMCA::apply_mask(fltarray &Data,fltarray &Mask,fltarray &DataRec) + void GMCA::apply_mask(dblarray &Data,dblarray &Mask,dblarray &DataRec) { int i,k,s; for (i=0;i 7 - void force_disjsources(fltarray &TabSource); + void force_disjsources(dblarray &TabSource); // Force the sources to have disjoint supports for KThrd > 7 - void apply_mat(fltarray & TabData,fltarray & Mat, fltarray & TabRes); + void apply_mat(dblarray & TabData,dblarray & Mat, dblarray & TabRes); // Apply the matrix to the data - void mat_wt1d_trans(fltarray &Mat); + void mat_wt1d_trans(dblarray &Mat); // 1D wavelet transform along the column of the matrix MixingMat - void mat_wt1d_recons(fltarray &Mat); + void mat_wt1d_recons(dblarray &Mat); // 1D wavelet reconstruction along the column of the matrix MixingMat - void recons_sources(fltarray &DataIn, fltarray &EstSources); + void recons_sources(dblarray &DataIn, dblarray &EstSources); // reconstruct the estimated sources from the mixing matrix estimate and the input data - void init_TabSource(fltarray &TabSource,int & NbrSamples); + void init_TabSource(dblarray &TabSource,int & NbrSamples); // Initialize TabSource - void ortho_hypcube(fltarray &TabSource); - // Orthogonalize the hyperspectral source cubes + void ortho_hypcube(dblarray &TabSource); + // Orthogonalize the hyperspectral source cubes - void positive_cube_constraint(fltarray &Data); + void positive_cube_constraint(dblarray &Data); // Constrain a cube to be positive void positive_mixmat_constraint(); @@ -200,27 +200,27 @@ class GMCA { void PCAMixMat(); // Constrain the subspace spaned by Mixing Mat - void recons_data(fltarray &TabCannels,fltarray &Data); + void recons_data(dblarray &TabCannels,dblarray &Data); // To reconstruct the data - void reform_data(fltarray &TabSources,fltarray &Data); + void reform_data(dblarray &TabSources,dblarray &Data); // To reconstruct the data - void transform_data(fltarray &DataRec,fltarray &TabCannels); + void transform_data(dblarray &DataRec,dblarray &TabCannels); - void apply_mask(fltarray &Data,fltarray &Mask,fltarray &DataRec); + void apply_mask(dblarray &Data,dblarray &Mask,dblarray &DataRec); - void run_gmca(fltarray &TabCannels, fltarray & TabSource); - // run the gmca method + void run_gmca(dblarray &TabCannels, dblarray & TabSource); + // run the gmca method - virtual void run(fltarray &TabCannels, fltarray & TabSource, fltarray & MixingMat, fltarray & InvMixingMat) {}; + virtual void run(dblarray &TabCannels, dblarray & TabSource, dblarray & MixingMat, dblarray & InvMixingMat) {}; // run the gmca method calling run_gmca - virtual void transrecons_sources(fltarray &TabVect,fltarray &Recdata) {}; + virtual void transrecons_sources(dblarray &TabVect,dblarray &Recdata) {}; - virtual void transform_sources(fltarray &Data,fltarray &TabVect) {}; + virtual void transform_sources(dblarray &Data,dblarray &TabVect) {}; - virtual void inpainting_run(fltarray &TabCannels, fltarray & InpData){}; + virtual void inpainting_run(dblarray &TabCannels, dblarray & InpData){}; // run the inpainting method calling run_gmca virtual ~GMCA(){}; diff --git a/src/sparse/libsparse2d/MR1D1D.cc b/src/sparse/libsparse2d/MR1D1D.cc index bc2171a..317e0dc 100644 --- a/src/sparse/libsparse2d/MR1D1D.cc +++ b/src/sparse/libsparse2d/MR1D1D.cc @@ -417,6 +417,56 @@ void MR1D1D::transform_to_vectarray(fltarray &Data, fltarray &TabVect) fltarray Vect(Ny); intarray TabPosBand(NbrBandX); + // 2D wt transform per frame + for (y=0; y < Ny; y++) + { + int Pix=0; + for (i=0; i < Nx; i++) Frame(i) = Data(i,y); + WT_x.transform(Frame); + for (b=0; b < NbrBandX; b++) + { + for (i=0; i < WT_x.size_scale_np (b); i++) TabVect(Pix++,y) = WT_x(b,i); + if (y == 0) TabPosBand(b) = Pix; + } + } +#ifdef DB_MR1D1D + cout << "1D_Y " << NbrBandY << endl; +#endif + // 1D wt + if (NbrBandY >= 2) + { + for (b=0; b < NbrBandX; b++) + for (i=0; i < WT_x.size_scale_np (b); i++) + { + for (y=0; y< Ny; y++) Vect(y) = TabBand[b](i,y); TabVect(TabPosBand(b)+i*WT_x.size_scale_np(b)+i, y); + WT_y.transform(Vect); + y = 0; + for (int b1=0; b1 < NbrBandY; b1++) + for (int p=0; p < WT_y.size_scale_np (b1); p++) TabVect(TabPosBand(b)+i* WT_x.size_scale_np(b)+i,y) = WT_y(b1,p); + } + } +#ifdef DB_MR1D1D + cout << "END TRANS " << endl; +#endif +} + +/****************************************************************************/ + +void MR1D1D::transform_to_vectdblarray(dblarray &Data, dblarray &TabVect) +{ + int Nsx = nbr_coef_x(); + int Nsy = nbr_coef_y(); +#ifdef DB_MR1D1D + cout << "ALLOC transform_to_vectarray " << Nsx << " " << Nsy << endl; +#endif + TabVect.alloc(nbr_coef_x(), nbr_coef_y()); + int x,b,y,i; + Nx = Data.nx(); + Ny = Data.ny(); + fltarray Frame(Nx); + fltarray Vect(Ny); + intarray TabPosBand(NbrBandX); + // 2D wt transform per frame for (y=0; y < Ny; y++) { diff --git a/src/sparse/libsparse2d/MR1D1D.h b/src/sparse/libsparse2d/MR1D1D.h index 1dc9735..07dcbd5 100755 --- a/src/sparse/libsparse2d/MR1D1D.h +++ b/src/sparse/libsparse2d/MR1D1D.h @@ -90,8 +90,10 @@ class MR1D1D { void read(char *Name); // Read a transformation from a file - + void transform_to_vectarray(fltarray &Data, fltarray &TabVect); + + void transform_to_vectdblarray(dblarray &Data, dblarray &TabVect); // Apply the 2D-1D transfrom, but the transformation is not stored in TabBand, and // reconstruction is possible // get_band, put_band, read, write routines cannot be used diff --git a/src/sparse/libsparse3d/MR2D1D.cc b/src/sparse/libsparse3d/MR2D1D.cc index e9d1f74..b042c91 100644 --- a/src/sparse/libsparse3d/MR2D1D.cc +++ b/src/sparse/libsparse3d/MR2D1D.cc @@ -432,7 +432,57 @@ void MR2D1D::transform_to_vectarray(fltarray &Data, fltarray &TabVect) { int Nsx = nbr_coef_2d(); int Nsy = nbr_coef_1d(); - cout << "ALLOC " << Nsx << " " << Nsy << endl; + // cout << "ALLOC " << Nsx << " " << Nsy << endl; + TabVect.alloc(nbr_coef_2d(), nbr_coef_1d()); + int i,j,b,z; + Nx = Data.nx(); + Ny = Data.ny(); + Nz = Data.nz(); + Ifloat Frame(Ny, Nx); + fltarray Vect(Nz); + intarray TabPosBand(nbr_band_2d()); + + // 2D wt transform per frame + for (z=0; z < Nz; z++) + { + int Pix=0; + for (i=0; i < Ny; i++) + for (j=0; j < Nx; j++) Frame(i,j) = Data(j,i,z); + WT2D.transform(Frame); + for (b=0; b < nbr_band_2d(); b++) + { + for (i=0; i < WT2D.size_band_nl(b); i++) + for (j=0; j < WT2D.size_band_nc(b); j++) TabVect(Pix++,z) = WT2D(b,i,j); + if (z == 0) TabPosBand(b) = Pix; + } + } + // cout << "1D " << NbrBand1D << endl; + // 1D wt + if (NbrBand1D >= 2) + { + for (b=0; b < nbr_band_2d(); b++) + for (i=0; i < WT2D.size_band_nl(b); i++) + for (j=0; j < WT2D.size_band_nc(b); j++) + { + for (z=0; z < Nz; z++) Vect(z) = TabVect(TabPosBand(b)+i*WT2D.size_band_nc(b)+j,z); // TabBand[b](j,i,z); + WT1D.transform(Vect); + z = 0; + for (int b1=0; b1 < NbrBand1D; b1++) + { + for (int p=0; p < WT1D.size_scale_np (b1); p++) TabVect(TabPosBand(b)+i*WT2D.size_band_nc(b)+j,z) = WT1D(b1,p); + } + } + } + // cout << "END TRANS " << endl; +} + +/****************************************************************************/ + +void MR2D1D::transform_to_vectdblarray(dblarray &Data, dblarray &TabVect) +{ + int Nsx = nbr_coef_2d(); + int Nsy = nbr_coef_1d(); + // cout << "ALLOC " << Nsx << " " << Nsy << endl; TabVect.alloc(nbr_coef_2d(), nbr_coef_1d()); int i,j,b,z; Nx = Data.nx(); @@ -456,8 +506,8 @@ void MR2D1D::transform_to_vectarray(fltarray &Data, fltarray &TabVect) if (z == 0) TabPosBand(b) = Pix; } } - cout << "1D " << NbrBand1D << endl; - // 1D wt + // cout << "1D " << NbrBand1D << endl; + // 1D wt if (NbrBand1D >= 2) { for (b=0; b < nbr_band_2d(); b++) @@ -473,7 +523,7 @@ void MR2D1D::transform_to_vectarray(fltarray &Data, fltarray &TabVect) } } } - cout << "END TRANS " << endl; +// cout << "END TRANS " << endl; } diff --git a/src/sparse/libsparse3d/MR2D1D.h b/src/sparse/libsparse3d/MR2D1D.h index 938fbca..bf212ce 100755 --- a/src/sparse/libsparse3d/MR2D1D.h +++ b/src/sparse/libsparse3d/MR2D1D.h @@ -95,6 +95,7 @@ class MR2D1D { void read(char *Name); // Read a transformation from a file + void transform_to_vectdblarray(dblarray &Data, dblarray &TabVect); void transform_to_vectarray(fltarray &Data, fltarray &TabVect); // Apply the 2D-1D transfrom, but the transformation is not stored in TabBand, and // reconstruction is possible diff --git a/src/sparse/libtools/DefMath.cc b/src/sparse/libtools/DefMath.cc index b08f767..6794deb 100644 --- a/src/sparse/libtools/DefMath.cc +++ b/src/sparse/libtools/DefMath.cc @@ -282,6 +282,22 @@ float get_sigma_mad(float *Image, int N) /***************************************************************************/ +double get_sigma_mad(double *Image, int N) +{ + double Noise, Med; + int i; + dblarray Buff(N); + + for (i = 0; i < N; i++) Buff(i) = Image[i]; + Med = get_median(Buff.buffer(), N); + + for (i = 0; i < N; i++) Buff(i) = ABS(Image[i] - Med); + Noise = get_median(Buff.buffer(), N); + return (Noise/0.6745); +} + +/***************************************************************************/ + float get_sigma_clip (float *Data, int N, int Nit, Bool Average_Non_Null, Bool UseBadPixel, float BadPVal) { @@ -423,7 +439,36 @@ double curtosis(float *Dat, int N) else Curt = 0.; return Curt; } +/****************************************************************************/ +double curtosis(double *Dat, int N) +{ + double Curt; + double x1,x2,x3,x4,Sigma; + int i; + x1=x2=x3=x4=0.; + for (i=0; i < N; i++) + { + x1 += Dat[i]; + x2 += pow((double) Dat[i], (double) 2.); + x3 += pow((double) Dat[i], (double) 3.); + x4 += pow((double) Dat[i], (double) 4.); + } + x1 /= (double) N; + x2 /= (double) N; + x3 /= (double) N; + x4 /= (double) N; + Sigma = x2 - x1*x1; + if (Sigma > 0.) + { + double x1_2 = x1*x1; + double x1_4 = x1_2*x1_2; + Sigma = sqrt(Sigma); + Curt = 1. / pow(Sigma,(double) 4.) * (x4 -4*x1*x3 + 6.*x2*x1_2 -3.*x1_4 ) - 3.; + } + else Curt = 0.; + return Curt; +} /****************************************************************************/ // old version diff --git a/src/sparse/libtools/DefMath.h b/src/sparse/libtools/DefMath.h index 2ee1020..ecae17a 100644 --- a/src/sparse/libtools/DefMath.h +++ b/src/sparse/libtools/DefMath.h @@ -238,11 +238,14 @@ float get_random(); double b3_spline (double x); double entropy (float *Data, int Npix, float StepHisto=1.); float get_sigma_mad(float *Data, int N); -float get_sigma_clip(float *Data, int N, int Nit=3, Bool Average_Non_Null=True, +double get_sigma_mad(double *Data, int N); + +float get_sigma_clip(float *Data, int N, int Nit=3, Bool Average_Non_Null=True, Bool UseBadPixel=False, float BadPVal=0.); double skewness(float *Dat, int N); double curtosis(float *Dat, int N); -void moment4(float *Dat, int N, double &Mean, double &Pow, +double curtosis(double *Dat, int N); +void moment4(float *Dat, int N, double &Mean, double &Pow, double &M3, double & M4, float & Min, float & Max, bool not_centered=false); void moment4(double *Dat, int N, double &Mean, double &Sigma, double &Skew, double & Curt, double & Min, double & Max); diff --git a/src/sparse/libtools/OptMedian.cc b/src/sparse/libtools/OptMedian.cc index 030a2c1..4576b13 100644 --- a/src/sparse/libtools/OptMedian.cc +++ b/src/sparse/libtools/OptMedian.cc @@ -27,3 +27,7 @@ #undef elem_type #define elem_type float #include "GenMedian.h" + +#undef elem_type +#define elem_type double +#include "GenMedian.h" diff --git a/src/sparse/libtools/OptMedian.h b/src/sparse/libtools/OptMedian.h index 8128fdc..7de30df 100644 --- a/src/sparse/libtools/OptMedian.h +++ b/src/sparse/libtools/OptMedian.h @@ -38,10 +38,13 @@ float opt_med7(float *p); float opt_med9(float *p); float kth_smallest(float a[], int n, int k); float get_median(float a[], int n); +double get_median(double a[], int n); + float abs_kth_smallest(float a[], int n, int k); float get_abs_median(float a[], int n); int hmedian(int *ra, int n); float hmedian(float *ra, int n); +double hmedian(double *ra, int n); #endif diff --git a/src/sparse/libtools/TempArray.h b/src/sparse/libtools/TempArray.h index 9bd0bc5..7a1e5a4 100644 --- a/src/sparse/libtools/TempArray.h +++ b/src/sparse/libtools/TempArray.h @@ -1211,7 +1211,7 @@ double to_array::sigma () const { ad_val = po_Buffer[i] - ao_moy; ad_sigma += ad_val*ad_val; } - if ((ad_sigma /= i_NbElem) > 1e-07) ad_sigma = sqrt (ad_sigma); + if ((ad_sigma /= i_NbElem) > 0) ad_sigma = sqrt (ad_sigma); else ad_sigma = 0.; return ad_sigma; }