Xmipp  v3.23.11-Nereus
Public Member Functions | Public Attributes | List of all members
ProgResLocalFilter Class Reference

#include <resolution_localfilter.h>

Inheritance diagram for ProgResLocalFilter:
Inheritance graph
[legend]
Collaboration diagram for ProgResLocalFilter:
Collaboration graph
[legend]

Public Member Functions

void defineParams ()
 
void readParams ()
 
void produceSideInfo ()
 
void amplitudeMonogenicSignal3D (MultidimArray< std::complex< double > > &myfftV, double freq, double freqH, double freqL, MultidimArray< double > &amplitude, int count, FileName fnDebug)
 
void postProcessingLocalResolutions (MultidimArray< double > &resolutionVol, std::vector< double > &list, MultidimArray< double > &resolutionChimera, double &cut_value, MultidimArray< int > &pMask, double &resolutionThreshold)
 
void resolution2eval (int &count_res, double step, double &resolution, double &last_resolution, double &freq, double &freqL, int &last_fourier_idx, bool &continueIter, bool &breakIter, bool &doNextIteration)
 
void run ()
 
- Public Member Functions inherited from XmippProgram
const char * getParam (const char *param, int arg=0)
 
const char * getParam (const char *param, const char *subparam, int arg=0)
 
int getIntParam (const char *param, int arg=0)
 
int getIntParam (const char *param, const char *subparam, int arg=0)
 
double getDoubleParam (const char *param, int arg=0)
 
double getDoubleParam (const char *param, const char *subparam, int arg=0)
 
float getFloatParam (const char *param, int arg=0)
 
float getFloatParam (const char *param, const char *subparam, int arg=0)
 
void getListParam (const char *param, StringVector &list)
 
int getCountParam (const char *param)
 
bool checkParam (const char *param)
 
bool existsParam (const char *param)
 
void addParamsLine (const String &line)
 
void addParamsLine (const char *line)
 
ParamDefgetParamDef (const char *param) const
 
virtual void quit (int exit_code=0) const
 
virtual int tryRun ()
 
void initProgress (size_t total, size_t stepBin=60)
 
void setProgress (size_t value=0)
 
void endProgress ()
 
void processDefaultComment (const char *param, const char *left)
 
void setDefaultComment (const char *param, const char *comment)
 
virtual void initComments ()
 
void setProgramName (const char *name)
 
void addUsageLine (const char *line, bool verbatim=false)
 
void clearUsage ()
 
void addExampleLine (const char *example, bool verbatim=true)
 
void addSeeAlsoLine (const char *seeAlso)
 
void addKeywords (const char *keywords)
 
const char * name () const
 
virtual void usage (int verb=0) const
 
virtual void usage (const String &param, int verb=2)
 
int version () const
 
virtual void show () const
 
virtual void read (int argc, const char **argv, bool reportErrors=true)
 
virtual void read (int argc, char **argv, bool reportErrors=true)
 
void read (const String &argumentsLine)
 
 XmippProgram ()
 
 XmippProgram (int argc, const char **argv)
 
virtual ~XmippProgram ()
 

Public Attributes

FileName fnOut
 
FileName fnVol
 
FileName fnRes
 
FileName fnMask
 
FileName fnchim
 
FileName fnSpatial
 
FileName fnMeanVol
 
FileName fnMaskOut
 
FileName fnMd
 
FileName fnFilt
 
double sampling
 
double minRes
 
double maxRes
 
double R
 
int NVoxelsOriginalMask
 
int Nvoxels
 
int nthrs
 
double freq_step
 
double trimBound
 
double significance
 
bool noiseOnlyInHalves
 
bool automaticMode
 
MultidimArray< double > iu
 
MultidimArray< std::complex< double > > fftV
 
FourierTransformer transformer_inv
 
Image< double > Vfiltered
 
Image< double > VresolutionFiltered
 
Image< double > resVol
 
Matrix1D< double > freq_fourier
 
double sigma
 
- Public Attributes inherited from XmippProgram
bool doRun
 
bool runWithoutArgs
 
int verbose
 Verbosity level. More...
 
int debug
 

Additional Inherited Members

- Protected Member Functions inherited from XmippProgram
void defineCommons ()
 
- Protected Attributes inherited from XmippProgram
int errorCode
 
ProgramDefprogDef
 Program definition and arguments parser. More...
 
std::map< String, CommentListdefaultComments
 
int argc
 Original command line arguments. More...
 
const char ** argv
 

Detailed Description

SSNR parameters.

Definition at line 48 of file resolution_localfilter.h.

Member Function Documentation

◆ amplitudeMonogenicSignal3D()

void ProgResLocalFilter::amplitudeMonogenicSignal3D ( MultidimArray< std::complex< double > > &  myfftV,
double  freq,
double  freqH,
double  freqL,
MultidimArray< double > &  amplitude,
int  count,
FileName  fnDebug 
)

◆ defineParams()

void ProgResLocalFilter::defineParams ( )
virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 44 of file resolution_localfilter.cpp.

45 {
46  addUsageLine("This function performs a local filter of a map/tomogram based on the local resolution values");
47  addParamsLine(" --vol <vol_file=\"\"> : Volume");
48  addParamsLine(" --resvol <vol_file=\"\"> : Resolution map");
49  addParamsLine(" -o <output=\"MGresolution.vol\"> : Local resolution volume (in Angstroms)");
50  addParamsLine(" --filteredMap <output=\"filteredMap.vol\"> : Local resolution volume filtered (in Angstroms)");
51  addParamsLine(" [--sampling_rate <s=1>] : Sampling rate (A/px)");
52  addParamsLine(" [--step <s=0.25>] : The resolution is computed at a number of frequencies between minimum and");
53  addParamsLine(" [--significance <s=0.95>] : The level of confidence for the hypothesis test.");
54  addParamsLine(" [--threads <s=4>] : Number of threads");
55 }
void addUsageLine(const char *line, bool verbatim=false)
void addParamsLine(const String &line)

◆ postProcessingLocalResolutions()

void ProgResLocalFilter::postProcessingLocalResolutions ( MultidimArray< double > &  resolutionVol,
std::vector< double > &  list,
MultidimArray< double > &  resolutionChimera,
double &  cut_value,
MultidimArray< int > &  pMask,
double &  resolutionThreshold 
)

◆ produceSideInfo()

void ProgResLocalFilter::produceSideInfo ( )

Definition at line 58 of file resolution_localfilter.cpp.

59 {
60  std::cout << "Starting..." << std::endl;
61  std::cout << " " << std::endl;
62  std::cout << "IMPORTANT: If the angular step of the tilt series is higher than 3 degrees" << std::endl;
63  std::cout << " then, the tomogram is not properly for MonoTomo. Despite this is not "<< std::endl;
64  std::cout << " desired, MonoTomo will try to compute the local resolution." << std::endl;
65  std::cout << " " << std::endl;
66 
67  Image<double> V;
68 
69  V.read(fnVol);
70  V().setXmippOrigin();
71 
72  std::cout << "Map read" << std::endl;
73  MultidimArray<double> &inputVol = V();
74 
75  int N_smoothing = 10;
76 
77  int siz_z = ZSIZE(inputVol)*0.5;
78  int siz_y = YSIZE(inputVol)*0.5;
79  int siz_x = XSIZE(inputVol)*0.5;
80 
81  //Smoothing the boundaries of the volume
82  int limit_distance_x = (siz_x-N_smoothing);
83  int limit_distance_y = (siz_y-N_smoothing);
84  int limit_distance_z = (siz_z-N_smoothing);
85 
86  double uz, uy, ux, uz2, u2, uz2y2;
87  long n=0;
88  n=0;
89  for(int k=0; k<ZSIZE(inputVol); ++k)
90  {
91  uz = (k - siz_z);
92  for(int i=0; i<YSIZE(inputVol); ++i)
93  {
94  uy = (i - siz_y);
95  for(int j=0; j<XSIZE(inputVol); ++j)
96  {
97  ux = (j - siz_x);
98 
99  if (abs(ux)>=limit_distance_x)
100  {
101  DIRECT_MULTIDIM_ELEM(inputVol, n) *= 0.5*(1+cos(PI*(limit_distance_x - abs(ux))/N_smoothing));
102  }
103  if (abs(uy)>=limit_distance_y)
104  {
105  DIRECT_MULTIDIM_ELEM(inputVol, n) *= 0.5*(1+cos(PI*(limit_distance_y - abs(uy))/N_smoothing));
106  }
107  if (abs(uz)>=limit_distance_z)
108  {
109  DIRECT_MULTIDIM_ELEM(inputVol, n) *= 0.5*(1+cos(PI*(limit_distance_z - abs(uz))/N_smoothing));
110  }
111  ++n;
112  }
113  }
114  }
115  std::cout << "No fallo" << std::endl;
116 
117  FourierTransformer transformer;
118  transformer.FourierTransform(inputVol, fftV);
119 
120  MultidimArray<double> inputVol2;
121  size_t Zdim, Ydim, Xdim, Ndim;
122  fftV.getDimensions(Xdim, Ydim, Zdim, Ndim);
123  inputVol2.resizeNoCopy(Ndim, Zdim, Ydim, Xdim);
124 
125 
127  DIRECT_MULTIDIM_ELEM(inputVol2, n) = fabs(DIRECT_MULTIDIM_ELEM(fftV, n));
128 
129 // Image<double> filteredvolume;
130 // filteredvolume = inputVol2;
131 // filteredvolume.write("fourier.vol");
132  iu.resizeNoCopy(Ndim, Zdim, Ydim, Xdim);
133 
134 
135  // Calculate frequency map
136  n=0;
137  for(size_t k=0; k<ZSIZE(fftV); ++k)
138  {
139  FFT_IDX2DIGFREQ(k,ZSIZE(inputVol),uz);
140  uz2=uz*uz;
141 
142  for(size_t i=0; i<YSIZE(fftV); ++i)
143  {
144  FFT_IDX2DIGFREQ(i,YSIZE(inputVol),uy);
145  uz2y2=uz2+uy*uy;
146 
147  for(size_t j=0; j<XSIZE(fftV); ++j)
148  {
149  FFT_IDX2DIGFREQ(j,XSIZE(inputVol),ux);
150  u2=uz2y2+ux*ux;
151  if ((k != 0) || (i != 0) || (j != 0))
152  DIRECT_MULTIDIM_ELEM(iu,n) = 1.0/sqrt(u2);
153  else
154  DIRECT_MULTIDIM_ELEM(iu,n) = 1e38;
155  ++n;
156  }
157  }
158  }
159  V.clear();
160 
161  size_t len;
162  if ( (ZSIZE(fftV) <= XSIZE(fftV)) && (ZSIZE(fftV) <= YSIZE(fftV)) )
163  len = ZSIZE(fftV);
164  if ( (YSIZE(fftV) <= ZSIZE(fftV)) && (YSIZE(fftV) <= XSIZE(fftV)) )
165  len = YSIZE(fftV);
166  else
167  len = XSIZE(fftV);
168 
169  freq_fourier.initZeros(len);
170 
171 
172  resVol.read(fnRes);
173  std::cout << "resolution read" << std::endl;
174  MultidimArray<double> &pResVol = resVol();
175  maxRes = -1;
176  minRes = 1e38;
177  double sumres = 0, sumres2 = 0;
178 
179  long counter = 0;
180 
182  {
183  double resVal = DIRECT_MULTIDIM_ELEM(pResVol, n);
184 
185  sumres += resVal;
186  sumres2 += resVal*resVal;
187 
188  if (resVal>=maxRes)
189  maxRes=resVal;
190  if (resVal<minRes)
191  minRes=resVal;
192  }
193  double mean = sumres/counter;
194  sigma = sumres2/counter-mean*mean;
195 
196 // double u;
197 //
198 // VEC_ELEM(freq_fourier,0) = 1e-38;
199 // for(size_t k=0; k<ZSIZE(fftV); ++k)
200 // {
201 // FFT_IDX2DIGFREQ(k,ZSIZE(pMask), u);
202 // VEC_ELEM(freq_fourier,k) = u;
203 // }
204 }
#define YSIZE(v)
Matrix1D< double > freq_fourier
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
MultidimArray< double > iu
void abs(Image< double > &op)
MultidimArray< std::complex< double > > fftV
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define XSIZE(v)
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define ZSIZE(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
void initZeros()
Definition: matrix1d.h:592
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
#define j
#define len
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
#define PI
Definition: tools.h:43
int * n
void clear()
Definition: xmipp_image.h:144
void getDimensions(size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim) const

◆ readParams()

void ProgResLocalFilter::readParams ( )
virtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 31 of file resolution_localfilter.cpp.

32 {
33  fnVol = getParam("--vol");
34  fnRes = getParam("--resvol");
35  fnOut = getParam("-o");
36  fnFilt = getParam("--filteredMap");
37  sampling = getDoubleParam("--sampling_rate");
38  freq_step = getDoubleParam("--step");
39  significance = getDoubleParam("--significance");
40  nthrs = getIntParam("--threads");
41 }
double getDoubleParam(const char *param, int arg=0)
const char * getParam(const char *param, int arg=0)
int getIntParam(const char *param, int arg=0)

◆ resolution2eval()

void ProgResLocalFilter::resolution2eval ( int &  count_res,
double  step,
double &  resolution,
double &  last_resolution,
double &  freq,
double &  freqL,
int &  last_fourier_idx,
bool &  continueIter,
bool &  breakIter,
bool &  doNextIteration 
)

◆ run()

void ProgResLocalFilter::run ( )
virtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 207 of file resolution_localfilter.cpp.

208 {
209  produceSideInfo();
210 
211 
212  //exit(0);
213  //Determining the frequency range;
214  int lowIdx, highIdx;
215  double lowestfreq, highestfreq, freqL, freqH, freq;
216  lowestfreq = sampling/maxRes;
217  highestfreq = sampling/minRes;
218 
219  DIGFREQ2FFT_IDX(lowestfreq, ZSIZE(resVol()), lowIdx);
220  DIGFREQ2FFT_IDX(highestfreq, ZSIZE(resVol()), highIdx);
221 
223 
224  MultidimArray<double> filteredMap, filtered_aux;
225  filtered_aux.resizeNoCopy(resVol());
226  filteredMap.initZeros(resVol());
227 
228 
229  for (int idx = lowIdx; idx < highIdx; idx++)
230  {
231  FFT_IDX2DIGFREQ((double)idx, ZSIZE(resVol()), freq);
232 
233  freqL = freq - 0.02;
234  if (freqL < 0)
235  freqL = 0.001;
236 
237  freqH = freq + 0.02;
238 
239  if (freqH > 0.5)
240  freqH = 0.5;
241 
242  fftVaux.initZeros(fftV);
243 
244  // Filter the input volume
245  long n=0;
246  double ideltal=PI/(freq-freqH);
247  double ideltah=PI/(freqL-freq);
248 
249 
251  {
252  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
253  double un=1.0/iun;
254  if (un <= freqH && freq <= un)
255  {
256  //double H=0.5*(1+cos((un-w1)*ideltal));
258  DIRECT_MULTIDIM_ELEM(fftVaux, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
259  }
260  if (freqL<=un && un<=freq)
261  {
262  //double H=0.5*(1+cos((un-w1)*ideltal));
264  DIRECT_MULTIDIM_ELEM(fftVaux, n) *= 0.5*(1+cos((un-freq)*ideltah));//H;
265  }
266  }
267 
268  transformer_inv.inverseFourierTransform(fftVaux, filtered_aux);
269  double stdres = sampling/sigma;
270 
272  {
273  double weight, sumweight = 0., digitalfreq;
274  //TODO: make the converstion in the produdesideinfo
275  digitalfreq = sampling/DIRECT_MULTIDIM_ELEM(resVol, n);
276 
277  weight = exp(-(digitalfreq - freq)*(digitalfreq - freq)/stdres);
278  DIRECT_MULTIDIM_ELEM(filteredMap, n) += weight*DIRECT_MULTIDIM_ELEM(filtered_aux, n);
279  sumweight += weight;
280  }
281  }
282 
283  Image<double> saveMap;
284  saveMap = filteredMap;
285  saveMap.write(fnOut);
286 
287 
288 }
FourierTransformer transformer_inv
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
void resizeNoCopy(const MultidimArray< T1 > &v)
MultidimArray< double > iu
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
#define DIGFREQ2FFT_IDX(freq, size, idx)
Definition: xmipp_fft.h:61
MultidimArray< std::complex< double > > fftV
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define ZSIZE(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
void initZeros(const MultidimArray< T1 > &op)
#define PI
Definition: tools.h:43
int * n

Member Data Documentation

◆ automaticMode

bool ProgResLocalFilter::automaticMode

Definition at line 65 of file resolution_localfilter.h.

◆ fftV

MultidimArray< std::complex<double> > ProgResLocalFilter::fftV

Definition at line 94 of file resolution_localfilter.h.

◆ fnchim

FileName ProgResLocalFilter::fnchim

Definition at line 52 of file resolution_localfilter.h.

◆ fnFilt

FileName ProgResLocalFilter::fnFilt

Definition at line 52 of file resolution_localfilter.h.

◆ fnMask

FileName ProgResLocalFilter::fnMask

Definition at line 52 of file resolution_localfilter.h.

◆ fnMaskOut

FileName ProgResLocalFilter::fnMaskOut

Definition at line 52 of file resolution_localfilter.h.

◆ fnMd

FileName ProgResLocalFilter::fnMd

Definition at line 52 of file resolution_localfilter.h.

◆ fnMeanVol

FileName ProgResLocalFilter::fnMeanVol

Definition at line 52 of file resolution_localfilter.h.

◆ fnOut

FileName ProgResLocalFilter::fnOut

Filenames

Definition at line 52 of file resolution_localfilter.h.

◆ fnRes

FileName ProgResLocalFilter::fnRes

Definition at line 52 of file resolution_localfilter.h.

◆ fnSpatial

FileName ProgResLocalFilter::fnSpatial

Definition at line 52 of file resolution_localfilter.h.

◆ fnVol

FileName ProgResLocalFilter::fnVol

Definition at line 52 of file resolution_localfilter.h.

◆ freq_fourier

Matrix1D<double> ProgResLocalFilter::freq_fourier

Definition at line 97 of file resolution_localfilter.h.

◆ freq_step

double ProgResLocalFilter::freq_step

Step in digital frequency

Definition at line 62 of file resolution_localfilter.h.

◆ iu

MultidimArray<double> ProgResLocalFilter::iu

Definition at line 93 of file resolution_localfilter.h.

◆ maxRes

double ProgResLocalFilter::maxRes

Definition at line 56 of file resolution_localfilter.h.

◆ minRes

double ProgResLocalFilter::minRes

Definition at line 56 of file resolution_localfilter.h.

◆ noiseOnlyInHalves

bool ProgResLocalFilter::noiseOnlyInHalves

The search for resolutions is linear or inverse

Definition at line 65 of file resolution_localfilter.h.

◆ nthrs

int ProgResLocalFilter::nthrs

Definition at line 59 of file resolution_localfilter.h.

◆ Nvoxels

int ProgResLocalFilter::Nvoxels

Definition at line 59 of file resolution_localfilter.h.

◆ NVoxelsOriginalMask

int ProgResLocalFilter::NVoxelsOriginalMask

Is the volume previously masked?

Definition at line 59 of file resolution_localfilter.h.

◆ R

double ProgResLocalFilter::R

Definition at line 56 of file resolution_localfilter.h.

◆ resVol

Image<double> ProgResLocalFilter::resVol

Definition at line 96 of file resolution_localfilter.h.

◆ sampling

double ProgResLocalFilter::sampling

sampling rate, minimum resolution, and maximum resolution

Definition at line 56 of file resolution_localfilter.h.

◆ sigma

double ProgResLocalFilter::sigma

Definition at line 98 of file resolution_localfilter.h.

◆ significance

double ProgResLocalFilter::significance

Definition at line 62 of file resolution_localfilter.h.

◆ transformer_inv

FourierTransformer ProgResLocalFilter::transformer_inv

Definition at line 95 of file resolution_localfilter.h.

◆ trimBound

double ProgResLocalFilter::trimBound

Definition at line 62 of file resolution_localfilter.h.

◆ Vfiltered

Image<double> ProgResLocalFilter::Vfiltered

Definition at line 96 of file resolution_localfilter.h.

◆ VresolutionFiltered

Image<double> ProgResLocalFilter::VresolutionFiltered

Definition at line 96 of file resolution_localfilter.h.


The documentation for this class was generated from the following files: