frm.cpp File Reference

#include "frm.h"
#include <core/geometry.h>
#include <core/transformations.h>

Include dependency graph for frm.cpp:

Go to the source code of this file.

Functions
void	alignVolumesFRM (PyObject *pFunc, const MultidimArray< double > &Iref, MultidimArray< double > &I, PyObject *Imask, double &rot, double &tilt, double &psi, double &x, double &y, double &z, double &score, Matrix2D< double > &A, int maxshift, double maxFreq, const MultidimArray< int > *mask)

Function Documentation

alignVolumesFRM()

void alignVolumesFRM	(	PyObject *	pFunc,
		const MultidimArray< double > &	Iref,
		MultidimArray< double > &	I,
		PyObject *	Imask,
		double &	rot,
		double &	tilt,
		double &	psi,
		double &	x,
		double &	y,
		double &	z,
		double &	score,
		Matrix2D< double > &	A,
		int	maxshift = 10,
		double	maxFreq = 0.25,
		const MultidimArray< int > *	mask = nullptr
	)

Align two volumes using FRM. The first argument is the pointer to the FRM python function. You may obtain it with getPointerToPythonFRMFunction()

Imask is a mask in Fourier space for I Maxshift is in pixels. MaxFreq is in digital frequency.

A is the transformation matrix that needs to be applied on I to fit Iref.

If apply is set, then I is substituted with the aligned volume.

Definition at line 35 of file frm.cpp.

{
    PyObject *pyIref = Python::convertToNumpy(Iref);
    PyObject *pyI = Python::convertToNumpy(I);
    PyObject *pyMask=Py_None;
    if (mask!=nullptr)
        pyMask=Python::convertToNumpy(*mask);

//#define DEBUG
#ifdef DEBUG
    Image<double> save;
    save()=I;
    save.write("PPPI.vol");
    save()=Iref;
    save.write("PPPIref.vol");
    std::cout << "Imask=" << Imask << std::endl;
    std::cout << "Press any key\n";
    char c; std::cin >> c;
#endif

    // Call frm
    int bandWidthSphericalHarmonics0=4;
    int bandWidthSphericalHarmonicsF=64;
    auto frequencyPixels=(int)(XSIZE(Iref)*maxFreq);
    PyObject *arglistfrm = Py_BuildValue("OOOO(ii)iiO", pyI, Imask, pyIref, Py_None,
            bandWidthSphericalHarmonics0, bandWidthSphericalHarmonicsF, frequencyPixels, maxshift, pyMask);
    PyObject *resultfrm = PyObject_CallObject(pFunc, arglistfrm);
    Py_DECREF(arglistfrm);
    Py_DECREF(pyIref);
    Py_DECREF(pyI);
    if (mask!=nullptr)
        Py_DECREF(pyMask);
    if (resultfrm!=nullptr)
    {
        PyObject *shift=PyTuple_GetItem(resultfrm,0);
        PyObject *euler=PyTuple_GetItem(resultfrm,1);
        score=PyFloat_AsDouble(PyTuple_GetItem(resultfrm,2));
        double xfrm=PyFloat_AsDouble(PyList_GetItem(shift,0))-XSIZE(Iref)/2;
        double yfrm=PyFloat_AsDouble(PyList_GetItem(shift,1))-YSIZE(Iref)/2;
        double zfrm=PyFloat_AsDouble(PyList_GetItem(shift,2))-ZSIZE(Iref)/2;
        double angz1=PyFloat_AsDouble(PyList_GetItem(euler,0));
        double angz2=PyFloat_AsDouble(PyList_GetItem(euler,1));
        double angx=PyFloat_AsDouble(PyList_GetItem(euler,2));
        Matrix2D<double> Efrm, E(3,3);
        Euler_anglesZXZ2matrix(-angz1, -angx, -angz2, Efrm); // -angles because FRM rotation definition
                                                             // is the opposite of Xmipp

        // Reorganize the matrix because the Z and X axes in the coordinate system are reversed with
        // respect to Xmipp
        // Exmipp=[0 0 1; 0 1 0; 1 0 0]*Efrm*[0 0 1; 0 1 0; 1 0 0]
        MAT_ELEM(E,0,0)=MAT_ELEM(Efrm, 2, 2); MAT_ELEM(E,0,1)=MAT_ELEM(Efrm, 2, 1); MAT_ELEM(E,0,2)=MAT_ELEM(Efrm, 2, 0);
        MAT_ELEM(E,1,0)=MAT_ELEM(Efrm, 1, 2); MAT_ELEM(E,1,1)=MAT_ELEM(Efrm, 1, 1); MAT_ELEM(E,1,2)=MAT_ELEM(Efrm, 1, 0);
        MAT_ELEM(E,2,0)=MAT_ELEM(Efrm, 0, 2); MAT_ELEM(E,2,1)=MAT_ELEM(Efrm, 0, 1); MAT_ELEM(E,2,2)=MAT_ELEM(Efrm, 0, 0);
        E=E.inv();
        Euler_matrix2angles(E,rot,tilt,psi);
        Py_DECREF(resultfrm);

        x=-zfrm;
        y=-yfrm;
        z=-xfrm;

        // Apply
        Euler_angles2matrix(rot, tilt, psi, A, true);
        Matrix2D<double> Aaux;
        Matrix1D<double> r(3);
        XX(r)=x;
        YY(r)=y;
        ZZ(r)=z;
        translation3DMatrix(r,Aaux);
        A=A*Aaux;
#ifdef DEBUG
        std::cout << "Result: " << rot << " " << tilt << " " << psi << " " << x << " " << y << " " << z << " -> " << score << std::endl;
#endif
    }
    else
    {
        x=y=z=rot=tilt=psi=score=0;
        A.initIdentity(4);

        std::cout << "No result\n";
        std::cout << Python::print_traceback();
    }
}