SymList Class Reference

#include <symmetries.h>

Collaboration diagram for SymList:

Public Member Functions
	SymList ()
	SymList (const FileName &fn_sym, double accuracy=SYM_ACCURACY)
bool	isSymmetryGroup (FileName fn_sym, int &pgGroup, int &pgOrder)
void	fillSymmetryClass (const FileName &symmetry, int pgGroup, int pgOrder, std::vector< std::string > &fileContent)
void	getMatrices (int i, Matrix2D< double > &L, Matrix2D< double > &R, bool homogeneous=true) const
void	setMatrices (int i, const Matrix2D< double > &L, const Matrix2D< double > &R)
void	getShift (int i, Matrix1D< double > &shift) const
void	setShift (int i, const Matrix1D< double > &shift)
void	addShift (const Matrix1D< double > &shift)
int	readSymmetryFile (FileName fn_sym, double accuracy=SYM_ACCURACY)
void	addMatrices (const Matrix2D< double > &L, const Matrix2D< double > &R, int chain_length)
void	computeSubgroup (double accuracy=SYM_ACCURACY)
int	symsNo () const
int	spaceGroup () const
int	trueSymsNo () const
int	crystallographicSpaceGroup (double mag_a, double mag_b, double ang_a2b_deg) const
double	nonRedundantProjectionSphere (int pgGroup, int pgOrder)
double	computeDistance (double rot1, double tilt1, double psi1, double &rot2, double &tilt2, double &psi2, bool projdir_mode, bool check_mirrors, bool object_rotation=false, bool write_mirrors=true)
void	computeDistance (MetaData &md, bool projdir_mode, bool check_mirrors, bool object_rotation=false)
void	breakSymmetry (double rot1, double tilt1, double psi1, double &rot2, double &tilt2, double &psi2)

Public Attributes
Matrix2D< double >	__L
Matrix2D< double >	__R
Matrix2D< double >	__shift
Matrix1D< int >	__chain_length
int	true_symNo
int	__sym_elements

Detailed Description

Symmetry List class. Internally the symmetry list class is implemented as a single 2D matrix, where every 4 rows (remember that in 3D the geometrical transformation matrices are 4x4) comprise a symmetry matrix. Access, and ways to modify the symmetry list are supplied. Remind that any symmetry is expressed in terms of two matrices L and R, so that any Euler matrix must be transformed by L*Euler*R resulting into a new perspective of the volume which is equivalent to the original one.

The typical use of the symmetry lists is to read the symmetry file, and do nothing else but reading matrices from it.

The symmetry file format is

#This is a comment
# The following line is a 6-fold rotational symmetry axis along Z-axis.
# The fold is the number of times that the volume can be rotated along
# the symmetry axis giving the same view from different view points.
# the structure for the rotational axis is
# rot_axis      <fold> <X0> <Y0> <Z0>
# mirror_plane         <X0> <Y0> <Z0>
rot_axis      6 0 0 1
mirror_plane    0 0 1

Definition at line 138 of file symmetries.h.

Constructor & Destructor Documentation

SymList() [1/2]

SymList::SymList ( )

inline

Create an empty list. The 2D matrices are 0x0. \ Ex: SymList SL;

Definition at line 162 of file symmetries.h.

    {
        __sym_elements = true_symNo = space_group = 0;
    }

SymList() [2/2]

SymList::SymList ( const FileName &  fn_sym, double  accuracy = SYM_ACCURACY  )

inline

Create Symmetry List from a Symmetry file. All the subgroup elements are computed automatically. \ Ex: SymList SL("sym.txt");

Definition at line 170 of file symmetries.h.

    {
        readSymmetryFile(fn_sym, accuracy);
    }

Member Function Documentation

addMatrices()

void SymList::addMatrices	(	const Matrix2D< double > &	L,
		const Matrix2D< double > &	R,
		int	chain_length
	)

Add symmetry matrices to the symmetry list. The given matrix must specify a point of view equivalent to the actual point of view. The matrices are added to the subgroup generator but the subgroup is not updated, you must do it manually using computeSubgroup. What is more, the subgroup after the insertion is corrupted.

The chain length is the number of single matrices multiplication of which the inserted one is compound.

Definition at line 418 of file symmetries.cpp.

{
    if (MAT_XSIZE(L) != 4 || MAT_YSIZE(L) != 4 || MAT_XSIZE(R) != 4 || MAT_YSIZE(R) != 4)
        REPORT_ERROR(ERR_MATRIX_SIZE, "SymList::add_matrix: Transformation matrix is not 4x4");
    if (trueSymsNo() == symsNo())
    {
        __L.resize(MAT_YSIZE(__L) + 4, 4);
        __R.resize(MAT_YSIZE(__R) + 4, 4);
        __chain_length.resize(__chain_length.size() + 1);
    }

    setMatrices(true_symNo, L, R);
    __chain_length(__chain_length.size() - 1) = chain_length;
    true_symNo++;
}

addShift()

void SymList::addShift

(

const Matrix1D< double > &

shift

)

Add shift. Add a shift vector to the shift matrix. An exception is thrown if the input vector is not a 3x1 vector.

Definition at line 408 of file symmetries.cpp.

{
    if (shift.size() != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "SymList::add_shift: Shift vector is not 3x1");
    int i = MAT_YSIZE(__shift);
    __shift.resize(i + 1, 3);
    setShift(i, shift);
}

breakSymmetry()

void SymList::breakSymmetry	(	double	rot1,
		double	tilt1,
		double	psi1,
		double &	rot2,
		double &	tilt2,
		double &	psi2
	)

Return equivalent set of angles for a given symmetry

Definition at line 1253 of file symmetries.cpp.

{
    Matrix2D<double> E1;
    Euler_angles2matrix(rot1, tilt1, psi1, E1, true);
    static bool doRandomize=true;
    Matrix2D<double>  L(3, 3), R(3, 3);  // A matrix from the list

    int i;
    if (doRandomize)
    {
        srand ( time(NULL) );
        doRandomize=false;
    }
    int symOrder = symsNo()+1;
    //std::cerr << "DEBUG_ROB: symOrder: " << symOrder << std::endl;
    i = rand() % symOrder;//59+1
    //std::cerr << "DEBUG_ROB: i: " << i << std::endl;
    if (i < symOrder-1)
    {
        getMatrices(i, L, R);
        //std::cerr  << R << std::endl;
        Euler_matrix2angles(E1 * R, rot2, tilt2, psi2);
    }
    else
        {
        //std::cerr << "else" <<std::endl;
        rot2=rot1; tilt2=tilt1;psi2=psi1;
        }
//    if (rot2==0)
//:     std::cerr << "rot2  is zero " << i << R << L << std::endl;
}

computeDistance() [1/2]

double SymList::computeDistance	(	double	rot1,
		double	tilt1,
		double	psi1,
		double &	rot2,
		double &	tilt2,
		double &	psi2,
		bool	projdir_mode,
		bool	check_mirrors,
		bool	object_rotation = false,
		bool	write_mirrors = true
	)

Check symmetries. Given two sets of angles, modify set2 to be as close to set1 as possible considering the symmetries. Distances are measured over the sphere. If the projdir_mode is set to true, then the similarity is measured only using the projection direction. The function returns the minimum distance between the two angle sets after the symmetry is considered. If check_mirrors is set, up-down corrections are also considered. Object_rotation controls the way that symmetry is applied (LR if object_rotation=false or RL if object_rotation=true). Normally it is set to false.

Definition at line 1188 of file symmetries.cpp.

{
    Matrix2D<double> E1, E2;
    Euler_angles2matrix(rot1, tilt1, psi1, E1, false);

    int imax = symsNo() + 1;
    Matrix2D<double>  L(3, 3), R(3, 3);  // A matrix from the list
    double best_ang_dist = 3600;
    double best_rot2=0, best_tilt2=0, best_psi2=0;

    for (int i = 0; i < imax; i++)
    {
        double rot2p, tilt2p, psi2p;
        if (i == 0)
        {
            rot2p = rot2;
            tilt2p = tilt2;
            psi2p = psi2;
        }
        else
        {
            getMatrices(i - 1, L, R, false);
            if (object_rotation)
                Euler_apply_transf(R, L, rot2, tilt2, psi2, rot2p, tilt2p, psi2p);
            else
                Euler_apply_transf(L, R, rot2, tilt2, psi2, rot2p, tilt2p, psi2p);
        }

        double ang_dist = Euler_distanceBetweenAngleSets_fast(E1,rot2p, tilt2p, psi2p,
                          projdir_mode, E2);

        if (ang_dist < best_ang_dist)
        {
            best_rot2 = rot2p;
            best_tilt2 = tilt2p;
            best_psi2 = psi2p;
            best_ang_dist = ang_dist;
        }

        if (check_mirrors)
        {
            double rot2pm, tilt2pm, psi2pm;
            Euler_mirrorY(rot2p, tilt2p, psi2p, rot2pm, tilt2pm, psi2pm);
            double ang_dist_mirror = Euler_distanceBetweenAngleSets_fast(E1,
                                     rot2pm, tilt2pm, psi2pm,projdir_mode, E2);

            if (ang_dist_mirror < best_ang_dist)
            {
                best_rot2 = write_mirrors ? rot2pm : rot2p;
                best_tilt2 = write_mirrors ? tilt2pm : tilt2p;
                best_psi2 = write_mirrors ? psi2pm : psi2p;
                best_ang_dist = ang_dist_mirror;
            }

        }
    }
    rot2 = best_rot2;
    tilt2 = best_tilt2;
    psi2 = best_psi2;
    return best_ang_dist;
}

computeDistance() [2/2]

void SymList::computeDistance	(	MetaData &	md,
		bool	projdir_mode,
		bool	check_mirrors,
		bool	object_rotation = false
	)

auxiliary function that calls double computeDistance(double rot1, double tilt1, double psi1, double &rot2, double &tilt2, double &psi2, bool projdir_mode, bool check_mirrors, bool object_rotation=false); from a metadata. Input metadata _angleRot _angleRot2 _angleTilt _angleTilt2 _anglePsi _anglePsi2 _anglePsiDiff _image

output metadata

_angleRot _angleRot2 _angleRotDiff _angleTilt _angleTilt2 _angleTiltDiff _anglePsi _anglePsi2 _anglePsiDiff _angleDiff _image

Definition at line 1157 of file symmetries.cpp.

{
    double rot1, tilt1, psi1;
    double rot2, tilt2, psi2;
    double angDistance;
    for (const auto& row : md)
    {
        row.getValue(MDL_ANGLE_ROT, rot1);
        row.getValue(MDL_ANGLE_ROT2, rot2);

        row.getValue(MDL_ANGLE_TILT, tilt1);
        row.getValue(MDL_ANGLE_TILT2, tilt2);

        row.getValue(MDL_ANGLE_PSI, psi1);
        row.getValue(MDL_ANGLE_PSI2, psi2);

        angDistance=computeDistance( rot1,  tilt1,  psi1,
                                     rot2,  tilt2,  psi2,
                                     projdir_mode,  check_mirrors,
                                     object_rotation);

        md.setValue(MDL_ANGLE_ROT_DIFF,rot1 - rot2, row.id());
        md.setValue(MDL_ANGLE_TILT_DIFF,tilt1 - tilt2, row.id());
        md.setValue(MDL_ANGLE_PSI_DIFF,psi1 - psi2, row.id());
        md.setValue(MDL_ANGLE_DIFF,angDistance, row.id());
    }

}

computeSubgroup()

void SymList::computeSubgroup

(

double

accuracy = SYM_ACCURACY

)

Compute subgroup for this structure. After adding or setting a matrix, the subgroup information is lost, you must recalculate it using this function. The different matrices are multiplied until no more different matrices are produced. The accuracy is used in order to compare when two matrix elements are the same.

So far, all the shifts associated to generated matrices are set to 0

Definition at line 467 of file symmetries.cpp.

{
    Matrix2D<double> I(4, 4);
    I.initIdentity();
    Matrix2D<double> L1(4, 4), R1(4, 4), L2(4, 4), R2(4, 4), newL(4, 4), newR(4, 4),identity(4,4);
    Matrix2D<int>    tried(true_symNo, true_symNo);
    Matrix1D<double> shift(3);
    shift.initZeros();
    int i, j;
    int new_chain_length;
    identity.initIdentity();
    while (found_not_tried(tried, i, j, true_symNo))
    {
        tried(i, j) = 1;

        // Form new symmetry matrices
        // if (__chain_length(i)+__chain_length(j)>__sym_elements+2) continue;

        getMatrices(i, L1, R1);
        getMatrices(j, L2, R2);
        newL = L1 * L2;
        newR = R1 * R2;
        new_chain_length = __chain_length(i) + __chain_length(j);

        //if (newL.isIdentity() && newR.isIdentity()) continue;
        //rounding error make newR different from identity
        if (newL.equal(identity, accuracy) &&
            newR.equal(identity, accuracy))
            continue;

        // Try to find it in current ones
        bool found;
        found = false;
        for (int l = 0; l < symsNo(); l++)
        {
            getMatrices(l, L1, R1);
            if (newL.equal(L1, accuracy) && newR.equal(R1, accuracy))
            {
                found = true;
                break;
            }
        }

        if (!found)
        {
            //#define DEBUG
#ifdef DEBUG
            static int kjhg=0;
            /* std::cout << "Matrix size " << tried.Xdim() << " "
             << "trying " << i << " " << j << " "
             << "chain length=" << new_chain_length << std::endl;
             std::cout << "Result R Sh\n" << newR << shift;
             */
            //std::cerr << "shift" << __shift <<std::endl;
            std::cerr << "newR "  << kjhg++ << "\n" << newR <<std::endl;
#endif

            addMatrices(newL, newR, new_chain_length);
            addShift(shift);
            tried.resize(MAT_YSIZE(tried) + 1, MAT_XSIZE(tried) + 1);
        }
    }
#ifdef DEBUG
    std::cerr << "__R" << __R <<std::endl;
#endif
#undef DEBUG
}

crystallographicSpaceGroup()

int SymList::crystallographicSpaceGroup	(	double	mag_a,
		double	mag_b,
		double	ang_a2b_deg
	)		const

Guess Crystallographic space group. Return the http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-getgen number. So far it has only been implemented for P1 (1), P2122 & P2212 (17), P4 (75), P4212 (90) and P6 (168).

Mag_a and Mag_b are the crystal vector magnitude.

Guess Crystallographic space group. Return the group number http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-getgen. So far it has only been implemented for P1 (1), P2_122 (17) , P22_12, P4 (75), P4212 (90) and P6 (168)

Definition at line 541 of file symmetries.cpp.

{

    switch (space_group)
    {
    case sym_undefined:
    case sym_P1:
        return(space_group);
    case sym_P4:
        if (fabs((mag_a - mag_b)) > XMIPP_EQUAL_ACCURACY ||
            fabs(ang_a2b_deg - 90) > XMIPP_EQUAL_ACCURACY)
            std::cerr << "\nWARNING: P42 but mag_a != mag_b\n"
            << " or ang_a2b !=90" << std::endl;
        return(space_group);
        break;
    case sym_P2_122:
        if (fabs((mag_a - mag_b)) > XMIPP_EQUAL_ACCURACY ||
            fabs(ang_a2b_deg - 90) > XMIPP_EQUAL_ACCURACY)
            std::cerr << "\nWARNING: P2_122 but mag_a != mag_b\n"
            << " or ang_a2b !=90" << std::endl;
        return(space_group);
        break;
    case sym_P22_12:
        if (fabs((mag_a - mag_b)) > XMIPP_EQUAL_ACCURACY ||
            fabs(ang_a2b_deg - 90) > XMIPP_EQUAL_ACCURACY)
            std::cerr << "\nWARNING: P22_12 but mag_a != mag_b\n"
            << " or ang_a2b !=90" << std::endl;
        return(space_group);
        break;
    case sym_P42_12:
        if (fabs((mag_a - mag_b)) > XMIPP_EQUAL_ACCURACY ||
            fabs(ang_a2b_deg - 90) > XMIPP_EQUAL_ACCURACY)
            std::cerr << "\nWARNING: P42_12 but mag_a != mag_b\n"
            << " or ang_a2b !=90" << std::endl;
        return(space_group);
        break;
    case sym_P6:
        if (fabs((mag_a - mag_b)) > XMIPP_EQUAL_ACCURACY ||
            fabs(ang_a2b_deg - 120.) > XMIPP_EQUAL_ACCURACY)
        {
            std::cerr << "\nWARNING: marked as P6 but mag_a != mag_b\n"
            << "or ang_a2b !=120" << std::endl;
            std::cerr << "\nWARNING: P1 is assumed\n";
            return(sym_P1);
        }
        else
            return(space_group);
        break;
    default:
        std::cerr << "\n Congratulations: you have found a bug in the\n"
        << "routine crystallographic_space_group or\n"
        << "You have called to this rotuine BEFORE reading\n"
        << "the symmetry info" << std::endl;
        exit(0);
        break;
    }//switch(space_group)  end

}//crystallographicSpaceGroup end

fillSymmetryClass()

void SymList::fillSymmetryClass	(	const FileName &	symmetry,
		int	pgGroup,
		int	pgOrder,
		std::vector< std::string > &	fileContent
	)

fill fileContect with symmetry information

Definition at line 879 of file symmetries.cpp.

{
    std::ostringstream line1;
    std::ostringstream line2;
    std::ostringstream line3;
    std::ostringstream line4;
    if (pgGroup == pg_CN)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
    }
    else if (pgGroup == pg_CI)
    {
        line1 << "inversion ";
    }
    else if (pgGroup == pg_CS)
    {
        line1 << "mirror_plane 0 0 1";
    }
    else if (pgGroup == pg_CNV)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
        line2 << "mirror_plane 0 1 0";
    }
    else if (pgGroup == pg_CNH)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
        line2 << "mirror_plane 0 0 1";
    }
    else if (pgGroup == pg_SN)
    {
        int order = pgOrder / 2;
        if(2*order != pgOrder)
        {
            std::cerr << "ERROR: order for SN group must be even" << std::endl;
            exit(0);
        }
        line1 << "rot_axis " << order << " 0 0 1";
        line2 << "inversion ";
    }
    else if (pgGroup == pg_DN)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
        line2 << "rot_axis " << "2" << " 1 0 0";
    }
    else if (pgGroup == pg_DNV)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
        line2 << "rot_axis " << "2" << " 1 0 0";
        line3 << "mirror_plane 1 0 0";
    }
    else if (pgGroup == pg_DNH)
    {
        line1 << "rot_axis " << pgOrder << " 0 0 1";
        line2 << "rot_axis " << "2" << " 1 0 0";
        line3 << "mirror_plane 0 0 1";
    }
    else if (pgGroup == pg_T)
    {
        line1 << "rot_axis " << "3" << "  0. 0. 1.";
        line2 << "rot_axis " << "2" << " 0. 0.816496 0.577350";
    }
    else if (pgGroup == pg_TD)
    {
        line1 << "rot_axis " << "3" << "  0. 0. 1.";
        line2 << "rot_axis " << "2" << " 0. 0.816496 0.577350";
        line3 << "mirror_plane 1.4142136 2.4494897 0.0000000";
    }
    else if (pgGroup == pg_TH)
    {
        line1 << "rot_axis " << "3" << "  0. 0. 1.";
        line2 << "rot_axis " << "2" << " 0. -0.816496 -0.577350";
        line3 << "inversion";
    }
    else if (pgGroup == pg_O)
    {
        line1 << "rot_axis " << "3" << "  .5773502  .5773502 .5773502";
        line2 << "rot_axis " << "4" << " 0 0 1";
    }
    else if (pgGroup == pg_OH)
    {
        line1 << "rot_axis " << "3" << "  .5773502  .5773502 .5773502";
        line2 << "rot_axis " << "4" << " 0 0 1";
        line3 << "mirror_plane 0 1 1";
    }
    else if (pgGroup == pg_I || pgGroup == pg_I2)
    {
        line1 << "rot_axis 2  0 0 1";
        line2 << "rot_axis 5  0.525731114  0 0.850650807";
        line3 << "rot_axis 3  0 0.356822076 0.934172364";
    }
    else if (pgGroup == pg_I1)
    {
        line1 << "rot_axis 2  1      0        0";
        line2 << "rot_axis 5 0.85065080702670 0 -0.5257311142635";
        line3 << "rot_axis 3 0.9341723640 0.3568220765 0";
    }
    else if (pgGroup == pg_I3)
    {
        line1 << "rot_axis 2  -0.5257311143 0 0.8506508070";
        line3 << "rot_axis 5  0. 0. 1.";
        line2 << "rot_axis 3  -0.4911234778630044, 0.3568220764705179, 0.7946544753759428";
    }
    else if (pgGroup == pg_I4)
    {
        line1 << "rot_axis 2  0.5257311143 0 0.8506508070";
        line3 << "rot_axis 5  0.8944271932547096 0 0.4472135909903704";
        line2 << "rot_axis 3  0.4911234778630044 0.3568220764705179 0.7946544753759428";
    }
    else if (pgGroup == pg_I5)
    {
        std::cerr << "ERROR: Symmetry pg_I5 not implemented" << std::endl;
        exit(0);
    }
    else if (pgGroup == pg_IH || pgGroup == pg_I2H)
    {
        line1 << "rot_axis 2  0 0 1";
        line2 << "rot_axis 5  0.525731114  0 0.850650807";
        line3 << "rot_axis 3  0 0.356822076 0.934172364";
        line4 << "mirror_plane 1 0 0";
    }
    else if (pgGroup == pg_I1H)
    {
        line1 << "rot_axis 2  1      0        0";
        line2 << "rot_axis 5 0.85065080702670 0 -0.5257311142635";
        line3 << "rot_axis 3 0.9341723640 0.3568220765 0";
        line4 << "mirror_plane 0 0 -1";
    }
    else if (pgGroup == pg_I3H)
    {
        line1 << "rot_axis 2  -0.5257311143 0 0.8506508070";
        line3 << "rot_axis 5  0. 0. 1.";
        line2 << "rot_axis 3  -0.4911234778630044, 0.3568220764705179, 0.7946544753759428";
        line4 << "mirror_plane 0.850650807 0  0.525731114";
    }
    else if (pgGroup == pg_I4H)
    {
        line1 << "rot_axis 2  0.5257311143 0 0.8506508070";
        line3 << "rot_axis 5  0.8944271932547096 0 0.4472135909903704";
        line2 << "rot_axis 3  0.4911234778630044 0.3568220764705179 0.7946544753759428";
        line4 << "mirror_plane 0.850650807 0 -0.525731114";
    }
    else if (pgGroup == pg_I5H)
    {
        std::cerr << "ERROR: Symmetry pg_I5H not implemented" << std::endl;
        exit(0);
    }
    else
    {
        std::cerr << "ERROR: Symmetry " << symmetry  << "is not known" << std::endl;
        exit(0);
    }
    if (line1.str().size()>0)
        fileContent.push_back(line1.str());
    if (line2.str().size()>0)
        fileContent.push_back(line2.str());
    if (line3.str().size()>0)
        fileContent.push_back(line3.str());
    if (line4.str().size()>0)
        fileContent.push_back(line4.str());
    //#define DEBUG5
#ifdef DEBUG5

    for (int n=0; n<fileContent.size(); n++)
        std::cerr << fileContent[n] << std::endl;
    std::cerr << "fileContent.size()" << fileContent.size() << std::endl;
#endif
    #undef DEBUG5
}

getMatrices()

void SymList::getMatrices	(	int	i,
		Matrix2D< double > &	L,
		Matrix2D< double > &	R,
		bool	homogeneous = true
	)		const

Get matrices from the symmetry list. The number of matrices inside the list is given by symsNo. This function return the 4x4 (homogeneous=true) or 3x3 (homogeneous=false) transformation matrices associated to the one in the list which occupies the position 'i'. The matrix numbering within the list starts at 0. The output transformation matrices is given as a pointer to gain speed. \ Ex:

for (i=0; i<SL.symsNo; i++) {
    SL.getMatrices(i,L,R);
    ...
}

Definition at line 348 of file symmetries.cpp.

{
    int k, kp, l;
    if (homogeneous)
    {
        L.initZeros(4, 4);
        R.initZeros(4, 4);
        for (k = 4 * i, kp=0; k < 4*i + 4; k++, kp++)
            for (l = 0; l < 4; l++)
            {
                dMij(L,kp, l) = dMij(__L,k, l);
                dMij(R,kp, l) = dMij(__R,k, l);
            }
    }
    else
    {
        L.initZeros(3, 3);
        R.initZeros(3, 3);
        for (k = 4 * i, kp=0; k < 4*i + 3; k++, kp++)
            for (l = 0; l < 3; l++)
            {
                dMij(L,kp, l) = dMij(__L,k, l);
                dMij(R,kp, l) = dMij(__R,k, l);
            }
    }
}

getShift()

void SymList::getShift	(	int	i,
		Matrix1D< double > &	shift
	)		const

Get shift. Returns the shift associated to a certain symmetry.

Definition at line 391 of file symmetries.cpp.

{
    shift.resize(3);
    XX(shift) = __shift(i, 0);
    YY(shift) = __shift(i, 1);
    ZZ(shift) = __shift(i, 2);
}

isSymmetryGroup()

bool SymList::isSymmetryGroup	(	FileName	fn_sym,
		int &	pgGroup,
		int &	pgOrder
	)

translate string fn_sym to symmetry group, return false is translation is not possible. See http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/Symmetry for details. It also fill the symmetry information

Definition at line 601 of file symmetries.cpp.

{
    char G1,G2,G3='\0',G4;
    char auxChar[3];
    //each case check length, check first letter, second, is number
    //Non a point group

    //remove path
    FileName fn_sym_tmp;
    fn_sym_tmp=fn_sym.removeDirectories();
    int mySize=fn_sym_tmp.size();
    bool return_true;
    return_true=false;
    auxChar[2]='\0';
    //size maybe 4 because n maybe a 2 digit number
    if(mySize>4 || mySize<1)
    {
        pgGroup=-1;
        pgOrder=-1;
        return false;
    }
    //get the group character by character
    G1=toupper((fn_sym_tmp.c_str())[0]);
    G2=toupper((fn_sym_tmp.c_str())[1]);
    if (mySize > 2)
    {
        G3=toupper((fn_sym_tmp.c_str())[2]);
        if(mySize > 3)
            G4=toupper((fn_sym.c_str())[3]);
    }
    else
        G4='\0';
    //CN
    if (mySize==2 && G1=='C' && isdigit(G2))
    {
        pgGroup=pg_CN;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    if (mySize==3 && G1=='C' && isdigit(G2) && isdigit(G3))
    {
        pgGroup=pg_CN;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //CI
    else if (mySize==2 && G1=='C' && G2=='I')
    {
        pgGroup=pg_CI;
        pgOrder=-1;
        return_true=true;
    }
    //CS
    else if (mySize==2 && G1=='C' && G2=='S')
    {
        pgGroup=pg_CS;
        pgOrder=-1;
        return_true=true;
    }
    //CNH
    else if (mySize==3 && G1=='C' && isdigit(G2) && G3=='H')
    {
        pgGroup=pg_CNH;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    else if (mySize==4 && G1=='C' && isdigit(G2) && isdigit(G3) && G4=='H')
    {
        pgGroup=pg_CNH;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //CNV
    else if (mySize==3 && G1=='C' && isdigit(G2) && G3=='V')
    {
        pgGroup=pg_CNV;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    else if (mySize==4 && G1=='C' && isdigit(G2) && isdigit(G3) && G4=='V')
    {
        pgGroup=pg_CNV;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //SN
    else if (mySize==2 && G1=='S' && isdigit(G2) )
    {
        pgGroup=pg_SN;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    else if (mySize==3 && G1=='S' && isdigit(G2) && isdigit(G3) )
    {
        pgGroup=pg_SN;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //DN
    else if (mySize==2 && G1=='D' && isdigit(G2) )
    {
        pgGroup=pg_DN;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    if (mySize==3 && G1=='D' && isdigit(G2) && isdigit(G3))
    {
        pgGroup=pg_DN;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //DNV
    else if (mySize==3 && G1=='D' && isdigit(G2) && G3=='V')
    {
        pgGroup=pg_DNV;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    else if (mySize==4 && G1=='D' && isdigit(G2) && isdigit(G3) && G4=='V')
    {
        pgGroup=pg_DNV;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //DNH
    else if (mySize==3 && G1=='D' && isdigit(G2) && G3=='H')
    {
        pgGroup=pg_DNH;
        pgOrder=int(G2)-48;
        return_true=true;
    }
    else if (mySize==4 && G1=='D' && isdigit(G2) && isdigit(G3) && G4=='H')
    {
        pgGroup=pg_DNH;
        auxChar[0]=G2;
        auxChar[1]=G3;
        pgOrder=atoi(auxChar);
        return_true=true;
    }
    //T
    else if (mySize==1 && G1=='T')
    {
        pgGroup=pg_T;
        pgOrder=-1;
        return_true=true;
    }
    //TD
    else if (mySize==2 && G1=='T' && G2=='D')
    {
        pgGroup=pg_TD;
        pgOrder=-1;
        return_true=true;
    }
    //TH
    else if (mySize==2 && G1=='T' && G2=='H')
    {
        pgGroup=pg_TH;
        pgOrder=-1;
        return_true=true;
    }
    //O
    else if (mySize==1 && G1=='O')
    {
        pgGroup=pg_O;
        pgOrder=-1;
        return_true=true;
    }
    //OH
    else if (mySize==2 && G1=='O'&& G2=='H')
    {
        pgGroup=pg_OH;
        pgOrder=-1;
        return_true=true;
    }
    //I
    else if (mySize==1 && G1=='I')
    {
        pgGroup=pg_I;
        pgOrder=-1;
        return_true=true;
    }
    //I1
    else if (mySize==2 && G1=='I'&& G2=='1')
    {
        pgGroup=pg_I1;
        pgOrder=-1;
        return_true=true;
    }
    //I2
    else if (mySize==2 && G1=='I'&& G2=='2')
    {
        pgGroup=pg_I2;
        pgOrder=-1;
        return_true=true;
    }
    //I3
    else if (mySize==2 && G1=='I'&& G2=='3')
    {
        pgGroup=pg_I3;
        pgOrder=-1;
        return_true=true;
    }
    //I4
    else if (mySize==2 && G1=='I'&& G2=='4')
    {
        pgGroup=pg_I4;
        pgOrder=-1;
        return_true=true;
    }
    //I5
    else if (mySize==2 && G1=='I'&& G2=='5')
    {
        pgGroup=pg_I5;
        pgOrder=-1;
        return_true=true;
    }
    //IH
    else if (mySize==2 && G1=='I'&& G2=='H')
    {
        pgGroup=pg_IH;
        pgOrder=-1;
        return_true=true;
    }
    //I1H
    else if (mySize==3 && G1=='I'&& G2=='1'&& G3=='H')
    {
        pgGroup=pg_I1H;
        pgOrder=-1;
        return_true=true;
    }
    //I2H
    else if (mySize==3 && G1=='I'&& G2=='2'&& G3=='H')
    {
        pgGroup=pg_I2H;
        pgOrder=-1;
        return_true=true;
    }
    //I3H
    else if (mySize==3 && G1=='I'&& G2=='3'&& G3=='H')
    {
        pgGroup=pg_I3H;
        pgOrder=-1;
        return_true=true;
    }
    //I4H
    else if (mySize==3 && G1=='I'&& G2=='4'&& G3=='H')
    {
        pgGroup=pg_I4H;
        pgOrder=-1;
        return_true=true;
    }
    //I5H
    else if (mySize==3 && G1=='I'&& G2=='5'&& G3=='H')
    {
        pgGroup=pg_I5H;
        pgOrder=-1;
        return_true=true;
    }
    //#define DEBUG7
#ifdef DEBUG7
    std::cerr << "pgGroup" << pgGroup << " pgOrder " << pgOrder << std::endl;
#endif
#undef DEBUG7

    return return_true;
}

nonRedundantProjectionSphere()

double SymList::nonRedundantProjectionSphere	(	int	pgGroup,
		int	pgOrder
	)

Return the area of the non redundant part of the projection sphere

Definition at line 1048 of file symmetries.cpp.

{
    if (pgGroup == pg_CN)
    {
        return 4.*PI/pgOrder;
    }
    else if (pgGroup == pg_CI)
    {
        return 4.*PI/2.;
    }
    else if (pgGroup == pg_CS)
    {
        return 4.*PI/2.;
    }
    else if (pgGroup == pg_CNV)
    {
        return 4.*PI/pgOrder/2;
    }
    else if (pgGroup == pg_CNH)
    {
        return 4.*PI/pgOrder/2;
    }
    else if (pgGroup == pg_SN)
    {
        return 4.*PI/pgOrder;
    }
    else if (pgGroup == pg_DN)
    {
        return 4.*PI/pgOrder/2;
    }
    else if (pgGroup == pg_DNV)
    {
        return 4.*PI/pgOrder/4;
    }
    else if (pgGroup == pg_DNH)
    {
        return 4.*PI/pgOrder/4;
    }
    else if (pgGroup == pg_T)
    {
        return 4.*PI/12;
    }
    else if (pgGroup == pg_TD)
    {
        return 4.*PI/24;
    }
    else if (pgGroup == pg_TH)
    {
        return 4.*PI/24;
    }
    else if (pgGroup == pg_O)
    {
        return 4.*PI/24;
    }
    else if (pgGroup == pg_OH)
    {
        return 4.*PI/48;
    }
    else if (pgGroup == pg_I || pgGroup == pg_I2)
    {
        return 4.*PI/60;
    }
    else if (pgGroup == pg_I1)
    {
        return 4.*PI/60;
    }
    else if (pgGroup == pg_I3)
    {
        return 4.*PI/60;
    }
    else if (pgGroup == pg_I4)
    {
        return 4.*PI/60;
    }
    else if (pgGroup == pg_I5)
    {
        return 4.*PI/60;
    }
    else if (pgGroup == pg_IH || pgGroup == pg_I2H)
    {
        return 4.*PI/120;
    }
    else if (pgGroup == pg_I1H)
    {
        return 4.*PI/120;
    }
    else if (pgGroup == pg_I3H)
    {
        return 4.*PI/120;
    }
    else if (pgGroup == pg_I4H)
    {
        return 4.*PI/120;
    }
    else if (pgGroup == pg_I5H)
    {
        return 4.*PI/120;
    }
    else
    {
        std::cerr << "ERROR: Symmetry group, order=" << pgGroup
        << " "
        <<  pgOrder
        << "is not known"
        << std::endl;
        exit(0);
    }
}

readSymmetryFile()

int SymList::readSymmetryFile	(	FileName	fn_sym,
		double	accuracy = SYM_ACCURACY
	)

Read a symmetry file into a symmetry list. The former symmetry list is overwritten with the new one. All the subgroup members are added to the list. If the accuracy is negative then the subgroup is not generated. return symmetry group \ Ex: SL.readSymmetryFile("sym.txt");

Definition at line 33 of file symmetries.cpp.

{
    int i, j;
    FILE *fpoii;
    char line[80];
    char *auxstr;
    double ang_incr, rot_ang;
    int  fold;
    Matrix2D<double> L(4, 4), R(4, 4);
    Matrix1D<double> axis(3), shift(3);
    int pgGroup = sym_undefined;
    int pgOrder;
    std::vector<std::string> fileContent;

    //check if reserved word

    // Open file ---------------------------------------------------------
    if ((fpoii = fopen(fn_sym.c_str(), "r")) == NULL)
    {
        //check if reserved word and return group and order
        if (isSymmetryGroup(fn_sym, pgGroup, pgOrder))
        {
            fillSymmetryClass(fn_sym, pgGroup, pgOrder,fileContent);
        }
        else
            REPORT_ERROR(ERR_IO_NOTOPEN, (std::string)"SymList::read_sym_file:Can't open file: "
                         + " or do not recognize symmetry group" + fn_sym);
    }
    else
    {
        while (fgets(line, 79, fpoii) != NULL)
        {
            if (line[0] == ';' || line[0] == '#' || line[0] == '\0')
                continue;
            fileContent.push_back(line);
        }
        fclose(fpoii);
    }

    //reset space_group
    space_group = 0;
    // Count the number of symmetries ------------------------------------
    true_symNo = 0;
    // count number of axis and mirror planes. It will help to identify
    // the crystallographic symmetry

    int no_axis, no_mirror_planes, no_inversion_points;
    no_axis = no_mirror_planes = no_inversion_points = 0;

    for (size_t n=0; n<fileContent.size(); n++)
    {
        strcpy(line,fileContent[n].c_str());
        auxstr = firstToken(line);
        if (auxstr == NULL)
        {
            std::cout << line;
            std::cout << "Wrong line in symmetry file, the line is skipped\n";
            continue;
        }
        if (strcmp(auxstr, "rot_axis") == 0)
        {
            auxstr = nextToken();
            fold = textToInteger(auxstr);
            true_symNo += (fold - 1);
            no_axis++;
        }
        else if (strcmp(auxstr, "mirror_plane") == 0)
        {
            true_symNo++;
            no_mirror_planes++;
        }
        else if (strcmp(auxstr, "inversion") == 0)
        {
            true_symNo += 1;
            no_inversion_points = 1;
        }
        else if (strcmp(auxstr, "P4212") == 0)
            true_symNo += 7;
        else if (strcmp(auxstr, "P2_122") == 0)
            true_symNo += 3;
        else if (strcmp(auxstr, "P22_12") == 0)
            true_symNo += 3;
    }
    // Ask for memory
    __L.resize(4*true_symNo, 4);
    __R.resize(4*true_symNo, 4);
    __shift.resize(true_symNo, 3);
    __chain_length.resize(true_symNo);
    __chain_length.initConstant(1);

    // Read symmetry parameters
    i = 0;
    shift.initZeros();

    for (size_t n=0; n<fileContent.size(); n++)
    {
        strcpy(line,fileContent[n].c_str());
        auxstr = firstToken(line);
        // Rotational axis ---------------------------------------------------
        if (strcmp(auxstr, "rot_axis") == 0)
        {
            auxstr = nextToken();
            fold = textToInteger(auxstr);
            auxstr = nextToken();
            XX(axis) = textToFloat(auxstr);
            auxstr = nextToken();
            YY(axis) = textToFloat(auxstr);
            auxstr = nextToken();
            ZZ(axis) = textToFloat(auxstr);
            ang_incr = 360. / fold;
            L.initIdentity();
            for (j = 1, rot_ang = ang_incr; j < fold; j++, rot_ang += ang_incr)
            {
                rotation3DMatrix(rot_ang, axis, R);
                setShift(i, shift);
                setMatrices(i++, L, R.transpose());
            }
            __sym_elements++;
            // inversion ------------------------------------------------------
        }
        else if (strcmp(auxstr, "inversion") == 0)
        {
            L.initIdentity();
            L(2, 2) = -1;
            R.initIdentity();
            R(0, 0) = -1.;
            R(1, 1) = -1.;
            R(2, 2) = -1.;
            setShift(i, shift);
            setMatrices(i++, L, R);
            __sym_elements++;
            // mirror plane -------------------------------------------------------------
        }
        else if (strcmp(auxstr, "mirror_plane") == 0)
        {
            auxstr = nextToken();
            XX(axis) = textToFloat(auxstr);
            auxstr = nextToken();
            YY(axis) = textToFloat(auxstr);
            auxstr = nextToken();
            ZZ(axis) = textToFloat(auxstr);
            L.initIdentity();
            L(2, 2) = -1;
            Matrix2D<double> A;
            alignWithZ(axis,A);
            A = A.transpose();
            R = A * L * A.inv();
            setShift(i, shift);
            L.initIdentity();
            setMatrices(i++, L, R);
            __sym_elements++;
            // P4212 -------------------------------------------------------------
        }
        else if (strcmp(auxstr, "P4212") == 0)
        {
            space_group = sym_P42_12;
            accuracy = -1; // Do not compute subgroup
            L.initIdentity();

            // With 0 shift
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = R(1, 1) = -1;
            R(2, 2) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(2, 2) = -1;
            R(0, 1) = R(1, 0) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(2, 2) = R(0, 1) = R(1, 0) = -1;
            setShift(i, shift);
            setMatrices(i++, L, R);

            // With 1/2 shift
            VECTOR_R3(shift, 0.5, 0.5, 0);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 1) = -1;
            R(1, 0) = R(2, 2) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(1, 0) = -1;
            R(0, 1) = R(2, 2) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = R(2, 2) = -1;
            R(1, 1) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(1, 1) = R(2, 2) = -1;
            R(0, 0) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);

            __sym_elements++;
        }
        else if (strcmp(auxstr, "P2_122") == 0)
        {
            space_group = sym_P2_122;
            accuracy = -1; // Do not compute subgroup
            L.initIdentity();

            // With 0 shift
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = -1;
            R(1, 1) = -1;
            R(2, 2) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);

            // With 1/2 shift
            VECTOR_R3(shift, 0.5, 0.0, 0.0);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = -1;
            R(1, 1) = 1;
            R(2, 2) = -1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = 1;
            R(1, 1) = -1;
            R(2, 2) = -1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            __sym_elements++;
        }
        else if (strcmp(auxstr, "P22_12") == 0)
        {
            space_group = sym_P22_12;
            accuracy = -1; // Do not compute subgroup
            L.initIdentity();

            // With 0 shift
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = -1;
            R(1, 1) = -1;
            R(2, 2) = 1;
            setShift(i, shift);
            setMatrices(i++, L, R);

            // With 1/2 shift
            VECTOR_R3(shift, 0.0, 0.5, 0.0);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = 1;
            R(1, 1) = -1;
            R(2, 2) = -1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            R.initZeros();
            R(3, 3) = 1;
            R(0, 0) = -1;
            R(1, 1) = 1;
            R(2, 2) = -1;
            setShift(i, shift);
            setMatrices(i++, L, R);
            __sym_elements++;
        }
    }

    if (accuracy > 0)
        computeSubgroup(accuracy);

    //possible crystallographic symmetry
    if (no_axis == 0 && no_mirror_planes == 0 && no_inversion_points == 0 &&
        true_symNo == 7 && space_group == sym_P42_12)
        space_group = sym_P42_12;
    else if (no_axis == 0 && no_mirror_planes == 0 && no_inversion_points == 0 &&
             true_symNo == 3 && space_group == sym_P2_122)
        space_group = sym_P2_122;
    else if (no_axis == 0 && no_mirror_planes == 0 && no_inversion_points == 0 &&
             true_symNo == 3 && space_group == sym_P22_12)
        space_group = sym_P22_12;
    // P4 and P6
    else if (no_axis == 1 && no_mirror_planes == 0 && no_inversion_points == 0 &&
             fabs(R(2, 2) - 1.) < XMIPP_EQUAL_ACCURACY &&
             fabs(R(0, 0) - R(1, 1)) < XMIPP_EQUAL_ACCURACY &&
             fabs(R(0, 1) + R(1, 0)) < XMIPP_EQUAL_ACCURACY)
    {
        switch (true_symNo)
        {
        case(5):
                        space_group = sym_P6;
            break;
        case(3):
                        space_group = sym_P4;
            break;
        default:
            space_group = sym_undefined;
            break;
        }//switch end
    }//end else if (no_axis==1 && no_mirror_planes== 0
    else if (no_axis == 0 && no_inversion_points == 0 && no_mirror_planes == 0)
        space_group = sym_P1;
    else
        space_group = sym_undefined;
    return pgGroup;
}

setMatrices()

void SymList::setMatrices	(	int	i,
		const Matrix2D< double > &	L,
		const Matrix2D< double > &	R
	)

Set a couple of matrices in the symmetry list. The number of matrices inside the list is given by symsNo. This function sets the 4x4 transformation matrices associated to the one in the list which occupies the position 'i'. The matrix numbering within the list starts at 0. \ Ex:

for (i=0; i<SL.symsNo; i++) {
    SL.set_matrix(i,L,R);
    ...
}

Definition at line 378 of file symmetries.cpp.

{
    int k, l;
    for (k = 4 * i; k < 4*i + 4; k++)
        for (l = 0; l < 4; l++)
        {
            __L(k, l) = L(k - 4 * i, l);
            __R(k, l) = R(k - 4 * i, l);
        }
}

setShift()

void SymList::setShift	(	int	i,
		const Matrix1D< double > &	shift
	)

Set shift. Set the shift associated to a certain symmetry.

Definition at line 399 of file symmetries.cpp.

{
    if (shift.size() != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "SymList::add_shift: Shift vector is not 3x1");
    __shift(i, 0) = XX(shift);
    __shift(i, 1) = YY(shift);
    __shift(i, 2) = ZZ(shift);
}

spaceGroup()

int SymList::spaceGroup ( ) const

inline

Return space group

Definition at line 275 of file symmetries.h.

    {
        return space_group;
    }

symsNo()

int SymList::symsNo ( ) const

inline

Number of symmetry matrices inside the structure. This is the number of all the matrices inside the subgroup. \ Ex:

for (i=0; i<SL.symsNo; i++) {
    SL.get_matrix(i,A);
    ...
}

Definition at line 268 of file symmetries.h.

    {
        return MAT_YSIZE(__L) / 4;
    }

trueSymsNo()

int SymList::trueSymsNo ( ) const

inline

Number of symmetry matrices which generated the structure. This is the number of the matrices which generated the structure, notice that it should be always less or equal to the total number of matrices in the subgroup.

Definition at line 283 of file symmetries.h.

    {
        return true_symNo;
    }

Member Data Documentation

__chain_length

Matrix1D<int> SymList::__chain_length

Definition at line 148 of file symmetries.h.

__L

Matrix2D<double> SymList::__L

Definition at line 146 of file symmetries.h.

__R

Matrix2D<double> SymList::__R

Definition at line 146 of file symmetries.h.

__shift

Matrix2D<double> SymList::__shift

Definition at line 147 of file symmetries.h.

__sym_elements

int SymList::__sym_elements

Definition at line 156 of file symmetries.h.

true_symNo

int SymList::true_symNo

Definition at line 153 of file symmetries.h.

---

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

• xmippCore/core/symmetries.h
• xmippCore/core/symmetries.cpp