#include <resolution_fso.h>

Inheritance diagram for ProgFSO:

Collaboration diagram for ProgFSO:

Public Member Functions
void	defineParams ()

void	readParams ()

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 ()

Additional Inherited Members
Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

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

Definition at line 38 of file resolution_fso.h.

Member Function Documentation

◆ defineParams()

void ProgFSO::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 36 of file resolution_fso.cpp.

 {
     addUsageLine("Calculate Fourier Shell Occupancy - FSO curve - via directional FSC measurements.");
     addUsageLine("Following outputs are generated:");
     addUsageLine("  1) FSO curve");
     addUsageLine("  2) Global resolution from FSO and FSC");
     addUsageLine("  3) 3DFSC");
     addUsageLine("  4) Anisotropic filter");
     addUsageLine("Reference: J.L. Vilas, H.D. Tagare, XXXXX (2021)");
     addUsageLine("+* Fourier Shell Occupancy (FSO)", true);
     addUsageLine("+ The Fourier Shell Occupancy Curve can be obtained from a set of directional FSC (see below).");
     addUsageLine("+ To do that, two half maps are used to determine the Global FSC at threshold 0.143. Then, the ratio between the number");
     addUsageLine("+ of directions with resolution higher (better) than the Global resolution and the total number of measured directions is");
     addUsageLine("+ calculated at different frequencies (resolutions). Note that this ratio is between 0 (resolution of all directions is worse than the global FSC)");
     addUsageLine("+ resolution than the global FSC)  and 1 (all directions present better resolution than the FSC) at a given resolution.");
     addUsageLine("+ In the particular case, FSO curve takes the value of 0.5 (FSO=0.5), then half of the directions are better, and.");
     addUsageLine("+ the other half are worse than the FSC, this situation occurs at the resoltuion of the map. It means the FSO = 0.5 at the ");
     addUsageLine("+ FSC resolution. A map is isotropic if all directional resolution are similar, and anisotropic is there are significant resolution values along");
     addUsageLine("+ different directions. Thus, when the FSO presents a sharp cliff, a step-like function, the map will be isotropic.");
     addUsageLine("+ In contrast, when the OFSC shows a slope the map will be anisotropic. The lesser slope the higher resolution isotropy.");
     addUsageLine("+ ");
     addUsageLine("+* Directional Fourier Shell Correlation (dFSC)", true);
     addUsageLine("+ This program estimates the directional FSC between two half maps along all posible directions on the projection sphere.");
     addUsageLine("+ The directionality is measured by means of conical-like filters in Fourier Space. To avoid possible Gibbs effects ");
     addUsageLine("+ the filters are gaussian functions with their respective maxima along the filtering direction. A set of 321 directions ");
     addUsageLine("+ is used to cover the projection sphere, computing for each direction the directional FSC at 0.143 between the two half maps.");
     addUsageLine("+ The result is a set of 321 FSC curves (321 is the number of analyzed directions, densely covering the projection sphere).");
     addUsageLine("+ The 3DFSC is then obtained from all curves by interpolation. Note that as well as it occurs with global FSC, the directional FSC is mask dependent.");
     addUsageLine(" ");
     addUsageLine("+* Resolution Distribution and 3DFSC", true);
     addUsageLine("+ The directional-FSC, dFSC is estimated along 321 directions on the projection sphere. For each direction the corresponding");
     addUsageLine("+ resolution is determined. Thus, it is possible to determine the resolution distribution on the projection sphere.");
     addUsageLine("+ This distribution is saved in the output metadata named resolutionDistribution.xmd. Also by means of all dFSC, the 3DFSC");
     addUsageLine("+ is calculated and saved as 3dFSC.mrc, which gives an idea about the information distributionin Fourier Space.");
     addUsageLine(" ");
     addUsageLine(" ");
     addSeeAlsoLine("resolution_fsc");
 
     addParamsLine("   --half1 <input_file>               : Input Half map 1");
     addParamsLine("   --half2 <input_file>               : Input Half map 2");
 
     addParamsLine("   [-o <output_folder=\"\">]          : Folder where the results will be stored.");
 
     addParamsLine("   [--sampling <Ts=1>]                : (Optical) Pixel size (Angstrom). If it is not provided by default will be 1 A/px.");
     addParamsLine("   [--mask <input_file=\"\">]         : (Optional) Smooth mask to remove noise. If it is not provided, the computation will be carried out without mask.");
 
     addParamsLine("   [--anglecone <ang_con=17>]         : (Optional) Angle Cone (angle between the axis and the  generatrix) for estimating the directional FSC");
     addParamsLine("   [--threshold <thrs=0.143>]         : (Optional) Threshold for the FSC/directionalFSC estimation ");
 
     addParamsLine("   [--threedfsc_filter]               : (Optional) Put this flag to estimate the 3DFSC, and apply it as low pass filter to obtain a directionally filtered map. It mean to apply an anisotropic filter.");
     
     addParamsLine("   [--threads <Nthreads=1>]           : (Optional) Number of threads to be used");
 
     addExampleLine("Resolution of two half maps half1.mrc and half2.mrc with a sampling rate of 2 A/px", false);
     addExampleLine("xmipp_resolution_fso --half1 half1.mrc --half2 half2.mrc --sampling_rate 2 ");
     addExampleLine("Resolution of two half maps half1.mrc and half2.mrc with a sampling rate of 2 A/px and a mask mask.mrc", false);
     addExampleLine("xmipp_resolution_fso --half1 half1.mrc --half2 half2.mrc --mask mask.mrc --sampling_rate 2 ");
 }

◆ readParams()

void ProgFSO::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 95 of file resolution_fso.cpp.

 {
     fnhalf1 = getParam("--half1");
     fnhalf2 = getParam("--half2");
     fnOut = getParam("-o");
 
     sampling = getDoubleParam("--sampling");
     fnmask = getParam("--mask");
     ang_con = getDoubleParam("--anglecone");
     thrs = getDoubleParam("--threshold");
     do_3dfsc_filter = checkParam("--threedfsc_filter");
     
     Nthreads = getIntParam("--threads");
 }

◆ run()

void ProgFSO::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 1309 of file resolution_fso.cpp.

     {
         std::cout << "Starting ... " << std::endl << std::endl;
         
         MultidimArray<double> half1, half2;
         MultidimArray<double> &phalf1 = half1, &phalf2 = half2;
 
         //This read the data and applies a fftshift to in the next step compute the fft
         prepareData(half1, half2);
 
         //Computing the FFT
         FourierTransformer transformer2(FFTW_BACKWARD), transformer1(FFTW_BACKWARD);
         transformer1.setThreadsNumber(Nthreads);
         transformer2.setThreadsNumber(Nthreads);
 
         transformer1.FourierTransform(half1, FT1, false);
         transformer2.FourierTransform(half2, FT2, false);
 
         // Defining frequencies freq_fourier_x,y,z and freqMap
         // The number of frequencies in each shell freqElem is determined
         defineFrequencies(FT1, phalf1);
 
         half1.clear();
         half2.clear();
 
         // Storing the shell of both maps as vectors global
         // The global FSC is also computed
         MultidimArray<double> freq;
         arrangeFSC_and_fscGlobal(sampling, thrs, freq);
 
         FT2.clear();
 
         // Generating the set of directions to be analyzed
         // And converting the cone angle to radians
         generateDirections(angles, true);
         ang_con = ang_con*PI/180.0;
 
         // Preparing the metadata for storing the FSO
         MultidimArray<double> resDirFSC(angles.mdimx);//, aniParam;
         //aniParam.initZeros(xvoldim/2+1);
 
         // Computing directional FSC and 3DFSC
         // Create one copy for each thread
         std::vector<MultidimArray<float>> threeD_FSCs(Nthreads);
         std::vector<MultidimArray<float>> normalizationMaps(Nthreads);
         std::vector<MultidimArray<float>> aniParams(Nthreads);
         for (auto & ma : threeD_FSCs) {
             ma.initZeros(real_z1z2);
         }
         for (auto & ma : normalizationMaps) {
             ma.initZeros(real_z1z2);
         }
         for (auto & ma : aniParams) {
             ma.initZeros(xvoldim/2+1);
         }
 
         // Binham Test
         // We will have as many vectors as threads
         // Each thread considers a vector with the matrices of the Bingham test
         // instantiate vector object of type vector<Matrix2D<float>>
         std::vector<std::vector<Matrix2D<double>>> isotropyMatrices;
     
         // resize the vector to `freq.nzyxdim` elements of type vector<float>, each having size `C`
         isotropyMatrices.resize(Nthreads, std::vector<Matrix2D<double>>(freq.nzyxdim));
 
         // Each component of the vector is a 0-matrix (null matrix)
         for (auto & im : isotropyMatrices)
         {
             for (auto & m : im)
                 m.initZeros(3,3);
         }
 
         ctpl::thread_pool threadPool;
         threadPool.resize(Nthreads);
         auto futures = std::vector<std::future<void>>();
         futures.reserve(angles.mdimx);
         
         for (size_t k = 0; k<angles.mdimx; k++)
         {
             futures.emplace_back(threadPool.push(
                 [k, this, &resDirFSC, &threeD_FSCs, &normalizationMaps, &aniParams, &isotropyMatrices](int thrId){
                 float rot  = MAT_ELEM(angles, 0, k);
                 float tilt = MAT_ELEM(angles, 1, k);
 
                 double resInterp = -1;  
                 MultidimArray<float> fsc;
 
                 auto &threeD_FSC = threeD_FSCs.at(thrId);
                 auto &normalizationMap = normalizationMaps.at(thrId);
                 auto &freqMat = isotropyMatrices.at(thrId);
                 // Estimating the direction FSC along the direction given by rot and tilt
                 fscDir_fast(fsc, rot, tilt, threeD_FSC, normalizationMap, thrs, resInterp, freqMat, k);
 
                 printf ("Direction %zu/%zu -> %.2f A \n", k, angles.mdimx, resInterp);
 
                 dAi(resDirFSC, k) = resInterp;
                 
                 // Updating the FSO curve
                 anistropyParameter(fsc, aniParams.at(thrId), thrs);
             }));
         }
 
         // wait till done
         for (auto &f : futures) {
             f.get();
         }
 
         std::cout << "----- Directional resolution estimated -----" <<  std::endl <<  std::endl;
         std::cout << "Preparing results ..." <<  std::endl;
 
         // merge results from multiple threads
         for (size_t i = 1; i < aniParams.size(); ++i) {
             aniParams.at(0) += aniParams.at(i);
         }
 
         // For the Bingham Test
         // merge results from multiple threads
         for (size_t i = 0; i < freq.nzyxdim; ++i) 
         {
             for (size_t j = 1; j<Nthreads; ++j)
             {
                 isotropyMatrices.at(0)[i] += isotropyMatrices.at(j)[i];
             }
             isotropyMatrices.at(0)[i] /= aniParams.at(0)[i];
         }
 
         //auto &freqMat = isotropyMatrices.at(0);
         auto &aniParam = aniParams.at(0);
         MultidimArray<double> isotropyMatrix;
         isotropyMatrix.initZeros(freq.nzyxdim);
 
         for (size_t i = 0; i< freq.nzyxdim; ++i)
         {
             double trT2 = 0;
             if (aniParams.at(0)[i] >0)
             {
                 // Let us define T = 1/n*Sum(wi*xi*xi) => Tr(T^2) = x*x + y*y + z*z;
                 //This is the Bingham Test (1/2)(p-1)*(p-2)*n*Sum(Tr(T^2) - 1/p)
                 //std::cout << isotropyMatrices.at(0)[i] << aniParams.at(0)[i] << std::endl;
                 Matrix2D<double> T2;
                 Matrix2D<double> T;
                 T = isotropyMatrices.at(0)[i];//freqMat[i]/DIRECT_MULTIDIM_ELEM(aniParam, i);
                 T2 = T*T;
                 trT2 = T2.trace();
                 double pdim = 3;
                 dAi(isotropyMatrix, i) = 0.5*pdim*(pdim+2)*(2*aniParams.at(0)[i])*(trT2 - 1./pdim);
                 // The factor 2 is because we need to compute the point in the whole sphere, but currently
                 // we are measuing half of the sphere due to the symmetry of the problem
                 // S = 0.5 * p * (p+2) * n * (trace(T2) - 1/p); Where p is the
                 // dimension of the sphere (p=3) and n the number of points to
                 // determine if they are uniformly or not. T = x*x';
 
                 // The hypohtesis test considers p = 0.05.
                 // p=0.05; nu = 5; x = chi2inv(p,nu);
             }
         }
 
 
         // ANISOTROPY CURVE
         aniParams.at(0) /= (double) angles.mdimx;
         for (size_t k = 0; k<5; k++)
             DIRECT_MULTIDIM_ELEM(aniParams.at(0), k) = 1.0;
         MetaDataVec mdani;
         saveAnisotropyToMetadata(mdani, freq, aniParams.at(0), isotropyMatrix);
 
         
         // DIRECTIONAL RESOLUTION DISTRIBUTION ON THE PROJECTION SPHERE
         FileName fn;
         fn = fnOut+"/Resolution_Distribution.xmd";
         
         resolutionDistribution(resDirFSC, fn);
 
 
         // 3DFSC and ANISOTROPIC FILTER
         if (do_3dfsc_filter)
         {
             // HALF 3DFSC MAP
             MultidimArray<double> d3_FSCMap;
             MultidimArray<double> d3_aniFilter;
             d3_FSCMap.resizeNoCopy(FT1);
             d3_FSCMap.initConstant(0);
             d3_aniFilter.resizeNoCopy(FT1);
             d3_aniFilter.initConstant(0);
 
             // merge results from multiple threads
             for (size_t i = 1; i < Nthreads; ++i) {
                 threeD_FSCs.at(0) += threeD_FSCs.at(i);
                 normalizationMaps.at(0) += normalizationMaps.at(i);
             }
             auto &threeD_FSC = threeD_FSCs.at(0);
             auto &normalizationMap = normalizationMaps.at(0);
 
             FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(threeD_FSC)
             {
                     double value = DIRECT_MULTIDIM_ELEM(threeD_FSC, n) /= DIRECT_MULTIDIM_ELEM(normalizationMap, n);
                     if (std::isnan(value))
                         value = 1.0;
 
                     size_t idx = DIRECT_MULTIDIM_ELEM(arr2indx, n);
                     DIRECT_MULTIDIM_ELEM(d3_FSCMap, idx) = value;
                     
                     if (DIRECT_MULTIDIM_ELEM(threeD_FSC, n)> thrs)//&& (DIRECT_MULTIDIM_ELEM(aniFilter, n) <1))
                     {
                         DIRECT_MULTIDIM_ELEM(d3_aniFilter, idx) = 1; //;DIRECT_MULTIDIM_ELEM(aniFilter, n) = 1;
                     }
                     else
                     {
                         DIRECT_MULTIDIM_ELEM(d3_aniFilter, idx) = value;//DIRECT_MULTIDIM_ELEM(aniFilter, n);
                     }
                     
             }
 
             // This code fix the empty line line in Fourier space
             size_t auxVal;
             auxVal = YSIZE(d3_FSCMap)/2;
 
             size_t j = 0;
             for(size_t i=(auxVal+1); i<YSIZE(d3_FSCMap); ++i)
         {
                 for(size_t k=0; k<ZSIZE(d3_FSCMap); ++k)
                 {
                     DIRECT_A3D_ELEM(d3_FSCMap,k,i,0) = DIRECT_A3D_ELEM(d3_FSCMap,k,i,1);
                     DIRECT_A3D_ELEM(d3_aniFilter,k,i,0) = DIRECT_A3D_ELEM(d3_aniFilter,k,i,1);
                 }
             }
             size_t i=0;
             for(size_t k=0; k<ZSIZE(d3_FSCMap); ++k)
             {
                 DIRECT_A3D_ELEM(d3_FSCMap,k,0,0) = DIRECT_A3D_ELEM(d3_FSCMap,k,0, 1);
                 DIRECT_A3D_ELEM(d3_aniFilter,k,0,0) = DIRECT_A3D_ELEM(d3_aniFilter,k,0,1);
             }
 
 
             Image<double> img;
             img() = d3_FSCMap;
             img.write("threeD_FSC.mrc");
 
             double sigma = 3;
 
             realGaussianFilter(d3_aniFilter, sigma);
 
             // DIRECTIONAL FILTERED MAP
             MultidimArray<double> filteredMap;
             //directionalFilter(FT1, d3_aniFilter, filteredMap, xvoldim, yvoldim, zvoldim);
             directionalFilterHalves(FT1, d3_aniFilter);
 
             //FULL 3DFSC MAP
             fn = fnOut+"/3dFSC.mrc";
             createFullFourier(d3_FSCMap, fn, xvoldim, yvoldim, zvoldim);
         }
 
 
         
         std::cout << "-------------Finished-------------" << std::endl;
 }

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

xmipp/libraries/reconstruction/resolution_fso.h
xmipp/libraries/reconstruction/resolution_fso.cpp

Public Member Functions

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ defineParams()

◆ readParams()

◆ run()