common_lines.cpp File Reference

#include <fstream>
#include "common_lines.h"
#include "data/mask.h"
#include "core/metadata_extension.h"
#include "data/fourier_filter.h"
#include "reconstruction/radon.h"
#include "core/linear_system_helper.h"

Include dependency graph for common_lines.cpp:

Go to the source code of this file.

Classes
struct	ThreadPrepareImages
struct	ThreadCompareImages

Macros
#define	q0   dMi(q, 0)
#define	q1   dMi(q, 1)
#define	q2   dMi(q, 2)
#define	q3   dMi(q, 3)
#define	AXIS(var)   DVector var(3); var.initZeros()
#define	cl(i, j)   dMij(clMatrix, k##i, k##j)
#define	SET_POINT(x, y)   dMij(syncMatrix, x+2*k1, y+2*k2) = dMij(R, x, y)
#define	SORTED_INDEX(i)   dMi(indexes, i)
#define	SORTED_ELEM(M, i)   dMi(M, SORTED_INDEX(i))

Functions
void *	threadPrepareImages (void *args)
void	commonLineTwoImages (std::vector< MultidimArray< std::complex< double > > > &RTFsi, std::vector< MultidimArray< double > > &RTsi, int idxi, std::vector< MultidimArray< std::complex< double > > > &RTFsj, std::vector< MultidimArray< double > > &RTsj, int idxj, ProgCommonLine *parent, CommonLine &result, FourierTransformer &transformer)
void *	threadCompareImages (void *args)
void	randomQuaternions (int k, DMatrix &quaternions)
void	saveMatrix (const char *fn, DMatrix &matrix)
void	quaternionToMatrix (const DVector &q, DMatrix &rotMatrix)
void	quaternionCommonLines (const DMatrix &quaternions, CommonLineInfo &clInfo)
void	commonlineMatrixCheat (const DMatrix &quaternions, size_t nRays, DMatrix &clMatrix, DMatrix &clCorr)
void	anglesRotationMatrix (size_t nRays, int clI, int clJ, const DVector &Q1, const DVector &Q2, DMatrix &R)
int	tripletRotationMatrix (const DMatrix &clMatrix, size_t nRays, int k1, int k2, int k3, DMatrix &R)
void	putRotationMatrix (const DMatrix &R, int k1, int k2, DMatrix &syncMatrix)
void	computeSyncMatrix (const DMatrix &clMatrix, size_t nRays, DMatrix &sMatrix)
void	rotationsFromSyncMatrix (const DMatrix &sMatrix)

Variables
constexpr long double	EPS = 1.0e-13
constexpr int	MAX_COND = 1000

Macro Definition Documentation

AXIS

#define AXIS

(

var

)

DVector var(3); var.initZeros()

Definition at line 631 of file common_lines.cpp.

cl

#define cl	(		i,
			j
	)		dMij(clMatrix, k##i, k##j)

Definition at line 830 of file common_lines.cpp.

q0

#define q0   dMi(q, 0)

q1

#define q1   dMi(q, 1)

q2

#define q2   dMi(q, 2)

q3

#define q3   dMi(q, 3)

SET_POINT

#define SET_POINT	(		x,
			y
	)		dMij(syncMatrix, x+2*k1, y+2*k2) = dMij(R, x, y)

SORTED_ELEM

#define SORTED_ELEM	(		M,
			i
	)		dMi(M, SORTED_INDEX(i))

SORTED_INDEX

#define SORTED_INDEX

(

i

)

dMi(indexes, i)

Function Documentation

anglesRotationMatrix()

void anglesRotationMatrix	(	size_t	nRays,
		int	clI,
		int	clJ,
		const DVector &	Q1,
		const DVector &	Q2,
		DMatrix &	R
	)

Definition at line 793 of file common_lines.cpp.

{

    double alpha1 = TWOPI * clI / nRays;
    double alpha2 = TWOPI * clJ / nRays;
    DVector N1(3);
    vectorProduct(Q1, Q2, N1);
    N1.selfNormalize();

    DMatrix U(3, 3);

    //sincos(alpha1, &dMij(U, 1, 0), &dMij(U, 0, 0));
    dMij(U, 1, 0) = sin(alpha1);
    dMij(U, 0, 0) = cos(alpha1);
    //sincos(alpha2, &dMij(U, 1, 1), &dMij(U, 0, 1));
    dMij(U, 1, 1) = sin(alpha2);
    dMij(U, 0, 1) = cos(alpha2);
    dMij(U, 2, 2) = 1.;

    if (U.det3x3() < 0)
        dMij(U, 2, 2) = -1.; //% Make sure we have a right-handed system)

    DMatrix Q(3, 3);
    Q.setCol(0, Q1);
    Q.setCol(1, Q2);
    Q.setCol(2, N1);

    U.selfInverse();
    DMatrix T = Q * U;
    T.svd(U, N1, Q);
    R = U * Q.transpose();
}//function

commonLineTwoImages()

void commonLineTwoImages	(	std::vector< MultidimArray< std::complex< double > > > &	RTFsi,
		std::vector< MultidimArray< double > > &	RTsi,
		int	idxi,
		std::vector< MultidimArray< std::complex< double > > > &	RTFsj,
		std::vector< MultidimArray< double > > &	RTsj,
		int	idxj,
		ProgCommonLine *	parent,
		CommonLine &	result,
		FourierTransformer &	transformer
	)

Definition at line 247 of file common_lines.cpp.

                                                                                     {
    MultidimArray<std::complex<double> > &RTFi = RTFsi[idxi];
    MultidimArray<std::complex<double> > &RTFj = RTFsj[idxj];
    MultidimArray<double> &RTi = RTsi[idxi];
    MultidimArray<double> &RTj = RTsj[idxj];

    result.distanceij = 1e60;
    result.angj = -1;
    MultidimArray<std::complex<double> > lineFi, lineFj;
    MultidimArray<double> linei, linej;
    MultidimArray<double> correlationFunction;
    int jmax;
    for (size_t ii = 0; ii < YSIZE(RTFi) / 2 + 1; ii++) {
        lineFi.aliasRow(RTFi,ii);
        linei.aliasRow(RTi,ii);

        for (size_t jj=0; jj<YSIZE(RTFj); jj++)
        {
            lineFj.aliasRow(RTFj,jj);
            linej.aliasRow(RTj,jj);

            // Compute distance between the two lines
            fast_correlation_vector(lineFi, lineFj, correlationFunction, transformer);
            correlationFunction.maxIndex(jmax);
            double distance=1-A1D_ELEM(correlationFunction,jmax);

            // Check if this is the best match
            if (distance<result.distanceij)
            {
                result.distanceij=distance;
                result.angi=ii;
                result.angj=jj;
                result.jmax=jmax;
            }
        }
    }
    result.angi *= parent->stepAng;
    result.angj *= parent->stepAng;
}

putRotationMatrix()

void putRotationMatrix	(	const DMatrix &	R,
		int	k1,
		int	k2,
		DMatrix &	syncMatrix
	)

Helper function to set the 2x2 upper left corner of a rotation matrix into a bigger matrix indexed by images

Definition at line 882 of file common_lines.cpp.

{
#define SET_POINT(x, y) dMij(syncMatrix, x+2*k1, y+2*k2) = dMij(R, x, y)
    SET_POINT(0, 0);
    SET_POINT(0, 1);
    SET_POINT(1, 0);
    SET_POINT(1, 1);
}

threadCompareImages()

void* threadCompareImages

(

void *

args

)

Definition at line 306 of file common_lines.cpp.

{
    auto * master = (ThreadCompareImages *) args;
    ProgCommonLine * parent = master->parent;

    int blockIsize = master->RTFsi->size();
    int blockJsize = master->RTFsj->size();
    FourierTransformer transformer;
    MultidimArray<double> linei;
    linei.resize(parent->Xdim);
    transformer.setReal(linei);
    for (int i = 0; i < blockIsize; i++) {
        long int ii = parent->Nblock * master->i + i;
        for (int j = 0; j < blockJsize; j++) {
            // Check if this two images have to be compared
            long int jj = parent->Nblock * master->j + j;
            if (ii >= jj)
                continue;
            if ((ii * blockJsize + jj + 1) % parent->Nthr != master->myThreadID)
                continue;

            // Effectively compare the two images
            long int idx_ij = ii * parent->Nimg + jj;
            commonLineTwoImages(
                    *(master->RTFsi), *(master->RTsi), i,
                    *(master->RTFsj), *(master->RTsj), j,
                    parent, parent->CLmatrix[idx_ij], transformer);

            // Compute the symmetric element
            long int idx_ji = jj * parent->Nimg + ii;
            parent->CLmatrix[idx_ji].distanceij = parent->CLmatrix[idx_ij].distanceij;
            parent->CLmatrix[idx_ji].angi = parent->CLmatrix[idx_ij].angj;
            parent->CLmatrix[idx_ji].angj = parent->CLmatrix[idx_ij].angi;
            parent->CLmatrix[idx_ji].jmax = -parent->CLmatrix[idx_ij].jmax;
        }
    }
    return nullptr;
}

threadPrepareImages()

void* threadPrepareImages

(

void *

args

)

Definition at line 106 of file common_lines.cpp.

{
    auto * master = (ThreadPrepareImages *) args;
    ProgCommonLine * parent = master->parent;
    MetaDataVec SFi = *(master->SFi);
    size_t Ydim, Xdim, Zdim, Ndim;
    getImageSize(SFi, Xdim, Ydim, Zdim, Ndim);

    MultidimArray<int> mask;
    mask.resize(Ydim, Xdim);
    mask.setXmippOrigin();
    BinaryCircularMask(mask, Xdim / 2, OUTSIDE_MASK);
    auto NInsideMask = (int)(XSIZE(mask) * YSIZE(mask) - mask.sum());

    FourierFilter Filter;
    Filter.w1 = -1;
    if (parent->lpf > 0 && parent->hpf > 0)
    {
        Filter.FilterBand = BANDPASS;
        Filter.w1 = parent->lpf;
        Filter.w2 = parent->hpf;
    }
    else if (parent->lpf > 0)
    {
        Filter.FilterBand = LOWPASS;
        Filter.w1 = parent->lpf;
    }
    else if (parent->hpf > 0)
    {
        Filter.FilterBand = HIGHPASS;
        Filter.w1 = parent->hpf;
    }
    Filter.raised_w = Filter.w1 / 3;

    int ii = 0;
    bool first = true;
    Image<double> I;
    FileName fnImg;
    MultidimArray<double> RT, linei(Xdim);
    FourierTransformer transformer;
    transformer.setReal(linei);
    MultidimArray<std::complex<double> >&mlineiFourier=transformer.fFourier;
    MultidimArray<std::complex<double> > RTFourier;
    for (size_t objId : SFi.ids())
    {
        if ((ii + 1) % parent->Nthr == master->myThreadID) {
            I.readApplyGeo(SFi, objId);
            I().setXmippOrigin();
            MultidimArray<double> &mI = I();

            // Bandpass filter images
            if (Filter.w1 > 0)
            {
                if (first)
                {
                    Filter.generateMask(mI);
                    first = false;
                }
                Filter.applyMaskSpace(mI);
            }

            // Cut the image outside the largest circle
            // And compute the DC value inside the mask
            double meanInside = 0;
            FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
                if (A2D_ELEM(mask,i,j))A2D_ELEM(mI,i,j)=0;
                else
                meanInside+=A2D_ELEM(mI,i,j);
            meanInside /= NInsideMask;

            // Subtract the mean inside the circle
            FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
            if (!A2D_ELEM(mask,i,j))
                A2D_ELEM(mI,i,j)-=meanInside;

            // Compute the Radon transform
            Radon_Transform(I(), parent->stepAng, RT);

            // Normalize each line in the Radon Transform so that the
            // multiplication of any two lines is actually their correlation index
            RTFourier.resize(YSIZE(RT), XSIZE(mlineiFourier));
            for (size_t i = 0; i < YSIZE(RT); i++) {
                memcpy(&(DIRECT_A1D_ELEM(linei,0)),&DIRECT_A2D_ELEM(RT,i,0),XSIZE(linei)*sizeof(double));
                linei.statisticsAdjust(0.,1.);
                transformer.FourierTransform();
                transformer.fReal->initZeros();
                transformer.inverseFourierTransform();
                memcpy(&DIRECT_A2D_ELEM(RTFourier,i,0),&(DIRECT_A1D_ELEM(mlineiFourier,0)),XSIZE(mlineiFourier)*sizeof(std::complex<double>));

                memcpy(&DIRECT_A2D_ELEM(RT,i,0),&(DIRECT_A1D_ELEM(linei,0)),XSIZE(linei)*sizeof(double));
            }

            (*(master->blockRTFs))[ii] = RTFourier;
            (*(master->blockRTs))[ii] = RT;
        }
        ii++;
    }
    return nullptr;
}

Variable Documentation

EPS

constexpr long double EPS = 1.0e-13

Definition at line 828 of file common_lines.cpp.

MAX_COND

constexpr int MAX_COND = 1000

Definition at line 829 of file common_lines.cpp.