angular_projection_matching.h File Reference

#include "core/metadata_db.h"
#include "core/multidim_array.h"
#include "core/xmipp_program.h"
#include "core/xmipp_threads.h"
#include "data/polar.h"
#include "data/sampling.h"

Include dependency graph for angular_projection_matching.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
class	Image< T >
struct	structThreadRotationallyAlignOneImage
class	ProgAngularProjectionMatching

Functions
void *	threadRotationallyAlignOneImage (void *data)

Variables
constexpr int	MY_OUPUT_SIZE = 10

Function Documentation

threadRotationallyAlignOneImage()

void* threadRotationallyAlignOneImage

(

void *

data

)

Definition at line 530 of file angular_projection_matching.cpp.

{
    auto * thread_data = (structThreadRotationallyAlignOneImage *) data;

    // Variables from above
    size_t thread_id = thread_data->thread_id;
    ProgAngularProjectionMatching *prm = thread_data->prm;
    size_t thread_num = prm->threads;
    MultidimArray<double> *img = thread_data->img;
    size_t &this_image = thread_data->this_image;


    int *opt_refno = thread_data->opt_refno;
    double *opt_psi = thread_data->opt_psi;
    bool *opt_flip = thread_data->opt_flip;
    double *maxcorr = thread_data->maxcorr;

    // Local variables
    MultidimArray<double>       Maux;
    MultidimArray<double>       ang, corr;
    MultidimArray<int>          indxCorr;
    size_t                      myinit, myfinal, myincr;
    int                         refno;
    bool                        done_once=false;
    double                      mean, stddev;
    Polar<double>               P;
    Polar<std::complex <double> > fP,fPm;
    RotationalCorrelationAux    rotAux;
    Polar_fftw_plans            local_plans;
    size_t                         imgno = this_image - FIRST_IMAGE;

#ifdef TIMING

    TimeStamp t0,t1,t2;
    time_config();
    annotate_time(&t0);
    annotate_time(&t2);
#endif

    produceSplineCoefficients(xmipp_transformation::BSPLINE3,Maux,*img);
    // Precalculate polar transform of each translation
    // This loop is also threaded
    myinit = thread_id;
    myincr = thread_num;
    for (size_t itrans = myinit; itrans < prm->nr_trans; itrans+=myincr)
    {
        P.getPolarFromCartesianBSpline(Maux,prm->Ri,prm->Ro,3,
                                       (double)prm->search5d_xoff[itrans],
                                       (double)prm->search5d_yoff[itrans]);
        P.computeAverageAndStddev(mean,stddev);
        P -= mean; // for normalized cross-correlation coefficient
        if (itrans == myinit)
            P.calculateFftwPlans(local_plans);
        fourierTransformRings(P,prm->fP_img[itrans],local_plans,false);
        fourierTransformRings(P,prm->fPm_img[itrans],local_plans,true);
        prm->stddev_img[itrans] = stddev;
        done_once=true;
    }
    // If thread did not have to do any itrans, initialize fftw plans
    if (!done_once)
    {
        P.getPolarFromCartesianBSpline(Maux,prm->Ri,prm->Ro);
        P.calculateFftwPlans(local_plans);
    }
    // Prepare FFTW plan for rotational correlation
    corr.resize(P.getSampleNoOuterRing());
    rotAux.local_transformer.setReal(corr);
    rotAux.local_transformer.FourierTransform();

    // All threads have to wait until the itrans loop is done
    barrier_wait(&(prm->thread_barrier));

#ifdef TIMING

    float prepare_img = elapsed_time(t0);
    float get_refs = 0.;
    annotate_time(&t0);
#endif

    //pthread_mutex_lock(  &debug_mutex );
    // Switch the order of looping through the references every time.
    // That way, in case max_nr_refs_in_memory<total_nr_refs
    // the references read in memory for the previous image
    // will still be there when processing the next image
    // Every thread processes a part of the references
    if (prm->loop_forward_refs)
    {
        myinit = 0;
        myfinal = prm->mysampling.my_neighbors[imgno].size();
        myincr = +1;
    }
    else
    {
        myinit = prm->mysampling.my_neighbors[imgno].size() - 1;
        myfinal = -1;
        myincr = -1;
    }
    // Loop over all relevant "neighbours" (i.e. directions within the search range)
    //for (int i = myinit; i != myfinal; i+=myincr)
    for (size_t i = myinit; i != myfinal; i += myincr)
    {
        if (i%thread_num == thread_id)
        {

#ifdef DEBUG_THREADS
            pthread_mutex_lock(  &debug_mutex );
            std::cerr<<" thread_id= "<<thread_id<<" i= "<<i<<" "<<myinit<<" "<<myfinal<<" "<<myincr<<std::endl;
            pthread_mutex_unlock(  &debug_mutex );
#endif

#ifdef TIMING

            annotate_time(&t1);
#endif
            // Get pointer to the current reference image
#ifdef DEBUG

            if(prm->mysampling.my_neighbors[imgno][i]==58)
            {
                std::cerr << "XXXXpointer_allrefs2refsinmemXXXXXX" <<std::endl;
                for (std::vector<int>::iterator i = prm->
                                                    pointer_allrefs2refsinmem.begin();
                     i != prm->pointer_allrefs2refsinmem.end();
                     ++i)
                    std::cerr << *i << std::endl;
                std::cerr << "XXXXpointer_refsinmem2allrefsXXXXXX" <<std::endl;
                for (std::vector<int>
                     ::iterator i = prm->pointer_refsinmem2allrefs.begin();
                     i != prm->pointer_refsinmem2allrefs.end();
                     ++i)
                    std::cerr << *i << std::endl;
                std::cerr <<std::endl;

            }
#endif

            refno = prm->pointer_allrefs2refsinmem[prm->mysampling.my_neighbors[imgno][i]];
            if (refno == -1)
            {
                // Reference is not stored in memory (anymore): (re-)read from disc
                prm->getCurrentReference(prm->mysampling.my_neighbors[imgno][i],local_plans);
                refno = prm->pointer_allrefs2refsinmem[prm->mysampling.my_neighbors[imgno][i]];
            }


#ifdef TIMING
            get_refs += elapsed_time(t1);
#endif
            //#define DEBUG
#ifdef DEBUG

            std::cerr << "imgno " << imgno <<std::endl;
            std::cerr<<"Got refno= "<<refno
            <<" pointer= "<<prm->mysampling.my_neighbors[imgno][i]<<std::endl;
#endif

            // Loop over all 5D-search translations
            for (size_t itrans = 0; itrans < prm->nr_trans; itrans++)
            {
#ifdef DEBUG

                std::cerr<< "prm->stddev_ref[refno], prm->stddev_img[itrans]: " <<
                prm->stddev_ref[refno] << " " <<
                prm->stddev_img[itrans];
#endif

                MultidimArray<double>  allCorr, allAng;
                allCorr.resizeNoCopy(XSIZE(corr)*2);
                allAng.resizeNoCopy(XSIZE(corr)*2);

                // A. Check straight image
                rotationalCorrelation(prm->fP_img[itrans],
                                      prm->fP_ref[refno],
                                      ang,rotAux);
                corr /= prm->stddev_ref[refno] * prm->stddev_img[itrans]; // for normalized ccf


                memcpy(&dAi(allCorr,0),&dAi(corr,0),XSIZE(corr)*sizeof(double));
                memcpy(&dAi(allAng,0),&dAi(ang,0),XSIZE(ang)*sizeof(double));

                rotationalCorrelation(prm->fPm_img[itrans],prm->fP_ref[refno],ang,rotAux);
                corr /= prm->stddev_ref[refno] * prm->stddev_img[itrans]; // for normalized ccf
                memcpy(&dAi(allCorr,XSIZE(corr)),&dAi(corr,0),XSIZE(corr)*sizeof(double));
                memcpy(&dAi(allAng,XSIZE(corr)),&dAi(ang,0),XSIZE(corr)*sizeof(double));

                size_t nIter = XMIPP_MIN(thread_data->numOrientations,corr.getSize());
                double bestLastCorr = 99e99;
                for (size_t n = 0; n < nIter; n++)
                {
                    for (size_t k = 0; k < XSIZE(allCorr); k++)
                    {
                        if ( (DIRECT_A1D_ELEM(allCorr,k)> maxcorr[n]) && (DIRECT_A1D_ELEM(allCorr,k)<bestLastCorr))
                        {
                            maxcorr[n] = DIRECT_A1D_ELEM(allCorr,k);
                            opt_psi[n] = DIRECT_A1D_ELEM(allAng,k);
                            //FIXME not sure about FIRST_IMAGE
                            opt_refno[n] = prm->mysampling.my_neighbors[imgno][i];/*+FIRST_IMAGE;*/
                            if ( k >= XSIZE(corr))
                                opt_flip[n] = true;
                            else
                                opt_flip[n] = false;
                        }
                    }

                    bestLastCorr = maxcorr[n];
                }


#ifdef DEBUG
                std::cerr<<"mirror: corr "<<maxcorr;
                if (opt_flip)
                    std::cerr<<"**";
                std::cerr<<std::endl;
#endif
#undef DEBUG

            }
            //#define DEBUG
#ifdef DEBUG
            std::cerr << "DEBUG_ROB, imgno:" << imgno << std::endl;
            std::cerr << "DEBUG_ROB, i:" << i << std::endl;
            std::cerr << "DEBUG_ROB, prm->mysampling.my_neighbors[imgno][i]:" << prm->mysampling.my_neighbors[imgno][i] << std::endl;
            std::cerr<<"straight: corr "<<maxcorr<<std::endl;
#endif
#undef DEBUG

        }
    }


#ifdef TIMING
    float all_rot_align = elapsed_time(t0);
    float total_rot = elapsed_time(t2);
    std::cerr<<" rotal%% "<<total_rot
    <<" => prep: "<<prepare_img
    <<" all_refs: "<<all_rot_align
    <<" (of which "<<get_refs
    <<" to get "<< prm->mysampling.my_neighbors[imgno].size()
    <<" refs for imgno "<<imgno<<" )"
    <<std::endl;
#endif
    //pthread_mutex_unlock(  &debug_mutex );
    //std::cerr << "DEBUG_JM: threadRotationallyAlignOneImage END" <<std::endl;
    return nullptr;
}

Variable Documentation

MY_OUPUT_SIZE

constexpr int MY_OUPUT_SIZE = 10

Definition at line 40 of file angular_projection_matching.h.