ctf_estimate_from_micrograph.cpp File Reference

#include "ctf_estimate_from_micrograph.h"
#include "ctf_enhance_psd.h"
#include "core/xmipp_fftw.h"
#include "core/metadata_extension.h"
#include "core/transformations.h"
#include "core/xmipp_threads.h"
#include "core/xmipp_image_generic.h"
#include "data/basic_pca.h"
#include "data/normalize.h"
#include "data/numerical_tools.h"

Include dependency graph for ctf_estimate_from_micrograph.cpp:

Go to the source code of this file.

Classes
class	ThreadFastEstimateEnhancedPSDParams

Functions
void	threadFastEstimateEnhancedPSD (ThreadArgument &thArg)
void	fastEstimateEnhancedPSD (const FileName &fnMicrograph, double downsampling, MultidimArray< double > &enhancedPSD, int numberOfThreads)

Function Documentation

threadFastEstimateEnhancedPSD()

void threadFastEstimateEnhancedPSD

(

ThreadArgument &

thArg

)

Definition at line 867 of file ctf_estimate_from_micrograph.cpp.

{
    auto *args = (ThreadFastEstimateEnhancedPSDParams*) thArg.workClass;
    int Nthreads = thArg.getNumberOfThreads();
    int id = thArg.thread_id;
    ImageGeneric &I = *(args->I);
    const MultidimArrayGeneric& mI = I();
    size_t IXdim, IYdim, IZdim;
    I.getDimensions(IXdim, IYdim, IZdim);
    MultidimArray<double> &pieceSmoother = *(args->pieceSmoother);
    MultidimArray<int> &pieceMask = *(args->pieceMask);
    MultidimArray<double> localPSD, piece;
    MultidimArray<std::complex<double> > Periodogram;
    piece.initZeros(pieceMask);
    localPSD.initZeros(*(args->PSD));

    FourierTransformer transformer;
    transformer.setReal(piece);

    int pieceNumber = 0;
    int Nprocessed = 0;
    double pieceDim2 = XSIZE(piece) * XSIZE(piece);
    for (size_t i = 0; i < (IYdim - YSIZE(piece)); i+=YSIZE(piece))
        for (size_t j = 0; j < (IXdim - XSIZE(piece)); j+=XSIZE(piece), pieceNumber++)
        {
            if ((pieceNumber + 1) % Nthreads != id)
                continue;
            Nprocessed++;

            // Extract micrograph piece ..........................................
            for (size_t k = 0; k < YSIZE(piece); k++)
                for (size_t l = 0; l < XSIZE(piece); l++)
                    DIRECT_A2D_ELEM(piece, k, l)= mI(i+k, j+l);
            piece.statisticsAdjust(0., 1.);
            normalize_ramp(piece, &pieceMask);
            piece *= pieceSmoother;

            // Estimate the power spectrum .......................................
            transformer.FourierTransform();
            transformer.getCompleteFourier(Periodogram);
            FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(localPSD)
            {
                const auto *ptr = (double*) &DIRECT_MULTIDIM_ELEM(Periodogram, n);
                double re=*ptr;
                double im=*(ptr+1);
                double magnitude2=re*re+im*im;
                DIRECT_MULTIDIM_ELEM(localPSD,n)+=magnitude2*pieceDim2;
            }
        }

    // Gather results
    args->mutex->lock();
    args->Nprocessed += Nprocessed;
    *(args->PSD) += localPSD;
    args->mutex->unlock();
}