volume_initial_simulated_annealing.cpp File Reference

#include "volume_initial_simulated_annealing.h"
#include <data/filters.h>

Include dependency graph for volume_initial_simulated_annealing.cpp:

Go to the source code of this file.

Functions
void	alignSingleImage (size_t nImg, ProgVolumeInitialSimulatedAnnealing &prm, MetaData &mdReconstruction, double &newCorr, double &improvementFraction)
void	threadAlignSubset (ThreadArgument &thArg)

Function Documentation

alignSingleImage()

void alignSingleImage	(	size_t	nImg,
		ProgVolumeInitialSimulatedAnnealing &	prm,
		MetaData &	mdReconstruction,
		double &	newCorr,
		double &	improvementFraction
	)

Definition at line 90 of file volume_initial_simulated_annealing.cpp.

{
    MultidimArray<double> mGalleryProjection, mCurrentImage, mCurrentImageAligned;
    Matrix2D<double> M;

    mCurrentImage.aliasImageInStack(prm.inputImages(),nImg);
    mCurrentImage.setXmippOrigin();

    MultidimArray<double> allCorrs;
    allCorrs.resizeNoCopy(prm.mdGallery.size());

    size_t improvementCount=0;
    double oldCorr=prm.mdInp[nImg].maxcc;
#ifdef DEBUG
    FileName fnImg;
    prm.mdIn.getValue(MDL_IMAGE,fnImg,id);
    std::cout << "Image: " << fnImg << " oldCorr=" << oldCorr << std::endl;
#endif
    std::vector< Matrix2D<double> > allM;
    for (size_t nGallery=0; nGallery<XSIZE(allCorrs); ++nGallery)
    {
        mCurrentImageAligned=mCurrentImage;
        mGalleryProjection.aliasImageInStack(prm.gallery(),nGallery);
        mGalleryProjection.setXmippOrigin();
        double corr=alignImagesConsideringMirrors(mGalleryProjection,mCurrentImageAligned,M,xmipp_transformation::DONT_WRAP);
#ifdef DEBUG
        mdGallery.getValue(MDL_MAXCC,corr,__iter.objId);
#endif
        A1D_ELEM(allCorrs,nGallery)=corr;
        allM.push_back(M);
        if (corr>oldCorr)
        {
            ++improvementCount;
#ifdef DEBUG
        prm.mdGallery.getValue(MDL_IMAGE,fnImg,__iter.objId);
        std::cout << "   Matching Gallery: " << fnImg << " " << corr << std::endl;
#endif
        }
    }
    improvementFraction=(double)improvementCount/NSIZE(prm.gallery());
#ifdef DEBUG
        mdGallery.write("PPPgallery.xmd");
        std::cout << "   oldcorr=" << oldCorr << std::endl;
#endif

    MultidimArray<double> scaledCorrs;
    allCorrs.cumlativeDensityFunction(scaledCorrs);
    double bestCorr, bestAngleRot, bestAngleTilt, bestAnglePsi, bestShiftX, bestShiftY, bestWeight;
    bool bestFlip;
    bestCorr=-1;
    double correctionFactor=1;
    double currentT=prm.T0*(1.0-((double)prm.iter)/prm.NiterRandom);
    double icurrentT=1.0/currentT;
    for (size_t nGallery=0; nGallery<XSIZE(allCorrs); ++nGallery)
    {
        bool getThis=(A1D_ELEM(allCorrs,nGallery)>oldCorr) || (improvementCount==0 && A1D_ELEM(scaledCorrs,nGallery)>=0.98);
        if (A1D_ELEM(allCorrs,nGallery)<oldCorr && !getThis && prm.iter<=prm.NiterRandom)
        {
            // Simulated annealing
            double diff=oldCorr-A1D_ELEM(allCorrs,nGallery);
            double p=rnd_unif();
            getThis=(p<exp(-diff*icurrentT));
        }
        if (getThis)
        {
            bool flip;
            double scale, shiftX, shiftY, anglePsi;
            transformationMatrix2Parameters2D(allM[nGallery],flip,scale,shiftX,shiftY,anglePsi);
            anglePsi*=-1;
            double weight=A1D_ELEM(scaledCorrs,nGallery)*correctionFactor;
            double corr=A1D_ELEM(allCorrs,nGallery);
            double angleRot=prm.mdGallery[nGallery].rot;
            double angleTilt=prm.mdGallery[nGallery].tilt;
#ifdef DEBUG
            fnImg=prm.mdGallery[nGallery].fnImg;
            std::cout << "   Getting Gallery: " << fnImg << " corr=" << corr << " cdf=" << weight << " rot=" << angleRot
                      << " tilt=" << angleTilt << std::endl;
#endif

            if (prm.iter<=prm.NiterRandom)
            {
                size_t recId=mdReconstruction.addObject();
                FileName fnImg;
                fnImg=prm.mdInp[nImg].fnImg;
                mdReconstruction.setValue(MDL_IMAGE,fnImg,recId);
                mdReconstruction.setValue(MDL_ENABLED,1,recId);
                mdReconstruction.setValue(MDL_MAXCC,corr,recId);
                mdReconstruction.setValue(MDL_ANGLE_ROT,angleRot,recId);
                mdReconstruction.setValue(MDL_ANGLE_TILT,angleTilt,recId);
                mdReconstruction.setValue(MDL_ANGLE_PSI,anglePsi,recId);
                mdReconstruction.setValue(MDL_SHIFT_X,shiftX,recId);
                mdReconstruction.setValue(MDL_SHIFT_Y,shiftY,recId);
                mdReconstruction.setValue(MDL_FLIP,flip,recId);
                mdReconstruction.setValue(MDL_WEIGHT,weight,recId);
            }
            if (corr>bestCorr)
            {
                bestCorr=corr;
                bestAngleRot=angleRot;
                bestAngleTilt=angleTilt;
                bestAnglePsi=anglePsi;
                bestShiftX=shiftX;
                bestShiftY=shiftY;
                bestFlip=flip;
                bestWeight=weight;
            }
        }

        ++nGallery;
    }

    if (prm.iter>prm.NiterRandom)
    {
        size_t recId=mdReconstruction.addObject();
        FileName fnImg;
        fnImg=prm.mdInp[nImg].fnImg;
        mdReconstruction.setValue(MDL_IMAGE,fnImg,recId);
        mdReconstruction.setValue(MDL_ENABLED,1,recId);
        mdReconstruction.setValue(MDL_MAXCC,bestCorr,recId);
        mdReconstruction.setValue(MDL_ANGLE_ROT,bestAngleRot,recId);
        mdReconstruction.setValue(MDL_ANGLE_TILT,bestAngleTilt,recId);
        mdReconstruction.setValue(MDL_ANGLE_PSI,bestAnglePsi,recId);
        mdReconstruction.setValue(MDL_SHIFT_X,bestShiftX,recId);
        mdReconstruction.setValue(MDL_SHIFT_Y,bestShiftY,recId);
        mdReconstruction.setValue(MDL_FLIP,bestFlip,recId);
        mdReconstruction.setValue(MDL_WEIGHT,bestWeight,recId);
    }
    newCorr=bestCorr;
#ifdef DEBUG
    mdReconstruction.write("PPPreconstruction.xmd");
    std::cout << "Press any key" << std::endl;
    char c;
    std::cin >> c;
#endif
}

threadAlignSubset()

void threadAlignSubset

(

ThreadArgument &

thArg

)

Definition at line 228 of file volume_initial_simulated_annealing.cpp.

{
    ProgVolumeInitialSimulatedAnnealing &prm=*((ProgVolumeInitialSimulatedAnnealing *)thArg.workClass);
    ThreadVolumeInitialAlignment results=prm.threadResults[thArg.thread_id];

    results.sumCorr=results.sumImprovement=0.0;
    results.mdReconstruction.clear();
    auto nMax=(int)prm.mdInp.size();
    for (int nImg=0; nImg<nMax; ++nImg)
    {
        if ((nImg+1)%prm.Nthr==thArg.thread_id)
        {
            double corr, improvement;
            alignSingleImage(nImg, prm, results.mdReconstruction, corr, improvement);
            results.sumCorr+=corr;
            results.sumImprovement+=improvement;
            prm.mdInp[nImg].maxcc=corr;
        }

        if (thArg.thread_id==0)
            progress_bar(nImg+1);
    }

    // Update all MAXCC
    prm.mutexMaxCC.lock();

    prm.mdReconstruction.unionAll(results.mdReconstruction);
    prm.sumCorr+=results.sumCorr;
    prm.sumImprovement+=results.sumImprovement;

    prm.mutexMaxCC.unlock();
}