CrystalARTRecons Class Reference

#include <art_crystal.h>

Inheritance diagram for CrystalARTRecons:

Collaboration diagram for CrystalARTRecons:

Public Member Functions
void	readParams (XmippProgram *program)
void	print (std::ostream &o) const
	std::cout << crystal_prm; More...
void	preProcess (GridVolume &vol_basis0, int level=FULL, int rank=-1)
void	singleStep (GridVolume &vol_in, GridVolume *vol_out, Projection &theo_proj, Projection &read_proj, int sym_no, Projection &diff_proj, Projection &corr_proj, Projection &alig_proj, double &mean_error, int numIMG, double lambda, int act_proj, const FileName &fn_ctf, const MultidimArray< int > *maskPtr, bool refine)
void	postProcess (GridVolume &vol_basis)
void	applySymmetry (GridVolume &vol_in, GridVolume *vol_out, int grid_type)
Public Member Functions inherited from ARTReconsBase
virtual	~ARTReconsBase ()
void	initHistory (const GridVolume &vol_basis0)
void	iterations (GridVolume &vol_basis, int rank=-1)

Static Public Member Functions
static void	defineParams (XmippProgram *program, const char *prefix=nullptr, const char *comment=nullptr)
Static Public Member Functions inherited from ARTReconsBase
static void	defineParams (XmippProgram *program, bool mpiMode=false)

Public Attributes
double	a_mag
double	b_mag
double	ang_a2b_deg
	angle from a to b (degrees) More...
double	ang_x2a_deg
	angle from x axis to a (degrees) More...
Matrix1D< double >	avox
	First lattice vector (voxel units) More...
Matrix1D< double >	bvox
	Second lattice vector (voxel units) More...
bool	fill_space
	Fill space, repeat unit cell all over the space. More...
Matrix1D< double >	a
	First lattice vector (BCC units) More...
Matrix1D< double >	b
	Second lattice vector (BCC units) More...
Matrix1D< double >	aint
	First lattice vector approximated to integer numbers (BCC units) More...
Matrix1D< double >	bint
	Second lattice vector approximated to integer numbers (BCC units) More...
Matrix1D< double >	ai
	ai=aint/2 as double numbers More...
Matrix1D< double >	bi
	bi=aint/2 as double numbers More...
int	space_group
	space_group More...
MultidimArray< int >	unit_cell_mask
Public Attributes inherited from ARTReconsBase
BasicARTParameters	artPrm

Detailed Description

ART+crystal parameters. Here only those specific parameters for crystals are found, the rest of parameters common with normal ART should be looked up in BasicARTParameters

Definition at line 43 of file art_crystal.h.

Member Function Documentation

applySymmetry()

void CrystalARTRecons::applySymmetry ( GridVolume &  vol_in, GridVolume *  vol_out, int  grid_type  )

virtual

Force the trial volume to be symmetric.

Reimplemented from ARTReconsBase.

Definition at line 388 of file art_crystal.cpp.

{
    symmetrizeCrystalVolume(vol_in,aint,bint,space_group,unit_cell_mask,grid_type);
}

defineParams()

void CrystalARTRecons::defineParams ( XmippProgram *  program, const char *  prefix = nullptr, const char *  comment = nullptr  )

static

Definition at line 36 of file art_crystal.cpp.

{
    std::string line = "  [--crystal ]   : Crystal mode activation";

    if(nullptr != prefix)
        line = prefix + line;
    if (nullptr != comment)
        line = line + " : " + comment;

    program->addParamsLine(line);
    program->addParamsLine("   [--mag_a <mag>]         : Magnitude of the first  lattice vector");
    program->addParamsLine("   [--mag_b <mag>]         : Magnitude of the second lattice vector");
    program->addParamsLine("   [--ang_a2b_deg <angle>] : Angle from vector a to vector b");
    program->addParamsLine("   [--ang_x2a_deg <angle=0>] : Angle from vector x to vector a");
    program->addParamsLine("   [--fill_space]          : Repeat unit cell all over the space (pixels)");

}

postProcess()

void CrystalARTRecons::postProcess ( GridVolume &  vol_basis )

virtual

Finish iterations. For WLS: delete residual images Else: do nothing.

Reimplemented from ARTReconsBase.

Definition at line 382 of file art_crystal.cpp.

{
    if (fill_space)
        expandToFillSpace(artPrm, *this, vol_basis);
}

preProcess()

void CrystalARTRecons::preProcess ( GridVolume &  vol_basis0, int  level = FULL, int  rank = -1  )

virtual

Produce Initial and Side information for ART+crystal. This function computes from the ART+crystal parameters things like the integer<–>float lattice vector matrices, the unit cell mask. The reconstructed volume size can be modified such that the unit cell mask fits.

Reimplemented from ARTReconsBase.

Definition at line 114 of file art_crystal.cpp.

{
    // Lattice vectors in BCC units
    a = avox / artPrm.grid_relative_size;
    b = bvox / artPrm.grid_relative_size;

    // Compute space_group
    if (artPrm.fn_sym != "")
        space_group = artPrm.SL.crystallographicSpaceGroup(a_mag, b_mag,
                      ang_a2b_deg);
    else
        space_group = sym_P1;

    // Integer lattice vectors ----------------------------------------------
    Matrix2D<double> D;
    computeIntegerLattice(a, b, a_mag / artPrm.grid_relative_size,
                          b_mag / artPrm.grid_relative_size,
                          ang_a2b_deg, aint, bint, D,
                          space_group);
    // Define two more useful lattice vectors
    ai = aint / 2;
    bi = bint / 2;

    // Set general deformation pointer to this matrix
    artPrm.D = new Matrix2D<double>;
    *(artPrm.D) = D;
    artPrm.D->resize(3, 3);
    MAT_ELEM(*(artPrm.D), 2, 2) = 1;
    artPrm.Dinv = new Matrix2D<double>;
    *(artPrm.Dinv) = artPrm.D->inv();
    artPrm.basis.setD(artPrm.D);

    // Unit cell mask within volume -----------------------------------------
    // Compute the 4 parallelogram corners
    Matrix1D<double> c1 = ai + bi;
    Matrix1D<double> c2 = -ai + bi;
    Matrix1D<double> c3 = -c1;
    Matrix1D<double> c4 = -c2;
    Matrix1D<double> c1c3 = c1 - c3; // These extra variables are needed because
    Matrix1D<double> c2c4 = c2 - c4; // the compiler messes up

    // Resize unit cell mask
    // The unit mask is a little bigger to avoid the possibility of losing
    // any basis due to a too tight mask.
    int mask_xdim = CEIL(XMIPP_MAX(ABS(XX(c1c3)), ABS(XX(c2c4)))) + 3;
    int mask_ydim = CEIL(XMIPP_MAX(ABS(YY(c1c3)), ABS(YY(c2c4)))) + 3;
    unit_cell_mask.initZeros(mask_ydim, mask_xdim);
    unit_cell_mask.setXmippOrigin();

    // Resize the reconstructed volume
    Matrix1D<double> r1(3), r2(3);
    for (size_t n = 0; n < vol_basis0.VolumesNo(); n++)
    {
        Image<double> &V = vol_basis0(n);
        ZZ(r1) = XMIPP_MIN(ZZ(r1), STARTINGZ(V()));
        ZZ(r2) = XMIPP_MAX(ZZ(r2), FINISHINGZ(V()));
    }
    XX(r1) = STARTINGX(unit_cell_mask);
    YY(r1) = STARTINGY(unit_cell_mask);
    r1 *= artPrm.grid_relative_size;
    XX(r2) = FINISHINGX(unit_cell_mask);
    YY(r2) = FINISHINGY(unit_cell_mask);
    r2 *= artPrm.grid_relative_size;
    vol_basis0.resize(r1, r2);

    // Fill the unit cell mask
    Matrix1D<double> r(2);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(unit_cell_mask)
    {
        // Position vector of actual BCC element been considered
        YY(r) = i;
        XX(r) = j;

        // Vectors from r to each corner
        Matrix1D<double> c1r = c1 - r;
        Matrix1D<double> c2r = c2 - r;
        Matrix1D<double> c3r = c3 - r;
        Matrix1D<double> c4r = c4 - r;

        // Product of each of these vectors with tha parallelogram borders
        int sgn[4];
        sgn[0] = SGN0_VEC_PRODUCT(c1r, -ai);
        sgn[1] = SGN0_VEC_PRODUCT(c2r, -bi);
        sgn[2] = SGN0_VEC_PRODUCT(c3r, ai);
        sgn[3] = SGN0_VEC_PRODUCT(c4r, bi);

#ifdef DEBUG_A_LOT

        std::cout << "(x,y)=(" << XX(r) << "," << YY(r) << ") " << sgn[0] << sgn[1]
        << sgn[2] << sgn[3] << std::endl;
#endif

        // Now check if point is inside
        int inside_table[4][4] = {
                                     {1, 1, 1, 1},
                                     {0, 1, 1, 1},
                                     {1, 0, 1, 1},
                                     // 1,1,0,1, Take this side out
                                     // 1,1,1,0, Take this side out
                                     {0, 0, 1, 1}
                                     // 1,0,0,1, and their three Vertex
                                     // 1,1,0,0,
                                     // 0,1,1,0
                                 };
#define INSIDE_ACCORDING_TO(n) \
    sgn[0]==inside_table[n][0] && sgn[1]==inside_table[n][1] && \
    sgn[2]==inside_table[n][2] && sgn[3]==inside_table[n][3]

        for (int n = 0; n < 4; n++)
            if (INSIDE_ACCORDING_TO(n))
            {
                unit_cell_mask(i, j) = 1;
#ifdef DEBUG_A_LOT

                std::cout << "    Inside\n";
#endif

                break;
            }
    }

    // Now calculate side info of basel class
    ARTReconsBase::preProcess(vol_basis0);

    // Show all parameters --------------------------------------------------
#ifdef DEBUG
    std::cout << "avox= " << avox.transpose() << std::endl;
    std::cout << "bvox= " << bvox.transpose() << std::endl;
    std::cout << "grid_relative_size= " << artPrm.grid_relative_size << std::endl;
    std::cout << "a=    " << a.transpose()    << std::endl;
    std::cout << "b=    " << b.transpose()    << std::endl;
    std::cout << "aint= " << aint.transpose() << std::endl;
    std::cout << "bint= " << bint.transpose() << std::endl;
    std::cout << "ai=   " << ai.transpose()   << std::endl;
    std::cout << "bi=   " << bi.transpose()   << std::endl;
    std::cout << "D=    " << D;
    std::cout << "Check that a=D*aint " << (D*aint).transpose() << std::endl;
    std::cout << "Check that b=D*bint " << (D*bint).transpose() << std::endl;
    std::cout << "Symmetry group: " <<  space_group;
    ImageXmipp I;
    I = unit_cell_mask;
    I.write("unit_cell_mask.xmp");
    std::cout << "unit_cell_mask shape=";
    unit_cell_mask.printShape();
    std::cout << std::endl;
#endif
}

print()

void CrystalARTRecons::print ( std::ostream &  o ) const

virtual

std::cout << crystal_prm;

Reimplemented from ARTReconsBase.

Definition at line 89 of file art_crystal.cpp.

{
    o << "Crystal information ------------------------------------------" << std::endl;
    o << "Lattice vector a: " << a.transpose() << std::endl;
    o << "Lattice vector b: " << b.transpose() << std::endl;
    o << "mag_a: " << a_mag << std::endl;
    o << "mag_b: " << b_mag << std::endl;
    o << "ang_a2b_deg: " << ang_a2b_deg << std::endl;
    o << "ang_x2a_deg: " << ang_x2a_deg << std::endl;
    o << "Fill space: " << fill_space << std::endl;
    o << "Symmetry group: " <<  space_group << std::endl;
}

readParams()

void CrystalARTRecons::readParams ( XmippProgram *  program )

virtual

Read special parameters from command line.

Reimplemented from ARTReconsBase.

Definition at line 54 of file art_crystal.cpp.

{
    ARTReconsBase::readParams(program);
    artPrm.is_crystal = true;
    a_mag = program->getDoubleParam("--mag_a");
    a_mag /= artPrm.sampling;
    b_mag = program->getDoubleParam("--mag_b");
    b_mag /= artPrm.sampling;
    ang_a2b_deg = program->getDoubleParam("--ang_a2b_deg");
    ang_x2a_deg = program->getDoubleParam("--ang_x2a_deg");
    fill_space = program->checkParam("--fill_space");
    avox.resize(2);
    bvox.resize(2);
    //DISCLAMER: I know this 0 and 90 degrees cases are rather silly but
    //when debugging is so good that 90 is 90 and not 90.000001
    //NOTE:ang_x2a_deg is applied ONLY in the final volume
    //when moving from basis to voxels

    XX(avox) = a_mag;
    YY(avox) = 0.;

    if (ang_a2b_deg == 90.)
    {
        XX(bvox) = 0;
        YY(bvox) = b_mag;
    }
    else
    {
        XX(bvox) = b_mag * COSD(ang_a2b_deg);
        YY(bvox) = b_mag * SIND(ang_a2b_deg);
    }

}

singleStep()

void CrystalARTRecons::singleStep ( GridVolume &  vol_in, GridVolume *  vol_out, Projection &  theo_proj, Projection &  read_proj, int  sym_no, Projection &  diff_proj, Projection &  corr_proj, Projection &  alig_proj, double &  mean_error, int  numIMG, double  lambda, int  act_proj, const FileName &  fn_ctf, const MultidimArray< int > *  maskPtr, bool  refine  )

virtual

Run a single step of ART. An ART iteration is compound of as many steps as projections, this function runs a single step of the process. In ART the pointer to the output volume must point to the same vol_in, while in SIRT it should point to a second volume. The read projection must be provided to the algorithm but the rest of projections are output by the routine. The mean error is also an output. numIMG is a normalizing factor to be used in SIRT, if you are running pure ART then this factor should be 1.

The symmetry matrix from which the view is derived must be given in sym_no.

Reimplemented from ARTReconsBase.

Definition at line 267 of file art_crystal.cpp.

{
    // Only works for blob volumes .............................................
    if (artPrm.basis.type != Basis::blobs)
        REPORT_ERROR(ERR_VALUE_INCORRECT,
                     "This function only works with blob volumes");

    // Compute lattice vectors to be used ......................................
    Matrix1D<double> aint, bint, shift;
    aint.resize(2);
    bint.resize(2);
    shift.resize(3);
    shift.initZeros();

    symmetrizeCrystalVectors(aint, bint, shift, this->space_group, sym_no,
                               this->aint, this->bint);
    // Project structure .......................................................
    // The correction image is reused in this call to store the normalizing
    // projection, ie, the projection of an all-1 volume

    project_Crystal_Volume(vol_in, artPrm.basis, theo_proj,
                           corr_proj, YSIZE(read_proj()), XSIZE(read_proj()),
                           read_proj.rot(), read_proj.tilt(), read_proj.psi(), shift,
                           aint, bint, *(artPrm.D), *(artPrm.Dinv), this->unit_cell_mask,
                           FORWARD, artPrm.eq_mode);
    double shift_X, shift_Y;

    //   #define DEBUG_SHIFT
#ifdef DEBUG_SHIFT

    Matrix2D<double> A(3, 3);
    A.initIdentity();
    dMij(A, 0, 2) =  8;
    dMij(A, 1, 2) =  -5;
    std::cout << "A" << A;
    //move read_proj
    FOR_ALL_DIRECT_ELEMENTS_IN_MATRIX2D(theo_proj())
    dMij(theo_proj(), i, j) = dMij(read_proj(), i, j);
    applyGeometry(IMGMATRIX(read_proj), A, IMGMATRIX(theo_proj), IS_NOT_INV, WRAP);
#endif
#undef DEBUG_SHIFT

    if (artPrm.ref_trans_after != -1    &&
        imagen_no > artPrm.ref_trans_after && imagen_no != 0)
    {
        CorrelationAux aux;
        calculate_and_find_correlation_max_proj(read_proj, theo_proj,
                                                alig_proj,
                                                shift_X, shift_Y,
                                                artPrm.ref_trans_step,
                                                artPrm.ref_trans_after,
                                                imagen_no, aux);

        // Apply correction
        Matrix2D<double> Correction(3, 3);
        alig_proj().resize(read_proj());
        Correction.initIdentity();

        dMij(Correction, 0, 2) =  - shift_X;
        dMij(Correction, 1, 2) =  - shift_Y;
        //copy theo_proj to a temporal matrix
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(theo_proj())
        dAij(alig_proj(), i, j) = dAij(read_proj(), i, j);
        applyGeometry(xmipp_transformation::LINEAR, IMGMATRIX(alig_proj), IMGMATRIX(read_proj), Correction,
                      xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
    }
    // Now compute differences .................................................
    double applied_lambda = lambda / numIMG; // In ART mode, numIMG=1
    mean_error = 0;
    diff_proj().resize(read_proj());

    FOR_ALL_ELEMENTS_IN_ARRAY2D(IMGMATRIX(read_proj))
    {
        // Compute difference image and error
        IMGPIXEL(diff_proj, i, j) = IMGPIXEL(read_proj, i, j) - IMGPIXEL(theo_proj, i, j);
        mean_error += IMGPIXEL(diff_proj, i, j) * IMGPIXEL(diff_proj, i, j);

        // Compute the correction image
        if (ABS(IMGPIXEL(corr_proj, i, j)) < 1)
            IMGPIXEL(corr_proj, i, j) = SGN(IMGPIXEL(corr_proj, i, j));
        IMGPIXEL(corr_proj, i, j) =
            applied_lambda * IMGPIXEL(diff_proj, i, j) / IMGPIXEL(corr_proj, i, j);
    }
    mean_error /= XSIZE(diff_proj()) * YSIZE(diff_proj());

    // Backprojection of correction plane ......................................
    project_Crystal_Volume(*vol_out, artPrm.basis, theo_proj,
                           corr_proj, YSIZE(read_proj()), XSIZE(read_proj()),
                           read_proj.rot(), read_proj.tilt(), read_proj.psi(), shift,
                           aint, bint, *(artPrm.D), *(artPrm.Dinv), this->unit_cell_mask,
                           BACKWARD, artPrm.eq_mode);
}

Member Data Documentation

a

Matrix1D<double> CrystalARTRecons::a

First lattice vector (BCC units)

Definition at line 66 of file art_crystal.h.

a_mag

double CrystalARTRecons::a_mag

First lattice vector module (user supplies it in arbitrary units but it is divided by sampling as soon as it is initialized)

Definition at line 49 of file art_crystal.h.

ai

Matrix1D<double> CrystalARTRecons::ai

ai=aint/2 as double numbers

Definition at line 74 of file art_crystal.h.

aint

Matrix1D<double> CrystalARTRecons::aint

First lattice vector approximated to integer numbers (BCC units)

Definition at line 70 of file art_crystal.h.

ang_a2b_deg

double CrystalARTRecons::ang_a2b_deg

angle from a to b (degrees)

Definition at line 54 of file art_crystal.h.

ang_x2a_deg

double CrystalARTRecons::ang_x2a_deg

angle from x axis to a (degrees)

Definition at line 56 of file art_crystal.h.

avox

Matrix1D<double> CrystalARTRecons::avox

First lattice vector (voxel units)

Definition at line 58 of file art_crystal.h.

b

Matrix1D<double> CrystalARTRecons::b

Second lattice vector (BCC units)

Definition at line 68 of file art_crystal.h.

b_mag

double CrystalARTRecons::b_mag

Second lattice vector module (user supplies it in arbitrary units but it is divided by sampling as soon as it is initialized)

Definition at line 52 of file art_crystal.h.

bi

Matrix1D<double> CrystalARTRecons::bi

bi=aint/2 as double numbers

Definition at line 76 of file art_crystal.h.

bint

Matrix1D<double> CrystalARTRecons::bint

Second lattice vector approximated to integer numbers (BCC units)

Definition at line 72 of file art_crystal.h.

bvox

Matrix1D<double> CrystalARTRecons::bvox

Second lattice vector (voxel units)

Definition at line 60 of file art_crystal.h.

fill_space

bool CrystalARTRecons::fill_space

Fill space, repeat unit cell all over the space.

Definition at line 62 of file art_crystal.h.

space_group

int CrystalARTRecons::space_group

space_group

Definition at line 78 of file art_crystal.h.

unit_cell_mask

MultidimArray<int> CrystalARTRecons::unit_cell_mask

Unit cell mask. This mask is defined in the BCC space and it represent a parallelogram defined by points (-a-b)/2, (-a+b)/2, (a-b)/2, (a+b)/2. The reconstruction will be only performed for the basis inside this mask. Be careful that this is a 2D mask for a 3D reconstruction.

Definition at line 85 of file art_crystal.h.

---

The documentation for this class was generated from the following files:

• xmipp/libraries/reconstruction/art_crystal.h
• xmipp/libraries/reconstruction/art_crystal.cpp