Xmipp  v3.23.11-Nereus
Public Member Functions | Public Attributes | List of all members
ProgResDir Class Reference
Libraries » Reconstruction libraryLibraries » | Basic data library » Resolution

#include <resolution_directional.h>

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

Public Member Functions

void defineParams ()
 
void readParams ()
 
void produceSideInfo ()
 
void amplitudeMonogenicSignal3D_fast (const MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, double wH, MultidimArray< double > &amplitude, int count, int dir, FileName fnDebug)
 
void defineCone (MultidimArray< std::complex< double > > &myfftV, MultidimArray< std::complex< double > > &conefilter, double rot, double tilt)
 
void diagSymMatrix3x3 (Matrix2D< double > A, Matrix1D< double > &eigenvalues, Matrix2D< double > &P)
 
void resolution2eval_ (int &fourier_idx, double min_step, double &resolution, double &last_resolution, int &last_fourier_idx, double &freq, double &freqL, double &freqH, bool &continueIter, bool &breakIter, bool &doNextIteration)
 
double firstMonoResEstimation (MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, MultidimArray< double > &amplitude)
 
void generateGridProjectionMatching (Matrix2D< double > &angles)
 
void removeOutliers (Matrix2D< double > &resolutionMat)
 
void ellipsoidFitting (Matrix2D< double > &resolutionMat, Matrix2D< double > &axis)
 
void radialAzimuthalResolution (Matrix2D< double > &resolutionMat, MultidimArray< int > &pmask, MultidimArray< double > &radial, MultidimArray< double > &azimuthal, MultidimArray< double > &lowestResolution, MultidimArray< double > &highestResolution, MultidimArray< double > &doaResolution_1, MultidimArray< double > &doaResolution_2, double &radial_Thr, double &azimuthal_Thr, MetaDataVec &mdprefDirs)
 
void radialAverageInMask (MultidimArray< int > &mask, MultidimArray< double > &inputVol_1, MultidimArray< double > &inputVol_2, MultidimArray< double > &inputVol_3, MultidimArray< double > &inputVol_4, MultidimArray< double > &inputVol_5, MetaDataVec &md)
 
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)
 
ParamDef * getParamDef (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 fnMask
 
FileName fnradial
 
FileName fnazimuthal
 
FileName fnMonoRes
 
FileName fnRadialAvG
 
FileName fnMDThr
 
FileName fnLowestResolution
 
FileName fnHighestResolution
 
FileName fnprefMin
 
FileName fnZscore
 
FileName fnDoa1
 
FileName fnDoa2
 
double sampling
 
double minRes
 
double maxRes
 
double R
 
double N_directions
 
double Rparticle
 
double res_step
 
int NVoxelsOriginalMask
 
int Nvoxels
 
int Nthr
 
double N_freq
 
double significance
 
bool checkellipsoids
 
bool fastCompute
 
MultidimArray< std::complex< double > > fftVRiesz
 
MultidimArray< std::complex< double > > fftVRiesz_aux
 
MultidimArray< double > iu
 
MultidimArray< double > VRiesz
 
FourierTransformer transformer_inv
 
MultidimArray< std::complex< double > > fftV
 
MultidimArray< std::complex< double > > conefilter
 
Matrix2D< double > angles
 
Matrix2D< double > resolutionMatrix
 
Matrix2D< double > maskMatrix
 
Matrix2D< double > trigProducts
 
Matrix1D< double > freq_fourier
 
Image< int > mask
 
int N_smoothing
 
- 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
 
ProgramDef * progDef
 Program definition and arguments parser. More...
 
std::map< String, CommentList > defaultComments
 
int argc
 Original command line arguments. More...
 
const char ** argv
 

Detailed Description

SSNR parameters.

Definition at line 51 of file resolution_directional.h.

Member Function Documentation

◆ amplitudeMonogenicSignal3D_fast()

void ProgResDir::amplitudeMonogenicSignal3D_fast ( const MultidimArray< std::complex< double > > &  myfftV,
double  w1,
double  w1l,
double  wH,
MultidimArray< double > &  amplitude,
int  count,
int  dir,
FileName  fnDebug 
)

Definition at line 347 of file resolution_directional.cpp.

349 {
350  fftVRiesz.initZeros(myfftV);
351 // MultidimArray<double> coneVol;
352 // coneVol.initZeros(myfftV);
353  fftVRiesz_aux.initZeros(myfftV);
354  std::complex<double> J(0,1);
355 
356  // Filter the input volume and add it to amplitude
357  long n=0;
358  double ideltal=PI/(freq-freqH);
359 
360  for(size_t k=0; k<ZSIZE(myfftV); ++k)
361  {
362  for(size_t i=0; i<YSIZE(myfftV); ++i)
363  {
364  for(size_t j=0; j<XSIZE(myfftV); ++j)
365  {
366  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
367 // double iun = *ptriun;
368  double un=1.0/iun;
369  if (freqH<=un && un<=freq)
370  {
371  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
372 // DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= DIRECT_MULTIDIM_ELEM(conefilter, n);
373  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
374  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
375  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
376  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
377 // DIRECT_MULTIDIM_ELEM(coneVol, n) = DIRECT_MULTIDIM_ELEM(conefilter, n);
378  } else if (un>freq)
379  {
380  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
381 // DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= DIRECT_MULTIDIM_ELEM(conefilter, n);
382  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
383  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
384  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
385 // DIRECT_MULTIDIM_ELEM(coneVol, n) = real(DIRECT_MULTIDIM_ELEM(myfftV, n)*conj(DIRECT_MULTIDIM_ELEM(myfftV, n)));
386  }
387  ++n;
388  }
389  }
390  }
391 
392 
393 // #ifdef DEBUG_DIR
396 // Image<double> direction;
397 // direction = coneVol;
398 // direction.write(formatString("cone_%i_%i.vol", dir+1, count));
400 // #endif
401 
402  transformer_inv.inverseFourierTransform(fftVRiesz, amplitude);
403 
404 // #ifdef DEBUG_DIR
405 // Image<double> filteredvolume;
406 // filteredvolume = VRiesz;
407 // filteredvolume.write(formatString("Volumen_filtrado_%i_%i.vol", dir+1,count));
408 // #endif
409 
410 
411 // amplitude.resizeNoCopy(VRiesz);
412  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
413  DIRECT_MULTIDIM_ELEM(amplitude,n) *= DIRECT_MULTIDIM_ELEM(amplitude,n);
414 
415 
416  // Calculate first component of Riesz vector
417  double ux;
418  n=0;
419  for(size_t k=0; k<ZSIZE(myfftV); ++k)
420  {
421  for(size_t i=0; i<YSIZE(myfftV); ++i)
422  {
423  for(size_t j=0; j<XSIZE(myfftV); ++j)
424  {
425  ux = VEC_ELEM(freq_fourier,j);
426  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = ux*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
427  ++n;
428  }
429  }
430  }
431 
432  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
433 
434  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
435  {
436  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
437  }
438 
439  // Calculate second and third component of Riesz vector
440  n=0;
441  double uy, uz;
442  for(size_t k=0; k<ZSIZE(myfftV); ++k)
443  {
444  uz = VEC_ELEM(freq_fourier,k);
445  for(size_t i=0; i<YSIZE(myfftV); ++i)
446  {
447  uy = VEC_ELEM(freq_fourier,i);
448  for(size_t j=0; j<XSIZE(myfftV); ++j)
449  {
450  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = uz*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
451  DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = uy*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
452  ++n;
453  }
454  }
455  }
456 
457  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
458 
459  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
460  DIRECT_MULTIDIM_ELEM(amplitude,n) += DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
461 
462  transformer_inv.inverseFourierTransform(fftVRiesz_aux, VRiesz);
463 
464 
465 // amplitude.setXmippOrigin();
466  int z_size = ZSIZE(amplitude);
467  int siz = z_size*0.5;
468 
469  double limit_radius = (siz-N_smoothing);
470  n=0;
471  for(int k=0; k<z_size; ++k)
472  {
473  uz = (k - siz);
474  uz *= uz;
475  for(int i=0; i<z_size; ++i)
476  {
477  uy = (i - siz);
478  uy *= uy;
479  for(int j=0; j<z_size; ++j)
480  {
481  ux = (j - siz);
482  ux *= ux;
483  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
484  DIRECT_MULTIDIM_ELEM(amplitude,n)=sqrt(DIRECT_MULTIDIM_ELEM(amplitude,n));
485  double radius = sqrt(ux + uy + uz);
486  if ((radius>=limit_radius) && (radius<=siz))
487  DIRECT_MULTIDIM_ELEM(amplitude, n) *= 0.5*(1+cos(PI*(limit_radius-radius)/N_smoothing));
488  else if (radius>siz)
489  DIRECT_MULTIDIM_ELEM(amplitude, n) = 0;
490  ++n;
491  }
492  }
493  }
494 
495  //TODO: change (k - z_size*0.5)
496 
497 // #ifdef MONO_AMPLITUDE
498 // Image<double> saveImg2;
499 // saveImg2 = amplitude;
500 // if (fnDebug.c_str() != "")
501 // {
502 // FileName iternumber = formatString("smoothed_volume_%i_%i.vol", dir+1, count);
503 // saveImg2.write(fnDebug+iternumber);
504 // }
505 // saveImg2.clear();
506 // #endif
507 
508 
509  transformer_inv.FourierTransform(amplitude, fftVRiesz, false);
510 
511  double raised_w = PI/(freqL-freq);
512 
513  n=0;
514  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(fftVRiesz)
515  {
516  double un=1.0/DIRECT_MULTIDIM_ELEM(iu,n);
517  if (freqL>=un && un>=freq)
518  {
519  DIRECT_MULTIDIM_ELEM(fftVRiesz,n) *= 0.5*(1 + cos(raised_w*(un-freq)));
520  }
521  else
522  {
523  if (un>freqL)
524  {
525  DIRECT_MULTIDIM_ELEM(fftVRiesz,n) = 0;
526  }
527  }
528  }
529 
530  transformer_inv.inverseFourierTransform();
531 
532 // #ifdef MONO_AMPLITUDE
533 
534 // if (fnDebug.c_str() != "")
535 // {
536 // Image<double> saveImg2;
537 // saveImg2 = amplitude;
538 // FileName iternumber = formatString("_Filtered_Amplitude_%i_%i.vol", dir+1, count);
539 // saveImg2.write(fnDebug+iternumber);
540 // }
541 // #endif // DEBUG
542 }
ProgResDir::transformer_inv
FourierTransformer transformer_inv
Definition: resolution_directional.h:129
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgResDir::fftVRiesz
MultidimArray< std::complex< double > > fftVRiesz
Definition: resolution_directional.h:127
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
FourierTransformer::FourierTransform
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
ProgResDir::freq_fourier
Matrix1D< double > freq_fourier
Definition: resolution_directional.h:132
j
#define j
Definition: numerical_recipes.cpp:2493
ProgResDir::fftVRiesz_aux
MultidimArray< std::complex< double > > fftVRiesz_aux
Definition: resolution_directional.h:127
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
PI
#define PI
Definition: tools.h:43
n
int * n
Definition: numerical_recipes.cpp:2229

◆ defineCone()

void ProgResDir::defineCone ( MultidimArray< std::complex< double > > &  myfftV,
MultidimArray< std::complex< double > > &  conefilter,
double  rot,
double  tilt 
)

Definition at line 545 of file resolution_directional.cpp.

547 {
548 // conefilter.initZeros(myfftV);
549  conefilter = myfftV;
550  // Filter the input volume and add it to amplitude
551 
552 // MultidimArray<double> conetest;
553 // conetest.initZeros(myfftV);
554 // #ifdef DEBUG_DIR
555 // MultidimArray<double> coneVol;
556 // coneVol.initZeros(iu);
557 // #endif
558 
559  double x_dir, y_dir, z_dir;
560 
561  x_dir = sin(tilt*PI/180)*cos(rot*PI/180);
562  y_dir = sin(tilt*PI/180)*sin(rot*PI/180);
563  z_dir = cos(tilt*PI/180);
564 
565 // double ang_con = 10*PI/180;
566  double ang_con = 15*PI/180;
567 
568  double uz, uy, ux;
569  long n = 0;
570  for(size_t k=0; k<ZSIZE(myfftV); ++k)
571  {
572  uz = VEC_ELEM(freq_fourier,k);
573  uz *= z_dir;
574  for(size_t i=0; i<YSIZE(myfftV); ++i)
575  {
576  uy = VEC_ELEM(freq_fourier,i);
577  uy *= y_dir;
578  for(size_t j=0; j<XSIZE(myfftV); ++j)
579  {
580  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
581  ux = VEC_ELEM(freq_fourier,j);
582  ux *= x_dir;
583 
584  //BE CAREFULL with the order
585  //double dotproduct = (uy*x_dir + ux*y_dir + uz*z_dir)*iun;
586  iun *= (ux + uy + uz);
587  double acosine = acos(fabs(iun));
588 // DIRECT_MULTIDIM_ELEM(conetest, n) = log(real(conj(DIRECT_MULTIDIM_ELEM(myfftV, n))*DIRECT_MULTIDIM_ELEM(myfftV, n)));
589  if (acosine>ang_con)
590  {
591  DIRECT_MULTIDIM_ELEM(conefilter, n) = 0;
592 // DIRECT_MULTIDIM_ELEM(conetest, n) = 0;
593  }
594 /*
595  //4822.53 mean a smoothed cone angle of 20 degrees
596  double arg_exp = acosine*acosine*acosine*acosine*4822.53;
597  DIRECT_MULTIDIM_ELEM(conefilter, n) *= exp(-arg_exp);
598 // DIRECT_MULTIDIM_ELEM(conetest, n) = exp(-arg_exp);
599  */
600  ++n;
601  }
602  }
603  }
604 //
605 // Image<double> saveImg2;
606 // saveImg2 = conetest;
607 // saveImg2.write("cono.vol");
608 
609 }
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
ProgResDir::freq_fourier
Matrix1D< double > freq_fourier
Definition: resolution_directional.h:132
j
#define j
Definition: numerical_recipes.cpp:2493
PI
#define PI
Definition: tools.h:43
n
int * n
Definition: numerical_recipes.cpp:2229

◆ defineParams()

void ProgResDir::defineParams ( )
virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 60 of file resolution_directional.cpp.

61 {
62  addUsageLine("This function determines the local resolution of a map");
63  addParamsLine(" --vol <vol_file=\"\"> : Input volume");
64  addParamsLine(" --mask <vol_file=\"\"> : Mask defining the macromolecule");
65  addParamsLine(" -o <output=\"MGresolution.vol\"> : Local resolution volume (in Angstroms)");
66  addParamsLine(" [--sampling_rate <s=1>] : Sampling rate (A/px)");
67  addParamsLine(" [--resStep <s=0.5>] : Resolution step (precision) in A");
68  addParamsLine(" [--volumeRadius <s=100>] : This parameter determines the radius of a sphere where the volume is");
69  addParamsLine(" [--significance <s=0.95>] : The level of confidence for the hypothesis test.");
70  addParamsLine(" --radialRes <vol_file=\"\"> : Output radial resolution map");
71  addParamsLine(" --azimuthalRes <vol_file=\"\"> : Output azimuthal resolution map");
72  addParamsLine(" --highestResolutionVol <vol_file=\"\"> : Output highest resolution map");
73  addParamsLine(" --lowestResolutionVol <vol_file=\"\"> : Output lowest resolution map");
74  addParamsLine(" --doa1 <vol_file=\"\"> : Output doa1 map");
75  addParamsLine(" --doa2 <vol_file=\"\"> : Output doa2 map");
76  addParamsLine(" --radialAzimuthalThresholds <vol_file=\"\"> : Radial and azimuthal threshold for representation resolution maps");
77  addParamsLine(" --radialAvG <md_file=\"\"> : Radial Average of higesht, lowest, azimuthal, and radial, and local resolution map");
78  addParamsLine(" --monores <vol_file=\"\"> : Local resolution map");
79  addParamsLine(" --prefMin <vol_file=\"\"> : Metadata of highest resolution per direction");
80  addParamsLine(" [--threads <s=4>] : Number of threads");
81  addParamsLine(" --zScoremap <vol_file=\"\"> : Local zScore map, voxel with zscore higher than 3 are weird");
82  addParamsLine(" [--fast] : Fast computation");
83 }
XmippProgram::addUsageLine
void addUsageLine(const char *line, bool verbatim=false)
Definition: xmipp_program.cpp:282
XmippProgram::addParamsLine
void addParamsLine(const String &line)
Definition: xmipp_program.cpp:305

◆ diagSymMatrix3x3()

void ProgResDir::diagSymMatrix3x3 ( Matrix2D< double >  A,
Matrix1D< double > &  eigenvalues,
Matrix2D< double > &  P 
)

Definition at line 611 of file resolution_directional.cpp.

613 {
614  Matrix2D<double> B;
615  B.initZeros(3,3);
616 
617  MAT_ELEM(B,0,0) = 1;
618  MAT_ELEM(B,1,1) = 1;
619  MAT_ELEM(B,2,2) = 1;
620 
621  generalizedEigs(A, B, eigenvalues, eigenvectors);
622 }
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
generalizedEigs
void generalizedEigs(const Matrix2D< double > &A, const Matrix2D< double > &B, Matrix1D< double > &D, Matrix2D< double > &P)
Definition: matrix2d.cpp:267
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626

◆ ellipsoidFitting()

void ProgResDir::ellipsoidFitting ( Matrix2D< double > &  resolutionMat,
Matrix2D< double > &  axis 
)

Definition at line 828 of file resolution_directional.cpp.

830 {
831 
832  std::cout << "FITTIG" << std::endl;
833  double x, y, z, a, b, c, resolution, rot, tilt;
834  int numberdirections = angles.mdimx;
835  std::vector<double> list_distances;
836 
837  //std::cout << "xrows = " << xrows << std::endl;
838 
839  //MAT_ELEM(resolutionMat, direccion, resolucion)
840 
841  Matrix2D<double> ellipMat;
842  int dimMatrix = 0;
843  size_t mycounter;
844  Matrix2D<double> pseudoinv, quadricMatrix;
845  Matrix1D<double> onesVector, leastSquares, residuals, residualssorted;
846  Matrix2D<double> eigenvectors;
847  Matrix1D<double> eigenvalues;
848  //rows 1 2 3 (length axis a b c- where a is the smallest one)
849  //rows 4 5 6 x y z coordinates of the first eigenvector
850  //rows 7 8 9 x y z coordinates of the second eigenvector
851  //rows 10 11 12 x y z coordinates of the third eigenvector
852  axis.initZeros(12, NVoxelsOriginalMask);
853 
854  quadricMatrix.initZeros(3,3);
855  for (int k = 0; k<NVoxelsOriginalMask; ++k)
856  {
857  dimMatrix = 0;
858  for (int i = 0; i<numberdirections; ++i)
859  {
860  if (MAT_ELEM(resolutionMat, i, k) > 0)
861  ++dimMatrix;
862  }
863 
864  ellipMat.initZeros(dimMatrix, 6);
865  mycounter = 0; //It is required to store the matrix ellipMat
866  for (int i = 0; i<numberdirections; ++i)
867  {
868  resolution = MAT_ELEM(resolutionMat, i, k);
869 
870  if (resolution>0)
871  {
872  x = resolution*MAT_ELEM(trigProducts, 0, i);
873  y = resolution*MAT_ELEM(trigProducts, 1, i);
874  z = resolution*MAT_ELEM(trigProducts, 2, i);
875 
876  MAT_ELEM(ellipMat, mycounter, 0) = x*x;
877  MAT_ELEM(ellipMat, mycounter, 1) = y*y;
878  MAT_ELEM(ellipMat, mycounter, 2) = z*z;
879  MAT_ELEM(ellipMat, mycounter, 3) = 2*x*y;
880  MAT_ELEM(ellipMat, mycounter, 4) = 2*x*z;
881  MAT_ELEM(ellipMat, mycounter, 5) = 2*y*z;
882  ++mycounter;
883  }
884  }
885 
886 
887  ellipMat.inv(pseudoinv);
888 
889  onesVector.initConstant(mycounter, 1.0);
890  leastSquares = pseudoinv*onesVector;
891 
892  //Removing outliers
893  residuals = ellipMat*leastSquares - onesVector;
894  residualssorted = residuals.sort();
895 
896  double threshold_plus = VEC_ELEM(residualssorted, size_t(residualssorted.size()*significance));
897  double threshold_minus = VEC_ELEM(residualssorted, size_t(residualssorted.size()*(1.0-significance)));
898 
899  mycounter = 0;
900  size_t ellipsoidcounter = 0;
901  for (int i = 0; i<numberdirections; ++i)
902  {
903  resolution = MAT_ELEM(resolutionMat, i, k);
904 
905  if (resolution>0)
906  {
907  if ( (VEC_ELEM(residuals, mycounter) > threshold_plus) ||
908  (VEC_ELEM(residuals, mycounter) < threshold_minus) )
909  {
910  MAT_ELEM(resolutionMat, i, k) = -1;
911  }
912  else
913  {
914  ++ellipsoidcounter;
915  }
916  ++mycounter;
917  }
918  }
919 
920  ellipMat.initZeros(ellipsoidcounter, 6);
921  mycounter = 0; //It is required to store the matrix ellipMat
922  for (int i = 0; i<numberdirections; ++i)
923  {
924  resolution = MAT_ELEM(resolutionMat, i, k);
925 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
926 // {
927 // std::cout << k << " " << resolution << " " << MAT_ELEM(trigProducts, 0, i)
928 // << " " << MAT_ELEM(trigProducts, 1, i)
929 // << " " << MAT_ELEM(trigProducts, 2, i) << ";" << std::endl;
930 // }
931 
932  if (resolution>0)
933  {
934  x = resolution*MAT_ELEM(trigProducts, 0, i);
935  y = resolution*MAT_ELEM(trigProducts, 1, i);
936  z = resolution*MAT_ELEM(trigProducts, 2, i);
937 
938  MAT_ELEM(ellipMat, mycounter, 0) = x*x;
939  MAT_ELEM(ellipMat, mycounter, 1) = y*y;
940  MAT_ELEM(ellipMat, mycounter, 2) = z*z;
941  MAT_ELEM(ellipMat, mycounter, 3) = 2*x*y;
942  MAT_ELEM(ellipMat, mycounter, 4) = 2*x*z;
943  MAT_ELEM(ellipMat, mycounter, 5) = 2*y*z;
944  ++mycounter;
945  }
946  }
947 
948  // defining ellipsoid
949 
950  ellipMat.inv(pseudoinv);
951 
952  onesVector.initConstant(mycounter, 1.0);
953  leastSquares = pseudoinv*onesVector;
954 
955 
956  MAT_ELEM(quadricMatrix, 0, 0) = VEC_ELEM(leastSquares, 0);
957  MAT_ELEM(quadricMatrix, 0, 1) = VEC_ELEM(leastSquares, 3);
958  MAT_ELEM(quadricMatrix, 0, 2) = VEC_ELEM(leastSquares, 4);
959  MAT_ELEM(quadricMatrix, 1, 0) = VEC_ELEM(leastSquares, 3);
960  MAT_ELEM(quadricMatrix, 1, 1) = VEC_ELEM(leastSquares, 1);
961  MAT_ELEM(quadricMatrix, 1, 2) = VEC_ELEM(leastSquares, 5);
962  MAT_ELEM(quadricMatrix, 2, 0) = VEC_ELEM(leastSquares, 4);
963  MAT_ELEM(quadricMatrix, 2, 1) = VEC_ELEM(leastSquares, 5);
964  MAT_ELEM(quadricMatrix, 2, 2) = VEC_ELEM(leastSquares, 2);
965 
966  diagSymMatrix3x3(quadricMatrix, eigenvalues, eigenvectors);
967 
968  if (VEC_ELEM(eigenvalues, 0)<0)
969  {//This is de the vectorial product
970  VEC_ELEM(eigenvalues, 0) = VEC_ELEM(eigenvalues, 1);
971  MAT_ELEM(axis,3, k) = MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,2,2)-
972  MAT_ELEM(eigenvectors,2,1)*MAT_ELEM(eigenvectors,1,2);
973  MAT_ELEM(axis,4, k) = MAT_ELEM(eigenvectors,2,1)*MAT_ELEM(eigenvectors,0,2)-
974  MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,2,2);
975  MAT_ELEM(axis,5, k) = MAT_ELEM(eigenvectors,0,1)*MAT_ELEM(eigenvectors,1,2)-
976  MAT_ELEM(eigenvectors,1,1)*MAT_ELEM(eigenvectors,0,2);
977  }
978 
979  a = 1/sqrt(VEC_ELEM(eigenvalues, 0));
980  b = 1/sqrt(VEC_ELEM(eigenvalues, 1));
981  c = 1/sqrt(VEC_ELEM(eigenvalues, 2));
982 
983 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
984 // {
985 // std::cout << "a = " << a << std::endl;
986 // std::cout << "b = " << b << std::endl;
987 // std::cout << "c = " << c << std::endl;
988 // std::cout << "=--------------=" << std::endl;
989 // }
990 
991  MAT_ELEM(axis,0, k) = a;
992  MAT_ELEM(axis,1, k) = b;
993  MAT_ELEM(axis,2, k) = c;
994 
995  MAT_ELEM(axis,3, k) = MAT_ELEM(eigenvectors,0,0);
996  MAT_ELEM(axis,4, k) = MAT_ELEM(eigenvectors,1,0);
997  MAT_ELEM(axis,5, k) = MAT_ELEM(eigenvectors,2,0);
998 
999  MAT_ELEM(axis,6, k) = MAT_ELEM(eigenvectors,0,1);
1000  MAT_ELEM(axis,7, k) = MAT_ELEM(eigenvectors,1,1);
1001  MAT_ELEM(axis,8, k) = MAT_ELEM(eigenvectors,2,1);
1002 
1003  MAT_ELEM(axis,9, k) = MAT_ELEM(eigenvectors,0,2);
1004  MAT_ELEM(axis,10, k) = MAT_ELEM(eigenvectors,1,2);
1005  MAT_ELEM(axis,11, k) = MAT_ELEM(eigenvectors,2,2);
1006  }
1007 }
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
Matrix1D::size
size_t size() const
Definition: matrix1d.h:508
c
doublereal * c
Definition: numerical_recipes.cpp:2230
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
y
static double * y
Definition: numerical_recipes.cpp:8487
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
b
doublereal * b
Definition: numerical_recipes.cpp:2230
z
double z
Definition: numerical_recipes.cpp:8490
ProgResDir::diagSymMatrix3x3
void diagSymMatrix3x3(Matrix2D< double > A, Matrix1D< double > &eigenvalues, Matrix2D< double > &P)
Definition: resolution_directional.cpp:611
ProgResDir::trigProducts
Matrix2D< double > trigProducts
Definition: resolution_directional.h:131
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
Matrix1D::sort
Matrix1D< T > sort() const
Definition: matrix1d.cpp:850
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
a
doublereal * a
Definition: numerical_recipes.cpp:2230

◆ firstMonoResEstimation()

double ProgResDir::firstMonoResEstimation ( MultidimArray< std::complex< double > > &  myfftV,
double  w1,
double  w1l,
MultidimArray< double > &  amplitude 
)

Definition at line 1365 of file resolution_directional.cpp.

1367 {
1368  fftVRiesz.initZeros(myfftV);
1369  amplitude.resizeNoCopy(VRiesz);
1370  std::complex<double> J(0,1);
1371 
1372  // Filter the input volume and add it to amplitude
1373  long n=0;
1374  double iw=1.0/freq;
1375  double iwl=1.0/freqH;
1376  double ideltal=PI/(freq-freqH);
1377 
1378  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(myfftV)
1379  {
1380  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1381  double un=1.0/iun;
1382  if (freqH<=un && un<=freq)
1383  {
1384  //double H=0.5*(1+cos((un-w1)*ideltal));
1385  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1386  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
1387  } else if (un>freq)
1388  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1389  }
1390 
1391  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1392 
1393  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1394  DIRECT_MULTIDIM_ELEM(amplitude,n)=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1395 
1396  // Calculate first component of Riesz vector
1397  fftVRiesz.initZeros(myfftV);
1398  double uz, uy, ux;
1399  n=0;
1400 
1401  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1402  {
1403  for(size_t i=0; i<YSIZE(myfftV); ++i)
1404  {
1405  for(size_t j=0; j<XSIZE(myfftV); ++j)
1406  {
1407  ux = VEC_ELEM(freq_fourier,j);
1408  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1409  double un=1.0/iun;
1410  if (freqH<=un && un<=freq)
1411  {
1412  //double H=0.5*(1+cos((un-w1)*ideltal));
1413  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*ux*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1414  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *=0.5*(1+cos((un-w1)*ideltal));//H;
1415  //Next lines are an optimization of the commented ones
1416  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1417  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1418  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -ux*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1419  } else if (un>freq)
1420  {
1421  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-ux*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1422  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1423  }
1424  ++n;
1425  }
1426  }
1427  }
1428  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1429  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1430  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1431 
1432  // Calculate second component of Riesz vector
1433  fftVRiesz.initZeros(myfftV);
1434  n=0;
1435 
1436  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1437  {
1438  for(size_t i=0; i<YSIZE(myfftV); ++i)
1439  {
1440  uy = VEC_ELEM(freq_fourier,i);
1441  for(size_t j=0; j<XSIZE(myfftV); ++j)
1442  {
1443  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1444  double un=1.0/iun;
1445  if (freqH<=un && un<=freq)
1446  {
1447  //double H=0.5*(1+cos((un-w1)*ideltal));
1448  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*uy*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1449  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *=0.5*(1+cos((un-w1)*ideltal));//H;
1450  //Next lines are an optimization of the commented ones
1451  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1452  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1453  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -uy*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1454  } else if (un>freq)
1455  {
1456  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-uy*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1457  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1458  }
1459  ++n;
1460  }
1461  }
1462  }
1463  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1464  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1465  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1466 
1467  // Calculate third component of Riesz vector
1468  fftVRiesz.initZeros(myfftV);
1469  n=0;
1470  for(size_t k=0; k<ZSIZE(myfftV); ++k)
1471  {
1472  uz = VEC_ELEM(freq_fourier,k);
1473  for(size_t i=0; i<YSIZE(myfftV); ++i)
1474  {
1475  for(size_t j=0; j<XSIZE(myfftV); ++j)
1476  {
1477  double iun=DIRECT_MULTIDIM_ELEM(iu,n);
1478  double un=1.0/iun;
1479  if (freqH<=un && un<=freq)
1480  {
1481  //double H=0.5*(1+cos((un-w1)*ideltal));
1482  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-J*uz*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1483  //DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-w1)*ideltal));//H;
1484  //Next lines are an optimization of the commented ones
1485  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
1486  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1487  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= -uz*iun*0.5*(1+cos((un-freq)*ideltal));//H;
1488  } else if (un>freq)
1489  {
1490  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = (-uz*iun)*DIRECT_MULTIDIM_ELEM(myfftV, n);
1491  DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= J;
1492  }
1493  ++n;
1494  }
1495  }
1496  }
1497  transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
1498  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1499  {
1500  DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
1501  DIRECT_MULTIDIM_ELEM(amplitude,n)=sqrt(DIRECT_MULTIDIM_ELEM(amplitude,n));
1502  }
1503 
1504  // Low pass filter the monogenic amplitude
1505  // Prepare low pass filter
1506  FourierFilter lowPassFilter, FilterBand;
1507  lowPassFilter.FilterShape = RAISED_COSINE;
1508  lowPassFilter.raised_w = 0.01;
1509  lowPassFilter.do_generate_3dmask = false;
1510  lowPassFilter.FilterBand = LOWPASS;
1511  lowPassFilter.w1 = freq;
1512  amplitude.setXmippOrigin();
1513  lowPassFilter.applyMaskSpace(amplitude);
1514 
1515  double sumS=0, sumS2=0, sumN=0, sumN2=0, NN = 0, NS = 0;
1516  MultidimArray<int> &pMask = mask();
1517  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
1518  {
1519  double amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitude, n);
1520  if (DIRECT_MULTIDIM_ELEM(pMask, n)==0)
1521  {
1522  sumN += amplitudeValue;
1523  sumN2 += amplitudeValue*amplitudeValue;
1524  ++NN;
1525  }
1526  }
1527 
1528  double mean_noise = sumN/NN;
1529  return mean_noise;
1530 
1531 }
ProgResDir::transformer_inv
FourierTransformer transformer_inv
Definition: resolution_directional.h:129
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
FourierTransformer::inverseFourierTransform
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgResDir::fftVRiesz
MultidimArray< std::complex< double > > fftVRiesz
Definition: resolution_directional.h:127
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
RAISED_COSINE
#define RAISED_COSINE
Definition: fourier_filter.h:72
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
ProgResDir::freq_fourier
Matrix1D< double > freq_fourier
Definition: resolution_directional.h:132
j
#define j
Definition: numerical_recipes.cpp:2493
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
PI
#define PI
Definition: tools.h:43
n
int * n
Definition: numerical_recipes.cpp:2229
LOWPASS
#define LOWPASS
Definition: fourier_filter.h:77
FourierFilter::applyMaskSpace
void applyMaskSpace(MultidimArray< double > &v)
Definition: fourier_filter.cpp:710

◆ generateGridProjectionMatching()

void ProgResDir::generateGridProjectionMatching ( Matrix2D< double > &  angles)

Definition at line 205 of file resolution_directional.cpp.

206 {
207  if (fastCompute == false)
208  {
209  angles.initZeros(2,81);
210  MAT_ELEM(angles, 0, 0) = 0.000000; MAT_ELEM(angles, 1, 0) = 0.000000;
211  MAT_ELEM(angles, 0, 1) = 36.000000; MAT_ELEM(angles, 1, 1) = 15.858741;
212  MAT_ELEM(angles, 0, 2) = 36.000000; MAT_ELEM(angles, 1, 2) = 31.717482;
213  MAT_ELEM(angles, 0, 3) = 36.000000; MAT_ELEM(angles, 1, 3) = 47.576224;
214  MAT_ELEM(angles, 0, 4) = 36.000000; MAT_ELEM(angles, 1, 4) = 63.434965;
215  MAT_ELEM(angles, 0, 5) = 62.494295; MAT_ELEM(angles, 1, 5) = -76.558393;
216  MAT_ELEM(angles, 0, 6) = 54.000000; MAT_ELEM(angles, 1, 6) = 90.000000;
217  MAT_ELEM(angles, 0, 7) = 45.505705; MAT_ELEM(angles, 1, 7) = 76.558393;
218  MAT_ELEM(angles, 0, 8) = 108.000000; MAT_ELEM(angles, 1, 8) = 15.858741;
219  MAT_ELEM(angles, 0, 9) = 108.000000; MAT_ELEM(angles, 1, 9) = 31.717482;
220  MAT_ELEM(angles, 0, 10) = 108.000000; MAT_ELEM(angles, 1, 10) = 47.576224;
221  MAT_ELEM(angles, 0, 11) = 108.000000; MAT_ELEM(angles, 1, 11) = 63.434965;
222  MAT_ELEM(angles, 0, 12) = 134.494295; MAT_ELEM(angles, 1, 12) = -76.558393;
223  MAT_ELEM(angles, 0, 13) = 126.000000; MAT_ELEM(angles, 1, 13) = 90.000000;
224  MAT_ELEM(angles, 0, 14) = 117.505705; MAT_ELEM(angles, 1, 14) = 76.558393;
225  MAT_ELEM(angles, 0, 15) = 144.000000; MAT_ELEM(angles, 1, 15) = -15.858741;
226  MAT_ELEM(angles, 0, 16) = 144.000000; MAT_ELEM(angles, 1, 16) = -31.717482;
227  MAT_ELEM(angles, 0, 17) = 144.000000; MAT_ELEM(angles, 1, 17) = -47.576224;
228  MAT_ELEM(angles, 0, 18) = 144.000000; MAT_ELEM(angles, 1, 18) = -63.434965;
229  MAT_ELEM(angles, 0, 19) = 170.494295; MAT_ELEM(angles, 1, 19) = 76.558393;
230  MAT_ELEM(angles, 0, 20) = 162.000000; MAT_ELEM(angles, 1, 20) = 90.000000;
231  MAT_ELEM(angles, 0, 21) = 153.505705; MAT_ELEM(angles, 1, 21) = -76.558393;
232  MAT_ELEM(angles, 0, 22) = 72.000000; MAT_ELEM(angles, 1, 22) = -15.858741;
233  MAT_ELEM(angles, 0, 23) = 72.000000; MAT_ELEM(angles, 1, 23) = -31.717482;
234  MAT_ELEM(angles, 0, 24) = 72.000000; MAT_ELEM(angles, 1, 24) = -47.576224;
235  MAT_ELEM(angles, 0, 25) = 72.000000; MAT_ELEM(angles, 1, 25) = -63.434965;
236  MAT_ELEM(angles, 0, 26) = 98.494295; MAT_ELEM(angles, 1, 26) = 76.558393;
237  MAT_ELEM(angles, 0, 27) = 90.000000; MAT_ELEM(angles, 1, 27) = 90.000000;
238  MAT_ELEM(angles, 0, 28) = 81.505705; MAT_ELEM(angles, 1, 28) = -76.558393;
239  MAT_ELEM(angles, 0, 29) = 0.000000; MAT_ELEM(angles, 1, 29) = -15.858741;
240  MAT_ELEM(angles, 0, 30) = 0.000000; MAT_ELEM(angles, 1, 30) = -31.717482;
241  MAT_ELEM(angles, 0, 31) = 0.000000; MAT_ELEM(angles, 1, 31) = -47.576224;
242  MAT_ELEM(angles, 0, 32) = 0.000000; MAT_ELEM(angles, 1, 32) = -63.434965;
243  MAT_ELEM(angles, 0, 33) = 26.494295; MAT_ELEM(angles, 1, 33) = 76.558393;
244  MAT_ELEM(angles, 0, 34) = 18.000000; MAT_ELEM(angles, 1, 34) = 90.000000;
245  MAT_ELEM(angles, 0, 35) = 9.505705; MAT_ELEM(angles, 1, 35) = -76.558393;
246  MAT_ELEM(angles, 0, 36) = 12.811021; MAT_ELEM(angles, 1, 36) = 42.234673;
247  MAT_ELEM(angles, 0, 37) = 18.466996; MAT_ELEM(angles, 1, 37) = 59.620797;
248  MAT_ELEM(angles, 0, 38) = 0.000000; MAT_ELEM(angles, 1, 38) = 90.000000;
249  MAT_ELEM(angles, 0, 39) = 8.867209; MAT_ELEM(angles, 1, 39) = 75.219088;
250  MAT_ELEM(angles, 0, 40) = 72.000000; MAT_ELEM(angles, 1, 40) = 26.565058;
251  MAT_ELEM(angles, 0, 41) = 59.188979; MAT_ELEM(angles, 1, 41) = 42.234673;
252  MAT_ELEM(angles, 0, 42) = 84.811021; MAT_ELEM(angles, 1, 42) = 42.234673;
253  MAT_ELEM(angles, 0, 43) = 53.533003; MAT_ELEM(angles, 1, 43) = 59.620797;
254  MAT_ELEM(angles, 0, 44) = 72.000000; MAT_ELEM(angles, 1, 44) = 58.282544;
255  MAT_ELEM(angles, 0, 45) = 90.466996; MAT_ELEM(angles, 1, 45) = 59.620797;
256  MAT_ELEM(angles, 0, 46) = 72.000000; MAT_ELEM(angles, 1, 46) = 90.000000;
257  MAT_ELEM(angles, 0, 47) = 63.132791; MAT_ELEM(angles, 1, 47) = 75.219088;
258  MAT_ELEM(angles, 0, 48) = 80.867209; MAT_ELEM(angles, 1, 48) = 75.219088;
259  MAT_ELEM(angles, 0, 49) = 144.000000; MAT_ELEM(angles, 1, 49) = 26.565058;
260  MAT_ELEM(angles, 0, 50) = 131.188979; MAT_ELEM(angles, 1, 50) = 42.234673;
261  MAT_ELEM(angles, 0, 51) = 156.811021; MAT_ELEM(angles, 1, 51) = 42.234673;
262  MAT_ELEM(angles, 0, 52) = 125.533003; MAT_ELEM(angles, 1, 52) = 59.620797;
263  MAT_ELEM(angles, 0, 53) = 144.000000; MAT_ELEM(angles, 1, 53) = 58.282544;
264  MAT_ELEM(angles, 0, 54) = 162.466996; MAT_ELEM(angles, 1, 54) = 59.620797;
265  MAT_ELEM(angles, 0, 55) = 144.000000; MAT_ELEM(angles, 1, 55) = 90.000000;
266  MAT_ELEM(angles, 0, 56) = 135.132791; MAT_ELEM(angles, 1, 56) = 75.219088;
267  MAT_ELEM(angles, 0, 57) = 152.867209; MAT_ELEM(angles, 1, 57) = 75.219088;
268  MAT_ELEM(angles, 0, 58) = 180.000000; MAT_ELEM(angles, 1, 58) = -26.565058;
269  MAT_ELEM(angles, 0, 59) = 167.188979; MAT_ELEM(angles, 1, 59) = -42.234673;
270  MAT_ELEM(angles, 0, 60) = 180.000000; MAT_ELEM(angles, 1, 60) = -58.282544;
271  MAT_ELEM(angles, 0, 61) = 161.533003; MAT_ELEM(angles, 1, 61) = -59.620797;
272  MAT_ELEM(angles, 0, 62) = 171.132791; MAT_ELEM(angles, 1, 62) = -75.219088;
273  MAT_ELEM(angles, 0, 63) = 108.000000; MAT_ELEM(angles, 1, 63) = -26.565058;
274  MAT_ELEM(angles, 0, 64) = 120.811021; MAT_ELEM(angles, 1, 64) = -42.234673;
275  MAT_ELEM(angles, 0, 65) = 95.188979; MAT_ELEM(angles, 1, 65) = -42.234673;
276  MAT_ELEM(angles, 0, 66) = 126.466996; MAT_ELEM(angles, 1, 66) = -59.620797;
277  MAT_ELEM(angles, 0, 67) = 108.000000; MAT_ELEM(angles, 1, 67) = -58.282544;
278  MAT_ELEM(angles, 0, 68) = 89.533003; MAT_ELEM(angles, 1, 68) = -59.620797;
279  MAT_ELEM(angles, 0, 69) = 108.000000; MAT_ELEM(angles, 1, 69) = 90.000000;
280  MAT_ELEM(angles, 0, 70) = 116.867209; MAT_ELEM(angles, 1, 70) = -75.219088;
281  MAT_ELEM(angles, 0, 71) = 99.132791; MAT_ELEM(angles, 1, 71) = -75.219088;
282  MAT_ELEM(angles, 0, 72) = 36.000000; MAT_ELEM(angles, 1, 72) = -26.565058;
283  MAT_ELEM(angles, 0, 73) = 48.811021; MAT_ELEM(angles, 1, 73) = -42.234673;
284  MAT_ELEM(angles, 0, 74) = 23.188979; MAT_ELEM(angles, 1, 74) = -42.234673;
285  MAT_ELEM(angles, 0, 75) = 54.466996; MAT_ELEM(angles, 1, 75) = -59.620797;
286  MAT_ELEM(angles, 0, 76) = 36.000000; MAT_ELEM(angles, 1, 76) = -58.282544;
287  MAT_ELEM(angles, 0, 77) = 17.533003; MAT_ELEM(angles, 1, 77) = -59.620797;
288  MAT_ELEM(angles, 0, 78) = 36.000000; MAT_ELEM(angles, 1, 78) = 90.000000;
289  MAT_ELEM(angles, 0, 79) = 44.867209; MAT_ELEM(angles, 1, 79) = -75.219088;
290  MAT_ELEM(angles, 0, 80) = 27.132791; MAT_ELEM(angles, 1, 80) = -75.219088;
291  }
292  else
293  {
294  angles.initZeros(2,47);
295  MAT_ELEM(angles, 0,1) =0; MAT_ELEM(angles, 1,1) =0;
296  MAT_ELEM(angles, 0,2) =36; MAT_ELEM(angles, 1,2) =21.145;
297  MAT_ELEM(angles, 0,3) =36; MAT_ELEM(angles, 1,3) =42.29;
298  MAT_ELEM(angles, 0,4) =36; MAT_ELEM(angles, 1,4) =63.435;
299  MAT_ELEM(angles, 0,5) =59.6043; MAT_ELEM(angles, 1,5) =-81.0207;
300  MAT_ELEM(angles, 0,6) =48.3957; MAT_ELEM(angles, 1,6) =81.0207;
301  MAT_ELEM(angles, 0,7) =108; MAT_ELEM(angles, 1,7) =21.145;
302  MAT_ELEM(angles, 0,8) =108; MAT_ELEM(angles, 1,8) =42.29;
303  MAT_ELEM(angles, 0,9) =108; MAT_ELEM(angles, 1,9) =63.435;
304  MAT_ELEM(angles, 0,10) =131.6043; MAT_ELEM(angles, 1,10) =-81.0207;
305  MAT_ELEM(angles, 0,11) =120.3957; MAT_ELEM(angles, 1,11) =81.0207;
306  MAT_ELEM(angles, 0,12) =144; MAT_ELEM(angles, 1,12) =-21.145;
307  MAT_ELEM(angles, 0,13) =144; MAT_ELEM(angles, 1,13) =-42.29;
308  MAT_ELEM(angles, 0,14) =144; MAT_ELEM(angles, 1,14) =-63.435;
309  MAT_ELEM(angles, 0,15) =167.6043; MAT_ELEM(angles, 1,15) =81.0207;
310  MAT_ELEM(angles, 0,16) =156.3957; MAT_ELEM(angles, 1,16) =-81.0207;
311  MAT_ELEM(angles, 0,17) =72; MAT_ELEM(angles, 1,17) =-21.145;
312  MAT_ELEM(angles, 0,18) =72; MAT_ELEM(angles, 1,18) =-42.29;
313  MAT_ELEM(angles, 0,19) =72; MAT_ELEM(angles, 1,19) =-63.435;
314  MAT_ELEM(angles, 0,20) =95.6043; MAT_ELEM(angles, 1,20) =81.0207;
315  MAT_ELEM(angles, 0,21) =84.3957; MAT_ELEM(angles, 1,21) =-81.0207;
316  MAT_ELEM(angles, 0,22) =0; MAT_ELEM(angles, 1,22) =-21.145;
317  MAT_ELEM(angles, 0,23) =0; MAT_ELEM(angles, 1,23) =-42.29;
318  MAT_ELEM(angles, 0,24) =0; MAT_ELEM(angles, 1,24) =-63.435;
319  MAT_ELEM(angles, 0,25) =23.6043; MAT_ELEM(angles, 1,25) =81.0207;
320  MAT_ELEM(angles, 0,26) =12.3957; MAT_ELEM(angles, 1,26) =-81.0207;
321  MAT_ELEM(angles, 0,27) =12.3756; MAT_ELEM(angles, 1,27) =58.8818;
322  MAT_ELEM(angles, 0,28) =72; MAT_ELEM(angles, 1,28) =36.349;
323  MAT_ELEM(angles, 0,29) =59.6244; MAT_ELEM(angles, 1,29) =58.8818;
324  MAT_ELEM(angles, 0,30) =84.3756; MAT_ELEM(angles, 1,30) =58.8818;
325  MAT_ELEM(angles, 0,31) =72; MAT_ELEM(angles, 1,31) =80.2161;
326  MAT_ELEM(angles, 0,32) =144; MAT_ELEM(angles, 1,32) =36.349;
327  MAT_ELEM(angles, 0,33) =131.6244; MAT_ELEM(angles, 1,33) =58.8818;
328  MAT_ELEM(angles, 0,34) =156.3756; MAT_ELEM(angles, 1,34) =58.8818;
329  MAT_ELEM(angles, 0,35) =144; MAT_ELEM(angles, 1,35) =80.2161;
330  MAT_ELEM(angles, 0,36) =180; MAT_ELEM(angles, 1,36) =-36.349;
331  MAT_ELEM(angles, 0,37) =167.6244; MAT_ELEM(angles, 1,37) =-58.8818;
332  MAT_ELEM(angles, 0,38) =180; MAT_ELEM(angles, 1,38) =-80.2161;
333  MAT_ELEM(angles, 0,39) =108; MAT_ELEM(angles, 1,39) =-36.349;
334  MAT_ELEM(angles, 0,40) =120.3756; MAT_ELEM(angles, 1,40) =-58.8818;
335  MAT_ELEM(angles, 0,41) =95.6244; MAT_ELEM(angles, 1,41) =-58.8818;
336  MAT_ELEM(angles, 0,42) =108; MAT_ELEM(angles, 1,42) =-80.2161;
337  MAT_ELEM(angles, 0,43) =36; MAT_ELEM(angles, 1,43) =-36.349;
338  MAT_ELEM(angles, 0,44) =48.3756; MAT_ELEM(angles, 1,44) =-58.8818;
339  MAT_ELEM(angles, 0,45) =23.6244; MAT_ELEM(angles, 1,45) =-58.8818;
340  MAT_ELEM(angles, 0,46) =36; MAT_ELEM(angles, 1,46) =-80.2161;
341  }
342 
343 }
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626

◆ produceSideInfo()

void ProgResDir::produceSideInfo ( )

Definition at line 85 of file resolution_directional.cpp.

86 {
87 
88  std::cout << "Starting..." << std::endl;
89 
90  Image<double> V;
91  V.read(fnVol);
92 
93  V().setXmippOrigin();
94 
95  //Sweeping the projection sphere
96  std::cout << "Obtaining angular projections..." << std::endl;
97  generateGridProjectionMatching(angles);
98 
99  FourierTransformer transformer;
100 
101  MultidimArray<double> &inputVol = V();
102  VRiesz.resizeNoCopy(inputVol);
103  N_freq = ZSIZE(inputVol);
104  maxRes = ZSIZE(inputVol);
105  minRes = 2*sampling;
106 
107  transformer_inv.setThreadsNumber(Nthr);
108 
109  transformer.FourierTransform(inputVol, fftV);
110  iu.initZeros(fftV);
111 
112  // Frequency volume
113  double uz, uy, ux, uz2, u2, uz2y2;
114  long n=0;
115 
116  for(size_t k=0; k<ZSIZE(fftV); ++k)
117  {
118  FFT_IDX2DIGFREQ(k,ZSIZE(inputVol),uz);
119  uz2=uz*uz;
120 
121  for(size_t i=0; i<YSIZE(fftV); ++i)
122  {
123  FFT_IDX2DIGFREQ(i,YSIZE(inputVol),uy);
124  uz2y2=uz2+uy*uy;
125 
126  for(size_t j=0; j<XSIZE(fftV); ++j)
127  {
128  FFT_IDX2DIGFREQ(j,XSIZE(inputVol), ux);
129  u2=uz2y2+ux*ux;
130  if ((k != 0) || (i != 0) || (j != 0))
131  DIRECT_MULTIDIM_ELEM(iu,n) = 1.0/sqrt(u2);
132  else
133  DIRECT_MULTIDIM_ELEM(iu,n) = 1e38;
134  ++n;
135  }
136  }
137  }
138 
139  // Prepare mask
140  MultidimArray<int> &pMask=mask();
141 
142  if (fnMask != "")
143  {
144  mask.read(fnMask);
145  mask().setXmippOrigin();
146  }
147  else
148  {
149  std::cout << "Error: a mask ought to be provided" << std::endl;
150  exit(0);
151  }
152 
153  N_smoothing = 7;
154  NVoxelsOriginalMask = 0;
155  double radius = 0;
156  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
157  {
158  if (A3D_ELEM(pMask, k, i, j) == 1)
159  {
160  if ((k*k + i*i + j*j)>radius)
161  radius = k*k + i*i + j*j;
162  }
163 // std::cout << "i j k " << i << " " << j << " " << k << std::endl;
164 
165  if (A3D_ELEM(pMask, k, i, j) == 1)
166  ++NVoxelsOriginalMask;
167  if (i*i+j*j+k*k > (R-N_smoothing)*(R-N_smoothing))
168  A3D_ELEM(pMask, k, i, j) = -1;
169  }
170  Rparticle = round(sqrt(radius));
171  std::cout << "particle radius = " << Rparticle << std::endl;
172  size_t xrows = angles.mdimx;
173 
174  resolutionMatrix.initConstant(xrows, NVoxelsOriginalMask, maxRes);
175 
176 
177  #ifdef DEBUG_MASK
178  std::cout << "-------------DEBUG-----------" <<std::endl;
179  std::cout << "Next number ought to be the same than number of directions"
180  << std::endl;
181  std::cout << "number of angles" << xrows << std::endl;
182  std::cout << "---------END-DEBUG--------" <<std::endl;
183  mask.write("mask.vol");
184  #endif
185 
186  freq_fourier.initZeros(ZSIZE(inputVol));
187  int size = ZSIZE(inputVol);
188  maxRes = size;
189  minRes = 1;
190  V.clear();
191 
192 
193  double u;
194  int size_fourier(ZSIZE(fftV));
195 
196  VEC_ELEM(freq_fourier,0) = 1e-38;
197  for(size_t k=1; k<size_fourier; ++k)
198  {
199  FFT_IDX2DIGFREQ(k,size, u);
200  VEC_ELEM(freq_fourier,k) = u;
201  }
202 }
ProgResDir::transformer_inv
FourierTransformer transformer_inv
Definition: resolution_directional.h:129
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
FourierTransformer::setThreadsNumber
void setThreadsNumber(int tNumber)
Definition: xmipp_fftw.h:119
Matrix2D::initConstant
void initConstant(T val)
Definition: matrix2d.h:602
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ImageBase::write
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)
Definition: xmipp_image_base.cpp:292
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
ProgResDir::resolutionMatrix
Matrix2D< double > resolutionMatrix
Definition: resolution_directional.h:131
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
FourierTransformer::FourierTransform
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
FFT_IDX2DIGFREQ
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
Definition: cuda_gpu_reconstruct_fourier.cpp:382
ProgResDir::freq_fourier
Matrix1D< double > freq_fourier
Definition: resolution_directional.h:132
j
#define j
Definition: numerical_recipes.cpp:2493
alglib::round
int round(double x)
Definition: ap.cpp:7245
u
doublereal * u
Definition: numerical_recipes.cpp:2230
ImageBase::read
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
Definition: xmipp_image_base.cpp:119
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
ProgResDir::fftV
MultidimArray< std::complex< double > > fftV
Definition: resolution_directional.h:130
n
int * n
Definition: numerical_recipes.cpp:2229
Image::clear
void clear()
Definition: xmipp_image.h:144
ProgResDir::generateGridProjectionMatching
void generateGridProjectionMatching(Matrix2D< double > &angles)
Definition: resolution_directional.cpp:205

◆ radialAverageInMask()

void ProgResDir::radialAverageInMask ( MultidimArray< int > &  mask,
MultidimArray< double > &  inputVol_1,
MultidimArray< double > &  inputVol_2,
MultidimArray< double > &  inputVol_3,
MultidimArray< double > &  inputVol_4,
MultidimArray< double > &  inputVol_5,
MetaDataVec &  md 
)

Definition at line 1010 of file resolution_directional.cpp.

1014 {
1015  double u_inf, u_sup, u;
1016 
1017  int step = 1;
1018 
1019  double N, NN;
1020  MultidimArray<double> radialAvg(XSIZE(inputVol_1)*0.5);
1021  int uk, uj, ui;
1022 
1023  int siz = XSIZE(inputVol_1);
1024  size_t objId;
1025 
1026  inputVol_1.setXmippOrigin();
1027  inputVol_2.setXmippOrigin();
1028  inputVol_3.setXmippOrigin();
1029  inputVol_4.setXmippOrigin();
1030  inputVol_5.setXmippOrigin();
1031  MultidimArray<double> zVolume;
1032  zVolume.initZeros(inputVol_1);
1033  zVolume.setXmippOrigin();
1034 
1035  MultidimArray<int> &pMask = mask;
1036  pMask.setXmippOrigin();
1037 
1038  Matrix2D<double> std_mean_Radial_1, std_mean_Radial_2, std_mean_Radial_3,
1039  std_mean_Radial_4, std_mean_Radial_5;
1040 
1041  std_mean_Radial_1.initZeros(2, (size_t) siz*0.5 + 1);
1042  std_mean_Radial_2.initZeros(2, (size_t) siz*0.5 + 1);
1043  std_mean_Radial_3.initZeros(2, (size_t) siz*0.5 + 1);
1044  std_mean_Radial_4.initZeros(2, (size_t) siz*0.5 + 1);
1045  std_mean_Radial_5.initZeros(2, (size_t) siz*0.5 + 1);
1046 
1047  for(size_t kk=1; kk<siz*0.5; ++kk)
1048  {
1049  double cum_mean_1 = 0;
1050  double cum_mean_2 = 0;
1051  double cum_mean_3 = 0;
1052  double cum_mean_4 = 0;
1053  double cum_mean_5 = 0;
1054 
1055  double cum2_mean_1 = 0;
1056  double cum2_mean_2 = 0;
1057  double cum2_mean_3 = 0;
1058  double cum2_mean_4 = 0;
1059  double cum2_mean_5 = 0;
1060 
1061  N = 0;
1062  u_sup = kk + step;
1063  u_inf = kk - step;
1064 
1065  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
1066  {
1067  if (A3D_ELEM(pMask, k, i, j)>0)
1068  {
1069  u = sqrt(k*k + i*i + j*j);
1070  if ((u<u_sup) && (u>=u_inf))
1071  {
1072  double aux;
1073  aux = A3D_ELEM(inputVol_1, k, i, j);
1074  cum_mean_1 += aux;
1075  cum2_mean_1 += aux*aux;
1076  aux = A3D_ELEM(inputVol_2, k, i, j);
1077  cum_mean_2 += aux;
1078  cum2_mean_2 += aux*aux;
1079  aux = A3D_ELEM(inputVol_3, k, i, j);
1080  cum_mean_3 += aux;
1081  cum2_mean_3 += aux*aux;
1082  aux = A3D_ELEM(inputVol_4, k, i, j);
1083  cum_mean_4 += aux;
1084  cum2_mean_4 += aux*aux;
1085  aux = A3D_ELEM(inputVol_5, k, i, j);
1086  cum_mean_5 += aux;
1087  cum2_mean_5 += aux*aux;
1088 
1089  N = N + 1;
1090  }
1091  }
1092  }
1093 
1094  double sigma1, sigma2, sigma3, sigma4, sigma5;
1095 
1096 
1097 
1098 // if ((cum_mean_1==0) || (cum_mean_2==0) || (cum_mean_3==0) || (cum_mean_4==0) || (cum_mean_5==0))
1099 // {
1100 // md.setValue(MDL_IDX, kk, objId);
1101 // md.setValue(MDL_VOLUME_SCORE1, cum_mean_1, objId);
1102 // md.setValue(MDL_VOLUME_SCORE2, cum_mean_2, objId);
1103 // md.setValue(MDL_VOLUME_SCORE3, cum_mean_3, objId);
1104 // md.setValue(MDL_VOLUME_SCORE4, cum_mean_4, objId);
1105 // md.setValue(MDL_AVG, cum_mean_5, objId);
1106 // }
1107 // else
1108 // {
1109  if ((cum_mean_1>0) || (cum_mean_2>0) || (cum_mean_3>0) || (cum_mean_4>0) || (cum_mean_5>0))
1110  {
1111  objId = md.addObject();
1112  cum_mean_1 = (cum_mean_1/N);
1113  cum_mean_2 = (cum_mean_2/N);
1114  cum_mean_3 = (cum_mean_3/N);
1115  cum_mean_4 = (cum_mean_4/N);
1116  cum_mean_5 = (cum_mean_5/N);
1117 
1118  MAT_ELEM(std_mean_Radial_1,1, kk) = cum_mean_1;
1119  MAT_ELEM(std_mean_Radial_2,1, kk) = cum_mean_2;
1120  MAT_ELEM(std_mean_Radial_3,1, kk) = cum_mean_3;
1121  MAT_ELEM(std_mean_Radial_4,1, kk) = cum_mean_4;
1122  MAT_ELEM(std_mean_Radial_5,1, kk) = cum_mean_5;
1123 
1124  md.setValue(MDL_IDX, kk, objId);
1125  md.setValue(MDL_VOLUME_SCORE1, cum_mean_1, objId);
1126  md.setValue(MDL_VOLUME_SCORE2, cum_mean_2, objId);
1127  md.setValue(MDL_VOLUME_SCORE3, cum_mean_3, objId);
1128  md.setValue(MDL_VOLUME_SCORE4, cum_mean_4, objId);
1129  md.setValue(MDL_AVG, cum_mean_5, objId);
1130 
1131  MAT_ELEM(std_mean_Radial_1,0, kk) = sqrt(cum2_mean_1/N - cum_mean_1*cum_mean_1);
1132  MAT_ELEM(std_mean_Radial_2,0, kk) = sqrt(cum2_mean_2/N - cum_mean_2*cum_mean_2);
1133  MAT_ELEM(std_mean_Radial_3,0, kk) = sqrt(cum2_mean_3/N - cum_mean_3*cum_mean_3);
1134  MAT_ELEM(std_mean_Radial_4,0, kk) = sqrt(cum2_mean_4/N - cum_mean_4*cum_mean_4+0.001);
1135  MAT_ELEM(std_mean_Radial_5,0, kk) = sqrt(cum2_mean_5/N - cum_mean_5*cum_mean_5);
1136 
1137 // std::cout << "cum2_mean_4/(N) = " << cum2_mean_4/(N) << " " << "cum_mean_4*cum_mean_4) = " << cum_mean_4*cum_mean_4 << std::endl;
1138 // std::cout << "MAT_ELEM(std_mean_Radial_1,0, kk) = " << MAT_ELEM(std_mean_Radial_1,0, kk) << " " << MAT_ELEM(std_mean_Radial_2,0, kk) << " " << MAT_ELEM(std_mean_Radial_3,0, kk) << " " << MAT_ELEM(std_mean_Radial_4,0, kk) << " " << MAT_ELEM(std_mean_Radial_5,0, kk)<< std::endl;
1139 // std::cout << "MAT_ELEM(std_mean_Radial_1,1, kk) = " << MAT_ELEM(std_mean_Radial_1,1, kk) << " " << MAT_ELEM(std_mean_Radial_2,1, kk) << " " << MAT_ELEM(std_mean_Radial_3,1, kk) << " " << MAT_ELEM(std_mean_Radial_4,1, kk) << " " << MAT_ELEM(std_mean_Radial_5,1, kk)<< std::endl;
1140  }
1141 
1142 
1143  double lastz;
1144  FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
1145  {
1146  if (A3D_ELEM(pMask, k, i, j)>0)
1147  {
1148  lastz = -1;
1149  u = sqrt(k*k + i*i + j*j);
1150  if ((u<u_sup) && (u>=u_inf) && (MAT_ELEM(std_mean_Radial_1,1, kk)>0))
1151  {
1152  double z, aux;
1153  aux = abs((A3D_ELEM(inputVol_1, k, i, j) - MAT_ELEM(std_mean_Radial_1,1, kk))/MAT_ELEM(std_mean_Radial_1,0, kk)) + + 0.002;
1154  if (aux > lastz)
1155  lastz = aux;
1156  aux = abs((A3D_ELEM(inputVol_2, k, i, j) - MAT_ELEM(std_mean_Radial_2,1, kk))/MAT_ELEM(std_mean_Radial_2,0, kk)) + 0.002;
1157  if (aux > lastz)
1158  lastz = aux;
1159  aux = abs((A3D_ELEM(inputVol_3, k, i, j) - MAT_ELEM(std_mean_Radial_3,1, kk))/MAT_ELEM(std_mean_Radial_3,0, kk)) + 0.002;
1160  if (aux > lastz)
1161  lastz = aux;
1162 // aux = abs((A3D_ELEM(inputVol_4, k, i, j) - MAT_ELEM(std_mean_Radial_4,1, kk))/MAT_ELEM(std_mean_Radial_4,0, kk));
1163 // if (aux > lastz)
1164 // lastz = aux;
1165  aux = abs((A3D_ELEM(inputVol_5, k, i, j) - MAT_ELEM(std_mean_Radial_5,1, kk))/MAT_ELEM(std_mean_Radial_5,0, kk)) + 0.002;
1166  if (aux > lastz)
1167  lastz = aux;
1168 
1169  if (lastz>5.0) //This line considers zscores higher than 5 as 5(saturated at 5sigma)
1170  lastz = 5;
1171 
1172  A3D_ELEM(zVolume, k, i, j) = lastz;
1173  }
1174  }
1175  }
1176  }
1177  Image<double> zVolumesave;
1178  zVolumesave = zVolume;
1179  zVolumesave.write(fnZscore);
1180 }
radialAvg
void radialAvg(Image< double > &op)
Definition: image_operate.cpp:269
MDL_VOLUME_SCORE3
< Score 2 for volumes
Definition: metadata_label.h:478
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ImageBase::write
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)
Definition: xmipp_image_base.cpp:292
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
MetaDataVec::addObject
size_t addObject() override
Definition: metadata_vec.cpp:435
MDL_VOLUME_SCORE4
< Score 3 for volumes
Definition: metadata_label.h:479
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
z
double z
Definition: numerical_recipes.cpp:8490
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_VOLUME_SCORE1
Definition: metadata_label.h:476
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
MDL_IDX
Index within a list (size_t)
Definition: metadata_label.h:245
u
doublereal * u
Definition: numerical_recipes.cpp:2230
MDL_VOLUME_SCORE2
< Score 1 for volumes
Definition: metadata_label.h:477
MDL_AVG
average value (double)
Definition: metadata_label.h:70
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723

◆ radialAzimuthalResolution()

void ProgResDir::radialAzimuthalResolution ( Matrix2D< double > &  resolutionMat,
MultidimArray< int > &  pmask,
MultidimArray< double > &  radial,
MultidimArray< double > &  azimuthal,
MultidimArray< double > &  lowestResolution,
MultidimArray< double > &  highestResolution,
MultidimArray< double > &  doaResolution_1,
MultidimArray< double > &  doaResolution_2,
double &  radial_Thr,
double &  azimuthal_Thr,
MetaDataVec &  mdprefDirs 
)

Definition at line 1182 of file resolution_directional.cpp.

1193 {
1194 
1195  radial.initZeros(pmask);
1196  azimuthal.initZeros(pmask);
1197  lowestResolution.initZeros(pmask);
1198  highestResolution.initZeros(pmask);
1199  doaResolution_1.initZeros(pmask);
1200  doaResolution_2.initZeros(pmask);
1201 
1202  double radial_angle = 45*PI/180;
1203  double azimuthal_resolution = 0;
1204  double radial_resolution = 0;
1205  double azimuthal_angle = 70*PI/180;
1206  double resolution, dotproduct, x, y, z, iu, arcos;
1207  int xrows = angles.mdimx;
1208  int idx = 0;
1209 
1210  double count_radial = 0.0, count_azimuthal = 0.0;
1211 
1212  Matrix1D<int> PrefferredDirHist, resolutionMeanVector;
1213  PrefferredDirHist.initZeros(xrows);
1214  resolutionMeanVector.initZeros(xrows);
1215 
1216 // meanResolution.initZeros(pmask);
1217  size_t objId;
1218  FOR_ALL_ELEMENTS_IN_ARRAY3D(pmask)
1219  {
1220  //i defines the direction and k the voxel
1221  if (A3D_ELEM(pmask,k,i,j) > 0 )
1222  {
1223  iu = 1/sqrt(i*i + j*j + k*k);
1224  count_radial = 0;
1225  count_azimuthal = 0;
1226  std::vector<double> ResList;
1227 
1228  double lastRes = 100; //A non-sense value
1229 
1230  for (int ii = 0; ii<xrows; ++ii)
1231  {
1232  resolution = MAT_ELEM(resolutionMat, ii, idx);
1233 
1234  if (resolution>0)
1235  {
1236  ResList.push_back(resolution);
1237  x = MAT_ELEM(trigProducts, 0, ii);
1238  y = MAT_ELEM(trigProducts, 1, ii);
1239  z = MAT_ELEM(trigProducts, 2, ii);
1240 
1241  dotproduct = (x*i + y*j + z*k)*iu;
1242  arcos = acos(fabs(dotproduct));
1243  if (arcos>=azimuthal_angle)
1244  {
1245  count_azimuthal = count_azimuthal + 1;
1246  azimuthal_resolution += resolution;
1247  }
1248  if (arcos<=radial_angle)
1249  {
1250  count_radial = count_radial + 1;
1251  radial_resolution += resolution;
1252  }
1253 
1254  }
1255 
1256  }
1257 
1258 // std::cout << "count_radial = " << count_radial << std::endl;
1259 // std::cout << "count_azimuthal = " << count_azimuthal << std::endl;
1260 // std::cout << " " << std::endl;
1261 
1262 // A3D_ELEM(meanResolution,k,i,j) = meanRes[(size_t) floor(0.5*meanRes.size())];
1263  std::sort(ResList.begin(),ResList.end());
1264 
1265  double Mres, mres, medianResol, res75, res25;
1266 
1267  Mres = ResList[ (size_t) floor(0.95*ResList.size()) ];
1268  A3D_ELEM(lowestResolution,k,i,j) = Mres;
1269 
1270  mres = ResList[ (size_t) floor(0.05*ResList.size()) ];
1271 
1272  res75 = ResList[ (size_t) floor(0.83*ResList.size()) ];
1273  res25 = ResList[ (size_t) floor(0.17*ResList.size()) ];
1274 
1275  A3D_ELEM(highestResolution,k,i,j) = mres;
1276 
1277  A3D_ELEM(doaResolution_1,k,i,j) = 0.5*(res75 - res25);//( (Mres - mres)/(Mres + mres) );
1278 
1279  A3D_ELEM(doaResolution_2,k,i,j) = 0.5*( Mres + mres );
1280 
1281  ResList.clear();
1282 
1283  //Prefferred directions
1284 
1285  for (int ii = 0; ii<xrows; ++ii)
1286  {
1287  resolution = MAT_ELEM(resolutionMat, ii, idx);
1288 
1289  if (resolution>0)
1290  {
1291  if ((mres>(resolution-0.1)) && (mres<(resolution+0.1)))
1292  {
1293  VEC_ELEM(PrefferredDirHist,ii) += 1;
1294  VEC_ELEM(resolutionMeanVector,ii) += resolution;
1295  }
1296  }
1297 
1298  }
1299  ++idx;
1300  }
1301 
1302  if (count_radial<1)
1303  A3D_ELEM(radial,k,i,j) = A3D_ELEM(doaResolution_2,k,i,j);
1304  else
1305  A3D_ELEM(radial,k,i,j) = radial_resolution/count_radial;
1306 
1307  if (count_azimuthal<1)
1308  A3D_ELEM(azimuthal,k,i,j) = A3D_ELEM(doaResolution_2,k,i,j);
1309  else
1310  A3D_ELEM(azimuthal,k,i,j) = azimuthal_resolution/count_azimuthal;
1311 
1312 
1313  azimuthal_resolution = 0;
1314  radial_resolution = 0;
1315  }
1316 
1317 
1318  for (size_t ii = 0; ii<xrows; ++ii)
1319  {
1320  objId = mdprefDirs.addObject();
1321  size_t con;
1322  con = (size_t) VEC_ELEM(PrefferredDirHist,ii);
1323 // std::cout << ii << " " << con << std::endl;
1324  double rot = MAT_ELEM(angles, 0, ii);
1325  double tilt = MAT_ELEM(angles, 1, ii);
1326 
1327  if (tilt<0)
1328  {
1329  tilt = abs(tilt);
1330  rot = rot + 180;
1331  }
1332  double meanRes;
1333 
1334  meanRes = VEC_ELEM(resolutionMeanVector,ii)/((double) VEC_ELEM(PrefferredDirHist,ii));
1335 
1336  mdprefDirs.setValue(MDL_ANGLE_ROT, rot, objId);
1337  mdprefDirs.setValue(MDL_ANGLE_TILT, tilt, objId);
1338  mdprefDirs.setValue(MDL_WEIGHT, (double) con, objId);
1339  mdprefDirs.setValue(MDL_RESOLUTION_FREQ, meanRes, objId);
1340  mdprefDirs.setValue(MDL_X, (double) ii, objId);
1341  mdprefDirs.setValue(MDL_COUNT, con, objId);
1342  }
1343  mdprefDirs.write(fnprefMin);
1344 
1345  std::vector<double> radialList, azimuthalList;
1346 
1347  FOR_ALL_ELEMENTS_IN_ARRAY3D(radial)
1348  {
1349  //i defines the direction and k the voxel
1350  if (A3D_ELEM(pmask,k,i,j) > 0 )
1351  {
1352  radialList.push_back(A3D_ELEM(radial,k,i,j));
1353  azimuthalList.push_back(A3D_ELEM(azimuthal,k,i,j));
1354  }
1355  }
1356 
1357  std::sort(radialList.begin(),radialList.end());
1358  std::sort(azimuthalList.begin(),azimuthalList.end());
1359 
1360  radial_Thr = radialList[(size_t) floor(radialList.size()*0.9)];
1361  azimuthal_Thr = azimuthalList[(size_t) floor(azimuthalList.size()*0.9)];
1362 }
MDL_ANGLE_ROT
Rotation angle of an image (double,degrees)
Definition: metadata_label.h:55
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
MDL_ANGLE_TILT
Tilting angle of an image (double,degrees)
Definition: metadata_label.h:59
y
static double * y
Definition: numerical_recipes.cpp:8487
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
A3D_ELEM
#define A3D_ELEM(V, k, i, j)
Definition: multidim_array_base.h:253
FOR_ALL_ELEMENTS_IN_ARRAY3D
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
Definition: multidim_array_base.h:269
MDL_COUNT
Number of elements of a type (int) [this is a genereic type do not use to transfer information to ano...
Definition: metadata_label.h:92
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
MetaDataVec::addObject
size_t addObject() override
Definition: metadata_vec.cpp:435
z
double z
Definition: numerical_recipes.cpp:8490
ProgResDir::iu
MultidimArray< double > iu
Definition: resolution_directional.h:128
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_RESOLUTION_FREQ
Frequency in 1/A (double)
Definition: metadata_label.h:388
ProgResDir::trigProducts
Matrix2D< double > trigProducts
Definition: resolution_directional.h:131
MDL_X
X component (double)
Definition: metadata_label.h:495
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
PI
#define PI
Definition: tools.h:43
MDL_WEIGHT
< Score 4 for volumes
Definition: metadata_label.h:480

◆ readParams()

void ProgResDir::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 35 of file resolution_directional.cpp.

36 {
37  fnVol = getParam("--vol");
38  fnOut = getParam("-o");
39  fnMask = getParam("--mask");
40  sampling = getDoubleParam("--sampling_rate");
41  R = getDoubleParam("--volumeRadius");
42  significance = getDoubleParam("--significance");
43  res_step = getDoubleParam("--resStep");
44  fnradial = getParam("--radialRes");
45  fnazimuthal = getParam("--azimuthalRes");
46  fnRadialAvG = getParam("--radialAvG");
47  fnHighestResolution = getParam("--highestResolutionVol");
48  fnLowestResolution = getParam("--lowestResolutionVol");
49  fnDoa1 = getParam("--doa1");
50  fnDoa2 = getParam("--doa2");
51  fnMDThr = getParam("--radialAzimuthalThresholds");
52  fnMonoRes = getParam("--monores");
53  fnprefMin = getParam("--prefMin");
54  fnZscore = getParam("--zScoremap");
55  Nthr = getIntParam("--threads");
56  fastCompute = checkParam("--fast");
57 }
XmippProgram::getDoubleParam
double getDoubleParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:341
ProgResDir::fnLowestResolution
FileName fnLowestResolution
Definition: resolution_directional.h:55
XmippProgram::getParam
const char * getParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:321
ProgResDir::fnHighestResolution
FileName fnHighestResolution
Definition: resolution_directional.h:55
XmippProgram::checkParam
bool checkParam(const char *param)
Definition: xmipp_program.cpp:379
XmippProgram::getIntParam
int getIntParam(const char *param, int arg=0)
Definition: xmipp_program.cpp:331

◆ removeOutliers()

void ProgResDir::removeOutliers ( Matrix2D< double > &  resolutionMat)

Definition at line 719 of file resolution_directional.cpp.

720 {
721  std::cout << "Removing outliers..." << std::endl;
722 
723  double x1, y1, z1, x2, y2, z2, distance, resolution, sigma,
724  rot, tilt, threshold, sigma2, lastMinDistance;
725  double meandistance = 0, distance_2 = 0;
726  int numberdirections = angles.mdimx, N=0, count = 0;
727 
728  double criticalZ = icdf_gauss(significance);
729  double threshold_gauss;
730  double counter = 0;
731 
732  double ang = 20.0;
733 
734  Matrix2D<double> neigbour_dir;
735 
736  for (int k = 0; k<NVoxelsOriginalMask; ++k)
737  {
738  meandistance = 0;
739  sigma = 0;
740  distance_2 = 0;
741  counter = 0;
742 
743 // std::vector<double> neighbours;
744  neigbour_dir.initZeros(numberdirections, 2);
745  //Computing closest neighbours and its mean distance
746  for (int i = 0; i<numberdirections; ++i)
747  {
748 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
749 // {
750 // std::cout << k << " " << MAT_ELEM(resolutionMat, i, k) << " " << MAT_ELEM(trigProducts, 0, i) << " " <<
751 // MAT_ELEM(trigProducts, 1, i) << " " << MAT_ELEM(trigProducts, 2, i) << ";" << std::endl;
752 // }
753 
754  double resi = MAT_ELEM(resolutionMat, i, k);
755  if (resi>0)
756  {
757  x1 = MAT_ELEM(trigProducts, 0, i);
758  y1 = MAT_ELEM(trigProducts, 1, i);
759  z1 = MAT_ELEM(trigProducts, 2, i);
760 
761  for (int j = 0; j<numberdirections; ++j)
762  {
763  if (i != j)
764  {
765  x2 = MAT_ELEM(trigProducts, 0, j);
766  y2 = MAT_ELEM(trigProducts, 1, j);
767  z2 = MAT_ELEM(trigProducts, 2, j);
768  double resj = MAT_ELEM(resolutionMat, j, k);
769  if (resj>0)
770  {
771  distance = (180/PI)*acos(x1*x2 + y1*y2 + z1*z2);
772 
773  if (distance < ang)
774  {
775  x2 *= resj;
776  y2 *= resj;
777  z2 *= resj;
778  distance = sqrt( (resi*x1-x2)*(resi*x1-x2) + (resi*y1-y2)*(resi*y1-y2) + (resi*z1-z2)*(resi*z1-z2) );
779  MAT_ELEM(neigbour_dir, i, 0) += distance;
780  MAT_ELEM(neigbour_dir, i, 1) += 1;
781 // neighbours.push_back(distance);
782  meandistance += distance;
783  ++counter;
784  sigma += distance*distance;
785  }
786  }
787  }
788  }
789  }
790  }
791 
792  double thresholdDirection;
793 
794  meandistance= meandistance/counter;
795  sigma = sigma/counter - meandistance*meandistance;
796 
797 // std::sort(neighbours.begin(), neighbours.end());
798 // thresholdDirection = neighbours[size_t(neighbours.size()*significance)];
799  threshold_gauss = meandistance + criticalZ*sqrt(sigma);
800 
801 // neighbours.clear();
802 
803  //A direction is an outlier if is significative higher than overal distibution
804 
805  for (int i = 0; i<numberdirections; ++i)
806  {
807  double meandistance;
808  if (MAT_ELEM(neigbour_dir, i, 0)>0)
809  {
810  meandistance = MAT_ELEM(neigbour_dir, i, 0)/MAT_ELEM(neigbour_dir, i, 1);
811 
812  if ((meandistance>threshold_gauss) || MAT_ELEM(neigbour_dir, i, 0) <= 1)
813  {
814  MAT_ELEM(resolutionMat, i, k)=-1;
815  }
816  }
817 
818  }
819 
820 // if ((k == 201311) || (k == 201312) || (k == 283336) || (k == 324353) || (k == 324362) || (k == 324512))
821 // std::cout << "threshold_gauss--------------=" << threshold_gauss << std::endl;
822 
823  }
824 }
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
icdf_gauss
double icdf_gauss(double p)
Definition: xmipp_funcs.cpp:212
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
threshold
void threshold(double *phi, unsigned long nvox, double limit)
Definition: lib_vwk.cpp:524
j
#define j
Definition: numerical_recipes.cpp:2493
distance
TYPE distance(struct Point_T *p, struct Point_T *q)
Definition: point.cpp:28
ProgResDir::trigProducts
Matrix2D< double > trigProducts
Definition: resolution_directional.h:131
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
PI
#define PI
Definition: tools.h:43

◆ resolution2eval_()

void ProgResDir::resolution2eval_ ( int &  fourier_idx,
double  min_step,
double &  resolution,
double &  last_resolution,
int &  last_fourier_idx,
double &  freq,
double &  freqL,
double &  freqH,
bool &  continueIter,
bool &  breakIter,
bool &  doNextIteration 
)

Definition at line 624 of file resolution_directional.cpp.

629 {
630  int volsize = ZSIZE(VRiesz);
631 
632  FFT_IDX2DIGFREQ(fourier_idx, volsize, freq);
633 
634  resolution = sampling/freq;
635 // std::cout << "res = " << resolution << std::endl;
636 // std::cout << "min_step = " << min_step << std::endl;
637  if (resolution>8)
638  min_step =1;
639 
640 
641 
642  if ( fabs(resolution - last_resolution)<min_step )
643  {
644  freq = sampling/(last_resolution-min_step);
645  DIGFREQ2FFT_IDX(freq, volsize, fourier_idx);
646  FFT_IDX2DIGFREQ(fourier_idx, volsize, freq);
647 
648  if (fourier_idx == last_fourier_idx)
649  {
650  continueIter = true;
651  ++fourier_idx;
652  return;
653  }
654  }
655 
656  resolution = sampling/freq;
657  last_resolution = resolution;
658 
659  double step = 0.05*resolution;
660 
661  double resolution_L, resolution_H;
662 
663  if ( step < min_step)
664  {
665  resolution_L = resolution - min_step;
666  resolution_H = resolution + min_step;
667  }
668  else
669  {
670  resolution_L = 0.95*resolution;
671  resolution_H = 1.05*resolution;
672  }
673 
674  freqH = sampling/resolution_H;
675  freqL = sampling/resolution_L;
676 
677  if (freqH>0.5 || freqH<0)
678  freqH = 0.5;
679 
680  if (freqL>0.5 || freqL<0)
681  freqL = 0.5;
682  int fourier_idx_H, fourier_idx_L;
683 
684  DIGFREQ2FFT_IDX(freqH, volsize, fourier_idx_H);
685  DIGFREQ2FFT_IDX(freqL, volsize, fourier_idx_L);
686 
687  if (fourier_idx_H == fourier_idx)
688  fourier_idx_H = fourier_idx - 1;
689 
690  if (fourier_idx_L == fourier_idx)
691  fourier_idx_L = fourier_idx + 1;
692 
693  FFT_IDX2DIGFREQ(fourier_idx_H, volsize, freqH);
694  FFT_IDX2DIGFREQ(fourier_idx_L, volsize, freqL);
695 
696 // std::cout << "freq_H = " << freqH << std::endl;
697 // std::cout << "freq_L = " << freqL << std::endl;
698 
699  if (freq>0.49 || freq<0)
700  {
701  std::cout << "Nyquist limit reached" << std::endl;
702  breakIter = true;
703  doNextIteration = false;
704  return;
705  }
706  else
707  {
708  breakIter = false;
709  doNextIteration = true;
710  }
711 // std::cout << "resolution = " << resolution << " resolutionL = " <<
712 // sampling/(freqL) << " resolutionH = " << sampling/freqH
713 // << " las_res = " << last_resolution << std::endl;
714  last_fourier_idx = fourier_idx;
715  ++fourier_idx;
716 }
DIGFREQ2FFT_IDX
#define DIGFREQ2FFT_IDX(freq, size, idx)
Definition: xmipp_fft.h:61
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
FFT_IDX2DIGFREQ
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
Definition: cuda_gpu_reconstruct_fourier.cpp:382

◆ run()

void ProgResDir::run ( )
virtual

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

Reimplemented from XmippProgram.

Definition at line 1533 of file resolution_directional.cpp.

1534 {
1535  produceSideInfo();
1536  bool continueIter = false, breakIter = false;
1537  double criticalZ=icdf_gauss(significance);
1538 
1539  double step;
1540  step = res_step;
1541 
1542  std::cout << "Analyzing directions " << std::endl;
1543 
1544  double w = 0.0;
1545  double wH = 0.0;
1546  int volsize = ZSIZE(VRiesz);
1547 
1548  //Checking with MonoRes at 50A;
1549  int aux_idx;
1550 
1551  if (maxRes>18)
1552  {
1553  DIGFREQ2FFT_IDX(sampling/18, volsize, aux_idx);
1554 
1555  FFT_IDX2DIGFREQ(aux_idx, volsize, w);
1556  FFT_IDX2DIGFREQ(aux_idx+1, volsize, wH); //Frequency chosen for a first estimation
1557  }
1558  else
1559  {
1560  FFT_IDX2DIGFREQ(3, volsize, w);
1561  FFT_IDX2DIGFREQ(4, volsize, w);
1562  aux_idx = 3;
1563  }
1564  //std::cout << "fourier idx = " << aux_idx << std::endl;
1565  //std::cout << "Calling MonoRes core as a first estimation at " << sampling/w << "A." << std::endl;
1566 
1567  MultidimArray<double> amplitudeMS;
1568  double AvgNoise;
1569  AvgNoise = firstMonoResEstimation(fftV, w, wH, amplitudeMS)/9.0;
1570 
1571  N_directions=angles.mdimx;
1572 
1573  std::cout << "N_directions = " << N_directions << std::endl;
1574 
1575  double cone_angle = 45.0; //(degrees)
1576  cone_angle = PI*cone_angle/180;
1577 
1578  trigProducts.initZeros(3, N_directions);
1579 
1580  Image<double> outputResolution;
1581 
1582  for (size_t dir=0; dir<N_directions; dir++)
1583  {
1584  outputResolution().initZeros(VRiesz);
1585 // MultidimArray<double> &pOutputResolution = outputResolution();
1586  double freq, freqL, freqH, counter, resolution_2;
1587  MultidimArray<int> mask_aux = mask();
1588  MultidimArray<int> &pMask = mask_aux;
1589  std::vector<double> list;
1590  double resolution; //A huge value for achieving last_resolution < resolution
1591 
1592  double max_meanS = -1e38;
1593  double cut_value = 0.025;
1594 
1595  bool doNextIteration=true;
1596 
1597  int fourier_idx, last_fourier_idx = -1, iter = 0, fourier_idx_2;
1598  fourier_idx = aux_idx;
1599  int count_res = 0;
1600  double rot = MAT_ELEM(angles, 0, dir);
1601  double tilt = MAT_ELEM(angles, 1, dir);
1602  MAT_ELEM(trigProducts, 0, dir) = sin(tilt*PI/180)*cos(rot*PI/180);
1603  MAT_ELEM(trigProducts, 1, dir) = sin(tilt*PI/180)*sin(rot*PI/180);
1604  MAT_ELEM(trigProducts, 2, dir) = cos(tilt*PI/180);
1605  std::cout << "--------------NEW DIRECTION--------------" << std::endl;
1606  std::cout << "direction = " << dir+1 << " rot = " << rot << " tilt = " << tilt << std::endl;
1607 
1608 
1609  std::vector<float> noiseValues;
1610  FileName fnDebug;
1611  double last_resolution = 0;
1612 
1613  defineCone(fftV, conefilter, rot, tilt);
1614  maskMatrix.initConstant(1, NVoxelsOriginalMask, 1);
1615  do
1616  {
1617  continueIter = false;
1618  breakIter = false;
1619  //std::cout << "--------------Frequency--------------" << std::endl;
1620 
1621  resolution2eval_(fourier_idx, step,
1622  resolution, last_resolution, last_fourier_idx,
1623  freq, freqL, freqH,
1624  continueIter, breakIter, doNextIteration);
1625 
1626  if (breakIter)
1627  break;
1628 
1629  if (continueIter)
1630  continue;
1631 
1632  list.push_back(resolution);
1633 
1634  if (iter<2)
1635  resolution_2 = list[0];
1636  else
1637  resolution_2 = list[iter - 2];
1638 
1639  fnDebug = "Signal";
1640 
1641  amplitudeMonogenicSignal3D_fast(conefilter, freq, freqH, freqL, amplitudeMS, iter, dir, fnDebug);
1642 
1643  double sumS=0, sumS2=0, sumN=0, sumN2=0, NN = 0, NS = 0;
1644  noiseValues.clear();
1645 
1646 
1647  double x_dir = sin(tilt*PI/180)*cos(rot*PI/180);
1648  double y_dir = sin(tilt*PI/180)*sin(rot*PI/180);
1649  double z_dir = cos(tilt*PI/180);
1650 
1651  double uz, uy, ux;
1652 
1653  int n=0;
1654  int z_size = ZSIZE(amplitudeMS);
1655  int x_size = XSIZE(amplitudeMS);
1656  int y_size = YSIZE(amplitudeMS);
1657 
1658  size_t idx_mask;
1659  idx_mask = 0;
1660 
1661  double amplitudeValue;
1662 
1663  for(int k=0; k<z_size; ++k)
1664  {
1665  for(int i=0; i<y_size; ++i)
1666  {
1667  for(int j=0; j<x_size; ++j)
1668  {
1669  if (DIRECT_MULTIDIM_ELEM(pMask, n)>=1)
1670  {
1671  if (MAT_ELEM(maskMatrix, 0, idx_mask) >0)
1672  {
1673  amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitudeMS, n);
1674  sumS += amplitudeValue;
1675  ++NS;
1676  }
1677  ++idx_mask;
1678 
1679  }
1680  else if (DIRECT_MULTIDIM_ELEM(pMask, n)==0)
1681  {
1682  uz = (k - z_size*0.5);
1683  ux = (j - x_size*0.5);
1684  uy = (i - y_size*0.5);
1685 
1686  double rad = sqrt(ux*ux + uy*uy + uz*uz);
1687  double iun = 1/rad;
1688 
1689  //BE CAREFULL with the order
1690  double dotproduct = (uy*y_dir + ux*x_dir + uz*z_dir)*iun;
1691 
1692  double acosine = acos(dotproduct);
1693 
1694  //TODO: change efficiency the if condition by means of fabs(cos(angle))
1695  if (((acosine<cone_angle) || (acosine>(PI-cone_angle)) )
1696  && (rad>Rparticle))
1697  {
1698 // DIRECT_MULTIDIM_ELEM(coneVol, n) = 1;
1699  amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitudeMS, n);
1700  noiseValues.push_back((float) amplitudeValue);
1701  sumN += amplitudeValue;
1702  sumN2 += amplitudeValue*amplitudeValue;
1703  ++NN;
1704  }
1705  }
1706  ++n;
1707  }
1708  }
1709  }
1710 
1711  #ifdef DEBUG_DIR
1712  if (iter == 0)
1713  {
1714  Image<double> img;
1715 
1716  FileName iternumber;
1717  iternumber = formatString("cone_noise_%i_%i.vol", dir, iter);
1718  img = coneVol;
1719  img.write(iternumber);
1720  }
1721  #endif
1722 
1723  if ( (NS/(double) NVoxelsOriginalMask)<cut_value ) //when the 2.5% is reached then the iterative process stops
1724  {
1725  std::cout << "Search of resolutions stopped due to mask has been completed" << std::endl;
1726  doNextIteration =false;
1727  Nvoxels = 0;
1728 
1729  #ifdef DEBUG_MASK
1730  mask.write("partial_mask.vol");
1731  #endif
1732  }
1733  else
1734  {
1735  if (NS == 0)
1736  {
1737  std::cout << "There are no points to compute inside the mask" << std::endl;
1738  std::cout << "If the number of computed frequencies is low, perhaps the provided"
1739  "mask is not enough tight to the volume, in that case please try another mask" << std::endl;
1740  break;
1741  }
1742 
1743  double meanS=sumS/NS;
1744  // double sigma2S=sumS2/NS-meanS*meanS;
1745  double meanN=sumN/NN;
1746  double sigma2N=sumN2/NN-meanN*meanN;
1747 
1748  if (meanS>max_meanS)
1749  max_meanS = meanS;
1750 
1751  if (meanS<0.001*AvgNoise)//0001*max_meanS)
1752  {
1753  //std::cout << " meanS= " << meanS << " sigma2S= " << sigma2S << " NS = " << NS << std::endl;
1754  //std::cout << " meanN= " << meanN << " sigma2N= " << sigma2N << " NN= " << NN << std::endl;
1755  std::cout << "Search of resolutions stopped due to too low signal" << std::endl;
1756  std::cout << "\n"<< std::endl;
1757  doNextIteration = false;
1758  }
1759  else
1760  {
1761  // Check local resolution
1762  double thresholdNoise;
1763  //thresholdNoise = meanN+criticalZ*sqrt(sigma2N);
1764 
1765  std::sort(noiseValues.begin(),noiseValues.end());
1766  thresholdNoise = (double) noiseValues[size_t(noiseValues.size()*significance)];
1767 
1768  //std::cout << "thr="<< thresholdNoise << " " << meanN+criticalZ*sqrt(sigma2N) << " " << NN << std::endl;
1769  noiseValues.clear();
1770 
1771  std::cout << "Iteration = " << iter << ", Resolution= " << resolution << ", Signal = " << meanS << ", Noise = " << meanN << ", Threshold = " << thresholdNoise <<std::endl;
1772 
1773 
1774  size_t maskPos = 0;
1775  FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitudeMS)
1776  {
1777  if (DIRECT_MULTIDIM_ELEM(pMask, n)>=1)
1778  {
1779  if (MAT_ELEM(maskMatrix, 0, maskPos) >=1)
1780  {
1781  if (DIRECT_MULTIDIM_ELEM(amplitudeMS, n)>thresholdNoise)
1782  {
1783  MAT_ELEM(resolutionMatrix, dir, maskPos) = resolution;
1784  MAT_ELEM(maskMatrix, 0, maskPos) = 1;
1785  }
1786  else
1787  {
1788  MAT_ELEM(maskMatrix, 0, maskPos) += 1;
1789  if (MAT_ELEM(maskMatrix, 0, maskPos) >2)
1790  {
1791  MAT_ELEM(maskMatrix, 0, maskPos) = 0;
1792  MAT_ELEM(resolutionMatrix, dir, maskPos) = resolution_2;
1793  }
1794  }
1795  }
1796  ++maskPos;
1797  }
1798  }
1799 
1800  if (doNextIteration)
1801  if (resolution <= (minRes-0.001))
1802  doNextIteration = false;
1803  }
1804  }
1805  ++iter;
1806  last_resolution = resolution;
1807  }while(doNextIteration);
1808 
1809 // amplitudeMS.clear();
1810 // fftVRiesz.clear();
1811 
1812 // size_t maskPos=0;
1813 // Image<double> ResolutionVol;
1814 // MultidimArray<double> &pResolutionVol = ResolutionVol();
1815 //
1816 // pResolutionVol.initZeros(amplitudeMS);
1817 // FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(pResolutionVol)
1818 // {
1819 // if (DIRECT_MULTIDIM_ELEM(mask(), n) == 1)
1820 // {
1821 // double myres = MAT_ELEM(resolutionMatrix, dir, maskPos);
1822 // DIRECT_MULTIDIM_ELEM(pResolutionVol, n) = myres;
1825 // ++maskPos;
1826 // }
1827 // }
1828 // #ifdef DEBUG_DIR
1829 // Image<double> saveImg;
1830 // saveImg = pResolutionVol;
1831 // FileName fnres = formatString("resolution_dir_%i.vol", dir+1);
1832 // saveImg.write(fnres);
1833 // saveImg.clear();
1834 // #endif
1835 // pResolutionVol.clear();
1836  list.clear();
1837 
1838  std::cout << "----------------direction-finished----------------" << std::endl;
1839  }
1840 
1841 
1842 
1843 
1845 
1846  int maskPos = 0;
1847 
1848  //Remove outliers
1849  removeOutliers(resolutionMatrix);
1850 
1851  MultidimArray<double> radial, azimuthal, lowestResolution, highestResolution, doavol1, doavol2;
1852  MetaDataVec prefDir;
1853 
1854  double radialThr, azimuthalThr;
1855  radialAzimuthalResolution(resolutionMatrix, mask(), radial, azimuthal,
1856  lowestResolution, highestResolution, doavol1, doavol2, radialThr, azimuthalThr, prefDir);
1857 
1858  Image<double> imgdoa;
1859  imgdoa = radial;
1860  imgdoa.write(fnradial);
1861  imgdoa = azimuthal;
1862  imgdoa.write(fnazimuthal);
1863  imgdoa = lowestResolution;
1864  imgdoa.write(fnLowestResolution);
1865  imgdoa = highestResolution;
1866  imgdoa.write(fnHighestResolution);
1867  imgdoa = doavol1;
1868  imgdoa.write(fnDoa1);
1869  imgdoa = doavol2;
1870  imgdoa.write(fnDoa2);
1871 
1872  MetaDataVec mdRadialAzimuthalThr;
1873  size_t objIdx;
1874  objIdx = mdRadialAzimuthalThr.addObject();
1875  mdRadialAzimuthalThr.setValue(MDL_RESOLUTION_FREQ, radialThr, objIdx);
1876  mdRadialAzimuthalThr.setValue(MDL_RESOLUTION_FREQ2, azimuthalThr, objIdx);
1877 
1878  mdRadialAzimuthalThr.write(fnMDThr);
1879 
1880  //std::cout << "radial = " << radialThr << " azimuthal = " << azimuthalThr << std::endl;
1881  std::cout << "Calculating the radial and azimuthal resolution " << std::endl;
1882 
1883 
1884  MetaDataVec mdRadial, mdAvg, mdHighest, mdLowest;
1885 
1886  Image<double> monores;
1887  monores.read(fnMonoRes);
1888  MultidimArray<double> monoresVol;
1889  monoresVol = monores();
1890  radialAverageInMask(mask(), radial, azimuthal, highestResolution, lowestResolution, monoresVol, mdAvg);
1891  mdAvg.write(fnRadialAvG);
1892 }
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
Matrix2D::initConstant
void initConstant(T val)
Definition: matrix2d.h:602
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
ProgResDir::defineCone
void defineCone(MultidimArray< std::complex< double > > &myfftV, MultidimArray< std::complex< double > > &conefilter, double rot, double tilt)
Definition: resolution_directional.cpp:545
ImageBase::write
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)
Definition: xmipp_image_base.cpp:292
icdf_gauss
double icdf_gauss(double p)
Definition: xmipp_funcs.cpp:212
w
doublereal * w
Definition: numerical_recipes.cpp:2477
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
ProgResDir::radialAzimuthalResolution
void radialAzimuthalResolution(Matrix2D< double > &resolutionMat, MultidimArray< int > &pmask, MultidimArray< double > &radial, MultidimArray< double > &azimuthal, MultidimArray< double > &lowestResolution, MultidimArray< double > &highestResolution, MultidimArray< double > &doaResolution_1, MultidimArray< double > &doaResolution_2, double &radial_Thr, double &azimuthal_Thr, MetaDataVec &mdprefDirs)
Definition: resolution_directional.cpp:1182
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
DIGFREQ2FFT_IDX
#define DIGFREQ2FFT_IDX(freq, size, idx)
Definition: xmipp_fft.h:61
ProgResDir::angles
Matrix2D< double > angles
Definition: resolution_directional.h:131
ProgResDir::firstMonoResEstimation
double firstMonoResEstimation(MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, MultidimArray< double > &amplitude)
Definition: resolution_directional.cpp:1365
ProgResDir::resolutionMatrix
Matrix2D< double > resolutionMatrix
Definition: resolution_directional.h:131
i
#define i
Definition: numerical_recipes.cpp:2493
k
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
ProgResDir::radialAverageInMask
void radialAverageInMask(MultidimArray< int > &mask, MultidimArray< double > &inputVol_1, MultidimArray< double > &inputVol_2, MultidimArray< double > &inputVol_3, MultidimArray< double > &inputVol_4, MultidimArray< double > &inputVol_5, MetaDataVec &md)
Definition: resolution_directional.cpp:1010
ProgResDir::fnLowestResolution
FileName fnLowestResolution
Definition: resolution_directional.h:55
MAT_ELEM
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
ProgResDir::resolution2eval_
void resolution2eval_(int &fourier_idx, double min_step, double &resolution, double &last_resolution, int &last_fourier_idx, double &freq, double &freqL, double &freqH, bool &continueIter, bool &breakIter, bool &doNextIteration)
Definition: resolution_directional.cpp:624
MDL_RESOLUTION_FREQ2
Frequency in 1/A squared (double)
Definition: metadata_label.h:389
ProgResDir::amplitudeMonogenicSignal3D_fast
void amplitudeMonogenicSignal3D_fast(const MultidimArray< std::complex< double > > &myfftV, double w1, double w1l, double wH, MultidimArray< double > &amplitude, int count, int dir, FileName fnDebug)
Definition: resolution_directional.cpp:347
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
MetaDataVec::addObject
size_t addObject() override
Definition: metadata_vec.cpp:435
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
Definition: multidim_array_base.h:176
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
DIRECT_MULTIDIM_ELEM
#define DIRECT_MULTIDIM_ELEM(v, n)
Definition: multidim_array_base.h:161
ProgResDir::VRiesz
MultidimArray< double > VRiesz
Definition: resolution_directional.h:128
sort
void sort(struct DCEL_T *dcel)
Definition: sorting.cpp:18
FFT_IDX2DIGFREQ
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
Definition: cuda_gpu_reconstruct_fourier.cpp:382
j
#define j
Definition: numerical_recipes.cpp:2493
MDL_RESOLUTION_FREQ
Frequency in 1/A (double)
Definition: metadata_label.h:388
ProgResDir::fnHighestResolution
FileName fnHighestResolution
Definition: resolution_directional.h:55
ProgResDir::trigProducts
Matrix2D< double > trigProducts
Definition: resolution_directional.h:131
ProgResDir::conefilter
MultidimArray< std::complex< double > > conefilter
Definition: resolution_directional.h:130
formatString
String formatString(const char *format,...)
Definition: xmipp_strings.cpp:602
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
ProgResDir::removeOutliers
void removeOutliers(Matrix2D< double > &resolutionMat)
Definition: resolution_directional.cpp:719
ImageBase::read
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
Definition: xmipp_image_base.cpp:119
Matrix2D::mdimx
size_t mdimx
Definition: matrix2d.h:410
ProgResDir::maskMatrix
Matrix2D< double > maskMatrix
Definition: resolution_directional.h:131
ProgResDir::fftV
MultidimArray< std::complex< double > > fftV
Definition: resolution_directional.h:130
PI
#define PI
Definition: tools.h:43
n
int * n
Definition: numerical_recipes.cpp:2229

Member Data Documentation

◆ angles

Matrix2D<double> ProgResDir::angles

Definition at line 131 of file resolution_directional.h.

◆ checkellipsoids

bool ProgResDir::checkellipsoids

Analyze radial and azimuthal resolutoin

Definition at line 69 of file resolution_directional.h.

◆ conefilter

MultidimArray< std::complex<double> > ProgResDir::conefilter

Definition at line 130 of file resolution_directional.h.

◆ fastCompute

bool ProgResDir::fastCompute

Definition at line 69 of file resolution_directional.h.

◆ fftV

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

Definition at line 130 of file resolution_directional.h.

◆ fftVRiesz

MultidimArray< std::complex<double> > ProgResDir::fftVRiesz

Definition at line 127 of file resolution_directional.h.

◆ fftVRiesz_aux

MultidimArray< std::complex<double> > ProgResDir::fftVRiesz_aux

Definition at line 127 of file resolution_directional.h.

◆ fnazimuthal

FileName ProgResDir::fnazimuthal

Definition at line 55 of file resolution_directional.h.

◆ fnDoa1

FileName ProgResDir::fnDoa1

Definition at line 55 of file resolution_directional.h.

◆ fnDoa2

FileName ProgResDir::fnDoa2

Definition at line 55 of file resolution_directional.h.

◆ fnHighestResolution

FileName ProgResDir::fnHighestResolution

Definition at line 55 of file resolution_directional.h.

◆ fnLowestResolution

FileName ProgResDir::fnLowestResolution

Definition at line 55 of file resolution_directional.h.

◆ fnMask

FileName ProgResDir::fnMask

Definition at line 55 of file resolution_directional.h.

◆ fnMDThr

FileName ProgResDir::fnMDThr

Definition at line 55 of file resolution_directional.h.

◆ fnMonoRes

FileName ProgResDir::fnMonoRes

Definition at line 55 of file resolution_directional.h.

◆ fnOut

FileName ProgResDir::fnOut

Filenames

Definition at line 55 of file resolution_directional.h.

◆ fnprefMin

FileName ProgResDir::fnprefMin

Definition at line 55 of file resolution_directional.h.

◆ fnradial

FileName ProgResDir::fnradial

Definition at line 55 of file resolution_directional.h.

◆ fnRadialAvG

FileName ProgResDir::fnRadialAvG

Definition at line 55 of file resolution_directional.h.

◆ fnVol

FileName ProgResDir::fnVol

Definition at line 55 of file resolution_directional.h.

◆ fnZscore

FileName ProgResDir::fnZscore

Definition at line 55 of file resolution_directional.h.

◆ freq_fourier

Matrix1D<double> ProgResDir::freq_fourier

Definition at line 132 of file resolution_directional.h.

◆ iu

MultidimArray<double> ProgResDir::iu

Definition at line 128 of file resolution_directional.h.

◆ mask

Image<int> ProgResDir::mask

Definition at line 133 of file resolution_directional.h.

◆ maskMatrix

Matrix2D<double> ProgResDir::maskMatrix

Definition at line 131 of file resolution_directional.h.

◆ maxRes

double ProgResDir::maxRes

Definition at line 60 of file resolution_directional.h.

◆ minRes

double ProgResDir::minRes

Definition at line 60 of file resolution_directional.h.

◆ N_directions

double ProgResDir::N_directions

Definition at line 60 of file resolution_directional.h.

◆ N_freq

double ProgResDir::N_freq

Step in digital frequency

Definition at line 66 of file resolution_directional.h.

◆ N_smoothing

int ProgResDir::N_smoothing

Definition at line 134 of file resolution_directional.h.

◆ Nthr

int ProgResDir::Nthr

Definition at line 63 of file resolution_directional.h.

◆ Nvoxels

int ProgResDir::Nvoxels

Definition at line 63 of file resolution_directional.h.

◆ NVoxelsOriginalMask

int ProgResDir::NVoxelsOriginalMask

Is the volume previously masked?

Definition at line 63 of file resolution_directional.h.

◆ R

double ProgResDir::R

Definition at line 60 of file resolution_directional.h.

◆ res_step

double ProgResDir::res_step

Definition at line 60 of file resolution_directional.h.

◆ resolutionMatrix

Matrix2D<double> ProgResDir::resolutionMatrix

Definition at line 131 of file resolution_directional.h.

◆ Rparticle

double ProgResDir::Rparticle

Definition at line 60 of file resolution_directional.h.

◆ sampling

double ProgResDir::sampling

sampling rate, minimum resolution, and maximum resolution

Definition at line 60 of file resolution_directional.h.

◆ significance

double ProgResDir::significance

Definition at line 66 of file resolution_directional.h.

◆ transformer_inv

FourierTransformer ProgResDir::transformer_inv

Definition at line 129 of file resolution_directional.h.

◆ trigProducts

Matrix2D<double> ProgResDir::trigProducts

Definition at line 131 of file resolution_directional.h.

◆ VRiesz

MultidimArray<double> ProgResDir::VRiesz

Definition at line 128 of file resolution_directional.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.13