tomo_align_tilt_series.cpp File Reference

#include <fstream>
#include <queue>
#include "tomo_align_tilt_series.h"
#include "core/metadata_sql.h"
#include "core/xmipp_image.h"
#include "core/transformations.h"
#include "core/utils/memory_utils.h"
#include "data/fourier_filter.h"
#include "data/mask.h"
#include "data/morphology.h"
#include "data/numerical_tools.h"

Include dependency graph for tomo_align_tilt_series.cpp:

Go to the source code of this file.

Classes
class	AffineFitness
struct	ThreadComputeTransformParams
struct	ThreadGenerateLandmarkSetParams

Namespaces
	TomographAlignment

Macros
#define	DEBUG

Functions
void	generateMask (const MultidimArray< double > &I, MultidimArray< unsigned char > &mask, int patchSize)
void	setupAffineFitness (AffineFitness &fitness, const MultidimArray< double > &I1, const MultidimArray< double > &I2, int maxShift, bool isMirror, bool checkRotation, int pyramidLevel)
double	computeAffineTransformation (const MultidimArray< unsigned char > &I1, const MultidimArray< unsigned char > &I2, int maxShift, int maxIterDE, Matrix2D< double > &A12, Matrix2D< double > &A21, bool show, double thresholdAffine, bool globalAffine, bool isMirror, bool checkRotation, int pyramidLevel)
void *	threadComputeTransform (void *args)
void *	threadgenerateLandmarkSetGrid (void *args)
void *	threadgenerateLandmarkSetBlind (void *args)
void *	threadgenerateLandmarkSetCriticalPoints (void *args)
double	wrapperError (double *p, void *prm)
std::ostream &	operator<< (std::ostream &out, Alignment &alignment)

Variables
const ProgTomographAlignment *	TomographAlignment::global_prm

Macro Definition Documentation

DEBUG

#define DEBUG

Definition at line 2517 of file tomo_align_tilt_series.cpp.

Function Documentation

computeAffineTransformation()

double computeAffineTransformation	(	const MultidimArray< unsigned char > &	I1,
		const MultidimArray< unsigned char > &	I2,
		int	maxShift,
		int	maxIterDE,
		Matrix2D< double > &	A12,
		Matrix2D< double > &	A21,
		bool	show,
		double	thresholdAffine,
		bool	globalAffine,
		bool	isMirror,
		bool	checkRotation,
		int	pyramidLevel
	)

Definition at line 345 of file tomo_align_tilt_series.cpp.

{
    try
    {
        AffineFitness affy;
        MultidimArray<double> I1d, I2d;
        typeCast(I1, I1d);
        typeCast(I2, I2d);
        setupAffineFitness(affy, I1d, I2d, maxShift, isMirror, checkRotation,
                           pyramidLevel);

        // Return result
        double cost = 0.0;

        // Optimize with differential evolution
        Matrix1D<double> A(6);
        A(0)=A(3)=1;
        if (isMirror)
            A(3)*=-1;
        if (MAT_XSIZE(A12)==0)
        {
            if (globalAffine)
            {
                // Exhaustive search
                double bestCost=2;
                double stepX=(affy.maxAllowed(4)-affy.minAllowed(4))/40.0;
                double stepY=(affy.maxAllowed(5)-affy.minAllowed(5))/40.0;
                for (double shiftY=affy.minAllowed(5); shiftY<=affy.maxAllowed(5); shiftY+=stepY)
                    for (double shiftX=affy.minAllowed(4); shiftX<=affy.maxAllowed(4); shiftX+=stepX)
                    {
                        A(4)=shiftX;
                        A(5)=shiftY;
                        double cost=affy.affine_fitness_individual(MATRIX1D_ARRAY(A));
                        if (cost<bestCost)
                            bestCost=cost;
                    }
            }
            else
            {
                // Initialize with cross correlation
                double tx, ty;
                pthread_mutex_lock( &globalAffineMutex );
                CorrelationAux aux;
                if (!isMirror)
                    bestShift(I1d,I2d,tx,ty,aux);
                else
                {
                    MultidimArray<double> auxI2d=I2d;
                    auxI2d.selfReverseY();
                    STARTINGX(auxI2d)=STARTINGX(I2d);
                    STARTINGY(auxI2d)=STARTINGY(I2d);
                    bestShift(I1d,auxI2d,tx,ty,aux);
                    ty=-ty;
                }
                pthread_mutex_unlock( &globalAffineMutex );
                A(4)=-tx;
                A(5)=-ty;
            }

            // Optimize with Powell
            Matrix1D<double> steps(A);
            steps.initConstant(1);
            int iter;
            powellOptimizer(A, 1, VEC_XSIZE(A),
                            AffineFitness::Powell_affine_fitness_individual, &affy, 0.005,
                            cost, iter, steps, false);

            // Separate solution
            A12.initIdentity(3);
            A12(0,0)=A(0);
            A12(0,1)=A(1);
            A12(0,2)=A(4);
            A12(1,0)=A(2);
            A12(1,1)=A(3);
            A12(1,2)=A(5);

            A21=A12.inv();
        }

        if (show)
        {
            affy.showMode=true;
            Matrix1D<double> p(6);
            p(0)=A12(0,0);
            p(1)=A12(0,1);
            p(4)=A12(0,2);
            p(2)=A12(1,0);
            p(3)=A12(1,1);
            p(5)=A12(1,2);
            cost = affy.affine_fitness_individual(MATRIX1D_ARRAY(p));
            affy.showMode=false;
        }

        return cost;
    }
    catch (XmippError &XE)
    {
        std::cout << XE;
        exit(1);
    }
}

generateMask()

void generateMask	(	const MultidimArray< double > &	I,
		MultidimArray< unsigned char > &	mask,
		int	patchSize
	)

Definition at line 40 of file tomo_align_tilt_series.cpp.

{
    MultidimArray<double> dmask;
    dmask.initZeros(I);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(dmask)
    if (I(i,j)!=0)
        dmask(i,j)=1;

    MultidimArray<double> maskEroded;
    maskEroded.initZeros(dmask);
    erode2D(dmask,maskEroded,8,0,patchSize);
    typeCast(maskEroded,mask);
    mask.setXmippOrigin();

#ifdef DEBUG

    ImageXmipp save;
    save()=I;
    save.write("PPP.xmp");
    save()=maskEroded;
    save.write("PPPmask.xmp");
    std::cout << "Masks saved. Press any key\n";
    char c;
    std::cin >> c;
#endif
}

operator<<()

std::ostream& operator<<	(	std::ostream &	out,
		Alignment &	alignment
	)

Print an alignment

Definition at line 2879 of file tomo_align_tilt_series.cpp.

{
    out << "Alignment parameters ===========================================\n"
    << "rot=" << alignment.rot << " tilt=" << alignment.tilt << std::endl
    << "raxis=" << alignment.raxis.transpose() << std::endl;
    out << "Images ---------------------------------------------------------\n";
    for (int i=0; i<alignment.Nimg; i++)
        out << "Image " << i << " psi= " << alignment.psi(i)
        << " di= " << alignment.di[i].transpose()
        << " diaxis= " << alignment.diaxis[i].transpose()
        << " (" << alignment.prm->ni(i) << ")\n";
    out << "Landmarks ------------------------------------------------------\n";
    alignment.computeErrorForLandmarks();
    for (int j=0; j<alignment.Nlandmark; j++)
        out << "Landmark " << j << " rj= " << alignment.rj[j].transpose()
        << " " << alignment.errorLandmark(j) << std::endl;
    return out;
}

setupAffineFitness()

void setupAffineFitness	(	AffineFitness &	fitness,
		const MultidimArray< double > &	I1,
		const MultidimArray< double > &	I2,
		int	maxShift,
		bool	isMirror,
		bool	checkRotation,
		int	pyramidLevel
	)

Definition at line 251 of file tomo_align_tilt_series.cpp.

{
    fitness.checkRotation=checkRotation;

    // Set images
    int level=0;
    MultidimArray<double> I1aux=I1;
    MultidimArray<double> I2aux=I2;

    // Remove the borders
    int borderY=CEIL(0.1*YSIZE(I1));
    int borderX=CEIL(0.1*XSIZE(I1));
    I1aux.selfWindow(STARTINGY(I1)+borderY,STARTINGX(I1)+borderX,
                     FINISHINGY(I1)-borderY,FINISHINGX(I1)-borderX);
    I2aux.selfWindow(STARTINGY(I1)+borderY,STARTINGX(I1)+borderX,
                     FINISHINGY(I1)-borderY,FINISHINGX(I1)-borderX);

    MultidimArray<double> Mask1;
    Mask1.initZeros(I1aux);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I1aux)
    if (A2D_ELEM(I1,i,j)!=0)
        A2D_ELEM(Mask1,i,j)=1;

    MultidimArray<double> Mask2;
    Mask2.initZeros(I2aux);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I2aux)
    if (A2D_ELEM(I2,i,j)!=0)
        A2D_ELEM(Mask2,i,j)=1;

    do
    {
        // Push back the current images
        I1aux.setXmippOrigin();
        I2aux.setXmippOrigin();
        Mask1.setXmippOrigin();
        Mask2.setXmippOrigin();

        MultidimArray<double> *dummy=nullptr;
        dummy=new MultidimArray<double>;
        *dummy=I1aux;
        fitness.I1.push_back(dummy);
        dummy=new MultidimArray<double>;
        *dummy=I2aux;
        fitness.I2.push_back(dummy);
        dummy=new MultidimArray<double>;
        *dummy=Mask1;
        fitness.Mask1.push_back(dummy);
        dummy=new MultidimArray<double>;
        *dummy=Mask2;
        fitness.Mask2.push_back(dummy);

        // Prepare for next level
        level++;
        if (pyramidLevel>=level)
        {
            selfScaleToSize(xmipp_transformation::LINEAR,I1aux,YSIZE(I1aux)/2,XSIZE(I1aux)/2);
            selfScaleToSize(xmipp_transformation::LINEAR,I2aux,YSIZE(I2aux)/2,XSIZE(I2aux)/2);
            selfScaleToSize(xmipp_transformation::LINEAR,Mask1,YSIZE(Mask1)/2,XSIZE(Mask1)/2);
            selfScaleToSize(xmipp_transformation::LINEAR,Mask2,YSIZE(Mask2)/2,XSIZE(Mask2)/2);
            FOR_ALL_ELEMENTS_IN_ARRAY2D(Mask1)
            Mask1(i,j)=(Mask1(i,j)>0.5)? 1:0;
            FOR_ALL_ELEMENTS_IN_ARRAY2D(Mask2)
            Mask2(i,j)=(Mask2(i,j)>0.5)? 1:0;
        }
    }
    while (level<=pyramidLevel);

    // Set limits for the affine matrices
    // Order: 1->2: 4 affine params+2 translations
    // Order: 2->1: 4 affine params+2 translations
    fitness.minAllowed.resize(6);
    fitness.maxAllowed.resize(6);

    // Scale factors
    fitness.minAllowed(0)=fitness.minAllowed(3)=0.9;
    fitness.maxAllowed(0)=fitness.maxAllowed(3)=1.1;
    if (isMirror)
    {
        fitness.minAllowed(3)=-1.1;
        fitness.maxAllowed(3)=-0.9;
    }

    // Rotation factors
    fitness.minAllowed(1)=fitness.minAllowed(2)=-0.2;
    fitness.maxAllowed(1)=fitness.maxAllowed(2)= 0.2;

    // Shifts
    fitness.minAllowed(4)=fitness.minAllowed(5)=-maxShift;
    fitness.maxAllowed(4)=fitness.maxAllowed(5)= maxShift;
}

threadComputeTransform()

void* threadComputeTransform

(

void *

args

)

Definition at line 584 of file tomo_align_tilt_series.cpp.

{
    auto * master =
        (ThreadComputeTransformParams *) args;

    ProgTomographAlignment * parent = master->parent;
    int thread_id = master->myThreadID;
    int localnumThreads = parent->numThreads;
    bool isCapillar = parent->isCapillar;
    int Nimg = parent->Nimg;
    double maxShiftPercentage = parent->maxShiftPercentage;
    int maxIterDE = parent->maxIterDE;
    bool showAffine = parent->showAffine;
    double thresholdAffine = parent->thresholdAffine;
    std::vector < MultidimArray<unsigned char> *> & img = parent->img;
    std::vector< std::vector< Matrix2D<double> > > & affineTransformations = parent->affineTransformations;
    double globalAffine = parent->globalAffine;

    int maxShift=FLOOR(XSIZE(*img[0])*maxShiftPercentage);
    int initjj=1;
    if (isCapillar)
        initjj=0;

    initjj += thread_id;

    double cost;
    for (int jj=initjj; jj<Nimg; jj+= localnumThreads)
    {
        int jj_1;
        if (isCapillar)
            jj_1=intWRAP(jj-1,0,Nimg-1);
        else
            jj_1=jj-1;
        MultidimArray<unsigned char>& img_i=*img[jj_1];
        MultidimArray<unsigned char>& img_j=*img[jj];
        bool isMirror=(jj==0) && (jj_1==Nimg-1);

        if (MAT_XSIZE(affineTransformations[jj_1][jj])==0)
        {
            Matrix2D<double> Aij, Aji;
            cost = computeAffineTransformation(img_i, img_j, maxShift,
                                               maxIterDE, Aij, Aji,
                                               showAffine, thresholdAffine, globalAffine,
                                               isMirror,true, parent->pyramidLevel);
            parent->correlationList[jj]=1-cost;

            pthread_mutex_lock( &printingMutex );
            affineTransformations[jj_1][jj]=Aij;
            affineTransformations[jj][jj_1]=Aji;
            if (cost<1)
                parent->writeTransformations(
                    parent->fnRoot+"_transformations.txt");
            std::cout << "Cost for [" << jj_1 << "] - ["
            << jj << "] = " << cost << std::endl;
            parent->iteration++;
            pthread_mutex_unlock( &printingMutex );
        }
        else
        {
            Matrix2D<double> Aij;
            Aij=affineTransformations[jj_1][jj];

            MultidimArray<double> img_id, img_jd;
            typeCast(img_i, img_id);
            typeCast(img_j, img_jd);

            AffineFitness fitness;
            setupAffineFitness(fitness, img_id, img_jd, maxShift, isMirror,
                               false, parent->pyramidLevel);
            Matrix1D<double> p(6);
            p(0)=Aij(0,0);
            p(1)=Aij(0,1);
            p(4)=Aij(0,2);
            p(2)=Aij(1,0);
            p(3)=Aij(1,1);
            p(5)=Aij(1,2);
            cost = fitness.affine_fitness_individual(MATRIX1D_ARRAY(p));
            parent->correlationList[jj]=1-cost;

            pthread_mutex_lock( &printingMutex );
            std::cout << "Cost for [" << jj_1 << "] - ["
            << jj << "] = " << cost << std::endl;
            pthread_mutex_unlock( &printingMutex );
            parent->iteration++;
        }
    }

    return nullptr;
}

threadgenerateLandmarkSetBlind()

void* threadgenerateLandmarkSetBlind

(

void *

args

)

Definition at line 1216 of file tomo_align_tilt_series.cpp.

{
    auto * master =
        (ThreadGenerateLandmarkSetParams *) args;
    ProgTomographAlignment * parent = master->parent;
    int thread_id = master->myThreadID;
    int Nimg=parent->Nimg;
    int numThreads=parent->numThreads;
    const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
        parent->affineTransformations;
    int gridSamples=parent->gridSamples;

    auto deltaShift=(int)floor(XSIZE(*(parent->img)[0])/gridSamples);
    master->chainList=new std::vector<LandmarkChain>;
    Matrix1D<double> rii(3), rjj(3);
    ZZ(rii)=1;
    ZZ(rjj)=1;
    if (thread_id==0)
        init_progress_bar(gridSamples);
    int includedPoints=0;
    int maxSideLength=(parent->blindSeqLength-1)/2;
    for (int nx=thread_id; nx<gridSamples; nx+=numThreads)
    {
        XX(rii)=STARTINGX(*(parent->img)[0])+ROUND(deltaShift*(0.5+nx));
        for (int ny=0; ny<gridSamples; ny+=1)
        {
            YY(rii)=STARTINGY(*(parent->img)[0])+ROUND(deltaShift*(0.5+ny));
            for (int ii=0; ii<=Nimg-1; ++ii)
            {
                // Check if inside the mask
                if (!(*(parent->maskImg[ii]))((int)YY(rii),(int)XX(rii)))
                    continue;
                if (parent->isOutlier(ii))
                    continue;

                LandmarkChain chain;
                chain.clear();
                Landmark l;
                l.x=XX(rii);
                l.y=YY(rii);
                l.imgIdx=ii;
                chain.push_back(l);

                // Follow this landmark backwards
                bool acceptLandmark=true;
                int jj;
                int sideLength=0;
                if (parent->isCapillar)
                    jj=intWRAP(ii-1,0,Nimg-1);
                else
                    jj=ii-1;
                Matrix2D<double> Aij, Aji;
                Matrix1D<double> rcurrent=rii;
                while (jj>=0 && acceptLandmark && sideLength<maxSideLength &&
                       !parent->isOutlier(jj))
                {
                    // Compute the affine transformation between ii and jj
                    int jj_1;
                    if (parent->isCapillar)
                        jj_1=intWRAP(jj+1,0,Nimg-1);
                    else
                        jj_1=jj+1;
                    Aij=affineTransformations[jj][jj_1];
                    Aji=affineTransformations[jj_1][jj];
                    rjj=Aji*rcurrent;
                    auto iYYrjj=(int)YY(rjj);
                    auto iXXrjj=(int)XX(rjj);
                    if (!(*(parent->maskImg[jj])).outside(iYYrjj,iXXrjj))
                        acceptLandmark=(*(parent->maskImg[jj]))(iYYrjj,iXXrjj);
                    else
                        acceptLandmark=false;
                    if (acceptLandmark)
                    {
                        l.x=XX(rjj);
                        l.y=YY(rjj);
                        l.imgIdx=jj;
                        chain.push_back(l);
                        if (parent->isCapillar)
                            jj=intWRAP(jj-1,0,Nimg-1);
                        else
                            jj=jj-1;
                        rcurrent=rjj;
                        sideLength++;
                    }
                }

                // Follow this landmark forward
                acceptLandmark=true;
                if (parent->isCapillar)
                    jj=intWRAP(ii+1,0,Nimg-1);
                else
                    jj=ii+1;
                rcurrent=rii;
                sideLength=0;
                while (jj<Nimg && acceptLandmark && sideLength<maxSideLength &&
                       !parent->isOutlier(jj))
                {
                    // Compute the affine transformation between ii and jj
                    int jj_1;
                    if (parent->isCapillar)
                        jj_1=intWRAP(jj-1,0,Nimg-1);
                    else
                        jj_1=jj-1;
                    Aij=affineTransformations[jj_1][jj];
                    Aji=affineTransformations[jj][jj_1];
                    rjj=Aij*rcurrent;
                    auto iYYrjj=(int)YY(rjj);
                    auto iXXrjj=(int)XX(rjj);
                    if (!(*(parent->maskImg[jj])).outside(iYYrjj,iXXrjj))
                        acceptLandmark=(*(parent->maskImg[jj]))(iYYrjj,iXXrjj);
                    else
                        acceptLandmark=false;
                    if (acceptLandmark)
                    {
                        l.x=XX(rjj);
                        l.y=YY(rjj);
                        l.imgIdx=jj;
                        chain.push_back(l);
                        if (parent->isCapillar)
                            jj=intWRAP(jj+1,0,Nimg-1);
                        else
                            jj=jj+1;
                        rcurrent=rjj;
                        sideLength++;
                    }
                }

#ifdef DEBUG
                std::cout << "img=" << ii << " chain length="
                << chain.size() << " [" << jjleft
                << " - " << jjright << "]\n";
#endif

                master->chainList->push_back(chain);
                includedPoints+=chain.size();
            }
#ifdef DEBUG
            std::cout << "Point nx=" << nx << " ny=" << ny
            << " Number of points="
            << includedPoints
            << " Number of chains=" << master->chainList->size()
            << " ( " << ((double) includedPoints)/
            master->chainList->size() << " )\n";
#endif

        }
        if (thread_id==0)
            progress_bar(nx);
    }
    if (thread_id==0)
        progress_bar(gridSamples);
    return nullptr;
}

threadgenerateLandmarkSetCriticalPoints()

void* threadgenerateLandmarkSetCriticalPoints

(

void *

args

)

Definition at line 1371 of file tomo_align_tilt_series.cpp.

{
    auto * master =
        (ThreadGenerateLandmarkSetParams *) args;
    ProgTomographAlignment * parent = master->parent;
    int thread_id = master->myThreadID;
    int Nimg=parent->Nimg;
    int numThreads=parent->numThreads;
    const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
        parent->affineTransformations;

    master->chainList=new std::vector<LandmarkChain>;
    std::vector<LandmarkChain> candidateChainList;
    if (thread_id==0)
        init_progress_bar(Nimg);
    int halfSeqLength=parent->seqLength/2;

    // Design a mask for the dilation
    int radius=4;
    MultidimArray<int> mask;
    mask.resize(2*radius+1,2*radius+1);
    mask.setXmippOrigin();
    BinaryCircularMask(mask,4,OUTSIDE_MASK);

    Image<double> I;
    for (int ii=thread_id; ii<=Nimg-1; ii+=numThreads)
    {
        if (parent->isOutlier(ii))
            continue;
        I.read(parent->name_list[ii]);

        // Generate mask
        MultidimArray<unsigned char> largeMask;
        generateMask(I(),largeMask,
                     XMIPP_MAX(ROUND(parent->localSize*XSIZE(I()))/2,5));

        // Filter the image
        MultidimArray<double> Ifiltered;
        Ifiltered=I();
        Ifiltered.setXmippOrigin();
        FourierFilter FilterBP;
        FilterBP.FilterBand=BANDPASS;
        FilterBP.w1=2.0/XSIZE(Ifiltered);
        FilterBP.w2=128.0/XSIZE(Ifiltered);
        FilterBP.raised_w=1.0/XSIZE(Ifiltered);
        ;
        FilterBP.generateMask(Ifiltered);
        FilterBP.applyMaskSpace(Ifiltered);

        // Identify low valued points and perform dilation
        MultidimArray<double> Iaux=Ifiltered;
        Iaux.selfWindow(
            -ROUND(0.45*YSIZE(Ifiltered)),-ROUND(0.45*XSIZE(Ifiltered)),
            ROUND(0.45*YSIZE(Ifiltered)), ROUND(0.45*XSIZE(Ifiltered)));
        Histogram1D hist;
        compute_hist(Iaux, hist, 400);
        double th=hist.percentil(2);
        std::vector< Matrix1D<double> > Q;
        FOR_ALL_ELEMENTS_IN_ARRAY2D(Iaux)
        if (Ifiltered(i,j)<th && largeMask(i,j))
        {
            // Check if it is a local minimum
            bool localMinimum=true;
            int x0=j;
            int y0=i;
            FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
            {
                int x=x0+j;
                int y=y0+i;
                if (x>=STARTINGX(Iaux) && x<=FINISHINGX(Iaux) &&
                    y>=STARTINGY(Iaux) && y<=FINISHINGY(Iaux))
                    if (A2D_ELEM(Iaux,y,x)<A2D_ELEM(Iaux,y0,x0))
                    {
                        localMinimum=false;
                        break;
                    }
            }
            if (localMinimum)
                Q.push_back(vectorR2(x0,y0));
        }

        // Check that the list is not too long
        if (Q.size()>10*parent->Ncritical)
        {
            size_t qmax=Q.size();
            MultidimArray<double> minDistance;
            minDistance.resize(qmax);
            minDistance.initConstant(1e20);
            for (size_t q1=0; q1<qmax; q1++)
                for (size_t q2=q1+1; q2< qmax; q2++)
                {
                    double diffX=XX(Q[q1])-XX(Q[q2]);
                    double diffY=YY(Q[q1])-YY(Q[q2]);
                    double d12=diffX*diffX+diffY*diffY;
                    if (d12<minDistance(q1))
                        minDistance(q1)=d12;
                    if (d12<minDistance(q2))
                        minDistance(q2)=d12;
                }
            MultidimArray<int> idxDistanceSort;
            minDistance.indexSort(idxDistanceSort);
            std::vector< Matrix1D<double> > Qaux;
            int qlimit=XMIPP_MIN(10*(parent->Ncritical),qmax);
            for (int q=0; q<qlimit; q++)
                Qaux.push_back(Q[idxDistanceSort(qmax-1-q)-1]);
            Q.clear();
            Q=Qaux;
        }

        int qmax=Q.size();
        Matrix1D<double> rii(3), rjj(3);
        ZZ(rii)=1;
        ZZ(rjj)=1;
        MultidimArray<double> corrQ;
        corrQ.initZeros(qmax);
        for (int q=0; q<qmax; q++)
        {
            XX(rii)=XX(Q[q]);
            YY(rii)=YY(Q[q]);

            // Initiate a chain here
            LandmarkChain chain;
            chain.clear();
            Landmark l;
            l.x=XX(rii);
            l.y=YY(rii);
            l.imgIdx=ii;
            chain.push_back(l);

            // Follow this landmark backwards
            int jjmin=XMIPP_MAX(0,ii-halfSeqLength);
            int jjmax=ii-1;
            Matrix2D<double> Aij, Aji;
            Matrix1D<double> rcurrent=rii;
            for (int jj=jjmax; jj>=jjmin; --jj)
            {
                if (parent->isOutlier(jj))
                    break;

                // Compute the affine transformation between jj and jj+1
                int jj_1=jj+1;
                Aij=affineTransformations[jj][jj_1];
                Aji=affineTransformations[jj_1][jj];
                rjj=Aji*rcurrent;
                double corr;
                parent->refineLandmark(jj_1,jj,rcurrent,rjj, corr,true);
                l.x=XX(rjj);
                l.y=YY(rjj);
                l.imgIdx=jj;
                chain.push_back(l);
                rcurrent=rjj;
            }

            // Follow this landmark forwards
            jjmin=ii+1;
            jjmax=XMIPP_MIN(Nimg-1,ii+halfSeqLength);
            rcurrent=rii;
            for (int jj=jjmin; jj<=jjmax; ++jj)
            {
                if (parent->isOutlier(jj))
                    break;

                // Compute the affine transformation between jj-1 and jj
                int jj_1=jj-1;
                Aij=affineTransformations[jj_1][jj];
                Aji=affineTransformations[jj][jj_1];
                rjj=Aij*rcurrent;
                double corr;
                parent->refineLandmark(jj_1,jj,rcurrent,rjj,corr,true);
                l.x=XX(rjj);
                l.y=YY(rjj);
                l.imgIdx=jj;
                chain.push_back(l);
                rcurrent=rjj;
            }

            // Refine chain
            parent->refineChain(chain,corrQ(q));
            candidateChainList.push_back(chain);
        }
        if (thread_id==0)
            progress_bar(ii);

        // Sort all chains according to its correlation
        MultidimArray<int> idx;
        corrQ.indexSort(idx);
        int imax=XMIPP_MIN(XSIZE(idx),parent->Ncritical);
        for (int iq=0; iq<imax; iq++)
        {
            int q=idx(XSIZE(idx)-1-iq)-1;
            if (corrQ(q)>0.5)
            {
                master->chainList->push_back(candidateChainList[q]);
#ifdef DEBUG

                std::cout << "Corr " << iq << ": " << corrQ(q) << ":";
                for (int i=0; i<candidateChainList[q].size(); i++)
                    std::cout << candidateChainList[q][i].imgIdx << " ";
                std::cout << std::endl;
#endif

            }
            else
                std::cout << "It's recommended to reduce the number of critical points\n";
        }
        candidateChainList.clear();

#ifdef DEBUG

        Image<double> save;
        typeCast(*(parent->img[ii]),save());
        save.write("PPPoriginal.xmp");
        save()=Ifiltered;
        save.write("PPPfiltered.xmp");
        double minval=Ifiltered.computeMin();
        FOR_ALL_ELEMENTS_IN_MATRIX2D(Ifiltered)
        if (Ifiltered(i,j)>=th || !largeMask(i,j))
            save(i,j)=th;
        for (int q=0; q<Q.size(); q++)
            FOR_ALL_ELEMENTS_IN_MATRIX2D(mask)
            if (YY(Q[q])+i>=STARTINGY(Ifiltered) && YY(Q[q])+i<=FINISHINGY(Ifiltered) &&
                XX(Q[q])+j>=STARTINGX(Ifiltered) && XX(Q[q])+j<=FINISHINGX(Ifiltered))
                save(YY(Q[q])+i,XX(Q[q])+j)=minval;
        imax=XMIPP_MIN(XSIZE(idx),parent->Ncritical);
        for (int iq=0; iq<imax; iq++)
        {
            int q=idx(XSIZE(idx)-1-iq)-1;
            FOR_ALL_ELEMENTS_IN_MATRIX2D(mask)
            if (YY(Q[q])+i>=STARTINGY(Ifiltered) && YY(Q[q])+i<=FINISHINGY(Ifiltered) &&
                XX(Q[q])+j>=STARTINGX(Ifiltered) && XX(Q[q])+j<=FINISHINGX(Ifiltered))
                save(YY(Q[q])+i,XX(Q[q])+j)=(minval+th)/2;
        }

        save.write("PPPcritical.xmp");
        std::cout << "Number of critical points=" << Q.size() << std::endl;
        std::cout << "CorrQ stats:";
        corrQ.printStats();
        std::cout << std::endl;
        std::cout << "Press any key\n";
        char c;
        std::cin >> c;
#endif

    }
    if (thread_id==0)
        progress_bar(Nimg);
    return nullptr;
}

threadgenerateLandmarkSetGrid()

void* threadgenerateLandmarkSetGrid

(

void *

args

)

Definition at line 1047 of file tomo_align_tilt_series.cpp.

{
    auto * master =
        (ThreadGenerateLandmarkSetParams *) args;
    ProgTomographAlignment * parent = master->parent;
    int thread_id = master->myThreadID;
    int Nimg=parent->Nimg;
    int numThreads=parent->numThreads;
    const std::vector< std::vector< Matrix2D<double> > > &affineTransformations=
        parent->affineTransformations;
    int gridSamples=parent->gridSamples;

    auto deltaShift=(int)floor(XSIZE(*(parent->img)[0])/gridSamples);
    master->chainList=new std::vector<LandmarkChain>;
    Matrix1D<double> rii(3), rjj(3);
    ZZ(rii)=1;
    ZZ(rjj)=1;
    if (thread_id==0)
        init_progress_bar(gridSamples);
    int includedPoints=0;
    Matrix1D<int> visited(Nimg);
    for (int nx=thread_id; nx<gridSamples; nx+=numThreads)
    {
        XX(rii)=STARTINGX(*(parent->img)[0])+ROUND(deltaShift*(0.5+nx));
        for (int ny=0; ny<gridSamples; ny+=1)
        {
            YY(rii)=STARTINGY(*(parent->img)[0])+ROUND(deltaShift*(0.5+ny));
            for (int ii=0; ii<=Nimg-1; ++ii)
            {
                // Check if inside the mask
                if (!(*(parent->maskImg[ii]))((int)YY(rii),(int)XX(rii)))
                    continue;
                if (parent->isOutlier(ii))
                    continue;

                LandmarkChain chain;
                chain.clear();
                Landmark l;
                l.x=XX(rii);
                l.y=YY(rii);
                l.imgIdx=ii;
                chain.push_back(l);
                visited.initZeros();
                visited(ii)=1;

                // Follow this landmark backwards
                bool acceptLandmark=true;
                int jj;
                if (parent->isCapillar)
                    jj=intWRAP(ii-1,0,Nimg-1);
                else
                    jj=ii-1;
                Matrix2D<double> Aij, Aji;
                Matrix1D<double> rcurrent=rii;
                while (jj>=0 && acceptLandmark && !visited(jj) &&
                       !parent->isOutlier(jj))
                {
                    // Compute the affine transformation between ii and jj
                    int jj_1;
                    if (parent->isCapillar)
                        jj_1=intWRAP(jj+1,0,Nimg-1);
                    else
                        jj_1=jj+1;
                    visited(jj)=1;
                    Aij=affineTransformations[jj][jj_1];
                    Aji=affineTransformations[jj_1][jj];
                    rjj=Aji*rcurrent;
                    double corr;
                    acceptLandmark=parent->refineLandmark(jj_1,jj,rcurrent,rjj,
                                                          corr,true);
                    if (acceptLandmark)
                    {
                        l.x=XX(rjj);
                        l.y=YY(rjj);
                        l.imgIdx=jj;
                        chain.push_back(l);
                        if (parent->isCapillar)
                            jj=intWRAP(jj-1,0,Nimg-1);
                        else
                            jj=jj-1;
                        rcurrent=rjj;
                    }
                }

                // Follow this landmark forward
                acceptLandmark=true;
                if (parent->isCapillar)
                    jj=intWRAP(ii+1,0,Nimg-1);
                else
                    jj=ii+1;
                rcurrent=rii;
                while (jj<Nimg && acceptLandmark && !visited(jj) &&
                       !parent->isOutlier(jj))
                {
                    // Compute the affine transformation between ii and jj
                    int jj_1;
                    if (parent->isCapillar)
                        jj_1=intWRAP(jj-1,0,Nimg-1);
                    else
                        jj_1=jj-1;
                    visited(jj)=1;
                    Aij=affineTransformations[jj_1][jj];
                    Aji=affineTransformations[jj][jj_1];
                    rjj=Aij*rcurrent;
                    double corr;
                    acceptLandmark=parent->refineLandmark(jj_1,jj,rcurrent,rjj,
                                                          corr,true);
                    if (acceptLandmark)
                    {
                        l.x=XX(rjj);
                        l.y=YY(rjj);
                        l.imgIdx=jj;
                        chain.push_back(l);
                        if (parent->isCapillar)
                            jj=intWRAP(jj+1,0,Nimg-1);
                        else
                            jj=jj+1;
                        rcurrent=rjj;
                    }
                }

#ifdef DEBUG
                std::cout << "img=" << ii << " chain length="
                << chain.size() << " [" << jjleft
                << " - " << jjright << "]";
#endif

                if (chain.size()>parent->seqLength)
                {
                    double corrChain;
                    bool accepted=parent->refineChain(chain,corrChain);
                    if (accepted)
                    {
#ifdef DEBUG
                        std::cout << " Accepted with length= "
                        << chain.size()
                        << " [" << chain[0].imgIdx << " - "
                        << chain[chain.size()-1].imgIdx << "]= ";
                        for (int i=0; i<chain.size(); i++)
                            std::cout << chain[i].imgIdx << " ";
#endif

                        master->chainList->push_back(chain);
                        includedPoints+=chain.size();
                    }
                }
#ifdef DEBUG
                std::cout << std::endl;
#endif

            }
#ifdef DEBUG
            std::cout << "Point nx=" << nx << " ny=" << ny
            << " Number of points="
            << includedPoints
            << " Number of chains=" << master->chainList->size()
            << " ( " << ((double) includedPoints)/
            master->chainList->size() << " )\n";
#endif

        }
        if (thread_id==0)
            progress_bar(nx);
    }
    if (thread_id==0)
        progress_bar(gridSamples);
    return nullptr;
}

wrapperError()

double wrapperError	(	double *	p,
		void *	prm
	)

Definition at line 2508 of file tomo_align_tilt_series.cpp.

{
    Alignment alignment(TomographAlignment::global_prm);
    alignment=*(TomographAlignment::global_prm->bestPreviousAlignment);
    alignment.rot=p[1];
    alignment.tilt=p[2];
    return alignment.optimizeGivenAxisDirection();
}