From 120730e415d6f8e297db06c6c22027ccefaa8cbe Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 17:36:29 +0100 Subject: [PATCH 1/6] added optional sparse3d build --- CMakeLists.txt | 30 + src/libsparse3d/Atrou3D.cc | 648 ++++++++++++ src/libsparse3d/Atrou3D.h | 104 ++ src/libsparse3d/Atrou3DFil.cc | 180 ++++ src/libsparse3d/Atrou3DFil.h | 68 ++ src/libsparse3d/CUBE_Segment.cc | 487 +++++++++ src/libsparse3d/CUBE_Segment.h | 30 + src/libsparse3d/MR2D1D.cc | 566 +++++++++++ src/libsparse3d/MR2D1D.h | 139 +++ src/libsparse3d/MR3D_Obj.cc | 1338 +++++++++++++++++++++++++ src/libsparse3d/MR3D_Obj.h | 247 +++++ src/libsparse3d/Mr3d_FewEvent.cc | 1578 ++++++++++++++++++++++++++++++ src/libsparse3d/Mr3d_FewEvent.h | 144 +++ src/mr2d1d_stat.cc | 406 ++++++++ src/mr2d1d_trans.cc | 799 +++++++++++++++ src/mr3d_filter.cc | 474 +++++++++ src/mr3d_recons.cc | 156 +++ src/mr3d_stat.cc | 465 +++++++++ src/mr3d_trans.cc | 310 ++++++ 19 files changed, 8169 insertions(+) create mode 100644 src/libsparse3d/Atrou3D.cc create mode 100644 src/libsparse3d/Atrou3D.h create mode 100644 src/libsparse3d/Atrou3DFil.cc create mode 100644 src/libsparse3d/Atrou3DFil.h create mode 100755 src/libsparse3d/CUBE_Segment.cc create mode 100755 src/libsparse3d/CUBE_Segment.h create mode 100644 src/libsparse3d/MR2D1D.cc create mode 100755 src/libsparse3d/MR2D1D.h create mode 100755 src/libsparse3d/MR3D_Obj.cc create mode 100755 src/libsparse3d/MR3D_Obj.h create mode 100644 src/libsparse3d/Mr3d_FewEvent.cc create mode 100644 src/libsparse3d/Mr3d_FewEvent.h create mode 100644 src/mr2d1d_stat.cc create mode 100644 src/mr2d1d_trans.cc create mode 100644 src/mr3d_filter.cc create mode 100644 src/mr3d_recons.cc create mode 100644 src/mr3d_stat.cc create mode 100644 src/mr3d_trans.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d03bbe..6297cff 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -103,6 +103,19 @@ project(sparse2d) add_dependencies(sparse2d fftw3) endif(USE_FFTW) + option(SPARSE3D "Build Sparse3D library" ON) + if(SPARSE3D) + # Build sparse3d library + FILE(GLOB src_lib2 "${PROJECT_SOURCE_DIR}/src/libsparse3d/*.cc") + include_directories("${PROJECT_SOURCE_DIR}/src/libsparse3d") + add_library(sparse3d STATIC ${src_lib2}) + target_link_libraries(sparse3d ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) + if(USE_FFTW) + add_dependencies(sparse3d fftw3) + endif(USE_FFTW) + endif(SPARSE3D) + message(STATUS "Sparse3D Build: ${SPARSE3D}") + # Build mga2d library FILE(GLOB src_mgalib2 "${PROJECT_SOURCE_DIR}/src/libmga2d/*.cc") include_directories("${PROJECT_SOURCE_DIR}/src/libmga2d") @@ -115,10 +128,17 @@ project(sparse2d) # Compile and link executables set(BINMR2D mr_transform mr_recons mr_filter mr_deconv mr_info cur_contrast cur_deconv cur_filter cur_stat cur_trans) + if(SPARSE3D) + set(BINMR2D ${BINMR2D} mr3d_trans mr3d_filter mr3d_recons mr3d_stat + mr2d1d_trans) + endif(SPARSE3D) foreach(program ${BINMR2D}) add_executable(${program} ${PROJECT_SOURCE_DIR}/src/${program}.cc) target_link_libraries(${program} mga2d sparse2d sparse1d tools) + if(SPARSE3D) + target_link_libraries(${program} sparse3d) + endif(SPARSE3D) endforeach(program) # Install headers @@ -130,9 +150,16 @@ project(sparse2d) INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libtools/*.h") INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) + if(SPARSE3D) + FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libsparse3d/*.h") + INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) + endif(SPARSE3D) # Install library INSTALL(TARGETS sparse1d sparse2d mga2d tools DESTINATION lib) + if(SPARSE3D) + INSTALL(TARGETS sparse3d DESTINATION lib) + endif(SPARSE3D) # install executables INSTALL(TARGETS ${BINMR2D} DESTINATION bin) @@ -148,5 +175,8 @@ project(sparse2d) foreach(program ${UNIT_TESTS}) add_executable(${program} ${PROJECT_SOURCE_DIR}/tests/${program}.cpp) target_link_libraries(${program} mga2d sparse2d sparse1d tools) + if(SPARSE3D) + target_link_libraries(${program} sparse3d) + endif(SPARSE3D) add_test(${program} ${program}) endforeach(program) diff --git a/src/libsparse3d/Atrou3D.cc b/src/libsparse3d/Atrou3D.cc new file mode 100644 index 0000000..c5857ec --- /dev/null +++ b/src/libsparse3d/Atrou3D.cc @@ -0,0 +1,648 @@ +/******************************************************************************* +** +** DESCRIPTION Sub-Band decomposition with A trous algorithm +** ----------- +** +******************************************************************************/ + + +#include "Atrou3D.h" + + +/****************************************************************************/ +// A TROUS ALGORITHM 3D +/****************************************************************************/ + + +/****************************************************************************/ +// filtering by the filter h +// compute Vi space from Vi-1 +extern Bool Verbose; + +static float TabSigma_Atrou3D[6]= + {0.956546, 0.120909, 0.0354989, 0.0122210, 0.00442556, 0.00222700}; +static float TabSigma_Atrou3D_adj[6]= + {0.982969, 0.133002, 0.0399047, 0.0138409, 0.00506464, 0.00223828}; + +ATROUS_3D_WT::ATROUS_3D_WT () +{ + ModifiedAWT=False; + Adjoint=False; + Bord=I_MIRROR; + no_fine=false; +} + +void ATROUS_3D_WT::b3spline_filtering(fltarray & Old, fltarray & New, int s) +{ + int Nx = Old.nx(); + int Ny = Old.ny(); + int Nz = Old.nz(); + int i,j,k,Step; + double Coeff_h0 = 3. / 8.; + double Coeff_h1 = 1. / 4.; + double Coeff_h2 = 1. / 16.; + fltarray Buff1(Nx,Ny,Nz); + fltarray Buff2(Nx,Ny,Nz); + + Step = (int)(pow((double)2., (double) s) + 0.5); + //convolution in Z axe + for (i = 0; i < Nx; i ++) + for (j = 0; j < Ny; j ++) + for (k = 0; k < Nz; k ++) + { + double Val = Coeff_h0 * (double)get_pix(Old,i,j,k) + + Coeff_h1 * ((double)get_pix(Old,i,j,k-Step) + + (double)get_pix(Old,i,j,k+Step)) + + Coeff_h2 * ((double)get_pix(Old,i,j,k-2*Step) + + (double)get_pix(Old,i,j,k+2*Step)); + Buff1(i,j,k) = (float) Val; + } + + //convolution in Y axe + for (i = 0; i < Nx; i ++) + for (j = 0; j < Ny; j ++) + for (k = 0; k < Nz; k ++) + { + double Val = Coeff_h0 * (double)get_pix(Buff1,i,j,k) + + Coeff_h1 * ((double)get_pix(Buff1,i,j-Step,k) + + (double)get_pix(Buff1,i,j+Step,k)) + + Coeff_h2 * ((double)get_pix(Buff1,i,j-2*Step,k) + + (double)get_pix(Buff1,i,j+2*Step,k)); + Buff2(i,j,k) = (float) Val; + } + + + //convolution in X axe + for (i = 0; i < Nx; i ++) + for (j = 0; j < Ny; j ++) + for (k = 0; k < Nz; k ++) + { + New(i,j,k) = Coeff_h0 * (double)get_pix(Buff2,i,j,k) + + Coeff_h1 * ((double)get_pix(Buff2,i-Step,j,k) + + (double)get_pix(Buff2,i+Step,j,k)) + + Coeff_h2 * ((double)get_pix(Buff2,i-2*Step,j,k) + + (double)get_pix(Buff2,i+2*Step,j,k)); + } + +} + +/****************************************************************************/ +// A trous transformation +void ATROUS_3D_WT::transform(fltarray &Cube, fltarray * & Wavelet_Coef, int NbrScale) +{ + _NbrScale = NbrScale; + _Nx = Cube.nx(); + _Ny = Cube.ny(); + _Nz = Cube.nz(); + + Wavelet_Coef[0] = Cube; + fltarray Data_out; + if (ModifiedAWT == True) Data_out.alloc(Cube.nx(), Cube.ny(), Cube.nz()); + + //cout << "TRANS: " << Cube.min() << " " << Cube.max()<< endl; + for (int s = 0; s < NbrScale-1; s++) + { + //compute the signal in Vs+1 space: Signal(Vs+1)=h*Signal(Vs) + b3spline_filtering(Wavelet_Coef[s], Wavelet_Coef[s+1], s); + if (ModifiedAWT == True) + { + b3spline_filtering(Wavelet_Coef[s+1], Data_out, s); + Wavelet_Coef[s] -= Data_out; + } + else Wavelet_Coef[s] -= Wavelet_Coef[s+1]; + } + + if(no_fine) + Wavelet_Coef[0].init(0.0); + +} + +/****************************************************************************/ + +void ATROUS_3D_WT::recons(fltarray * & Wavelet_Coef, fltarray & Data_Out, + int NbrScale, Bool AddLastScale) +{ +// if(Verbose) cerr<<"ATROUS_3D_WT::recons(.,.,"<0) _NbrScale=NbrScale; + + Data_Out.alloc(Nx,Ny,Nz); + + //cout << "REC: " << Wavelet_Coef[_NbrScale-1].min() << " " << Wavelet_Coef[_NbrScale-1].max()<< endl; + if ((ModifiedAWT == False) && (Adjoint == False)) + { + int Last_Scale_Used= (AddLastScale == True) ? _NbrScale : _NbrScale-1; + Data_Out = Wavelet_Coef[0]; + for (s = 1; s < Last_Scale_Used; s++) Data_Out += Wavelet_Coef[s]; + } + else + { + fltarray temp(Nx, Ny, Nz); + if (AddLastScale == True) Data_Out= Wavelet_Coef[_NbrScale-1]; + else Data_Out.init(); + for (s=_NbrScale-2; s>= 0 ; s--) + { + if (ModifiedAWT == True) + { + //cout << "BAND: " << s << " " << Wavelet_Coef[s].min() << " " << Wavelet_Coef[s].max()<< endl; + b3spline_filtering (Data_Out, temp, s); + for (int i=0; i < Nx; i++) + for (int j=0; j < Ny; j++) + for (int k=0; k < Nz; k++) + Data_Out(i,j,k) = temp(i,j,k) + (Wavelet_Coef[s])(i,j,k); + } + else//Adjoint == True + { + b3spline_filtering (Data_Out, temp, s); + temp += Wavelet_Coef[s]; + b3spline_filtering (Wavelet_Coef[s], Data_Out, s); + Data_Out += temp; + } + } + } + +// if(Verbose) cerr<<"ATROUS_3D_WT::recons(.,.,"<<_NbrScale<<","<"<<(float(Nx*Ny*Nz)-cnt)/float(Nx*Ny*Nz)<0 )-1)*lvl; + if(Verbose) cerr<<" threshold scale "<"<<(float(Nx*Ny*Nz)-cnt)/float(Nx*Ny*Nz)< zero) + { + bool single=true; + + int ii=-1; + + while ((single==true) && (ii<2)) + { + int jj=-1; + while ((single==true) && (jj<2)) + { + int kk=-1; + while ((single==true) && (kk<2)) + { + if(!(ii==0 && jj==0 && kk==0)) + if((i+ii>0) && (j+jj>0) && (k+kk>0) + && (i+ii<_Nx-1) && (j+jj<_Ny-1) && (k+kk<_Nz-1) ) + if( abs((TabBandIn[s])(i+ii,j+jj,k+kk)) > zero ) single=false; + kk++; + } + jj++; + } + ii++; + } + if(single) (TabBandIn[s])(i,j,k) = 0; + } +// if(Verbose) cerr<<"...end atrou::clean_single"< &C) +{ + // coresponds to ATROUS_3D_WT::alloc + delete [] C[0]; +} + +void iwt3d_transform(fltarray &Data, vector< fltarray* > &vTabBand, int _NbrScale, bool modAWT) +{ +cerr<<"iwt3d_transform(fltarray &Data, vector< fltarray* > &vTabBand, "<<_NbrScale<<","<set_use_modAWT(modAWT); + + atrou->alloc(TabBand, Data.nx(), Data.ny(), Data.nz(), _NbrScale); + +// Calculus + atrou->transform(Data,TabBand,_NbrScale); +// atrou->normalize_self(TabBand,false); + +// vTabBand allocation + vTabBand.resize(_NbrScale); + for(int s=0;s<_NbrScale;s++) + vTabBand[s] = &TabBand[s]; + + delete atrou; + return ; +} + +void iwt3d_recons(vector< fltarray* > &vTabBand, fltarray &Data, bool modAWT, bool adjoint) +{ +// Number of scales + int _NbrScale = vTabBand.size(); + +// iwt allocation + fltarray* TabBand; + ATROUS_3D_WT *atrou = new ATROUS_3D_WT(); + atrou->set_use_adjoint(adjoint); + atrou->set_use_modAWT(modAWT); + atrou->set_nbr_sale(_NbrScale); + int nx=vTabBand[0]->nx(), ny=vTabBand[0]->ny(), nz=vTabBand[0]->nz(); + +// TabBand allocation + TabBand = new fltarray [_NbrScale]; + for(int s=0;s<_NbrScale;s++) + TabBand[s].alloc(vTabBand[s]->buffer(), nx, ny, nz); + +// Calculus +// atrou->normalize_self(TabBand,true); + atrou->recons(TabBand,Data,_NbrScale); +// atrou->normalize_self(TabBand,false); + + delete [] TabBand; + delete atrou; +} + +void iwt3d_threshold(vector< fltarray* > &vTabBand, bool modAWT, float threshold, filter_type FilterType) +{ +// Number of scales + int _NbrScale = vTabBand.size(); + +// iwt allocation + fltarray* TabBand; + ATROUS_3D_WT *atrou = new ATROUS_3D_WT(); + atrou->set_nbr_sale(_NbrScale); + atrou->set_use_modAWT(modAWT); + int nx=vTabBand[0]->nx(), ny=vTabBand[0]->ny(), nz=vTabBand[0]->nz(); + +// TabBand allocation + TabBand = new fltarray [_NbrScale]; + for(int s=0;s<_NbrScale;s++) + TabBand[s].alloc(vTabBand[s]->buffer(), nx, ny, nz); + +// Calculus + if(FilterType==FT_HARD || FilterType==FT_SOFT) atrou->threshold(TabBand, threshold, FilterType==FT_SOFT, true);// true=normalized + else cerr<<"Unknown filtering method"<set_use_adjoint(adjoint); + atrou->set_use_modAWT(modAWT); + atrou->alloc(TabBand, Data.nx(), Data.ny(), Data.nz(), _NbrScale); + +// Calculus + atrou->transform(Data,TabBand,_NbrScale); + if(FilterType==FT_HARD || FilterType==FT_SOFT) atrou->threshold(TabBand, threshold, FilterType==FT_SOFT, false);//false=!normalized + else cerr<<"Unknown filtering method"<recons(TabBand,Recons,_NbrScale); + + delete [] TabBand;// allocated by atrou->alloc -> ATROUS_3D_WT::alloc + delete atrou; +} + diff --git a/src/libsparse3d/Atrou3D.h b/src/libsparse3d/Atrou3D.h new file mode 100644 index 0000000..821c7da --- /dev/null +++ b/src/libsparse3d/Atrou3D.h @@ -0,0 +1,104 @@ +/****************************************************************************** +** Copyright (C) 2001 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 15/10/01 +** +** File: ATrou3D.h +** +******************************************************************************* +** +** DESCRIPTION 3D wavelet transform using the a trous algorithm +** ----------- +** +** +******************************************************************************/ + +#ifndef _Atrou3D_H_ +#define _Atrou3D_H_ + +#include +#include "GlobalInc.h" +#include "SB_Filter1D.h" +#include "IM_IO.h" + +/************************************************************************/ + +class ATROUS_3D_WT { + +private: + int _Nx; + int _Ny; + int _Nz; + int _NbrScale; + inline int index(int i, int N) {return get_index(i, N, Bord);} + + bool no_fine; + +public: + // constructor & destructor + ATROUS_3D_WT (); + ~ATROUS_3D_WT() {} + + Bool ModifiedAWT; // kept public for old functions, use set_use_modAWT(bool) instead + Bool Adjoint; // kept public for old functions, use set_use_modAWT(bool) instead + + type_border Bord; // Type of border + + // get a pixel value, taking into account the border. + float get_pix(fltarray & Cube, int x, int y, int z) + {return (Cube(index(x,Cube.nx()), index(y,Cube.ny()), index(z,Cube.nz())));} + + // allocate the memory space for a wavelet transform + void alloc(fltarray * & TabBand,int Nx,int Ny,int Nz, int NbrBand) + { + TabBand = new fltarray [NbrBand]; + for (int s = 0; s < NbrBand; s++) TabBand[s].alloc (Nx, Ny, Nz); + } + + // deallocate a wavelet transform. + void free(fltarray * & TabBand, int Nbr_Plan) + { if (Nbr_Plan != 0) delete [] TabBand;} + + // B3 spline spline filtering of a cube, taking into accounts + // the steps. + void b3spline_filtering(fltarray & Old, fltarray & New, int s); + + // 3D wavelet transform + void transform(fltarray & CubeIn, fltarray * & TabBandOut, int Nbr_Plan); + + // 3D wavelet reconstruction + // Nbr_Plan==0 means it's already saved in the transform + void recons(fltarray * & TabBandIn, fltarray & CubeOut, int Nbr_Plan=0, + Bool AddLastScale=True); + + void set_no_fine(bool nf) {no_fine=nf;} + + // added functions for mr3d_atrou + dblarray TabStat; + void set_use_modAWT(bool ad=true) {ModifiedAWT=(Bool)ad;} + void set_use_adjoint(bool ad=true) {Adjoint=(Bool)ad;} + void set_nbr_sale(int ns) {_NbrScale=ns;} + void normalize_self(fltarray *TabBand, bool inverse=false); + void threshold(fltarray * & TabBandIn, float thresh, bool soft=false, bool normalize=false); + void clean_single(fltarray * & TabBandIn, float Sigma); + void extract_stat(fltarray * & TabBandIn, char* Outname); + void read(char *Name, fltarray * &TabBand, bool *NormalizeInv); + void write(char *Name, fltarray * TabBand, bool Normalize); +}; + +void iwt3d_clear(vector< fltarray* > &C); +void iwt3d_transform(fltarray &Data, vector< fltarray* > &vTabBand, int NbrScale3D, bool modAWT); +void iwt3d_recons(vector< fltarray* > &vTabBand, fltarray &Data, bool modAWT, bool adjoint); +void iwt3d_threshold(vector< fltarray* > &vTabBand, bool modAWT, float threshold, filter_type FilterType); +void iwt3d_filter(fltarray &Data, fltarray &Recons, int NbrScale3D, bool modAWT, bool adjoint, float threshold, filter_type FilterType); + +/************************************************************************/ + +#endif diff --git a/src/libsparse3d/Atrou3DFil.cc b/src/libsparse3d/Atrou3DFil.cc new file mode 100644 index 0000000..e3ca77f --- /dev/null +++ b/src/libsparse3d/Atrou3DFil.cc @@ -0,0 +1,180 @@ +/******************************************************************************* +** +** DESCRIPTION tools for filtering using the A trous algorithm +** ----------- +** +******************************************************************************/ + +#include "Atrou3DFil.h" + + +//****************************************************************************** + +void ATROUS_3D_FIL::Filtering_Structure(intarray *& Support, int LastScale,Bool InfoBand, + Bool RemoveBorder, int FistScaleDetect) +{ + //remove the structure on the border of the cube + if((InfoBand == True)|| (RemoveBorder == True)) + { + int Nx = (Support[0]).nx(); + int Ny = (Support[0]).ny(); + int Nz = (Support[0]).nz(); + intarray seg(Nx, Ny,Nz); + for(int s=FistScaleDetect; s 0 and all neighboors are equal + // to zero, then Support[b](i,j,k) is set to zero. + void No_Single_Point(int Nx,int Ny, int Nz,intarray * & Support, + FewEventPoisson FEP,int Nbr_Scale, Bool Verbose=False); + + void threshold(intarray *& Support, int NbPlan); + void Filtering_Structure(intarray *& Support,int LastScale,Bool InfoBand, + Bool RemoveBorder, int FistScaleDetect); + void wttransform(fltarray & Cube, int Nbr_Plan, + Bool WriteRes,Bool InfoBand, char Name_Cube_Out[]); + void wtrecons(fltarray & Cube, int Nbr_Plan, + Bool InfoBand,Bool AddLastScale, Bool Adjoint); + fltarray *Wavelet_Coef; +}; + +/************************************************************************/ + +void free(fltarray * TabBand, int Nbr_Plan); +void free(intarray * TabBand, int Nbr_Plan); +void alloc_array (intarray * & Nb_Event, int Nx, int Ny, int Nz, int NbrBand); +void alloc_array (fltarray * & TabBand, int Nx, int Ny, int Nz, int NbrBand); +void init (fltarray * & TabBand, int NbrBand, float val=0.); + +#endif diff --git a/src/libsparse3d/CUBE_Segment.cc b/src/libsparse3d/CUBE_Segment.cc new file mode 100755 index 0000000..65574bc --- /dev/null +++ b/src/libsparse3d/CUBE_Segment.cc @@ -0,0 +1,487 @@ +/****************************************************************************** +** +** DESCRIPTION cube segmentation +** ----------- +** +** void cube_segment (fltarray &cube, fltarray &Segment, int &NLabel, float Level) +** +** +** cube : input cube we want to segment +** Segment: output segmented cube +** NLabel: Number of labels or regions +** Level: segmentation level +** +*******************************************************************************/ +#include "CUBE_Segment.h" +#include "IM_IO.h" +#include "IM3D_IO.h" + + +//---------------------------------------------------------- +// Pile +//---------------------------------------------------------- +Pile::Pile(void) +{ + data = NULL; + num = 0; + bloc = 0; +} + +//---------------------------------------------------------- +// Pile +//---------------------------------------------------------- +Pile::Pile (int n) +{ + data = NULL; + bloc = n; + if (data != NULL) data = (int *)realloc (data, bloc * sizeof (int)); + else data = (int *) malloc (bloc * sizeof (int)); + + num = 0; +} + +//---------------------------------------------------------- +// Pile :: destructeur +//---------------------------------------------------------- +Pile::~Pile (void) +{ + if (data) free (data); +} + +//---------------------------------------------------------- +// Pile :: max retourne le plus grand ellement +//---------------------------------------------------------- +int Pile::max (void) +{ + int max = -100000; + for (int i=0 ; i data[i]) ? max : data[i]; + return max; +} + +//---------------------------------------------------------- +// Pile :: allocation de la memoire de data +//---------------------------------------------------------- +void Pile::alloue (void) +{ + bloc += PILE_SIZE_BLOC; + if (data != NULL) data = (int *)realloc (data, bloc * sizeof (int)); + else data = (int *) malloc (bloc * sizeof (int)); +} + +//---------------------------------------------------------- +// Pile :: ajoute une donnee +//---------------------------------------------------------- +void Pile::add_data (int d) +{ + if (num >= bloc) + alloue (); + data[num++] = d; +} + +//---------------------------------------------------------- +// Pile :: trie les donnes +//---------------------------------------------------------- + +void Pile::trie (int max) +{ + int i, c; + Pile histo (max); + for (i=0 ; i S) + { + Bool NewLabel= True; + int UseLabel = 0; + TabNeighbour[0] = Segment(i-1, j-1, k-1, I_ZERO); + TabNeighbour[1] = Segment(i , j-1, k-1, I_ZERO); + TabNeighbour[2] = Segment(i+1, j-1, k-1, I_ZERO); + TabNeighbour[3] = Segment(i-1, j, k-1, I_ZERO); + TabNeighbour[4] = Segment(i, j, k-1, I_ZERO); + TabNeighbour[5] = Segment(i+1, j, k-1, I_ZERO); + TabNeighbour[6] = Segment(i-1, j+1, k-1, I_ZERO); + TabNeighbour[7] = Segment(i, j+1, k-1, I_ZERO); + TabNeighbour[8] = Segment(i+1, j+1, k-1, I_ZERO); + + TabNeighbour[9] = Segment(i-1, j-1, k, I_ZERO); + TabNeighbour[10] = Segment(i, j-1, k, I_ZERO); + TabNeighbour[11] = Segment(i+1, j-1, k, I_ZERO); + TabNeighbour[12] = Segment(i-1, j, k, I_ZERO); + + // Find the label to apply to the pixel + for (int l=0; l < Nvoisin; l++) + { + if (TabNeighbour[l] != 0) + { + int L = TabNeighbour[l]; + int ColNeighbour = P.data[L]; + NewLabel = False; + if (L > n) + { + cout << "Error: Segmentation problem ... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + // TabNeighbour[l] = P.data[L]; + if (UseLabel == 0) UseLabel = ColNeighbour; + else if (UseLabel > ColNeighbour) UseLabel = ColNeighbour; + } + } + // New label + if (NewLabel == True) + { + Segment(i,j,k) = n; + P.add_data (n++); + } + else + { + Segment(i,j,k) = UseLabel; + // neighbourhood label must have the same label + // and therefore may need to be changed + for (int l=0; l < Nvoisin; l++) + { + if (TabNeighbour[l] != 0) + { + int L = TabNeighbour[l]; + int ColNeighbour = P.data[L]; + if (ColNeighbour != UseLabel) + { + if (L > n) + { + cout << "Error: Segmentation problem .... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + P.remplace (P.data[L], UseLabel); + } + } + } + } + } + } + + if (CleanBord == True) + { + // cout << "CLEAN BORD" << endl; + // plan (*,*,0) et plan (*,*,Nz-1) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Ny ; j++) + { + if (Segment(i,j,0) != 0) P.remplace((int) Segment(i,j,0) ,0); + if (Segment(i,j,Nz-1) != 0) P.remplace((int) Segment(i,j,Nz-1) ,0); + } + // plan (*,0,*) et plan (*,Ny-1,*) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Nz ; j++) + { + if (Segment(i,0,j) != 0) P.remplace((int) Segment(i,0,j) ,0); + if (Segment(i,Ny-1,j) != 0) P.remplace((int) Segment(i,Ny-1,j),0); + } + // plan (0,*,*) et plan (Nx-1,*,*) + for (int i=0 ; i< Ny ; i++) + for (int j=0 ; j< Nz ; j++) + { + if (Segment(0,i,j) != 0) P.remplace((int) Segment(0,i,j), 0); + if (Segment(Nx-1,i,j) != 0) P.remplace((int) Segment(Nx-1,i,j), 0); + } + } // END CLEAN bord + + P.trie (n); + int FL = FirstLabel-1; + for (int i=0 ; i < Nx; i++) + for (int j=0 ; j < Ny; j++) + for (int k=0 ; k < Nz; k++) + { + int L = Segment(i,j,k); + if (L > n) + { + cout << "Error: Segmentation problem ..... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + Segment(i,j,k) = P.data[L]; + if (Segment(i,j,k) != 0) Segment(i,j,k) += FL; + } + nb_seg = P.max(); +} + + + +////---------------------------------------------------------- +////---------------------------------------------------------- +////---------------------------------------------------------- +////---------------------------------------------------------- +void cube_segment (intarray & cube, intarray & Segment, + int &nb_seg, float S, Bool CleanBord, int FirstLabel) +{ + int Nx=cube.nx(); + int Ny=cube.ny(); + int Nz=cube.nz(); + + + int TabNeighbour[10]; + float old_colour; + Pile P; + int n = 1; + P.add_data (0); + + if ((Segment.nx() != Nx) || (Segment.ny() != Ny)|| (Segment.nz() != Nz)) + Segment.resize(Nx,Ny,Nz); + Segment.init(); + + // First pixel cube(0, 0, 0) + if (ABS(cube(0,0,0)) > S){ + Segment(0,0,0) = 1; + P.add_data (n++); + } + + // X axe cube(1:Nx , 0, 0) + for (int i=1 ; i< Nx ; i++) + if (abs(cube(i,0,0)) > S){ + if (Segment(i-1,0,0) != 0) + Segment(i,0,0) = Segment(i-1,0,0); + else{ + Segment(i,0,0) = n; + P.add_data (n++); + } + } + + // Y axe cube(0, 1:Ny, 0) + for (int i=1 ; i< Ny ; i++) + if (abs(cube(0,i,0)) > S){ + if (ABS(Segment(0,i-1,0)) != 0 ) + Segment(0,i,0) = Segment(0,i-1,0); + else{ + Segment(0,i,0) = n; + P.add_data (n++); + } + } + + // Z axe cube(0, 0, 1:Nz) + for (int i=1 ; i< Nz ; i++) + if (abs(cube(0,0,i)) > S){ + if (ABS(Segment(0,0,i-1)) != 0 ) + Segment(0,0,i) = Segment(0,0,i-1); + else{ + Segment(0,0,i) = n; + P.add_data (n++); + } + } + + // image segmentation + for (int i=1; i S) + { + Bool NewLabel= True; + int UseLabel = 0; + + TabNeighbour[0] = Segment(i-1,j-1,k-1); + TabNeighbour[1] = Segment(i,j-1,k-1); + TabNeighbour[2] = (i+1 n) + { + cout << "Error: Segmentation problem ... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + TabNeighbour[l] = P.data[L]; + if (UseLabel == 0) UseLabel = TabNeighbour[l]; + else if (UseLabel > TabNeighbour[l]) UseLabel = TabNeighbour[l]; + } + } + + // New label + if (NewLabel == True) + { + Segment(i,j,k) = n; + P.add_data (n++); + } + else + { + Segment(i,j,k) = UseLabel; + // neighbourhood label must have the same label + // and therefore may need to be changed + for (int l=0; l < 10; l++) + { + if ((TabNeighbour[l] != 0) && (TabNeighbour[l] != UseLabel)) + { + int L = TabNeighbour[l]; + if (L > n) + { + cout << "Error: Segmentation problem .... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + P.remplace (P.data[L], UseLabel); + } + } + } + } + } + + if (CleanBord == True) + { + // plan (*,*,0) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Ny ; j++) + if (Segment(i,j,0) != 0) + { + old_colour = Segment(i,j,0); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + + } + // plan (*,*,Nz-1) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Ny ; j++) + if (Segment(i,j,Nz-1) != 0) + { + old_colour = Segment(i,j,Nz-1); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + + } + + // plan (*,0,*) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Nz ; j++) + if (Segment(i,0,j) != 0) + { + old_colour = Segment(i,0,j); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + + } + + // plan (*,Ny-1,*) + for (int i=0 ; i< Nx ; i++) + for (int j=0 ; j< Nz ; j++) + if (Segment(i,Ny-1,j) != 0) + { + old_colour = Segment(i,Ny-1,j); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + + } + + // plan (0,*,*) + for (int i=0 ; i< Ny ; i++) + for (int j=0 ; j< Nz ; j++) + if (Segment(0,i,j) != 0) + { + old_colour = Segment(0,i,j); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + } + + // plan (Nx-1,*,*) + for (int i=0 ; i< Ny ; i++) + for (int j=0 ; j< Nz ; j++) + if (Segment(Nx-1,i,j) != 0) + { + old_colour = Segment(Nx-1,i,j); + P.remplace((int)old_colour,0); + for(int a=0; a< Nx; a++) + for(int b=0; b< Ny; b++) + for(int c=0; c< Nz; c++) + if (Segment(a,b,c) == old_colour) Segment(a,b,c) = 0; + } + + } // END CLEAN bord + + P.trie (n); + int FL = FirstLabel-1; + for (int i=0 ; i < Nx; i++) + for (int j=0 ; j < Ny; j++) + for (int k=0 ; k < Nz; k++) + { + int L = Segment(i,j,k); + if (L > n) + { + cout << "Error: Segmentation problem ..... " << endl; + cout << " n = " << n << " TabNeighbour[l] = " << L << endl; + exit(-1); + } + Segment(i,j,k) = P.data[Segment(i,j,k)]; + if (Segment(i,j,k) != 0) Segment(i,j,k) += FL; + } + nb_seg = P.max(); + } diff --git a/src/libsparse3d/CUBE_Segment.h b/src/libsparse3d/CUBE_Segment.h new file mode 100755 index 0000000..ca36521 --- /dev/null +++ b/src/libsparse3d/CUBE_Segment.h @@ -0,0 +1,30 @@ +#ifndef _SEG3D_H_ +#define _SEG3D_H_ + +#include "IM_Obj.h" +#include "IM3D_IO.h" +#include "IM_IO.h" +#define PILE_SIZE_BLOC 100000 + + +void cube_segment (fltarray & cube, intarray & Segment, + int &nb_seg, float S, Bool CleanBord, int FirstLabel); +void cube_segment (intarray & cube, intarray & Segment, + int &nb_seg, float S, Bool CleanBord, int FirstLabel); + +class Pile { + public: + int * data; + int num; + int bloc; + Pile(void); + Pile (int n); + ~Pile (void); + int max (void); + void alloue (void); + void add_data (int d); + void trie (int max); + void remplace (int d, int f); +}; + +#endif diff --git a/src/libsparse3d/MR2D1D.cc b/src/libsparse3d/MR2D1D.cc new file mode 100644 index 0000000..431a56b --- /dev/null +++ b/src/libsparse3d/MR2D1D.cc @@ -0,0 +1,566 @@ +/****************************************************************************** +** Copyright (C) 2007 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 14/06/2007 +** +** File: MR2D1D.cc +** +******************************************************************************* +** +** DESCRIPTION multiresolution transform of cube, using 2D wavelet transform +** ----------- followed (or not) by a 1D wavelet transform +** +** +******************************************************************************/ + + +#include "MR2D1D.h" + + +/****************************************************************************/ + +static void mr_io_name (char *File_Name_In, char *File_Name_Out) +{ + int L; + + strcpy (File_Name_Out, File_Name_In); + + L = strlen (File_Name_In); + if ((L < 3) || (File_Name_In[L-1] != 'r') + || (File_Name_In[L-2] != 'm') + || (File_Name_In[L-3] != '.')) + { + strcat (File_Name_Out, ".mr"); + } +} + +/****************************************************************************/ + +/*--------------------------------------------------------------------------*/ +static void PrintError( int status) +{ + /*****************************************************/ + /* Print out cfitsio error messages and exit program */ + /*****************************************************/ + + char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG]; + + if (status) + fprintf(stderr, "\n*** Error occurred during program execution ***\n"); + + ffgerr(status, status_str); /* get the error status description */ + fprintf(stderr, "\nstatus = %d: %s\n", status, status_str); + + if ( ffgmsg(errmsg) ) /* get first message; null if stack is empty */ + { + fprintf(stderr, "\nError message stack:\n"); + fprintf(stderr, " %s\n", errmsg); + + while ( ffgmsg(errmsg) ) /* get remaining messages */ + fprintf(stderr, " %s\n", errmsg); + } + + exit( status ); /* terminate the program, returning error status */ +} + +/*--------------------------------------------------------------------------*/ + + +int MR2D1D::mr_io_fill_header( fitsfile *fptr) +{ + int status = 0; // this means OK for cfitsio !!! + /*****************************************************/ + /* write optional keyword to the header */ + /*****************************************************/ + +// cout << " DDDD " << (long) NbrBand2D << " " << (long)NbrBand1D << endl; + +if ( ffpkyj(fptr, "Nx", (long)Nx,"Nx of the input cube",&status)) + PrintError( status ); +if ( ffpkyj(fptr,"Ny",(long)Ny,"Ny of the input cube",&status)) + PrintError( status ); +if ( ffpkyj(fptr,"Nz",(long)Nz,"Nz of the input cube",&status)) + PrintError( status ); + if ( ffpkyj(fptr, "NSCALE2D", (long) NbrScale2D, "Number of bands 2D", &status)) + PrintError( status ); + if ( ffpkyj(fptr, "NSCALE1D", (long)NbrScale1D, "Number of bands 1D", &status)) + PrintError( status ); + if ( ffpkyj(fptr, "Type_Tra", (long) WT2D.Type_Transform, + (char*)StringTransform(WT2D.Type_Transform), &status)) + PrintError( status ); + return(status); +} + + +/****************************************************************************/ + +void MR2D1D::read(char *Name) +{ + char filename[256]; + fitsfile *fptr; /* pointer to the FITS file */ + int status=0, hdutype ; + long hdunum; + char comment[FLEN_COMMENT]; + int naxis; + long naxes[3]; + long mon_long; + int anynul = 0; + long nulval = 0; + long inc[3]; + void PrintError( int status); + long nelements = 0 ; // naxes[0] * naxes[1] in the image + //long fpixels[3]; + //long int lpixels[3]; + + // for multiresol + float *Ptr; + //int my_logical; // sais pas... + + mr_io_name (Name, filename); + + inc[0]=1; inc[1]=1; inc[2]=1; + +#if DEBUG_IO + cout << "Read in " << filename << endl; +#endif + + /* open the file */ + status = 0; /* initialize status before calling fitsio routines */ + if ( ffopen(&fptr, filename, (int) READONLY, &status) ) + PrintError( status ); + + hdunum = 1; /*read table */ + if ( ffmahd(fptr, hdunum, &hdutype, &status) ) /* move to the HDU */ + PrintError( status ); + + int simple, bitpix, extend; + long pcount, gcount; + if ( ffghpr(fptr, 3, &simple, &bitpix, &naxis, naxes, &pcount, + &gcount, &extend, &status)) /* move to the HDU */ + PrintError( status ); + + nelements = naxes[0]; + + if (ffgkyj(fptr,"Nx", &mon_long, comment, &status)) PrintError( status ); + int Nxi = (int) mon_long; + if (ffgkyj(fptr,"Ny", &mon_long, comment, &status)) PrintError( status ); + int Nyi = (int) mon_long; + if (ffgkyj(fptr,"Nz", &mon_long, comment, &status)) PrintError( status ); + int Nzi = (int) mon_long; + + if (ffgkyj(fptr,"NSCALE2D", &mon_long, comment, &status)) PrintError( status ); + int iNbrScale2D = (int) mon_long; + if (ffgkyj(fptr,"NSCALE1D", &mon_long, comment, &status)) PrintError( status ); + int iNbrScale1D = (int) mon_long; + + type_transform TT; + if (ffgkyj(fptr,"Type_Tra", &mon_long, comment, &status)) PrintError( status ); + else TT = (type_transform) mon_long; + + alloc(Nxi,Nyi,Nzi,TT, iNbrScale2D, iNbrScale1D); + + fltarray Tab(nelements); + //cout << " Nl = " << (TabCF_Band[0][1]).nl() << " Nc = " << TabCF_Band[0][1].nc() << endl; + + Ptr = Tab.buffer(); + if ( ffgpve(fptr, 1, 1, nelements, nulval, Ptr, &anynul, &status)) + PrintError( status ); + + int ind=2; + for (int s=0; s = size_band_nx(s2)) || + (j < 0) || (j >= size_band_ny(s2)) || + (k < 0) || (k >= size_band_nz(s2)) || + (s2 < 0) || (s2 >= nbr_band_2d()) ) + { + printf("Error: (s2,i,j,k)=(%d,%d,%d,%d), size band = (%d,%d,%d) , scale=(%d,%d) \n", s2,i,j,k,size_band_nx(s2), size_band_ny(s2), size_band_nz(s2), nbr_band_2d(), nbr_band_1d()); + exit(-1); + } +#endif + + return TabBand[s2] (i,j,k); +} + +/****************************************************************************/ + +fltarray MR2D1D::get_band(int s2, int s1) +{ + int Nxb = size_band_nx(s2,s1); + int Nyb = size_band_ny(s2,s1); + int Nzb = size_band_nz(s2,s1); + fltarray *Cube_Return = NULL; + + Cube_Return = new fltarray(Nxb, Nyb, Nzb); + for (int i=0; i < Nxb; i++) + for (int j=0; j < Nyb; j++) + for (int k=0; k < Nzb; k++) (*Cube_Return)(i,j,k) = (*this)(s2,s1,i,j,k); + + return (*Cube_Return); +} + +/****************************************************************************/ + +void MR2D1D::put_band(fltarray Band, int s2, int s1) +{ + int Nxb = size_band_nx(s2,s1); + int Nyb = size_band_ny(s2,s1); + int Nzb = size_band_nz(s2,s1); + + for (int i=0; i < Nxb; i++) + for (int j=0; i < Nyb; j++) + for (int k=0; k < Nzb; k++) (*this)(s2,s1,i,j,k) = Band(i,j,k); +} + +/****************************************************************************/ + +void MR2D1D::alloc (int iNx, int iNy, int iNz, type_transform Trans2D, int Ns2D, int Ns1D, Bool NoAlloc) +{ + Nx = iNx; + Ny = iNy; + Nz = iNz; + NbrScale2D = Ns2D; + NbrScale1D = Ns1D; + if (NbrScale1D < 2) + { + NbrScale1D = 1; + Apply1DTrans = False; + } + else Apply1DTrans = True; + ; + Norm = NORM_L2; + SB_Filter = F_MALLAT_7_9; + Bord = I_CONT; + U_Filter = DEF_UNDER_FILTER; + FilterAnaSynt *PtrFAS = NULL; + if ((Trans2D == TO_MALLAT) || (Trans2D == TO_UNDECIMATED_MALLAT)) + { + FAS.Verbose = Verbose; + FAS.alloc(SB_Filter); + PtrFAS = &FAS; + } + int NbrUndec = -1; /*number of undecimated scale */ + type_lift LiftingTrans = DEF_LIFT; + if (Trans2D == TO_LIFTING) WT2D.LiftingTrans = LiftingTrans; + WT2D.Border = Bord; + WT2D.Verbose = Verbose; + WT2D.alloc (Ny, Nx, Ns2D, Trans2D, PtrFAS, Norm, NbrUndec, U_Filter); + NbrBand2D = WT2D.nbr_band(); + + + Bool Rebin=False; + WT1D.U_Filter = U_Filter; + type_trans_1d Trans1D = TO1_MALLAT; + if (Apply1DTrans == True) + { + WT1D.alloc (Nz, Trans1D, Ns1D, PtrFAS, Norm, Rebin); + NbrBand1D = WT1D.nbr_band(); + } + else NbrBand1D = 1; + + if (NoAlloc == False) TabBand = new fltarray [NbrBand2D]; + else TabBand=NULL; + TabSizeBandNx.resize(NbrBand2D, NbrBand1D); + TabSizeBandNy.resize(NbrBand2D, NbrBand1D); + TabSizeBandNz.resize(NbrBand2D, NbrBand1D); + TabFirstPosBandNz.resize(NbrBand1D); + TabFirstPosBandNz(0) =0; + NbrCoef2D = 0; + for (int b=0; b < NbrBand2D; b++) + { + int Npz = (Apply1DTrans == True) ? WT1D.size_ima_np (): Nz; + if (NoAlloc == False) TabBand[b].alloc(WT2D.size_band_nc(b), WT2D.size_band_nl(b), Npz); + NbrCoef2D += WT2D.size_band_nc(b) * WT2D.size_band_nl(b); + for (int b1=0; b1 < NbrBand1D; b1++) + { + TabSizeBandNx(b,b1) = WT2D.size_band_nc(b); + TabSizeBandNy(b,b1) = WT2D.size_band_nl(b); + if (Apply1DTrans == True) TabSizeBandNz(b,b1) = WT1D.size_scale_np(b1); + else TabSizeBandNz(b,b1) = Nz; + } + } + NbrCoef1D = Nz; + if (Apply1DTrans == True) + for (int b1=1; b1 < NbrBand1D; b1++) TabFirstPosBandNz(b1) = TabFirstPosBandNz(b1-1) + TabSizeBandNz(0,b1-1); +} + +/****************************************************************************/ + +void MR2D1D::transform_to_vectarray(fltarray &Data, fltarray &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(); + 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 (fltarray &Data) +{ + int i,j,b,z; + Nx = Data.nx(); + Ny = Data.ny(); + Nz = Data.nz(); + Ifloat Frame(Ny, Nx); + fltarray Vect(Nz); + + // 2D wt transform per frame + for (z=0; z < Nz; z++) + { + 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 < NbrBand2D; b++) + { + for (i=0; i < WT2D.size_band_nl(b); i++) + for (j=0; j < WT2D.size_band_nc(b); j++) TabBand[b](j,i,z) = WT2D(b,i,j); + } + } + + // 1D wt + if (NbrBand1D >= 2) + { + for (b=0; b < NbrBand2D; 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) = 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++) TabBand[b](j,i,z++) = WT1D(b1,p); + } + } + } +} + +/****************************************************************************/ + + + +void MR2D1D::recons (fltarray &Data) +{ + if ((Data.nx() != Nx) || (Data.ny() != Ny) || (Data.nz() != Nz)) Data.resize(Nx, Ny, Nz); + int i,j,b,z; + Nx = Data.nx(); + Ny = Data.ny(); + Nz = Data.nz(); + Ifloat Frame(Ny, Nx); + fltarray Vect(Nz); + + // 1D wt + if (NbrBand1D >= 2) + { + for (b=0; b < NbrBand2D; 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) = TabBand[b](j,i,z); + z = 0; + for (int b1=0; b1 < NbrBand1D; b1++) + for (int p=0; p < WT1D.size_scale_np (b1); p++) WT1D(b1,p) = TabBand[b](j,i,z++); + Vect.init(); + WT1D.recons(Vect); + for (z=0; z < Nz; z++) TabBand[b](j,i,z) = Vect(z); + } + } + + // 2D wt + for (z=0; z < Nz; z++) + { + for (b=0; b < NbrBand2D; b++) + { + for (i=0; i < WT2D.size_band_nl(b); i++) + for (j=0; j < WT2D.size_band_nc(b); j++) WT2D(b,i,j) = TabBand[b](j,i,z); + } + WT2D.recons(Frame); + for (i=0; i < Ny; i++) + for (j=0; j < Nx; j++) Data(j,i,z) = Frame(i,j); + } +} + + /****************************************************************************/ + + diff --git a/src/libsparse3d/MR2D1D.h b/src/libsparse3d/MR2D1D.h new file mode 100755 index 0000000..153e4ec --- /dev/null +++ b/src/libsparse3d/MR2D1D.h @@ -0,0 +1,139 @@ +/****************************************************************************** +** Copyright (C) 2007 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 14/06/2007 +** +** File: MR2D1D.h +** +******************************************************************************* +** +** DESCRIPTION multiresolution transform of cube, using 2D wavelet transform +** ----------- followed (or not) by a 1D wavelet transform +** +** +******************************************************************************/ + + +#ifndef _MR2D1D_H_ +#define _MR2D1D_H_ + +#include "IM_Obj.h" +#include "IM_IO.h" +#include "SB_Filter.h" +#include "MR1D_Obj.h" +#include "MR_Obj.h" +#include "IM3D_IO.h" +#include "MR3D_Obj.h" +#include "DefFunc.h" + +#define TESTIND 1 + +/****************************************************************************/ + +class MR2D1D { + Bool Apply1DTrans; // If True, a 1D Wavelet transform is applied on the z-axis + int Nx, Ny, Nz; // input cube size + int NbrScale2D; // Number of scales of the 2D multiscale transform. + int NbrScale1D; // Number of scales of the 1D multiscale transform (NbrScale1D=1 if no WT1D is applied) + int NbrBand2D; // Number of bands in the 2D multiscale transform + int NbrBand1D; // Number of bands in the 1D multiscale transform + int NbrCoef2D; // Number of coefficients in the 2D multiscale transform + int NbrCoef1D; // Number of coefficients in the 1D multiscale transform + + public: // ATTENTION JE L'AI MIS ICI ET NON PLUS LOIN !! + + MultiResol WT2D; // 2D Multiresolution object + MR_1D WT1D; // 1D Multiresolution object + fltarray *TabBand; // Array of 3D bands + intarray TabFirstPosBandNz; + // TabBand[b2](Nx, Ny, TabFirstPosBandNz[b1]:TabFirstPosBandNz[b1]:TabSizeBandNz[b2, b1]-1) + // contains the band (b2,b1) b2=band 2D, b1=band 1D + intarray TabSizeBandNx; // TabSizeBandNx[b2, b1] = size along x-axis of the (b2,b1) b2=band 2D, b1=band 1D + intarray TabSizeBandNy; // TabSizeBandNy[b2, b1] = size along y-axis of the (b2,b1) b2=band 2D, b1=band 1D + intarray TabSizeBandNz; // TabSizeBandNz[b2, b1] = size along z-axis of the (b2,b1) b2=band 2D, b1=band 1D + + // WT 2D + sb_type_norm Norm; // Type or normalisation of the 2D transform + type_sb_filter SB_Filter; // Type of filter in case of 2D WT + type_border Bord; // Parameter for border management + type_undec_filter U_Filter; // Type of filter in case of undecimated WT + FilterAnaSynt FAS; // Filter bank object + + int mr_io_fill_header(fitsfile *fptr); + //public: + Bool Verbose; + MR2D1D (){ NbrBand2D=NbrBand1D=0;Verbose=False;} + + void alloc(int iNx, int iNy, int iNz, type_transform Trans2D, int Ns2D, int Ns1D=0, Bool NoAlloc=False); + // Allocate the class for a cube of size (iNx, iNy, iNz) using Ns2D scale in 2D and Ns1D scale in 1D + // If Ns1D < 2 then no wavelet transform is performed along z axis + // If NoAlloc=True, then TabBand is not allocated and ONLY routine transform_to_vectarray can be used + + void transform (fltarray &Data); + // Apply the 2D-1D wavelet transform. The transformation is stored in TabBand + + void recons (fltarray &Data); + // Apply the 2D-1D reconstruction + + fltarray get_band(int s2, int s1=0); + // Extract one band + + void put_band(fltarray Band, int s2, int s1=0); + // Insert one band + + void write(char *Name); + // Write the transformation to a file + + void read(char *Name); + // Read a transformation from a file + + 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 + // get_band, put_band, read, write routines cannot be used + + int nbr_band_2d () const { return NbrBand2D;} + // Return the number of bands of the 2D multiresolution transform + + int nbr_band_1d () const { return NbrBand1D;} + // Return the number of bands of the 1D wavelet transform + + int size_band_nx(int s2d, int s1d=0) const { return TabSizeBandNx(s2d, s1d);} + // Return the size of the band (s2d,s1d) along x-axis + + int size_band_ny(int s2d, int s1d=0) const { return TabSizeBandNy(s2d, s1d);} + // Return the size of the band (s2d,s1d) along y-axis + + int size_band_nz(int s2d, int s1d=0) const { return TabSizeBandNz(s2d, s1d);} + // Return the size of the band (s2d,s1d) along z-axis + + int nbr_coef_2d() const { return NbrCoef2D;} + // Return the number of coefficients of the 2D decomposition + + int nbr_coef_1d() const { return NbrCoef1D;} + // Return the number of coefficients of the 1D decomposition + + int nbr_pix_nx() const { return Nx;} + // Return the number of coefficients of the 1D decomposition + + int nbr_pix_ny() const { return Ny;} + // Return the number of coefficients of the 1D decomposition + + float & operator() (int s2, int s1, int i, int j, int k) const; + // Return one coefficient + + float & operator() (int s2, int i, int j, int k) const; + // Return one coefficient, for the case where no 1D WT is applied + + ~MR2D1D() { if (TabBand != NULL) delete [] TabBand; NbrBand1D=NbrBand2D=NbrScale1D=NbrScale2D=0;} +}; + +/****************************************************************************/ + #endif diff --git a/src/libsparse3d/MR3D_Obj.cc b/src/libsparse3d/MR3D_Obj.cc new file mode 100755 index 0000000..76953db --- /dev/null +++ b/src/libsparse3d/MR3D_Obj.cc @@ -0,0 +1,1338 @@ +/****************************************************************************** +** Copyright (C) 1999 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 04/09/99 +** +** File: mr3d_trans.cc +** +******************************************************************************* +** +** DESCRIPTION multiresolution transform of cube +** ----------- +** +** Usage: mr3d_trans options cube output +** +** Append : +** 07/2008 : Arnaud Woiselle +** undecimated transforms : PAVE_3D_WT +** +******************************************************************************/ + +#include "MR3D_Obj.h" + +/************************************************************************/ +// 3D sub-band decomposision: an image is transformed into eight sub-cubes + +void SubBand3D::transform3d (fltarray &Data) +{ +// cerr<<"SubBand3D::transform..."<transform(Nz, PtrHigh, PtrLow, PtrDet); + for (k=0; k< Nz2; k++) Data(i,j,k) = PtrLow[k]; + for (k=0; k< Nz/2; k++) Data(i,j,k+Nz2) = PtrDet[k]; + } + + delete [] PtrHigh; + delete [] PtrLow; + delete [] PtrDet; +// cerr<<"End SubBand3D::transform"< s z + // Horiz => x xz + // Diag => xy xyz + // Vert => y yz + // bands : 0 1 2 3 4 5 6 7 + // xyz : HHH HHL HLH HLL LHH LHL LLH LLL (high/low) + // DH DL HH HL VH VL SH SL + + // each z-vector is convolved by h and g + for (i = 0; i < Nx; i++) + for (j = 0; j < Ny; j++) + { + for (k=0; ktransform(Nz, PtrHigh, PtrLow, PtrDet, Step); + for (k=0; ktransform(Nz, PtrHigh, PtrLow, PtrDet, Step); + for (k=0; ktransform(Nz, PtrHigh, PtrLow, PtrDet, Step); + for (k=0; ktransform(Nz, PtrHigh, PtrLow, PtrDet, Step); + for (k=0; krecons(Nz,PtrLow, PtrDet, PtrHigh); + for (k=0; krecons(Nz, PtrLow, PtrDet, PtrHigh, Step); + for (k=0; krecons(Nz, PtrLow, PtrDet, PtrHigh, Step); + for (k=0; krecons(Nz, PtrLow, PtrDet, PtrHigh, Step); + for (k=0; krecons(Nz, PtrLow, PtrDet, PtrHigh, Step); + for (k=0; k= 0; s--) + { + Nxs = size_resol(s, Nx); + Nys = size_resol(s, Ny); + Nzs = size_resol(s, Nz); + + Nx_2 = (Nxs+1)/2; + Ny_2 = (Nys+1)/2; + Nz_2 = (Nzs+1)/2; + fltarray Data_Aux(Nxs,Nys,Nzs); + + // copy the transform part in Imag + for (i = 0; i < Nxs; i++) + for (j = 0; j < Nys; j++) + for (k = 0; k < Nzs; k++) Data_Aux(i,j,k) = Data(i,j,k); + + // inverse transform Ima_Aux + recons3d(Data_Aux); + + for (i = 0; i < Nxs; i++) + for (j = 0; j < Nys; j++) + for (k = 0; k < Nzs; k++) Data(i,j,k) = Data_Aux(i,j,k); + } + +} + +/************************************************************************/ +// PAVE_3D_WT OBJ +/****************************************************************************/ + +int PAVE_3D_WT::alloc (fltarray * & TabBand, int Nx, int Ny, int Nz, int NbrScale) +{ + int s,NbrBand_per_Resol = 7; + int NbrBand = NbrBand_per_Resol*(NbrScale-1)+1; + TabBand = new fltarray [NbrBand]; + for (s = 0; s < NbrBand; s++) + TabBand[s].alloc(Nx,Ny,Nz); + return NbrBand; +} + +/****************************************************************************/ + +void PAVE_3D_WT::free(fltarray *TabBand, int NbrScale) +{ + if (NbrScale != 0) delete [] TabBand; +} + +/****************************************************************************/ + +void PAVE_3D_WT::transform (fltarray &Cube, fltarray *TabTrans, int NbrScale) +{ +// cerr<<"PAVE_3D_WT::transform..."<= 0; s--) + { + int Step = POW2(s); + if (s == 0) + recons3d(TabTrans, Cube, Step); + else + recons3d(&(TabTrans[7*s]), TabTrans[7*s], Step); + } +// cerr<<"End PAVE_3D_WT::recons"<= nbr_band()) + || (i < 0) || (i >= size_band_nx(b)) + || (j < 0) || (j >= size_band_ny(b)) + || (k < 0) || (k >= size_band_nz(b))) + { + cout << "Error: band coefficient index ... " << endl; + cout << " Band number = " << b << " Nb = " << nbr_band() << endl; + cout << " X pos = " << i << " Nx = " << size_band_nx(b) << endl; + cout << " Y pos = " << j << " Ny = " << size_band_ny(b) << endl; + cout << " Z pos = " << k << " Nz = " << size_band_nz(b) << endl; + exit(-1); + } + if (Set_Transform == TRANS3_PAVE) + { + return TabBand[b](i,j,k); + } + else + { + int Indi = TabPosX[b] + i; + int Indj = TabPosY[b] + j; + int Indk = TabPosZ[b] + k; + + if ((Indi < 0) || (Indi >= Nx) || + (Indj < 0) || (Indj >= Ny) || (Indk < 0) || (Indk >= Nz)) + { + cout << "Error: band coefficient index ... " << endl; + cout << " Band number = " << b << " Nb = " << nbr_band() << endl; + cout << " Indi = " << Indi << " Nx = " << Nx << endl; + cout << " Indj = " << Indj << " Ny = " << Ny << endl; + cout << " Indk = " << Indk << " Nz = " << Nz << endl; + exit(-1); + } + return Data(Indi,Indj,Indk); + } +} + +/****************************************************************************/ + + +void set_size_band(int N, Bool FilterLow, int & P, int &Ns) +{ + int L = (N+1)/2; + + if (FilterLow == True) + { + Ns = L; + P = 0; + } + else + { + Ns = N/2; + P = L; + } +} + +/****************************************************************************/ + +void MR_3D::init() +{ + Nx=0;Ny=0;Nz=0; + Nbr_Plan=0; + Type_Transform = T3_UNDEFINED; + Set_Transform = S3_UNDEFINED; + Verbose=False; + FilterBank=NULL; + FilterBankAlloc=False; + LiftingTrans = DEF_LIFT; + SBFilter = DEF_SB_FILTER; + TypeNorm = DEF_SB_NORM; + TabBand=NULL; +} + +/****************************************************************************/ + +void MR_3D::alloc (int Npx, int Npy, int Npz, type_trans_3d T, int Nbr_Scale, + FilterAnaSynt *FAS, sb_type_norm Norm) +{ + extern type_3d_format IO_3D_Format; + int b,s; + + Nx = Npx; + Ny = Npy; + Nz = Npz; + Nbr_Plan = Nbr_Scale; + Type_Transform = T; + Set_Transform = SetTransform (T); + Border = DEFAULT_BORDER_3D; + + DataFormat = IO_3D_Format; + if (Set_Transform == TRANS3_PAVE) + { + Nbr_Band = Nbr_Plan; + TabSizeNx = new int [Nbr_Band]; + TabSizeNy = new int [Nbr_Band]; + TabSizeNz = new int [Nbr_Band]; + TabPosX = new int [Nbr_Band]; + TabPosY = new int [Nbr_Band]; + TabPosZ = new int [Nbr_Band]; + AT3D_WT.alloc(TabBand, Nx, Ny, Nz, Nbr_Plan); + for (b =0 ; b < Nbr_Band; b++) + { + TabSizeNx[b] = Nx; + TabSizeNy[b] = Ny; + TabSizeNz[b] = Nz; + TabPosX[b] = 0; + TabPosY[b] = 0; + TabPosZ[b] = 0; + } + } + else + { + FilterBank = FAS; + if (FAS != NULL) SBFilter = FAS->type_filter(); + else if (T == TO3_MALLAT) + { + if (FilterBank == NULL) + { + SBFilter = DEF_SB_FILTER; + TypeNorm = DEF_SB_NORM; + FilterBank = new FilterAnaSynt; + FilterBank->alloc(SBFilter); + FilterBankAlloc=True; + } + } + TypeNorm = Norm; + + // LiftingTrans = DEF_LIFT; + // SB_Filter = F_MALLAT_7_9; + // Norm = NORM_L1; + Nbr_Band = 7 * (Nbr_Plan-1) + 1; + TabSizeNx = new int [Nbr_Band]; + TabSizeNy = new int [Nbr_Band]; + TabSizeNz = new int [Nbr_Band]; + TabPosX = new int [Nbr_Band]; + TabPosY = new int [Nbr_Band]; + TabPosZ = new int [Nbr_Band]; + + Data.alloc (Nx,Ny,Nz); + TabBand = &Data; + if (Set_Transform == TRANS3_MALLAT) + { + int Lx = Nx; + int Ly = Ny; + int Lz = Nz; + + for (s=0; s < Nbr_Plan-1; s++) + { + int Ns,P; + + // Horizontal low + b = 7*s; + set_size_band(Lx, False, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, True, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, True, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // Vertical low + b = 7*s+1; + set_size_band(Lx, True, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, False, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, True, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // Diagonal low + b = 7*s+2; + set_size_band(Lx, False, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, False, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, True, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // Horizontal High + b = 7*s+3; + set_size_band(Lx, False, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, True, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, False, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // Vertical high + b = 7*s+4; + set_size_band(Lx, True, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, False, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, False, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // Diagonal high + b = 7*s+5; + set_size_band(Lx, False, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, False, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, False, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + // z-high, and x,y low + b = 7*s+6; + set_size_band(Lx, True, P, Ns); + TabPosX[b] = P; + TabSizeNx[b] = Ns; + set_size_band(Ly, True, P, Ns); + TabPosY[b] = P; + TabSizeNy[b] = Ns; + set_size_band(Lz, False, P, Ns); + TabPosZ[b] = P; + TabSizeNz[b] = Ns; + + Lx = (Lx+1) / 2; + Ly = (Ly+1) / 2; + Lz = (Lz+1) / 2; + } + b = Nbr_Band-1; + TabPosX[b] = 0; + TabPosY[b] = 0; + TabPosZ[b] = 0; + TabSizeNx[b] = Lx; + TabSizeNy[b] = Ly; + TabSizeNz[b] = Lz; + + }} +} + +/****************************************************************************/ + +void MR_3D::info_pos_band() +{ + int s; + + switch (Set_Transform) + { + case TRANS3_MALLAT: + cout << "Number of bands = " << Nbr_Band << endl; + for (s=0; s < Nbr_Band; s++) + { + cout << "Band " << s+1; + cout << " Pos = [" << TabPosX[s] << "," << TabPosY[s] << "," << TabPosZ[s] << "]"; + cout << " Size = [" << TabSizeNx[s] << "," << TabSizeNy[s] << "," << TabSizeNz[s] << "]"; + cout << " Band = [" << TabPosX[s] << ":" << TabPosX[s]+TabSizeNx[s]-1 << "," << + TabPosY[s] << ":" << TabPosY[s]+TabSizeNy[s]-1 << "," << + TabPosZ[s] << ":" << TabPosZ[s]+TabSizeNz[s]-1 << "]" << endl; + } + break; + case TRANS3_PAVE: + cout << "Number of bands = " << Nbr_Band << endl; + for (s=0; s < Nbr_Band; s++) + { + cout << "Band " << s+1; + cout << " Size = [" << TabSizeNx[s] << "," << TabSizeNy[s] << "," << TabSizeNz[s] << "]" << endl; + } + break; + default: cerr << "Error: bad transform ... " << endl; + exit(-1); + } +} + +/****************************************************************************/ + +void MR_3D::get_band(int b, fltarray &Band) +{ + int i,j,k; + int Nxb = size_band_nx(b); + int Nyb = size_band_ny(b); + int Nzb = size_band_nz(b); + + if (Band.n_elem() == 0) Band.alloc(Nxb,Nyb,Nzb); + else if ((Band.naxis() != 3) || (Band.nx() != Nxb) + || (Band.ny() != Nyb) || (Band.nz() != Nzb)) + { + Band.free(); + Band.alloc(Nxb,Nyb,Nzb); + } + for (i=0; i < Nxb; i++) + for (j=0; j < Nyb; j++) + for (k=0; k < Nzb; k++) Band(i,j,k) = (*this)(b,i,j,k); +} + +/****************************************************************************/ + +void MR_3D::insert_band(int b, fltarray &Band) +{ + int i,j,k; + int Nxb = size_band_nx(b); + int Nyb = size_band_ny(b); + int Nzb = size_band_nz(b); + + if ((Band.n_elem() == 0) || (Band.naxis() != 3) + || (Band.nx() != Nxb) + || (Band.ny() != Nyb) || (Band.nz() != Nzb)) + { + cerr << "Error: band to insert has not the correct dimensions ... " << endl; + exit(-1); + } + for (i=0; i < Nxb; i++) + for (j=0; j < Nyb; j++) + for (k=0; k < Nzb; k++) (*this)(b,i,j,k) = Band(i,j,k); +} + +/****************************************************************************/ + +void MR_3D::info_band(int b) +{ + int i,j,k; + double Min = (*this)(b,0,0,0); + double Max = Min; + double Flux = 0; + double Mean, Sigma=0; + int Nxb = size_band_nx(b); + int Nyb = size_band_ny(b); + int Nzb = size_band_nz(b); + long Np = Nxb*Nyb*Nzb; + + for (i=0; i < Nxb; i++) + for (j=0; j < Nyb; j++) + for (k=0; k < Nzb; k++) + { + double Coef = (*this)(b,i,j,k); + if (Coef < Min) Min = (*this)(b,i,j,k); + if (Coef > Max) Max = (*this)(b,i,j,k); + Flux += Coef; + } + Mean = Flux / (double) Np; + for (i=0; i < Nxb; i++) + for (j=0; j < Nyb; j++) + for (k=0; k < Nzb; k++) + { + double Coef = (*this)(b,i,j,k) - Mean; + Sigma += Coef*Coef; + } + Sigma /= (double) Np; + + cout << "Band " << b+1 << endl; + cout << " Pos = [" << TabPosX[b] << "," << TabPosY[b] << "," << TabPosZ[b] << "]"; + cout << " Size = [" << Nxb << "," << Nyb << "," << Nzb << "]"; + cout << " Band = [" << TabPosX[b] << ":" << TabPosX[b]+Nxb-1 << "," << + TabPosY[b] << ":" << TabPosY[b]+Nyb-1 << "," << + TabPosZ[b] << ":" << TabPosZ[b]+Nzb-1 << "]" << endl; + cout << " Min = " << Min << " Max = " << Max << endl; + cout << " Mean = " << Mean << " Flux = " << Flux << endl; + cout << " Sigma = " << Sigma << endl << endl; +} + + +/****************************************************************************/ + +MR_3D::~MR_3D() +{ + free (); +} + +/****************************************************************************/ + +void MR_3D::free () +{ + + Border = DEFAULT_BORDER_3D; + switch (Set_Transform) + { + case TRANS3_MALLAT: + Data.free(); + break; + case TRANS3_PAVE: + AT3D_WT.free(TabBand, Nbr_Plan); + break; + default: cerr << "Error: bad transform ... " << endl; + exit(-1); + } + Nbr_Plan = 0; + Nbr_Band = 0; + Nx = Ny = Nz = 0; + Set_Transform= S3_UNDEFINED; + Type_Transform=T3_UNDEFINED; + if (TabPosX != NULL) delete [] TabPosX; + if (TabSizeNx != NULL) delete [] TabSizeNx; + if (TabPosY != NULL) delete [] TabPosY; + if (TabSizeNy != NULL) delete [] TabSizeNy; + if (TabPosZ != NULL) delete [] TabPosZ; + if (TabSizeNz != NULL) delete [] TabSizeNz; + if ((FilterBankAlloc == True) && (FilterBank != NULL)) + { + delete FilterBank; + FilterBank = NULL; + } + init(); +} + +/****************************************************************************/ + +void MR_3D::transform (fltarray &Cube, type_border Bord) +{ + if (Nbr_Plan < 2) + { + cerr << "Error: Object not correctly allocated: Nbr_Plan must be > 2 "; + cerr << endl; + exit (-1); + } + switch (Type_Transform) + { + case TO3_ATROUS: + AT3D_WT.Bord = Bord; + AT3D_WT.transform(Cube, TabBand, Nbr_Plan); + break; + case TO3_MALLAT: + { + Data = Cube; + SubBandFilter WT1D(*FilterBank, TypeNorm); + Ortho_3D_WT WT3D(WT1D); + WT3D.transform(Data, Nbr_Plan); + } + break; + case TO3_LIFTING: + { + Data = Cube; + Lifting Clift1D(LiftingTrans); + Ortho_3D_WT WT3D(Clift1D); + WT3D.transform(Data, Nbr_Plan); + } + break; + default: + fprintf (stderr, "Error (proc. MR_3D_transform): Unknown transform\n"); + exit (-1); + break; + } +} + +/****************************************************************************/ + +void MR_3D::transform (fltarray &Cube) +{ + this->transform(Cube, Border); +} + +/****************************************************************************/ + +void MR_3D::recons (fltarray &Cube) +{ + this->recons(Cube, Border); +} + +/****************************************************************************/ + +void MR_3D::recons (fltarray &Cube, type_border Bord) +{ + Cube = Data; + switch (Type_Transform) + { + case TO3_ATROUS: + AT3D_WT.Bord = Bord; + AT3D_WT.recons(TabBand, Cube, Nbr_Plan); + break; + case TO3_MALLAT: + { + SubBandFilter WT1D(*FilterBank, TypeNorm); + Ortho_3D_WT WT3D(WT1D); + WT3D.recons(Cube, Nbr_Plan); + } + break; + case TO3_LIFTING: + { + Lifting Clift1D(LiftingTrans); + Ortho_3D_WT WT3D(Clift1D); + WT3D.recons(Cube, Nbr_Plan); + } + break; + default: + fprintf (stderr, "Error (proc. MR_3D_transform): Unknown transform\n"); + exit (-1); + break; + } +} + +/****************************************************************************/ + +static void mr3d_io_name (char *File_Name_In, char *File_Name_Out) +{ + int L; + + strcpy (File_Name_Out, File_Name_In); + + L = strlen (File_Name_In); + if ((L < 3) || (File_Name_In[L-1] != 'r') + || (File_Name_In[L-2] != 'm') + || (File_Name_In[L-3] != '.')) + { + strcat (File_Name_Out, ".mr"); + } +} + +/****************************************************************************/ +static void PrintError( int status) +{ + /*****************************************************/ + /* Print out cfitsio error messages and exit program */ + /*****************************************************/ + + char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG]; + + if (status) + fprintf(stderr, "\n*** Error occurred during program execution ***\n"); + + ffgerr(status, status_str); /* get the error status description */ + fprintf(stderr, "\nstatus = %d: %s\n", status, status_str); + + if ( ffgmsg(errmsg) ) /* get first message; null if stack is empty */ + { + fprintf(stderr, "\nError message stack:\n"); + fprintf(stderr, " %s\n", errmsg); + + while ( ffgmsg(errmsg) ) /* get remaining messages */ + fprintf(stderr, " %s\n", errmsg); + } + + exit( status ); /* terminate the program, returning error status */ +} + +/******************************************************************/ + +int MR_3D::mr_io_fill_header(fitsfile *fptr) +{ + int status = 0; // this means OK for cfitsio !!! + /*****************************************************/ + /* write optional keyword to the header */ + /*****************************************************/ + + +// if ( ffpkyj(fptr, "Nx", (long)Nx,"x-axis size",&status)) PrintError( status ); +// if ( ffpkyj(fptr,"Ny",(long)Ny,"y-axis size",&status)) PrintError( status ); +// if ( ffpkyj(fptr,"Nz",(long)Nz,"z-axis size",&status)) PrintError( status ); + if ( ffpkyj(fptr, (char*)"Nbr_Plan", (long)Nbr_Plan, (char*)"Number of scales", &status)) + PrintError( status ); + if (ffpkyj(fptr, (char*)"Set_Transform", (long)Set_Transform, + StringSet3D(Set_Transform), &status)) + PrintError( status ); + if ( ffpkyj(fptr, (char*)"Type_Transform", (long)Type_Transform, + StringTransf3D(Type_Transform), &status)) + PrintError( status ); + + if (Type_Transform == TO3_MALLAT) + { + if ( ffpkyj(fptr, (char*)"SBFilter", (long) SBFilter, (char*)"Type of filters", &status)) + PrintError( status ); + if ( ffpkyj(fptr, (char*)"NORM", (long) TypeNorm, (char*)"normalization", &status)) + PrintError( status ); + // if ( ffpkyj(fptr, "UNDEC", (long) NbrUndec , "nbr of undec. scales", &status)) + // PrintError( status ); + } + + + if (Type_Transform == TO3_LIFTING) + if ( ffpkyj(fptr, (char*)"LiftingTrans", (long) LiftingTrans, + (char*)StringLSTransform(LiftingTrans), &status)) + PrintError( status ); + + if ( ffpkyj(fptr, (char*)"DataFormat",(long) DataFormat,(char*)"Input data format", &status)) + PrintError( status ); + + if ( ffpkyj(fptr, (char*)"Border",(long)Border,(char*)"border type", &status)) + PrintError( status ); + + if (SBFilter == F_USER) + { + if (UserFilterFileName != NULL) + { + if ( ffpkys(fptr, (char*)"FilBank", UserFilterFileName ,(char*)"Filter", &status)) + PrintError( status ); + } + else if ( ffpkys(fptr, (char*)"FilBank", (char*)"-" ,(char*)"Filter", &status)) + PrintError( status ); + } + + return(status); +} + + +/****************************************************************************/ + +void MR_3D::write (char *Name) +/* new version with fits */ +{ + char filename[256]; + Ifloat Ima; + fitsfile *fptr; + int status; + // float *Ptr; + int b,simple; + int bitpix; + long naxis=0; + long naxes[4]; + long nelements; + long group = 1; /* group to write in the fits file, 1= first group */ + long firstpixel = 1; /* first pixel to write (begin with 1) */ + Ifloat Aux; + +/* we keep mr as extension even if its fits ! */ + mr3d_io_name (Name, filename); + +#if DEGUG_IO + cout << "Write on " << filename << endl; +#endif + + FILE *FEXIST = fopen(filename, "rb"); + if (FEXIST) + { + fclose(FEXIST); + remove(filename); /* Delete old file if it already exists */ + } + + status = 0; /* initialize status before calling fitsio routines */ + + /* open the file */ + if ( ffinit(&fptr, filename, &status) ) /* create the new FITS file */ + PrintError( status ); /* call PrintError if error occurs */ + +/* write the header */ + simple = True; + bitpix = -32; /* 32-bit real pixel values */ + long pcount = 0; /* no group parameters */ + long gcount = 1; /* only a single image/group */ + int extend = False; + switch (Set_Transform) + { + case TRANS3_MALLAT: + naxis = 3; + naxes[0] = Nx; + naxes[1] = Ny; + naxes[2] = Nz; + break; + case TRANS3_PAVE: + naxis = 4; + naxes[0] = Nx; + naxes[1] = Ny; + naxes[2] = Nz; + naxes[3] = Nbr_Plan; + break; + default: + fprintf (stderr, "Error in mr_io_write: bad Set_Transform ... \n"); + exit(-1); + break; + } + + // write first header part (parameters) + if ( ffphpr(fptr,simple,bitpix,naxis,naxes,pcount,gcount,extend,&status) ) + PrintError( status ); /* call PrintError if error occurs */ + + // write the header of the multiresolution file + status = mr_io_fill_header(fptr); + + // save the data + // long fpixels[3]; + // long lpixels[3]; + switch (Set_Transform) + { + case TRANS3_MALLAT: + // this works because everything is put in one image + nelements = naxes[0] * naxes[1] * naxes[2]; + if ( ffppre(fptr, group, firstpixel, nelements, Data.buffer(), + &status) ) + PrintError( status ); + break; + case TRANS3_PAVE: + nelements = naxes[0] * naxes[1] * naxes[2]; + for (b=0; b < Nbr_Plan; b++) + { + if ( ffppre(fptr, group, firstpixel, nelements, (TabBand[b]).buffer(), + &status) ) + PrintError( status ); + firstpixel += nelements; + } + break; + default: + fprintf (stderr, "Error in mr_io_write: bad Type_Transform ..\n\n"); + break; + } + + /* close the FITS file */ + if ( ffclos(fptr, &status) ) PrintError(status); +// cout << " end of write fits " << endl; + +} + +/****************************************************************************/ + +void MR_3D::read (char *Name) +{ + // for fits + char filename[256]; + fitsfile *fptr; /* pointer to the FITS file */ + int status=0, hdutype ; + long hdunum; + char comment[FLEN_COMMENT]; + int naxis; + long naxes[4]; + long mon_long; + float nulval = 0.; + int anynul = 0; + long inc[3]; + void PrintError( int status); + long nelements = 0 ; // naxes[0] * naxes[1] in the image + Ifloat Ima; + long firstpixel = 1; + + // for multiresol + float *Ptr; + + mr3d_io_name (Name, filename); + + inc[0]=1; inc[1]=1; inc[2]=1; + +#if DEBUG_IO + cout << "Read in " << filename << endl; +#endif + + /* open the file */ + status = 0; /* initialize status before calling fitsio routines */ + if ( ffopen(&fptr, filename, (int) READONLY, &status) ) + PrintError( status ); + // cout << " open the file " << endl; + + // ******* read the HEADER ******* + + hdunum = 1; /*read table */ + if ( ffmahd(fptr, hdunum, &hdutype, &status) ) /* move to the HDU */ + PrintError( status ); + + int simple, bitpix, extend; + long pcount, gcount; + if ( ffghpr(fptr, 3, &simple, &bitpix, &naxis, naxes, &pcount, + &gcount, &extend, &status)) /* move to the HDU */ + PrintError( status ); + + nelements = naxes[0] * naxes[1] * naxes[2]; +// cout << " begin to read " << endl; + /* read Number of lines, columns, plans */ + //if (ffgkyj(fptr,"Nx", &mon_long, comment, &status)) PrintError( status ); + //int Nxi = (int)mon_long; /* there is no function for reading int */ + //if (ffgkyj(fptr,"Ny", &mon_long, comment, &status)) PrintError( status ); + //int Nyi = (int)mon_long; + //if (ffgkyj(fptr,"Nz", &mon_long, comment, &status)) PrintError( status ); + //int Nzi = (int)mon_long; + int Nxi = (int)naxes[0]; + int Nyi = (int)naxes[1]; + int Nzi = (int)naxes[2]; + + if (ffgkyj(fptr,(char*)"Nbr_Plan", &mon_long, comment, &status)) PrintError( status ); + int NbrPlan = (int)mon_long; + + type_trans_3d TT; + if (ffgkyj(fptr,(char*)"Type_Transform", &mon_long, comment, &status)) + PrintError( status ); + else TT = (type_trans_3d) mon_long; + + FilterAnaSynt *PtrFAS = NULL; + if (TT == TO3_MALLAT) + { + if (ffgkyj(fptr,(char*)"SBFilter", &mon_long, comment, &status)) SBFilter = DEF_SB_FILTER; + else SBFilter = (type_sb_filter) mon_long; + + if (ffgkyj(fptr,(char*)"NORM", &mon_long, comment, &status)) TypeNorm = DEF_SB_NORM; + else TypeNorm = (sb_type_norm) mon_long; + + // if (ffgkyj(fptr,"UNDEC", &mon_long, comment, &status)) NbrUndecimatedScale = -1; + // else NU = (int) mon_long; + + if (SBFilter == F_USER) + { + char FBName[256]; + if (ffgkys(fptr,(char*)"FilBank", FBName, comment, &status)) + PrintError( status ); + if (FBName[0] == '-') UserFilterFileName = NULL; + else UserFilterFileName = strdup(FBName); + } + + PtrFAS = new FilterAnaSynt; + PtrFAS->Verbose = Verbose; + PtrFAS->alloc(SBFilter); + FilterBankAlloc=True; + } + alloc(Nxi,Nyi,Nzi, TT, NbrPlan, PtrFAS, TypeNorm); + + if (Type_Transform == TO3_LIFTING) + { + if (ffgkyj(fptr,(char*)"LiftingTrans", &mon_long, comment, &status)) + LiftingTrans = DEF_LIFT; + else LiftingTrans = (type_lift) mon_long; + } + + if (ffgkyj(fptr,(char*)"DataFormat", &mon_long, comment, &status)) + PrintError( status ); + DataFormat = (type_3d_format) mon_long; + + if (ffgkyj(fptr,(char*)"Border", &mon_long, comment, &status)) + PrintError( status ); + Border = (type_border)mon_long; + + +#if DEBUG_IO + cout << "Read in " << filename << endl; + cout << "Nx = " << Nx << endl; + cout << "Ny = " << Ny << endl; + cout << "Nz = " << Nz << endl; + cout << "Nbr_Plan = " << nbr_scale () << endl; + cout << "Type_Transform = " << Type_Transform << " " << + StringTransf3D(Type_Transform) << endl; + cout << "Set_Transform = " << Set_Transform << " " << + StringSet3D(Set_Transform) << endl; +#endif + + // ******* read the images ******* + + long fpixels[3]; + long lpixels[3]; + fpixels[0] = 1; + fpixels[1] = 1; + fpixels[2] = 1; + lpixels[0] = Nx; + lpixels[1] = Ny; + lpixels[2] = Nz; + + switch (Set_Transform) + { + case TRANS3_MALLAT: + Ptr = Data.buffer(); + if ( ffgpve(fptr, 1, 1, nelements, nulval,Ptr, &anynul, &status)) + PrintError( status ); + break; + case TRANS3_PAVE: + nelements = naxes[0] * naxes[1] * naxes[2]; + for (int b=0; b < Nbr_Plan; b++) + { + Ptr = (TabBand[b]).buffer(); + if ( ffgpve(fptr, 1, firstpixel, nelements, nulval,Ptr, &anynul, &status)) + PrintError( status ); + firstpixel += nelements; + } + break; + default: + fprintf (stderr, "Error in mr_io_read: bad Set_Transform .. \n"); + break; + } + +/* close the FITS file */ + if ( ffclos(fptr, &status) ) PrintError( status ); +// cout << " end of read fits file " << endl; +} + +/****************************************************************************/ + + diff --git a/src/libsparse3d/MR3D_Obj.h b/src/libsparse3d/MR3D_Obj.h new file mode 100755 index 0000000..4774428 --- /dev/null +++ b/src/libsparse3d/MR3D_Obj.h @@ -0,0 +1,247 @@ + +/****************************************************************************** +** Copyright (C) 1999 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: J.L. Starck +** +** Date: 19.8.99 +** +** File: MR3D_Obj.h +** +******************************************************************************* +** +** DESCRIPTION Class for 3D wavelet transform +** ----------- +** +******************************************************************************/ + +#ifndef _CMR3D_H_ +#define _CMR3D_H_ + +#include "IM_Obj.h" +#include "IM_IO.h" +#include "SB_Filter.h" +#include "IM3D_IO.h" +#include "Atrou3D.h" +#include "MR3D_Obj.h" + +#define NBR_TRANS_3D 3 +#define DEFAULT_BORDER_3D I_MIRROR +#define MAX_SCALE_1D 100 + +enum type_trans_3d {TO3_MALLAT,TO3_LIFTING, TO3_ATROUS, T3_UNDEFINED=-1}; + +enum set_trans_3d {TRANS3_MALLAT, TRANS3_PAVE, S3_UNDEFINED=-1}; +inline char * StringTransf3D (type_trans_3d type) +{ + switch (type) + { + case TO3_MALLAT: + return ((char*)"(bi-) orthogonal transform ");break; + break; + case TO3_LIFTING: + return ((char*)"(bi-) orthogonal transform via lifting sheme"); + break; + case TO3_ATROUS: + return ((char*)"A trous wavelet transform"); + break; + case T3_UNDEFINED: + return ((char*)"undefined transform");break; + } + return ((char*)"Error: bad type of transform"); +} + +inline char * StringSet3D (set_trans_3d type) +{ + switch (type) + { + case TRANS3_MALLAT: + return ((char*)"non-redundant"); + break; + case TRANS3_PAVE: + return ((char*)"redundant"); + break; + case S3_UNDEFINED: + return ((char*)"undefined type transform");break; + } + return ((char*)"Error: bad type of transform"); +} + +inline set_trans_3d which_set_is_trans3d(type_trans_3d type) +{ + switch (type) + { + case TO3_MALLAT: + case TO3_LIFTING: + return TRANS3_MALLAT; break; + case TO3_ATROUS: + return TRANS3_PAVE; break; + case T3_UNDEFINED: + return S3_UNDEFINED;break; + } + return S3_UNDEFINED; +} + +inline set_trans_3d SetTransform (type_trans_3d Transform) +{ + return which_set_is_trans3d(Transform); +} + + +/************************************************************************/ + +// 3D sub-band decomposision: an image is transformed into eight sub-cubes +class SubBand3D +{ + private: + SubBand1D *Ptr_SB1D; + public: + SubBand3D(SubBand1D &SB1D) {Ptr_SB1D = &SB1D;} + void transform3d(fltarray &Data); + void transform3d(fltarray &Data, fltarray* Trans, int Step); + void recons3d(fltarray &Data); + void recons3d(fltarray* Trans, fltarray &Data, int Step); + ~SubBand3D(){Ptr_SB1D = NULL;} +}; + +/*********************************************************************** +* Orthogonal 3D wavelet transform (decimated) +************************************************************************/ + +// Orthogonal 3D wavelet transform +class Ortho_3D_WT: public SubBand3D { + public: + Ortho_3D_WT(SubBand1D &SB1D):SubBand3D(SB1D) {}; + ~Ortho_3D_WT() {}; + void transform (fltarray &Cube, int Nbr_Plan); + void transform(fltarray &Cube_in, fltarray &Transf_Out, int Nbr_Plan); + void recons(fltarray &Cube, int Nbr_Plan); + void recons(fltarray &Transf_in, fltarray &Cube_Out, int Nbr_Plan); +}; + +/*********************************************************************** +* Orthogonal Undecimated 3D wavelet transform +************************************************************************/ + +class PAVE_3D_WT: public SubBand3D +{ + public: + PAVE_3D_WT(SubBand1D &SB1D):SubBand3D(SB1D) {}; + ~PAVE_3D_WT() {}; + int alloc (fltarray * & TabBand, int Nx, int Ny, int Nz, int NbrScale); + void free(fltarray *TabBand, int NbrScale); + void one_scale_transform (fltarray &Imag, fltarray *TabTrans, int Step, int Pos=0); + void transform (fltarray &Cube, fltarray *TabTrans, int NbrScale); + void recons(fltarray *TabTrans, fltarray &Cube, int NbrScale); +}; + +/************************************************************************/ + +class MR_3D { + ATROUS_3D_WT AT3D_WT; // Wavelet transform class for the + // 3D wavelet transform + fltarray Data; // buffer where the wavelet transform is stored + // only used for non redundant transform + fltarray *TabBand; // used for the redundant transform. + // otherwise TabBand = &Data + int Nbr_Plan; // number of scales + int Nbr_Band; // number of band + int Nx,Ny,Nz; // number of points of the input signal + + // for non redundant transform only + int *TabPosX; // TabPos and TabSize gives the position + int *TabPosY; // and the size of a given band. + int *TabPosZ; // used only for non-redundant transform + int *TabSizeNx; // (i.e. size(signal) = size(transform)) + int *TabSizeNy; + int *TabSizeNz; + + type_3d_format DataFormat; // input data format + int mr_io_fill_header(fitsfile *fptr); + // write to the MR file header information + + FilterAnaSynt *FilterBank; // pointer to the filter bank to use + // in the orthogonal wavelet transform + Bool FilterBankAlloc; // True if the filter bank is allocated + // by the class + void init(); + public: + Bool Verbose; + type_trans_3d Type_Transform; // type of transform + set_trans_3d Set_Transform; // class of transform + type_border Border; // type of border to user for the + // border intrpolation + type_sb_filter SBFilter; // Filter used in the subband + // decomposition (case of Mallat transform). + sb_type_norm TypeNorm; // type of normalization + // (in L1 or L2) + type_lift LiftingTrans; // type of lifting (in case of lifting + // transform + + // return a pointer to the filter bank + FilterAnaSynt * filter_bank() {return FilterBank;} + + MR_3D (){init();} + MR_3D (int Nx, int Ny, int Nz, type_trans_3d T, int Nbr_Scale) + { alloc(Nx,Ny,Nz,T,Nbr_Scale);} + + inline int nbr_scale () const {return Nbr_Plan;} + inline int nbr_band () const {return Nbr_Band;} + int size_cube_nx () const {return Nx;} + int size_cube_ny () const {return Ny;} + int size_cube_nz () const {return Nz;} + int size_band_nx (int b) const {return TabSizeNx[b];}; + int size_band_ny (int b) const {return TabSizeNy[b];}; + int size_band_nz (int b) const {return TabSizeNz[b];}; + void alloc (int Nx, int Ny, int Nz, type_trans_3d T, int Nbr_Scale, + FilterAnaSynt *FAS=NULL, sb_type_norm Norm=NORM_L1); + + // return the data set for non redundant transform + fltarray & cube() {return Data;} + + // return the data set for redundant transform + fltarray * & TBand() {return TabBand;} + + float & operator() (int b, int x, int y, int z) const; + + void transform (fltarray &Image, type_border Border); + // apply the wavelet transform to a 1D signal + // using a given border interpolation + + void transform (fltarray &Image); + // apply the wavelet transform to a 1D signal + // using the default border interpolation + + void recons (fltarray &Image, type_border Border); + // reconstruction of the 1D signal from its wavelet coefficient + // using a given border interpolation + + void recons (fltarray &Image); + // reconstruction of the 1D signal from its wavelet coefficient + // using the default border interpolation + + void get_band(int b, fltarray &Band); + // return in Band a band of the transform + + void insert_band(int b, fltarray &Band); + // insert a Band a band of the transform + + void free (); + // deallocate the object + + void read(char *FileName); // read the MR3D object from the disk + void write(char *FileName); // write the MR3D object to the disk + + void info_pos_band(); // print information about size and position + // of each band + void info_band(int b); // print statistical information + // about the band b + ~MR_3D(); +}; + +#endif diff --git a/src/libsparse3d/Mr3d_FewEvent.cc b/src/libsparse3d/Mr3d_FewEvent.cc new file mode 100644 index 0000000..e112fe8 --- /dev/null +++ b/src/libsparse3d/Mr3d_FewEvent.cc @@ -0,0 +1,1578 @@ +//****************************************************************************** +//** +//** DESCRIPTION Few Event Poisson Thresholding +//** ----------- +//** +//****************************************************************************** + + +#include "Mr3d_FewEvent.h" + +//****************************************************************************** +FewEventPoisson::FewEventPoisson ( Bool InitOk, int NAutoConv, float epsilon) { +// initialisation of datas + _Epsilon=epsilon; + _Param.alloc(2); + _InitOk = InitOk; + if(_InitOk == True){ + fits_read_fltarr(Name_Param, _Param); + _NbAutoConv = (int)(_Param(1)+0.5); + }else{ + _NbAutoConv= NAutoConv; + } + _HistoConv.alloc (_NbAutoConv+1, MAXHISTONBPOINT); + _HistoConv.init(0.); + _HistoBound.alloc (_NbAutoConv+1, 2); + _HistoBin.alloc (_NbAutoConv+1, 2); + _HistoDistrib.alloc (_NbAutoConv+1, MAXHISTONBPOINT); + _Threshold.alloc(_NbAutoConv+1, 2, ""); +} + + + +//****************************************************************************** +void FewEventPoisson::compute_distribution (Bool WriteAllInfo, Bool WriteHisto, + Bool Verbose, int param) { + // compute histogramms distributions and the values of threshold + + // compute Bspline Histogramm + if (Verbose==True) + cout << "Compute the histogram of the Wavelet ... " << endl; + bspline_histo_3D (WriteAllInfo, Verbose); + + // computes the auto-convolued histogramms + if (Verbose==True) cout << "Compute the autoconvolutions of the histogram ... " << endl; + histo_convolution (WriteAllInfo,Verbose,param); + + + // computes the distribution function of histogramms at each leve + histo_distribution (WriteAllInfo,Verbose); + + _Param(0)=param; + _Param(1)=_NbAutoConv; + + // writing debug result + if ((WriteHisto)||(WriteAllInfo==True)) { + fits_write_fltarr (Name_HistoConv, _HistoConv); + fits_write_fltarr (Name_HistoBound, _HistoBound); + fits_write_fltarr (Name_HistoBin, _HistoBin); + fits_write_fltarr(Name_HistoDistrib, _HistoDistrib); + fits_write_fltarr(Name_Param, _Param); + } + + _InitOk = True; +} + + + +//****************************************************************************** +void FewEventPoisson::bspline_histo_3D (Bool WriteAllInfo, Bool Verbose) { + // create the histogramm of the wavelet B3Spline + + // init Phi cubic bspline + // performs B(x) for x in [-2;+2] : + // B(x) =1/12 (|x-2|^3 - 4 |x-1|^3 + 6 |x|^3 - 4 |x+1|^3 + |x+2| ^ 3) + fltarray phi_bspline(BIN3D_2); + float x; + for (int i=-BIN3D/2;i<=BIN3D/2;i++) { + + x = 2 * ((float)i)/((float)BIN3D/2); + phi_bspline(i+BIN3D) = ( CUBE((ABS(x-2))) - 4 * CUBE((ABS(x-1))) + + 6 * CUBE((ABS(x))) - 4 * CUBE((ABS(x+1))) + + CUBE((ABS(x+2))) ) / 12; + } + + if (WriteAllInfo==True) { + fits_write_fltarr (Name_Bspline_, (phi_bspline)); + if (Verbose==True) cout << "File "<< Name_Bspline_ << " is created (Bspline)"<< endl; + } + + // performs Psi wavelet + // psi(x,y,z) = B(x,y,z) - 1/8 B(x/2,y/2,z/2) + // psi(x,y) = B(x) B(y) B(z) - 1/8 B(x/2) B(y/2) B(z/2) + + + float MaxPsi = 0; + float MinPsi = 0; + float val; + for (int i=-BIN3D;i<=BIN3D;i++) { + for (int j=-BIN3D;j<= BIN3D;j++) + for (int k=-BIN3D;k<= BIN3D;k++) { + val=phi_bspline(i+BIN3D) * phi_bspline(j+BIN3D) * phi_bspline(k+BIN3D) + - 0.125 * phi_bspline(i/2+BIN3D) * phi_bspline(j/2+BIN3D) * phi_bspline(k/2+BIN3D); + if (val >= MaxPsi) + MaxPsi = val; + if (val <= MinPsi) + MinPsi = val; + } + } + + // Fullfills the histogramm (first length = 1025 points) + // + int index; + int LocalSizeSignal = LENGTH_FIRST_HISTO; + for (int i=-BIN3D;i<=BIN3D;i++) { + for (int j=-BIN3D;j<=BIN3D;j++) { + for (int k=-BIN3D;k<=BIN3D;k++) { + // we only use (HSP-1)/2 point in _HistoConv tab + val=phi_bspline(i+BIN3D) * phi_bspline(j+BIN3D) * phi_bspline(k+BIN3D) + - 0.125 * phi_bspline(i/2+BIN3D) * phi_bspline(j/2+BIN3D) * phi_bspline(k/2+BIN3D); + index = (int) ( (val-MinPsi)*(LocalSizeSignal-1) + / (MaxPsi-MinPsi)); + _HistoConv(0,index)++; // first histogram, no convolution + } + } + } + + // bin + float LocalHistoBin = (MaxPsi-MinPsi) / (LocalSizeSignal-1); + + // set min, max, bin and nb points of first level (histogram only) + // + _HistoBound(0,0) = MinPsi; + _HistoBound(0,1) = MaxPsi; + _HistoBin (0,0) = LocalHistoBin; + _HistoBin (0,1) = LocalSizeSignal; + + float sum_check=0.0; + for(int i=0;i<_HistoBin(0,1);i++) + sum_check += _HistoConv(0,i)*_HistoBin (0,0); + // it must be a probability density -> normalisation by sum_check + for(int i=0;i<_HistoBin(0,1);i++) + _HistoConv(0,i) /= sum_check; + + + + // normalisation of the first histogram + + //Reduction to have a sigma = 1 and mean = 0 + //Dilatation on x axes + // the histogramme must have a mean=0, and the bin is divided by sigma (std) + // calculation of the mean and the sigma (std) + float mean=0.0; + float std=0.0; + for (int j=0; j<_HistoBin(0,1);j++) { + float real_value = ((float) j)*_HistoBin (0,0) + _HistoBound(0,0); + mean += _HistoConv(0,j)*_HistoBin (0,0)*real_value; + } + + for(int j=0; j<_HistoBin(0,1);j++){ + float real_value = ((float) j)*_HistoBin (0,0) + _HistoBound(0,0); + std += _HistoConv(0,j)*_HistoBin (0,0) + *(real_value-mean)*(real_value-mean); + } + + if (std < 0.0) std = 0.0; + else std = sqrt(std); + + //if (Verbose==True) cout << "Old sigma = " << std << " and old mean = " + // << mean << endl; + + //Update of min, max and bin + //offset of min and max and dilatation of 1/std + _HistoBound(0,0)=(_HistoBound(0,0)-mean)/std; + _HistoBound(0,1)=(_HistoBound(0,1)-mean)/std; + //dilatation of 1/std + _HistoBin (0,0)/=std; + + //we do not have a density of probability anymore becauseof the new bin + //Reduction to have a density of probability + // sum must be equal to 1 + sum_check=0.0; + for(int i=0;i<_HistoBin(0,1);i++) + sum_check += _HistoConv(0,i)*_HistoBin (0,0); + // it must be a probability density -> normalisation by sum_check + for(int i=0;i<_HistoBin(0,1);i++) + _HistoConv(0,i) /= sum_check; + + +#if DEBUG_FEW + if (Verbose==True) show_param ("End Convol Compute", 0, + _HistoBound(0,0), _HistoBound(0,1), + _HistoBin(0,0), _HistoBin(0,1)); + if (WriteAllInfo==True) { + char FileName[256]; + fltarray Courbe ((int)(_HistoBin(0,1))); + for (int i=0;i<_HistoBin(0,1);i++) Courbe(i)=_HistoConv(0,i); + sprintf (FileName, "Histo_%d.fits", 0); + fits_write_fltarr (FileName, Courbe); + if (Verbose==True) cout << "File "<< FileName << " is created .First Convolution."<< endl; + } +#endif + +} + + + +//****************************************************************************** +void FewEventPoisson::histo_convolution (Bool WriteAllInfo, Bool Verbose, int param) { +// compute all convolutions + + + // param gives the rythm of autoconvolution. + // during 2*(2^param) all autoconvolution are computed: + // if param=2, 2*(2^param)=8 + // 1 2 3 4 ... 8 + // S=S1 S*S=S2 S*S*S=S3 S*S*S*S=S4 ... S*S*S*S*S*S*S*S=S8 + // the next 2^param convolutions are autoconvolutions of the 2^param before: + // S9=S5*S5 S10=S6*S6 S11=S7*S7 S12=S8*S8 + // same thing for the 4 next... + // S13=S9*S9 S14=S10*S10 S15=S11*S11 S16=S12*S12 ... + + + + // length of first histogram h(0,*) + int HistoNbPoint= LENGTH_FIRST_HISTO; + int param_2= 1 << param;// param_2=2^param + + int HistoNbPointBegin = HistoNbPoint*(param_2*2-1); + if(param<0) param=0; + + // Reduced Histo used for computation + // init with h0 (_HistoConv(0,*)) + int MaxHistoNbPoint2=MAXHISTONBPOINT; + MaxHistoNbPoint2*=2; + fltarray ReducedHisto (MaxHistoNbPoint2); + fltarray ReducedHistoBegin (param_2*2-1,HistoNbPointBegin ); + fltarray ReducedHistoTransi (HistoNbPointBegin); + ReducedHisto.init(0.); + ReducedHistoBegin.init(0.); + ReducedHistoTransi.init(0.); + + for (int i=0; i=0)) + ReducedHistoBegin(nbConv,i) += ReducedHistoBegin(IndexFirstHisto,i-j) + *ReducedHistoBegin(IndexSecHisto,j); + + + _HistoBin (nbConv,1) = HistoNbPoint; + + // new_min = min_first_histo + min_sec_histo + _HistoBound(nbConv,0)=_HistoBound(IndexFirstHisto,0)+_HistoBound(IndexSecHisto,0); + // new_max = max_first_histo + max_sec_histo + _HistoBound(nbConv,1)=_HistoBound(IndexFirstHisto,1)+_HistoBound(IndexSecHisto,1); + + + //the bin do not change + _HistoBin (nbConv,0) = LocalBin; + if(Verbose==True){ + cout << "nbConv= " << nbConv << endl; + cout << "Min= " << _HistoBound(nbConv,0) <<" Max= " << _HistoBound(nbConv,1) << endl; + cout << "Bin= " << _HistoBin (nbConv,0) <<" Nb Points= " + << _HistoBin (nbConv,1) << endl< MAXHISTONBPOINT) we have to decimate + // we used a filter with kernel [0.25, 0.5, 0.25] (function shape_signal) + + shape_signal (ReducedHistoTransi); + HistoNbPoint = (HistoNbPoint-1)/2 + 1; + LocalBin = (_HistoBound(nbConv,1)-_HistoBound(nbConv,0)) / (HistoNbPoint-1); + _HistoBin(nbConv,1) = HistoNbPoint; + _HistoBin(nbConv,0) = LocalBin; + + if(Verbose==True){ + cout << "Histo reshaped:"<=0)) + ReducedHisto(i) += _HistoConv(nbOldConv,i-j) + *_HistoConv(nbOldConv,j); + + _HistoBound(nbConv,0) = 2.*_HistoBound(nbOldConv,0); + _HistoBound(nbConv,1) = 2.*_HistoBound(nbOldConv,1); + + + + _HistoBin (nbConv,1) = HistoNbPoint; + ////////LocalBin = _HistoBin(nbOldConv,0); !!!!!!!!!!!1 + LocalBin = (_HistoBound(nbConv,1)-_HistoBound(nbConv,0)) / (float)(HistoNbPoint); + _HistoBin(nbConv,0) = LocalBin; + + float sum=0.; + for (long i=0; i MAXHISTONBPOINT) we have to decimate + // we used a filter with kernel [0.25, 0.5, 0.25] (function shape_signal) + + shape_signal (ReducedHisto); + HistoNbPoint = (HistoNbPoint-1)/2 + 1; + LocalBin = (_HistoBound(nbConv,1)-_HistoBound(nbConv,0)) / (HistoNbPoint-1); + _HistoBin(nbConv,0) = LocalBin; + _HistoBin(nbConv,1) = HistoNbPoint; + + if(Verbose==True){ + cout << "Histo reshaped"< // " << text << endl; + cout << " -min:" << min << ", -max:" << max << "' bin:" << bin + << ", (NBP:" << nbechant << ")" << endl; +} + +//****************************************************************************** +void FewEventPoisson::histo_threshold (Bool WriteAllInfo, Bool Verbose) { +// compute the threshold max and the threshold min for every +// distribution laws. + + + for (int nbConv=0; nbConv<_NbAutoConv+1; nbConv++) { + + // serach indice of max threshold (RepFunc >= 1-eps) + // and min threshold (RepFunc <= eps) + // + long IndexMax=0; + long IndexMin=(long)(_HistoBin(nbConv,1)+0.5)-1; + + + // _HistoDistrib is a distribution of probability + // _HistoDistrib(min)=~0 _HistoDistrib(max)=~1 + // RepFuc = 1-eps in [IndexMax-1:IndexMax] + while ( (_HistoDistrib(nbConv,IndexMax) < 1-_Epsilon) + && (IndexMax<_HistoBin(nbConv,1)-1)) {IndexMax++;} + // RepFuc = eps in [IndexMin:IndexMin+1] + while ( (_HistoDistrib(nbConv,IndexMin) > _Epsilon) + && (IndexMin>0)) {IndexMin--;} + //cout << "Ind min:" << IndexMin <<", IndMax : " << IndexMax << endl; + // test + if ((IndexMax==0) || (IndexMin==_HistoBin(nbConv,1)-1)) { + cout << "ERROR : bad sampling of the distribution function " << endl; + exit (-1); + } + + // inter between IndexMax and IndexMax-1 + // X1-eps = Xmax + (1-eps-Ymax)/a, a=(Ymax-1 - Ymax)/(Xmax-1 - Xmax) + // + float LocalBin=_HistoBin(nbConv,0); + float Min=_HistoBound(nbConv,0); + + float a = _HistoConv(nbConv,IndexMax);//*Sigma; + + //float b = _HistoDistrib(nbConv+2,IndexMax-1) + // - a * _HistoDistrib(nbConv+1,IndexMax-1); + + //_Threshold (nbConv, 1) = _HistoDistrib(nbConv+1,IndexMax) + // + (1-_Epsilon-_HistoDistrib(nbConv+2,IndexMax))/a; + + _Threshold (nbConv, 1) = (IndexMax*LocalBin+Min)///Sigma + +(1-_Epsilon-_HistoDistrib(nbConv,IndexMax))/a; + + + //_Threshold (1, nbConv) = (1-_Epsilon-b)/a; + + + // inter between IndexMin and IndexMin+1 + // Xeps = Xmin + (eps-Ymin)/a, a=(Ymin+1 - Ymin)/(Xmin+1 - Xmin) + a = _HistoConv(nbConv,IndexMin+1);//*Sigma; + + _Threshold (nbConv, 0) = (IndexMin*LocalBin+Min)///Sigma + + (_Epsilon-_HistoDistrib(nbConv,IndexMin))/a; + + + + if (Verbose==True) { + cout << "Nbr of histogram autoconv.:" << nbConv << ", SeuilMin := " << _Threshold (nbConv, 0) + << ", SeuilMax := " << _Threshold (nbConv, 1) + << endl; + } + } + if (WriteAllInfo==True) { + io_write_ima_float (Name_Threshold, _Threshold); + if(Verbose==True) cout << "File "<< Name_Threshold << " is created (Threshold) ." << endl; + + } +} + +//****************************************************************************** + +void FewEventPoisson::find_threshold (Bool WriteAllInfo, Bool WriteHisto, + Bool Verbose,int param) { +// test if the thresholds are computed or read. + if(Verbose==True) cout << "Initalisation of all histograms" << endl; + if (_InitOk == False) { + if(Verbose==True) cout <<"Computing Wavelet Distribution"< 0) + { + SeuilMin *= sqrt((float)(NEventReal)) * sigma_wavelet / Alpha; + SeuilMax *= sqrt((float)(NEventReal)) * sigma_wavelet / Alpha; + } + else + { + SeuilMin *= sigma_wavelet / Alpha; + SeuilMax *= sigma_wavelet / Alpha; + } +} + + +//****************************************************************************** + +void FewEventPoisson::event_set_support( +// create the support with the threshold values and the wavelet coeff + + + fltarray *& Data_In, + int CurrentScale, + int FirstDectectScale, + int NEGFirstDectectScale, + int MinEventNumber, + Bool OnlyPositivDetect, + type_border Border, + //number of events + intarray *& Nb_Event, + //cube of support of the current scale + intarray *& Support, + Bool WriteAllInfo) +{ + int FDS = FirstDectectScale; + int NFDS = (NEGFirstDectectScale < 0) ? FirstDectectScale: NEGFirstDectectScale; + int Nx = Data_In[CurrentScale].nx(); + int Ny = Data_In[CurrentScale].ny(); + int Nz = Data_In[CurrentScale].nz(); + + for (int i=0;i=SeuilMax)) + { + + Support[CurrentScale](i,j,k) = VAL_SupOK; + + // + if (Nb_Event[CurrentScale](i,j,k) < MinEventNumber) + Support[CurrentScale](i,j,k) = VAL_SupMinEv; + + // + if ((OnlyPositivDetect) && (Data_In[CurrentScale](i,j,k)<0)) + Support[CurrentScale](i,j,k) = VAL_SupNull; + + // + if ((FDS > CurrentScale) && (Data_In[CurrentScale](i,j,k) >= 0)) + Support[CurrentScale](i,j,k) = VAL_SupFirstScale; + + if ((NFDS > CurrentScale) && (Data_In[CurrentScale](i,j,k) < 0)) + Support[CurrentScale](i,j,k) = VAL_SupFirstScale; + } + } +} + +//****************************************************************************** +int FewEventPoisson::get_pix(int i,int j,int k, intarray & DataIn, type_border Border){ +// give the value of a pixel depending to the border + int Nx = DataIn.nx(); + int Ny = DataIn.ny(); + int Nz = DataIn.nz(); + int ind_i=0, ind_j=0, ind_k=0; + switch(Border){ + case I_CONT: + ind_i=(i < 0 ? 0 : i); + if(ind_i>=Nx) ind_i=Nx-1; + ind_j=(j < 0 ? 0 : j); + if(ind_j>=Ny) ind_j=Ny-1; + ind_k=(k < 0 ? 0 : k); + if(ind_k>=Nz) ind_k=Nz-1; + break; + case I_MIRROR: + ind_i=(i < 0 ? -i : i); + if(ind_i>=Nx) ind_i=2*(Nx-1)-ind_i; + ind_j=(j < 0 ? -j : j); + if(ind_j>=Ny) ind_j=2*(Ny-1)-ind_j; + ind_k=(k < 0 ? -k : k); + if(ind_k>=Nz) ind_k=2*(Nz-1)-ind_k; + break; + case I_ZERO: + ind_i=(i < 0 ? 0 : i); + if(ind_i>=Nx) ind_i=0; + ind_j=(j < 0 ? 0 : j); + if(ind_j>=Ny) ind_j=0; + ind_k=(k < 0 ? 0 : k); + if(ind_k>=Nz) ind_k=0; + break; + case I_PERIOD: + ind_i=(i < 0 ? Nx+i : i); + if(ind_i>=Nx) ind_i=i-Nx; + ind_j=(j < 0 ? Ny+j : j); + if(ind_j>=Ny) ind_j=j-Ny; + ind_k=(k < 0 ? Nz+k : k); + if(ind_k>=Nz) ind_k=k-Nz; + break; + } + return( DataIn(ind_i, ind_j, ind_k)); +} +//****************************************************************************** +int FewEventPoisson::get_ev(int i,int j,int k,intarray & nb_event, type_border Border){ + //get the number of events is a box already computed depending to the border + + int Nx = nb_event.nx(); + int Ny = nb_event.ny(); + int Nz = nb_event.nz(); + int ind_i, ind_j, ind_k; + int Return=0; + switch(Border){ + case I_CONT: + cout << "I_CONT border is not implemented in case of no approximation "<=Nx) ind_i=2*Nx-1-ind_i; + ind_j=(j < 0 ? -j+1 : j); + if(ind_j>=Ny) ind_j=2*Ny-1-ind_j; + ind_k=(k < 0 ? -k+1 : k); + if(ind_k>=Nz) ind_k=2*Nz-1-ind_k; + Return=nb_event(ind_i, ind_j, ind_k); + break; + case I_ZERO: + if(i<0 || i>=Nx || j<0 || j>=Ny || k<0 || k>=Nz) + Return=0; + else + Return=nb_event(i, j, k); + break; + case I_PERIOD: + ind_i=(i < 0 ? Nx+i : i); + if(ind_i>=Nx) ind_i=i-Nx; + ind_j=(j < 0 ? Ny+j : j); + if(ind_j>=Ny) ind_j=j-Ny; + ind_k=(k < 0 ? Nz+k : k); + if(ind_k>=Nz) ind_k=k-Nz; + Return=nb_event(ind_i, ind_j, ind_k); + break; + } + return( Return); +} + + + + + +//****************************************************************************** +void FewEventPoisson::event_scale_3D_Approx ( + intarray & last_nb_event, + intarray & nb_event, + int Current_Scale, + type_border Border, + Bool WriteAllInfo) { +//compute the number of event in boxes as event_scale_3D but with the approxination +// of box = 2^n * 2^n * 2^n and not box=(2^n+1) * (2^n+1) * (2^n+1) +// the convention is nb_event[n](i,j,k)=sum( cur_i from i-2^(n+2) to i+2^(n+2)-1, +// cur_j from j-2^(n+2) to j+2^(n+2)-1, +// cur_k from k-2^(n+2) to k+2^(n+2)-1, +// last_nb_event(cur_i, cur_j, cur_k) + int Nx = last_nb_event.nx(); + int Ny = last_nb_event.ny(); + int Nz = last_nb_event.nz(); + if (Current_Scale==-2) + { + int * tab; + tab= new int[2]; + //*****************computation of the first point: (0, 0, 0)**************** + //(0,0,0) + nb_event(0,0,0)=0; + for (int cur_i=-1;cur_i<=0;cur_i++) + for (int cur_j=-1;cur_j<=0;cur_j++) + for (int cur_k=-1;cur_k<=0;cur_k++) + nb_event(0,0,0)+=get_pix(cur_i, cur_j, cur_k, last_nb_event, Border); + + //*****************computation of the first ligne: (1:Nx-1 ,0 ,0)********** + //(1:2,0,0) + for (int i=1;i<=2;i++) + { + tab[i-1]=0; + nb_event(i,0,0)=nb_event(i-1,0,0); + for (int cur_j=-1;cur_j<=0;cur_j++) + for (int cur_k=-1;cur_k<=0;cur_k++){ + tab[i-1]+=get_pix(i, cur_j, cur_k, last_nb_event, Border); + nb_event(i,0,0)-=get_pix(i-2, cur_j, cur_k, last_nb_event, Border); + + } + nb_event(i,0,0)+=tab[i-1]; + } + //(3:Nx-1,0,0) + for (int i=3;i=0 && j-ind>=0 && k-ind>=0) + total += nb_event[s-1](i-ind,j-ind,k-ind) + +nb_event[s-1](i-ind,j-ind,k+ind) + +nb_event[s-1](i-ind,j+ind,k-ind) + +nb_event[s-1](i-ind,j+ind,k+ind) + +nb_event[s-1](i+ind,j-ind,k-ind) + +nb_event[s-1](i+ind,j-ind,k+ind) + +nb_event[s-1](i+ind,j+ind,k-ind) + +nb_event[s-1](i+ind,j+ind,k+ind); + + for(int cur_i=i-ind_2; cur_i <= i+ind_2; cur_i++) + for(int cur_j=j-ind_2; cur_j <= j+ind_2; cur_j++) + total -= get_pix(cur_i,cur_j,k,data_event, Border); + + for(int cur_i=i-ind_2; cur_i <= i+ind_2; cur_i++) + for(int cur_k=k-ind_2; cur_k <= k+ind_2; cur_k++) + total -= get_pix(cur_i,j,cur_k,data_event, Border); + + for(int cur_j=j-ind_2; cur_j <= j+ind_2; cur_j++) + for(int cur_k=k-ind_2; cur_k <= k+ind_2; cur_k++) + total -= get_pix(i,cur_j,cur_k,data_event, Border); + + for(int cur_i=i-ind_2; cur_i <= i+ind_2; cur_i++) + total -= get_pix(cur_i,j,k,data_event, Border); + + for(int cur_j=j-ind_2; cur_j <= j+ind_2; cur_j++) + total -= get_pix(i,cur_j,k,data_event, Border); + + for(int cur_k=k-ind_2; cur_k <= k+ind_2; cur_k++) + total -= get_pix(i,j,cur_k,data_event, Border); + + total -= get_pix(i,j,k,data_event, Border); + nb_event[s](i,j,k) = total; + } + } + } + +}*/ diff --git a/src/libsparse3d/Mr3d_FewEvent.h b/src/libsparse3d/Mr3d_FewEvent.h new file mode 100644 index 0000000..e5add52 --- /dev/null +++ b/src/libsparse3d/Mr3d_FewEvent.h @@ -0,0 +1,144 @@ + + +#ifndef _MR3DFEWEVENT_H_ +#define _MR3DFEWEVENT_H_ + +#include "MR1D_Obj.h" +#include "MR_Obj.h" +#include "IM_IO.h" +#include "IM3D_IO.h" +#include "IM_Noise.h" +#include + + +// for wavelet psi +#define BIN3D 150 +#define BIN3D_2 2*BIN3D+1 + +//sigma wavelet 3D +#define SIGMA_BSPLINE_WAVELET 0.0102225 + +// length max of autoconvolution +#define MAXHISTONBPOINT 2048+1 + +// length of first histo (construct without convolution) +#define LENGTH_FIRST_HISTO 1024+1 + +// output names +#define Name_Bspline_ "Aba_bspline" +#define Name_Wavelet_ "Aba_wavelet" +#define Name_HistoConv "_Aba_histo" +#define Name_Threshold "_threshold" +#define Name_HistoBin "_histobin" +#define Name_HistoBound "_histobound" +#define Name_HistoDistrib "_histodistrib" +#define Name_Param "_param" + + +#define CUBE_FABS(x) (fabs(x) * fabs(x) * fabs(x)) +#define CUBE(x) ((x) * (x) * (x)) +#define DEBUG_FEW 1 + + +//definition of differents values of the support +#define VAL_SupNull 0 +#define VAL_SupOK 1 +#define VAL_SupDill 2 +#define VAL_SupEdge 3 +#define VAL_SupLastOK 9 +#define VAL_SupKill 10 +#define VAL_SupMinEv 11 +#define VAL_SupFirstScale 12 + + +class FewEventPoisson { + +private: + float _Epsilon; //epsilon for threshold + int _NbAutoConv; // number max of autoconv used + bool _InitOk; // initialisation flag + + // Computes the one dimensional bspline + // Computes the histogram of the wavelet in histo(0,*) + void bspline_histo_3D (Bool WriteAllInfo=False, Bool Verbose=False); + + // performs auto-convolution of the first histogram. + void histo_convolution (Bool WriteAllInfo=False, Bool Verbose=False, int param=0); + + // Normalizes the histograms + //void histo_normalisation (Bool WriteAllInfo=False, Bool Verbose=False); + + // Computes the repartition function + void histo_distribution (Bool WriteAllInfo=False, Bool Verbose=False); + + void shape_signal (fltarray &Signal1d); + + void show_param (char* text, float nbconv, float min, float max, + float bin, float nbechant); + + int get_ev(int i,int j,int k,intarray & nb_event, type_border Border); + +public: + fltarray _Param; // _Param(0) =value of rytme + // _Param(1) =value of _NbAutoConv + fltarray _HistoBound; // min and max reduced coefficient + fltarray _HistoBin; // bin and number of points + fltarray _HistoConv; // Histogram information on autoconv + fltarray _HistoDistrib; // Distribution of Histograms + Ifloat _Threshold; // Threshold(p, 0) = Min value threshold + // Threshold(p, 1) = Max value threshold + + // init attribute + FewEventPoisson(Bool InitOk, int AutoConv,float epsilon); + + + // Computes the histogram autoconvolutions, and the repartition function. + void compute_distribution (Bool WriteAllInfo, Bool WriteHisto, Bool + Verbose, int param); + + // set the Threshold array for a confidence interval + // given by Epsilon + void find_threshold (Bool WriteAllInfo,Bool WriteHisto, + Bool Verbose, int param); + + // set the Threshold array for a confidence interval + // given by Epsilon + void histo_threshold (Bool WriteAllInfo=False, Bool + Verbose=False); + + void get_threshold(int NbrEvent, int CurrentScale, float & SeuilMin, float & SeuilMax); + // Return the threshold min and max corresponding the number of events + // NbrEvent + + void event_set_support( fltarray *& Data_In, + int CurrentScale, + int FirstDectectScale, + int NEGFirstDectectScale, + int MinEvent, + Bool OnlyPosVal, + type_border Border, + intarray *& Nb_Event, + intarray *& Support, + Bool WriteAllInfo); + + + // give the value of a pixel depending to the border + int get_pix(int i,int j,int k, intarray & data, type_border Border); + + // compute nb_event with no approximation + void event_scale_3D (intarray & data_event, intarray *& nb_event, + int Current_Scale, type_border Border, Bool WriteAllInfo) ; + + // compute nb_event with no approximation + void event_scale_3D_Approx ( intarray & last_nb_event, + intarray & nb_event, + int Current_Scale, + type_border Border, + Bool WriteAllInfo) ; + + void get_sigma_band( fltarray & Nb_Event_In_Boxes,fltarray & tab_sigma, + float N_Sigma,int Nb_Scale); +}; + + +#endif diff --git a/src/mr2d1d_stat.cc b/src/mr2d1d_stat.cc new file mode 100644 index 0000000..ed61825 --- /dev/null +++ b/src/mr2d1d_stat.cc @@ -0,0 +1,406 @@ +/****************************************************************************** +** Copyright (C) 1999 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 04/09/99 +** +** File: mr3d_trans.cc +** +******************************************************************************* +** +** DESCRIPTION multiresolution transform of cube +** ----------- +** +** Usage: mr3d_trans options cube output +** +******************************************************************************/ + + + +#include "IM_Obj.h" +#include "IM_IO.h" +#include "SB_Filter.h" +#include "MR1D_Obj.h" +#include "IM3D_IO.h" +#include "MR3D_Obj.h" +#include "DefFunc.h" + +/****************************************************************************/ + +char Name_Cube_In[256]; +char Name_Out[256]; +char NameSurv[256]; + +int Nbr_Plan=9; +type_trans_3d Transform=TO3_MALLAT; + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char *const*argv, char *opts); + +sb_type_norm Norm = NORM_L2; +type_sb_filter SB_Filter = F_HAAR; // F_MALLAT_7_9; +type_lift LiftingTrans = DEF_LIFT; +Bool Verbose=False; +Bool InfoBand=False; +char FileSimu[180]; +Bool WithFileSimu=False; +int NbrCoeffMax = 10; +Bool GetMax = False; +float N_Sigma=3; // number of sigma (for the noise) +Bool Normalize=False; // normalize data in +Bool Survival=False; + + +type_border Bord = I_MIRROR; + +/*********************************************************************/ + +/*static int max_scale_number (int Nc) +{ + int Nmin = Nc; + int ScaleMax; + + ScaleMax=iround(log((float)Nmin/(float)MAX_SIZE_LAST_SCALE) / log(2.)+ 1.); + return (ScaleMax); +}*/ + +/*********************************************************************/ + +static void usage(char *argv[]) +{ + int i; + + fprintf(OUTMAN, "Usage: %s options image output\n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + +// fprintf(OUTMAN, "\n"); +// fprintf(OUTMAN, " [-t type_of_multiresolution_transform]\n"); +// for (i = 0; i < NBR_TRANS_3D; i++) +// fprintf(OUTMAN, " %d: %s \n",i+1, +// StringTransf3D((type_trans_3d)i)); +// fprintf(OUTMAN, " default is %s\n", StringTransf3D((type_trans_3d) Transform)); + + manline(); + fprintf(OUTMAN, " [-n number_of_scales]\n"); + fprintf(OUTMAN, " number of scales used in the multiresolution transform\n"); + manline(); + + sb_usage(SB_Filter); + manline(); + + fprintf(OUTMAN, " [-N]\n"); + fprintf(OUTMAN, " Normalize the data. Default is no. \n"); + manline(); + + verbose_usage(); + manline(); + vm_usage(); + manline(); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ + +static void transinit(int argc, char *argv[]) +{ + int c; + Bool OptL = False, Optf = False; +#ifdef LARGE_BUFF + int VMSSize=-1; + Bool OptZ = False; + char VMSName[1024] = ""; +#endif + + /* get options */ + while ((c = GetOpt(argc,argv,"VNs:m:n:S:Ll:T:vzZ:")) != -1) + { + switch (c) + { + case 'V': Survival = (Survival == True) ? False: True; break; + case 'N': Normalize = (Normalize == True) ? False: True; break; + case 'v': Verbose = True; break; + case 'L': Norm = NORM_L1; OptL = True; break; + case 'T': + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "Error: bad filter: %s\n", OptArg); + exit(-1); + + } + Optf = True; + SB_Filter = get_filter_bank(OptArg); + // SB_Filter = (type_sb_filter) (c); + break; + case 'n': + /* -n */ + if (sscanf(OptArg,"%d",&Nbr_Plan) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((Nbr_Plan <= 1) || (Nbr_Plan > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + fprintf(OUTMAN, "1 < Nbr Scales <= %d\n", MAX_SCALE_1D); + exit(-1); + } + break; + case 'S': + /* -n */ + if (sscanf(OptArg,"%s",FileSimu) != 1) + { + fprintf(OUTMAN, "bad file name: %s\n", OptArg); + exit(-1); + } + WithFileSimu = True; + break; +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': + usage(argv); + } + } + + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Cube_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Error: too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Dat; + int N,b; + /* Get command line arguments, open input file(s) if necessary */ + fitsstruct Header; + char Cmd[512]; + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + + lm_check(LIC_MR3); + transinit(argc, argv); + + if (Verbose == True) + { + cout << "Filename in = " << Name_Cube_In << endl; + cout << "Filename out = " << Name_Out << endl; + // cout << "Transform = " << StringTransf3D((type_trans_3d) Transform) << endl; + } + io_3d_read_data(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 (Normalize == True) + { + double Mean = Dat.mean(); + double Sigma = Dat.sigma(); + for (int i=0;i= 2) + { + for (b=0; b < NbrBand2D; 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) = 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++) TabBand[b](j,i,z++) = WT1D(b1,p); + } + } + } +} + +/****************************************************************************/ + + +/****************************************************************************/ + +void MR2D1D::recons (fltarray &Data) +{ + if ((Data.nx() != Nx) || (Data.ny() != Ny) || (Data.nz() != Nz)) Data.resize(Nx, Ny, Nz); + int i,j,b,z; + Nx = Data.nx(); + Ny = Data.ny(); + Nz = Data.nz(); + Ifloat Frame(Ny, Nx); + fltarray Vect(Nz); + + // 1D wt + if (NbrBand1D >= 2) + { + for (b=0; b < NbrBand2D; 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) = TabBand[b](j,i,z); + z = 0; + for (int b1=0; b1 < NbrBand1D; b1++) + for (int p=0; p < WT1D.size_scale_np (b1); p++) WT1D(b1,p) = TabBand[b](j,i,z++); + Vect.init(); + WT1D.recons(Vect); + for (z=0; z < Nz; z++) TabBand[b](j,i,z) = Vect(z); + } + } + + // 2D wt + for (z=0; z < Nz; z++) + { + for (b=0; b < NbrBand2D; b++) + { + for (i=0; i < WT2D.size_band_nl(b); i++) + for (j=0; j < WT2D.size_band_nc(b); j++) WT2D(b,i,j) = TabBand[b](j,i,z); + } + WT2D.recons(Frame); + for (i=0; i < Ny; i++) + for (j=0; j < Nx; j++) Data(j,i,z) = Frame(i,j); + } +} + + +/****************************************************************************/ + +char Name_Cube_In[256]; +char Name_Out[256]; + +int NbrScale2d = 5; +int Nbr_Plan=4; +type_transform Transform=TO_MALLAT; + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char *const*argv, char *opts); + +// sb_type_norm Norm = NORM_L2; +// type_sb_filter SB_Filter = F_HAAR; // F_MALLAT_7_9; +// type_lift LiftingTrans = DEF_LIFT; +Bool Verbose=False; +Bool GetMax = False; +float N_Sigma=3; // number of sigma (for the noise) +Bool Normalize=False; // normalize data in + +type_border Bord = I_MIRROR; +Bool Reverse = False; + +/*********************************************************************/ + +/*static int max_scale_number (int Nc) +{ + int Nmin = Nc; + int ScaleMax; + + ScaleMax=iround(log((float)Nmin/(float)MAX_SIZE_LAST_SCALE) / log(2.)+ 1.); + return (ScaleMax); +}*/ + +/*********************************************************************/ + +static void usage(char *argv[]) +{ + fprintf(OUTMAN, "Usage: %s options cube output\n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + + all_transform_usage(Transform); + manline(); + fprintf(OUTMAN, " [-n number_of_scales_2D]\n"); + fprintf(OUTMAN, " number of scales used in the 2D wavelet transform\n"); + manline(); + fprintf(OUTMAN, " [-N number_of_scales_2D]\n"); + fprintf(OUTMAN, " number of scales used in the 1D wavelet transform\n"); + manline(); + fprintf(OUTMAN, " [-M]\n"); + fprintf(OUTMAN, " Normalize the data. Default is no. \n"); + manline(); + manline(); + fprintf(OUTMAN, " [-r]\n"); + fprintf(OUTMAN, " Reconstruction. Default is no. \n"); + manline(); + verbose_usage(); + manline(); + vm_usage(); + manline(); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ + +static void transinit(int argc, char *argv[]) +{ + int c; +#ifdef LARGE_BUFF + int VMSSize=-1; + Bool OptZ = False; + char VMSName[1024] = ""; +#endif + + /* get options */ + while ((c = GetOpt(argc,argv,"rN:t:n:MvzZ:")) != -1) + { + switch (c) + { + case 't': + /* -d type of transform */ + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "bad type of multiresolution transform: %s\n", OptArg); + exit(-1); + + } + if ((c > 0) && (c <= NBR_TRANSFORM+1)) + Transform = (type_transform) (c-1); + else + { + fprintf(OUTMAN, "bad type of transform: %s\n", OptArg); + exit(-1); + } + break; + case 'r': Reverse = True; break; + case 'M': Normalize = (Normalize == True) ? False: True; break; + case 'v': Verbose = True; break; + case 'n': + /* -n */ + if (sscanf(OptArg,"%d",&NbrScale2d) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((NbrScale2d <= 1) || (NbrScale2d > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + break; + case 'N': + /* -n */ + if (sscanf(OptArg,"%d",&Nbr_Plan) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((Nbr_Plan <= 0) || (Nbr_Plan > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + break; +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': + usage(argv); + } + } + + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Cube_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Error: too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Dat; + /* Get command line arguments, open input file(s) if necessary */ + fitsstruct Header; + char Cmd[512]; + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + + lm_check(LIC_MR3); + transinit(argc, argv); + + + + if (Reverse == False) + { + if (Verbose == True) + { + cout << "Filename in = " << Name_Cube_In << endl; + cout << "Filename out = " << Name_Out << endl; + cout << "Transform = " << StringTransform((type_transform) Transform) << endl; + cout << "NbrScale2d = " << NbrScale2d<< endl; + cout << "NbrScale1d = " << Nbr_Plan<< endl; + } + + io_3d_read_data(Name_Cube_In, Dat, &Header); + + 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 (Normalize == True) + { + double Mean = Dat.mean(); + double Sigma = Dat.sigma(); + for (int i=0;i NBR_OK_TRANS3D_METHOD)) + { + fprintf(OUTMAN, "Error: bad transform type : %s\n", OptArg); + exit(-1); + } + Transform = (type_trans_3d) (TabMethod[c-1]); + break; + case 'v': Verbose = True; break; + case 'C': MAD = True; break; + case 'g': + /* -g */ + if (sscanf(OptArg,"%f",&Noise_Ima) != 1) + { + fprintf(OUTMAN, "Error: bad sigma noise: %s\n", OptArg); + exit(-1); + } + break; + case 's': + /* -s */ + if (sscanf(OptArg,"%f",&N_Sigma) != 1) + { + fprintf(OUTMAN, "Error: bad N_Sigma: %s\n", OptArg); + exit(-1); + } + if (N_Sigma <= 0.) N_Sigma = 3; + break; + case 'i': + /* -i < Number of iterations> */ + if (sscanf(OptArg,"%d",&Max_Iter) != 1) + { + fprintf(OUTMAN, "Error: bad Max_Iter: %s\n", OptArg); + exit(-1); + } + if (Max_Iter <= 0) Max_Iter = 0; + break; + case 'T': + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "Error: bad filter: %s\n", OptArg); + exit(-1); + + } + SB_Filter = get_filter_bank(OptArg); + // SB_Filter = (type_sb_filter) (c); + break; + case 'n': + /* -n */ + if (sscanf(OptArg,"%d",&Nbr_Plan) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((Nbr_Plan <= 1) || (Nbr_Plan > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + fprintf(OUTMAN, "1 < Nbr Scales <= %d\n", MAX_SCALE_1D); + exit(-1); + } + break; +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': + usage(argv); + } + } + + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Cube_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Cube_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Error: too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + +/*********************************************************************/ + +// static void mr_regul_cube_rec(MR_3D &MR_Step, fltarray &Result, int MaxIter, float LevelScale, float *TabLevel) +// { +// int Nx = Result.nx(); +// int Ny = Result.ny(); +// int Nz = Result.nz(); +// int Nw = Nx*Ny*Nz; +// int b,i,j,k,w,Iter=0; +// MR_3D MR_Iter (Nx, Ny, Nz, MR_Step.Type_Transform, MR_Step.nbr_scale()); +// fltarray Lap(Nx, Ny, Nz); +// MR_Iter.SB_Filter = MR_Step.SB_Filter; +// MR_Iter.Norm = MR_Step.Norm; +// MR_Iter.LiftingTrans = MR_Step.LiftingTrans; +// // MR_Iter.NormHaar = MR_Step.NormHaar; +// +// unsigned char *TabSup = new unsigned char [Nw]; +// float F; +// float F_old=0; +// float Delta_F=10;; +// float Level = LevelScale; +// +// // type_border Bord=I_CONT; +// float *TabSave = new float [Nw]; +// float Delta_w,LevelQ; +// +// w =0; +// for (b=0; b < MR_Step.nbr_band()-1; b++) +// for (i=0; i < MR_Step.size_band_nx(b)-1; i++) +// for (j=0; j < MR_Step.size_band_ny(b)-1; j++) +// for (k=0; k < MR_Step.size_band_nz(b)-1; k++) +// { +// // Save initial value of dequantized coefficients for range constraint +// TabSave[w]=MR_Step(b,i,j,k); +// +// // compute the multiresolution support +// TabSup[w] = (ABS(MR_Step(b,i,j,k)) > FLOAT_EPSILON) ? 1 : 0; +// w++; +// } +// +// while (Iter .001 )&&(IterLevel) && (TabSup[w] == 0)) +// MR_Step(b,i,j,k) = (MR_Step(b,i,j,k)>0.0)? Level:-Level; +// +// // range constraint on quantized-restored coefficients +// LevelQ = Level / 2.; +// Delta_w = MR_Step(b,i,j,k) - TabSave[w]; +// if ((ABS(Delta_w) > LevelQ/2.)&&(TabSup[w]==1)) +// MR_Step(b,i,j,k) = (Delta_w>0.0)? TabSave[w]+LevelQ:TabSave[w]-LevelQ; +// w++; +// } +// Delta_F = (F-F_old)/ w; +// +// // cerr << "Iter = " << Iter+1<< " F = " << F << " Delta_F =" << Delta_F*Delta_F << endl; +// F_old=F; +// Iter++; +// } +// MR_Step.recons(Result); +// // positivity constraint on reconstructed image +// // threshold(Result); +// +// delete TabSup; +// delete TabSave; +// } +// +/*********************************************************************/ + +int threshold_3d_band(MR_3D & MR_Data, int b, float Level) +{ + int i,j,k; + int Nxb = MR_Data.size_band_nx(b); + int Nyb = MR_Data.size_band_ny(b); + int Nzb = MR_Data.size_band_nz(b); + int Cpt = 0; + + for (i=0; i < Nxb; i++) + for (j=0; j < Nyb; j++) + for (k=0; k < Nzb; k++) + { + if (ABS(MR_Data(b,i,j,k)) < Level) MR_Data(b,i,j,k) = 0.; + else Cpt ++; + } + return Cpt; +} + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fitsstruct Header; + fltarray Dat; + int b,k; + char Cmd[512]; + /* Get command line arguments, open input file(s) if necessary */ + + extern softinfo Soft; + + Soft.mr3(); + Cmd[0] = '\0'; + for (k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + + lm_check(LIC_MR3); + transinit(argc, argv); + + io_3d_read_data(Name_Cube_In, Dat, &Header); + Header.origin = Cmd; + + int Nx = Dat.nx(); + int Ny = Dat.ny(); + int Nz = Dat.nz(); + if (Verbose == True) + { + cout << "Filename in = " << Name_Cube_In << endl; + cout << "Filename out = " << Name_Cube_Out << endl; + cout << "Nx = " << Dat.nx() << " Ny = " << Dat.ny() << " Nz = " << Dat.nz() << endl; + } + + MR_3D MR_Data; + MR_Data.Verbose=Verbose; + FilterAnaSynt FAS; + FilterAnaSynt *PtrFAS = NULL; + if (Transform == TO3_MALLAT) + { + FAS.Verbose = Verbose; + FAS.alloc(SB_Filter); + PtrFAS = &FAS; + } + MR_Data.alloc (Nx, Ny, Nz, Transform, Nbr_Plan, PtrFAS, Norm); + + if (Verbose == True) + { + cout << "Transform = " << StringTransf3D((type_trans_3d) Transform) << endl; + if (Transform == TO3_MALLAT) + { + cout << "Filter = " << StringSBFilter(MR_Data.SBFilter) << endl; + if (MR_Data.TypeNorm == NORM_L1) cout << "Norm L1 " << endl; + else cout << "Norm L2 " << endl; + } + cout << "Number of scales = " << MR_Data.nbr_scale() << endl; + if ((MAD == False) && (Noise_Ima > FLOAT_EPSILON)) + cout << "Sigma Noise = " << Noise_Ima << endl; + } + + if (Verbose == True) cout << "3D transform ... " << endl; + MR_Data.transform (Dat); + if (Verbose == True) cout << "Thresholding ... " << endl; + + if (Noise_Ima < FLOAT_EPSILON) + { + if (Verbose == True) cout << "Automatic noise estimation ... " << endl; + fltarray DatBand; + MR_Data.get_band(0, DatBand); + Noise_Ima = DatBand.sigma_clip(); + if (Verbose == True) cout << "Sigma noise = " << Noise_Ima << endl; + } + int Ndetect = 0; + for (b=0; b < MR_Data.nbr_band()-1; b++) + { + float Level = Noise_Ima*N_Sigma; + if (MAD == True) + { + int i,j,k; + int N1 = MR_Data.size_band_nx(b); + int N2 = MR_Data.size_band_ny(b); + int N3 = MR_Data.size_band_nz(b); + fltarray Buff(N1*N2*N3); + int Ind = 0; + for (i=0;i 0) && (c <= NBR_LIFT)) LiftingTrans = (type_lift) (c); + else + { + fprintf(OUTMAN, "Error: bad lifting method: %s\n", OptArg); + exit(-1); + } + break; + case 'T': + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "Error: bad filter: %s\n", OptArg); + exit(-1); + + } + Optf = True; + SB_Filter = get_filter_bank(OptArg); + // SB_Filter = (type_sb_filter) (c); + break; + case 't': + /* -d type of transform */ + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "bad type of multiresolution transform: %s\n", OptArg); + exit(-1); + + } + if ((c > 0) && (c <= NBR_TRANS_3D)) + Transform = (type_trans_3d) (c-1); + else + { + fprintf(OUTMAN, "bad type of transform: %s\n", OptArg); + exit(-1); + } + break; + case 'n': + /* -n */ + if (sscanf(OptArg,"%d",&Nbr_Plan) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((Nbr_Plan <= 1) || (Nbr_Plan > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + fprintf(OUTMAN, "1 < Nbr Scales <= %d\n", MAX_SCALE_1D); + exit(-1); + } + break; + case 'S': + /* -n */ + if (sscanf(OptArg,"%s",FileSimu) != 1) + { + fprintf(OUTMAN, "bad file name: %s\n", OptArg); + exit(-1); + } + WithFileSimu = True; + break; +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': + usage(argv); + } + } + + + if ((Transform != TO3_MALLAT) && ((OptL == True) || (Optf == True))) + { + fprintf(OUTMAN, "Error: option -f and -L are only valid with transforms 14 and 16 ... \n"); + exit(0); + } + if ((Transform != TO3_LIFTING) && ( LiftingTrans != DEF_LIFT)) + { + fprintf(OUTMAN, "Error: option -f and -L are only valid with transforms 14 and 16 ... \n"); + exit(0); + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Cube_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Error: too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Dat; + int N,b; + /* Get command line arguments, open input file(s) if necessary */ + fitsstruct Header; + char Cmd[512]; + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + + lm_check(LIC_MR3); + transinit(argc, argv); + + if (Verbose == True) + { + cout << "Filename in = " << Name_Cube_In << endl; + cout << "Filename out = " << Name_Out << endl; + cout << "Transform = " << StringTransf3D((type_trans_3d) Transform) << endl; + if (Transform == TO3_MALLAT) + { + cout << "Filter = " << StringSBFilter(SB_Filter) << endl; + if (Norm == NORM_L1) cout << "Norm L1 " << endl; + else cout << "Norm L2 " << endl; + } + if (Transform == TO3_LIFTING) + cout << "Lifting method = " << StringLSTransform(LiftingTrans) << endl; + } + io_3d_read_data(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 (Normalize == True) + { + double Mean = Dat.mean(); + double Sigma = Dat.sigma(); + for (int i=0;i 0) && (c <= NBR_LIFT)) LiftingTrans = (type_lift) (c); + else + { + fprintf(OUTMAN, "Error: bad lifting method: %s\n", OptArg); + exit(-1); + } + break; + case 'T': + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "Error: bad filter: %s\n", OptArg); + exit(-1); + + } + Optf = True; + SB_Filter = get_filter_bank(OptArg); + // SB_Filter = (type_sb_filter) (c); + break; + case 't': + /* -d type of transform */ + if (sscanf(OptArg,"%d",&c ) != 1) + { + fprintf(OUTMAN, "bad type of multiresolution transform: %s\n", OptArg); + exit(-1); + + } + if ((c > 0) && (c <= NBR_TRANS_3D)) + Transform = (type_trans_3d) (c-1); + else + { + fprintf(OUTMAN, "bad type of transform: %s\n", OptArg); + exit(-1); + } + break; + case 'n': + /* -n */ + if (sscanf(OptArg,"%d",&Nbr_Plan) != 1) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + exit(-1); + } + if ((Nbr_Plan <= 1) || (Nbr_Plan > MAX_SCALE_1D)) + { + fprintf(OUTMAN, "bad number of scales: %s\n", OptArg); + fprintf(OUTMAN, "1 < Nbr Scales <= %d\n", MAX_SCALE_1D); + exit(-1); + } + break; + case 'S': + /* -n */ + if (sscanf(OptArg,"%s",FileSimu) != 1) + { + fprintf(OUTMAN, "bad file name: %s\n", OptArg); + exit(-1); + } + WithFileSimu = True; + break; +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': + usage(argv); + } + } + + + if ((Transform != TO3_MALLAT) && ((OptL == True) || (Optf == True))) + { + fprintf(OUTMAN, "Error: option -f and -L are only valid with transforms 14 and 16 ... \n"); + exit(0); + } + if ((Transform != TO3_LIFTING) && ( LiftingTrans != DEF_LIFT)) + { + fprintf(OUTMAN, "Error: option -f and -L are only valid with transforms 14 and 16 ... \n"); + exit(0); + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Cube_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Cube_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Error: too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Dat; + int b; + /* Get command line arguments, open input file(s) if necessary */ + lm_check(LIC_MR3); + transinit(argc, argv); + + io_3d_read_data(Name_Cube_In, Dat); + + int Nx = Dat.nx(); + int Ny = Dat.ny(); + int Nz = Dat.nz(); + if (Verbose == True) + { + cout << "Filename in = " << Name_Cube_In << endl; + cout << "Filename out = " << Name_Cube_Out << endl; + cout << "Nx = " << Dat.nx() << " Ny = " << Dat.ny() << " Nz = " << Dat.nz() << endl; + } + + + MR_3D MR_Data; + MR_Data.Verbose=Verbose; + FilterAnaSynt FAS; + FilterAnaSynt *PtrFAS = NULL; + if (Transform == TO3_MALLAT) + { + FAS.Verbose = Verbose; + FAS.alloc(SB_Filter); + PtrFAS = &FAS; + } + MR_Data.alloc (Nx, Ny, Nz, Transform, Nbr_Plan, PtrFAS, Norm); + + if (Transform == TO3_LIFTING) MR_Data.LiftingTrans = LiftingTrans; + + if (Verbose == True) + { + cout << "Transform = " << StringTransf3D((type_trans_3d) Transform) << endl; + if (Transform == TO3_MALLAT) + { + cout << "Filter = " << StringSBFilter(MR_Data.SBFilter) << endl; + if (MR_Data.TypeNorm == NORM_L1) cout << "Norm L1 " << endl; + else cout << "Norm L2 " << endl; + } + + if (Transform == TO3_LIFTING) + cout << "Lifting method = " << StringLSTransform(MR_Data.LiftingTrans) << endl; + cout << "Number of scales = " << MR_Data.nbr_scale() << endl; + } + if (Verbose == True) MR_Data.info_pos_band(); + + if (Verbose == True) cout << "3D transform ... " << endl; + MR_Data.transform (Dat); + + if (Verbose == True) cout << "Write ... " << endl; + // io_3d_write_data(Name_Cube_Out, MR_Data.cube()); + MR_Data.write(Name_Cube_Out); + + if (InfoBand == True) + for (b=0; b < MR_Data.nbr_band(); b++) + MR_Data.info_band(b); + + exit(0); +} + From a5dcdfd92285c973fe0a4e50a0c624ddfbc1c56d Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 17:52:17 +0100 Subject: [PATCH 2/6] updated travis for gcc-7 --- .travis.yml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index e302964..917fc77 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,4 +1,6 @@ -language: C++ +language: cpp +compiler: gcc +dist: trusty # GitHub branch branches: @@ -11,6 +13,7 @@ sudo: false before_install: - sudo apt-get -qq update - sudo apt-get install -y libcfitsio3-dev + - sudo apt-get install -qq g++-7 # set up installation install: From e6b677ab49528785c6d06c63468f6c2af09fc03c Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 17:57:21 +0100 Subject: [PATCH 3/6] updated travis --- .travis.yml | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 917fc77..c4217fd 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,15 @@ +# System Set-Up +dist: trusty language: cpp compiler: gcc -dist: trusty +addons: + apt: + sources: + - ubuntu-toolchain-r-test + packages: + - gcc-7 + - g++-7 + - cmake # GitHub branch branches: @@ -13,7 +22,6 @@ sudo: false before_install: - sudo apt-get -qq update - sudo apt-get install -y libcfitsio3-dev - - sudo apt-get install -qq g++-7 # set up installation install: From cd48e729269a4bf7feba09f02bbf3912aa83ad5c Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 18:07:10 +0100 Subject: [PATCH 4/6] updated travis --- .travis.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.travis.yml b/.travis.yml index c4217fd..1207205 100644 --- a/.travis.yml +++ b/.travis.yml @@ -22,6 +22,8 @@ sudo: false before_install: - sudo apt-get -qq update - sudo apt-get install -y libcfitsio3-dev + - export CXX=g++-7 + - export CC=gcc-7 # set up installation install: From be3282ffc7a1c364ef7db86ab7da6eec52bc362d Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 18:38:56 +0100 Subject: [PATCH 5/6] updated travis and test cmake change --- .travis.yml | 11 ----------- CMakeLists.txt | 3 ++- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/.travis.yml b/.travis.yml index 1207205..94ec8f6 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,15 +1,6 @@ # System Set-Up -dist: trusty language: cpp compiler: gcc -addons: - apt: - sources: - - ubuntu-toolchain-r-test - packages: - - gcc-7 - - g++-7 - - cmake # GitHub branch branches: @@ -22,8 +13,6 @@ sudo: false before_install: - sudo apt-get -qq update - sudo apt-get install -y libcfitsio3-dev - - export CXX=g++-7 - - export CC=gcc-7 # set up installation install: diff --git a/CMakeLists.txt b/CMakeLists.txt index 6297cff..4ca4f28 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,7 +137,8 @@ project(sparse2d) add_executable(${program} ${PROJECT_SOURCE_DIR}/src/${program}.cc) target_link_libraries(${program} mga2d sparse2d sparse1d tools) if(SPARSE3D) - target_link_libraries(${program} sparse3d) + # target_link_libraries(${program} sparse3d) + target_link_libraries(${program} mga2d sparse3d sparse2d sparse1d tools) endif(SPARSE3D) endforeach(program) From b6d4cae22e3b83218d7981689239a1c0dc25260e Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 26 Feb 2018 18:42:35 +0100 Subject: [PATCH 6/6] minor cmake change --- CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4ca4f28..92a4482 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -137,7 +137,6 @@ project(sparse2d) add_executable(${program} ${PROJECT_SOURCE_DIR}/src/${program}.cc) target_link_libraries(${program} mga2d sparse2d sparse1d tools) if(SPARSE3D) - # target_link_libraries(${program} sparse3d) target_link_libraries(${program} mga2d sparse3d sparse2d sparse1d tools) endif(SPARSE3D) endforeach(program) @@ -177,7 +176,7 @@ project(sparse2d) add_executable(${program} ${PROJECT_SOURCE_DIR}/tests/${program}.cpp) target_link_libraries(${program} mga2d sparse2d sparse1d tools) if(SPARSE3D) - target_link_libraries(${program} sparse3d) + target_link_libraries(${program} mga2d sparse3d sparse2d sparse1d tools) endif(SPARSE3D) add_test(${program} ${program}) endforeach(program)