Sampling Class Reference

#include <sampling.h>

Collaboration diagram for Sampling:

Classes
struct	Vertex

Public Types
typedef std::vector< Vertex >	Vect_angles

Public Member Functions
	Sampling ()
bool	operator== (const Sampling &op) const
void	computeSamplingPoints (bool only_half_sphere=true, double max_tilt=180, double min_tilt=0)
void	fillEdge (const Matrix1D< double > &starting_point, const Matrix1D< double > &ending_point, std::vector< Matrix1D< double > > &edge_vector, bool FLAG_END)
void	fillDistance (const Matrix1D< double > &starting_point, const Matrix1D< double > &ending_point, std::vector< Matrix1D< double > > &edge_vector, int number, bool only_half_spheree, double min_z=-10., double max_z=+10.)
void	fillLRRepository (void)
void	setSampling (double sampling)
void	setNoise (double deviation, int my_seed=-1)
void	setNeighborhoodRadius (double neighborhood)
void	removeRedundantPoints (const int symmetry, int sym_order)
void	removeRedundantPointsExhaustive (const int symmetry, int sym_order, bool only_half_sphere, double max_ang)
int	sortFunc (const Matrix1D< double > &a, const Matrix1D< double > &b)
void	createSymFile (const FileName &simFp, int symmetry, int sym_order)
void	createAsymUnitFile (const FileName &docfilename)
void	computeNeighbors (bool only_winner=false)
void	saveSamplingFile (const FileName &fn_base, bool write_vectors=true, bool write_sampling_sphere=false)
void	readSamplingFile (const FileName &infilename, bool read_vectors=true, bool read_sampling_sphere=false)
void	removePointsFarAwayFromExperimentalData ()
void	findClosestSamplingPoint (const MetaData &DFi, const FileName &output_file_root)
void	findClosestSamplingPoint (const FileName &FnexperimentalImages, const FileName &output_file_root)
void	findClosestExperimentalPoint ()
void	fillExpDataProjectionDirectionByLR (const MetaData &DFi)
void	fillExpDataProjectionDirectionByLR (const FileName &FnexperimentalImages)

Public Attributes
Vect_angles	vertices_angles
std::vector< Matrix1D< double > >	vertices_vectors
double	sampling_rate_rad
double	sampling_noise
size_t	number_of_samples
size_t	numberSamplesAsymmetricUnit
double	neighborhood_radius_rad
double	cos_neighborhood_radius
std::vector< std::vector< size_t > >	my_neighbors
std::vector< std::vector< size_t > >	my_exp_img_per_sampling_point
std::vector< std::vector< double > >	my_cross_correlation
std::vector< Matrix1D< double > >	sampling_points_vector
std::vector< Matrix1D< double > >	sampling_points_angles
std::vector< Matrix2D< double > >	R_repository
std::vector< Matrix2D< double > >	L_repository
std::vector< Matrix1D< double > >	exp_data_projection_direction_by_L_R
std::vector< FileName >	exp_data_fileNames
std::vector< Matrix1D< double > >	no_redundant_sampling_points_vector
std::vector< Matrix1D< double > >	no_redundant_sampling_points_angles
std::vector< size_t >	no_redundant_sampling_points_index
int	verbose
FileName	symmetry_file
SymList	SL

Detailed Description

Routines with sampling the direction Sphere A triangular grid based on an icosahedron was first introduced in a meteorological model by Sadourny et al. (1968) and Williamson (1969). The approach outlined here, especially the code implementation, is based on the work of Baumgardner (1995). http://www.wmo.int/pages/prog/www/DPS/Icosah.pdf

Definition at line 46 of file sampling.h.

Member Typedef Documentation

Vect_angles

typedef std::vector<Vertex> Sampling::Vect_angles

Typename to contain a list of vertices

Definition at line 58 of file sampling.h.

Constructor & Destructor Documentation

Sampling()

Sampling::Sampling

(

)

Default constructor. sampling in degrees

Definition at line 32 of file sampling.cpp.

{
    Vertex aux;
    aux.rot =   -PI / 5.;
    aux.tilt = PI / 2. - cte_w  ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot =    PI / 5.;
    aux.tilt = PI / 2. - cte_w  ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = 3. * PI / 5.;
    aux.tilt = PI / 2. - cte_w  ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = 5. * PI / 5.;
    aux.tilt = PI / 2. - cte_w  ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = -3. * PI / 5.;
    aux.tilt = PI / 2. - cte_w  ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot =    0. / 5.;
    aux.tilt = -cte_w + PI / 2. ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = 2. * PI / 5.;
    aux.tilt = -cte_w + PI / 2. ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = 4. * PI / 5.;
    aux.tilt = -cte_w + PI / 2. ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = -4. * PI / 5.;
    aux.tilt = -cte_w + PI / 2. ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);
    aux.rot = -2. * PI / 5.;
    aux.tilt = -cte_w + PI / 2. ;
    aux.psi =  0.;
    vertices_angles.push_back(aux);

    vertices_vectors.push_back(vectorR3(0., 0., 1.));
    vertices_vectors.push_back(vectorR3(0.723606900230461, -0.525731185781806, 0.447213343087301));
    vertices_vectors.push_back(vectorR3(0.723606900230461, 0.525731185781806, 0.447213343087301));
    vertices_vectors.push_back(vectorR3(-0.276393239417711, 0.850650928976665, 0.447213343087301));
    vertices_vectors.push_back(vectorR3(-0.8944273172062, 0., 0.447213343087301));
    vertices_vectors.push_back(vectorR3(-0.276393239417711, -0.850650928976665, 0.447213343087301));
    vertices_vectors.push_back(vectorR3(0.8944273172062, 0., -0.447213343087301));
    vertices_vectors.push_back(vectorR3(0.276393242471372, 0.850650927984471, -0.447213343087301));
    vertices_vectors.push_back(vectorR3(-0.723606898343194, 0.525731188379405, -0.447213343087301));
    vertices_vectors.push_back(vectorR3(-0.723606898343194, -0.525731188379405, -0.447213343087301));
    vertices_vectors.push_back(vectorR3(0.276393242471372, -0.850650927984471, -0.447213343087301));
    vertices_vectors.push_back(vectorR3(0., 0., -1.));

    sampling_noise=0.0;
    cos_neighborhood_radius=-1.01;

    numberSamplesAsymmetricUnit=-1;
    exp_data_fileNames.clear();

    verbose=1;
    //#define DEBUG1
#ifdef  DEBUG1

    for (int i = 0;
         i < vertices_vectors.size();
         i++)
        std::cout  <<  vertices_vectors[].transpose()  << std::endl;
#endif
#undef DEBUG1
}

Member Function Documentation

computeNeighbors()

void Sampling::computeNeighbors

(

bool

only_winner = false

)

for each point i in the asymmetric sampling unit cell compute the neighbors inside the asymmetric unit cell, save not only the neighbors but the angle psi

Definition at line 1715 of file sampling.cpp.

{

    double my_dotProduct;
    double winner_dotProduct;
    Matrix2D<double>  L(4, 4), R(4, 4);
    std::vector<size_t>  aux_neighbors;
    bool new_reference=true;
    my_neighbors.clear();
#ifdef MYPSI
    std::vector<double> aux_neighbors_psi;
    my_neighbors_psi.clear();
#endif

    // calculate some sizes only once
    size_t exp_data_projection_direction_by_L_R_size = exp_data_projection_direction_by_L_R.size();
    size_t no_redundant_sampling_points_vector_size = no_redundant_sampling_points_vector.size();

    if (verbose)
    {
        std::cout << "Find valid sampling points based on the neighborhood" <<std::endl;
        init_progress_bar(exp_data_projection_direction_by_L_R_size);
    }
    size_t ratio = exp_data_projection_direction_by_L_R_size / 60;
    ratio = XMIPP_MAX(ratio, 1);

    for(size_t j = 0; j < exp_data_projection_direction_by_L_R_size;)
    {
        if ((j%ratio) == 0 && verbose)
            progress_bar(j);

#ifdef MYPSI

        aux_neighbors_psi.clear();
#endif

        aux_neighbors.clear();
        if (cos_neighborhood_radius <= -1.0)
        {
            aux_neighbors=no_redundant_sampling_points_index;
            j+=R_repository.size();
        }
        else
        {
            size_t * aux_neighborsArray = nullptr;
            for (size_t k = 0; k < R_repository.size(); k++,j++)
            {
                winner_dotProduct = -1.;
                for (size_t i = 0; i < no_redundant_sampling_points_vector_size; ++i)
                {
                    my_dotProduct = dotProduct(no_redundant_sampling_points_vector[i],
                                               exp_data_projection_direction_by_L_R[j]);

                    if (my_dotProduct > cos_neighborhood_radius)
                    {
                        if(aux_neighbors.size()==0)
                        {
                            aux_neighbors.push_back(no_redundant_sampling_points_index[i]);
                            winner_dotProduct=my_dotProduct;

    #ifdef MYPSI

                            aux_neighbors_psi.push_back(exp_data_projection_direction_by_L_R_psi[j]);
    #endif

                        }
                        else
                        {
                            new_reference = true;
                            if(only_winner)
                            {
                                if(winner_dotProduct<my_dotProduct)
                                {
                                    if(winner_dotProduct!=-1)
                                        aux_neighbors.pop_back();
    #ifdef MYPSI

                                    if(winner_dotProduct!=-1)
                                        aux_neighbors_psi.pop_back();
    #endif

                                    winner_dotProduct=my_dotProduct;
                                }
                                else
                                {
                                    new_reference=false;
                                }
                            }
                            else
                            {
                                //precalculate size saves time here but
                                //not in the whole loop
                                aux_neighborsArray =  &aux_neighbors[0];
                                size_t _size = aux_neighbors.size();
                                //std::cerr << "DEBUG_JM: no_redundant_sampling_points_index[i]: " << no_redundant_sampling_points_index[i] << std::endl;
                                // for (size_t kkk = 0; kkk < aux_neighbors.size(); ++kkk)
                                //   std::cerr << aux_neighbors[kkk] << " ";
                                for( size_t l=0;l<  _size;l++)
                                {
                                    //if (aux_neighbors[l]==i)
                                    if (aux_neighborsArray[l]==no_redundant_sampling_points_index[i])
                                    {
                                        new_reference=false;
                                        break;
                                    }
                                }
                                //std::cerr << "DEBUG_JM: new_reference: " << new_reference << std::endl;
                            }
                            if (new_reference)
                            {
                                //std::cerr << formatString("DEBUG_JM: j %lu k %lu i %lu ", j, k, i) << std::endl;
                                aux_neighbors.push_back(no_redundant_sampling_points_index[i]);

                                //for (size_t kkk = 0; kkk < aux_neighbors.size(); ++kkk)
                                //  std::cerr << aux_neighbors[kkk] << " ";
                                // std::cerr << std::endl;

    #ifdef MYPSI

                                aux_neighbors_psi.push_back(exp_data_projection_direction_by_L_R_psi[j]);
    #endif

                            }
                        }
                        //same sampling point should appear only once
                        //note that psi recorded here may be different from psi
                        //recorded in _closest_sampling_points because
                        //may refer to a different sampling point
                        //in fact every point is degenerated
                    }
                }//for i;
            }//for k
        }
        my_neighbors.push_back(aux_neighbors);
#ifdef MYPSI

        my_neighbors_psi.push_back(aux_neighbors_psi);
#endif

    }//for j
    if (verbose)
        progress_bar(exp_data_projection_direction_by_L_R_size);

    //#define DEBUG
#ifdef DEBUG

    for (int i=0;i< my_neighbors.size();i++)
        for (int j=0;j< my_neighbors[i].size();j++)
            std::cerr << "image:" << i << " "<< my_neighbors[i][j]<<std::endl;
    exit(1);
#endif
#undef DEBUG

//#define CHIMERA
#ifdef CHIMERA

    std::ofstream filestr;
    filestr.open ("compute_neighbors.bild");
    filestr    << ".color white"
    << std::endl << ".sphere 0 0 0 .95" << std::endl
    ;
    int exp_image=0;
    filestr    <<  ".color yellow" << std::endl
    <<  ".sphere "
    << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](0) << " "
    << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](1) << " "
    << exp_data_projection_direction_by_L_R[exp_image*R_repository.size()](2) << " "
    <<  " .021"
    << std::endl;
    for(int i=(exp_image*R_repository.size());
        i< (exp_image+1)*R_repository.size();
        i++)
    {
        filestr    <<  ".color red" << std::endl
        <<  ".sphere "
//        << exp_data_projection_direction_by_L_R[i]
        << exp_data_projection_direction_by_L_R[i](0) << " "
        << exp_data_projection_direction_by_L_R[i](1) << " "
        << exp_data_projection_direction_by_L_R[i](2) << " "
        <<  " .017"      << std::endl;
    }
    double blue;
    for(int i=0;
        i< my_neighbors[exp_image].size();
        i++)
    {
#ifdef MYPSI
        blue = (my_neighbors_psi[exp_image][i]+180.)/360.;
#else

        blue = 1.;
#endif

        filestr    <<  ".color 0 0 " << blue  << std::endl
        <<  ".sphere "
//      << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]]
        << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](0) << " "
        << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](1) << " "
        << no_redundant_sampling_points_vector[my_neighbors[exp_image][i]](2) << " "

        <<  " .019"      << std::endl;
        //std::cerr << my_neighbors_psi[exp_image][i] << std::endl;
    }
    filestr.close();

#endif
    #undef CHIMERA

}

computeSamplingPoints()

void Sampling::computeSamplingPoints	(	bool	only_half_sphere = true,
		double	max_tilt = 180,
		double	min_tilt = 0
	)

Compute edge sampling points if you are looking only for directtions set only_half_sphere = true

vector to decimate the triangles

vector to decimate the triangles

Definition at line 155 of file sampling.cpp.

{
    std::vector <Matrix1D<double> > edge_vector_start;
    std::vector <Matrix1D<double> > edge_vector_end;
    // I need 10 auxiliary vector for edges
    Matrix1D<double> starting_point, ending_point;
    double max_z;
    double min_z;
    sampling_points_angles.clear();
    sampling_points_vector.clear();

    max_z=cos(PI * max_tilt / 180.);
    min_z=cos(PI * min_tilt / 180.);
    if ( (min_tilt > max_tilt) ||
         (min_tilt <0 ) ||
         (max_tilt < 0) ||
         (max_tilt > 180.)
         )
    {//OK
        std::cerr << "tilt angles cannot be negative or min_tilt > max_tilt" << std::endl;
        std::cerr << "min_tilt=" << min_tilt << "max_tilt=" << max_tilt << std::endl;
        exit(0);
    }

    if (min_z>max_z)
    {
        double aux=min_z;
        min_z=max_z;
        max_z=aux;
    }

    //01a
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[1];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[6];
    ending_point = vertices_vectors[1];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //01b
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[2];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[6];
    ending_point = vertices_vectors[2];
    fillEdge(starting_point, ending_point, edge_vector_end, true);
    //02a
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[2];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[7];
    ending_point = vertices_vectors[2];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //02b
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[3];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[7];
    ending_point = vertices_vectors[3];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //03a
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[3];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[8];
    ending_point = vertices_vectors[3];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //03b
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[4];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[8];
    ending_point = vertices_vectors[4];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //04a
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[4];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[9];
    ending_point = vertices_vectors[4];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //04b
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[5];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[9];
    ending_point = vertices_vectors[5];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //05a
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[5];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[10];
    ending_point = vertices_vectors[5];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //05b
    starting_point = vertices_vectors[0];
    ending_point = vertices_vectors[1];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[10];
    ending_point = vertices_vectors[1];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //06a
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[10];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[5];
    ending_point = vertices_vectors[10];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //06b
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[9];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[5];
    ending_point = vertices_vectors[9];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //07a
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[9];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[4];
    ending_point = vertices_vectors[9];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //07b
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[8];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[4];
    ending_point = vertices_vectors[8];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //08a
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[8];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[3];
    ending_point = vertices_vectors[8];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //08b
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[7];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[3];
    ending_point = vertices_vectors[7];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //09a
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[7];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[2];
    ending_point = vertices_vectors[7];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //09b
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[6];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[2];
    ending_point = vertices_vectors[6];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //10a
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[6];
    fillEdge(starting_point, ending_point, edge_vector_start, false);
    starting_point = vertices_vectors[1];
    ending_point = vertices_vectors[6];
    fillEdge(starting_point, ending_point, edge_vector_start, true);
    //10b
    starting_point = vertices_vectors[11];
    ending_point = vertices_vectors[10];
    fillEdge(starting_point, ending_point, edge_vector_end, false);
    starting_point = vertices_vectors[1];
    ending_point = vertices_vectors[10];
    fillEdge(starting_point, ending_point, edge_vector_end, true);

    //#define DEBUG2
#ifdef  DEBUG2

    std::ofstream filestr1;
    filestr1.open ("computeSamplingPoints_1.bild");

    for (int i = 0;
         i < edge_vector_start.size();
         i++)
    {
        filestr1  <<  ".sphere "
        << edge_vector_start[i](0) << " "
        << edge_vector_start[i](1) << " "
        << edge_vector_start[i](2) << " "
        << " .025" << std::endl;
        filestr1  <<  ".sphere "
        << edge_vector_end[i](0) << " "
        << edge_vector_end[i](1) << " "
        << edge_vector_end[i](2) << " "
        << " .025" << std::endl;
    }
    filestr1.close();

    //std::cout  <<  ending_point.transpose()    << " 1.1 1.5 " << std::endl;
#endif
#undef DEBUG2
    // add main corners that are not part of the basic octahedrons
    {
        int i=11;
        if (   (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
               || ZZ(vertices_vectors[i]) < min_z
               || ZZ(vertices_vectors[i]) > max_z
           )
            //if (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
            ;
        else
            sampling_points_vector.push_back(vertices_vectors[i]);
    }
    {
        int i=0;
        if (   (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
               || ZZ(vertices_vectors[i]) < min_z
               || ZZ(vertices_vectors[i]) > max_z
           )
            //if (only_half_sphere && ZZ(vertices_vectors[i]) < 0.0)
            ;
        else
            sampling_points_vector.push_back(vertices_vectors[i]);
    }

    // add edges
    for (size_t i = 0;
         i < edge_vector_start.size();
         i++)
    {
        if (i < number_of_samples * 10 - 15)
        {
            if (   (only_half_sphere && ZZ(edge_vector_start[i]) < 0.0)
                   || ZZ(edge_vector_start[i]) < min_z
                   || ZZ(edge_vector_start[i]) > max_z
               )
                //if (only_half_sphere && ZZ(edge_vector_start[i]) < 0.0)
                continue;
            else
                sampling_points_vector.push_back(edge_vector_start[i]);
        }
        else
        {
            if (   (only_half_sphere && ZZ(edge_vector_end[i]) < 0.0)
                   || ZZ(edge_vector_end[i]) < min_z
                   || ZZ(edge_vector_end[i]) > max_z
               )
                //if (only_half_sphere && ZZ(edge_vector_end[i]) < 0.0)
                continue;
            else
                sampling_points_vector.push_back(edge_vector_end[i]);
        }
    }
    //#define DEBUG3
#ifdef  DEBUG3
    std::ofstream filestr;
    filestr.open ("computeSamplingPoints_2.bild");
    for (int i = 0;
         i < sampling_points_vector.size();
         i++)
    {
        filestr  <<  ".color 1 0 0" << std::endl
        <<  ".sphere "
        << sampling_points_vector[i](0) << " "
        << sampling_points_vector[i](1) << " "
        << sampling_points_vector[i](2) << " "
        << " .025" << std::endl;
        //sampling_points_vector[i];//.add(0.0, sampling_noise, "gaussian");
        sampling_points_vector[i].selfNormalize();
        filestr  <<  ".color 0 0 1" << std::endl
        <<  ".sphere "
        << sampling_points_vector[i](0) << " "
        << sampling_points_vector[i](1) << " "
        << sampling_points_vector[i](2) << " "
        <<" .027" << std::endl;
    }
    filestr.close();

#endif
#undef DEBUG3

    // add in between points
    int j = 0;
    bool j_flag = false;
    for (size_t i = 0;
         i < edge_vector_start.size();
         i++)
    {
        if ((j % (number_of_samples - 1)) == 0 && j != 0)
        {
            j = 0;
            j_flag = true;
        }
        if ((j % (number_of_samples - 2)) == 0 && j != 0  && j_flag == true)
        {
            j = 0;
            j_flag = false;
        }
        fillDistance(edge_vector_start[i],
                     edge_vector_end[i],
                     sampling_points_vector,
                     (j + 1) % number_of_samples,
                     only_half_sphere,min_z,max_z);
        j++;
    }
    //#define DEBUG3
#ifdef  DEBUG3
    std::ofstream filestr3;
    filestr3.open ("computeSamplingPoints_3.bild");
    for (int i = 0;
         i < sampling_points_vector.size();
         i++)
    {
        filestr3  <<  ".color 1 0 0" << std::endl
        <<  ".sphere "
        << sampling_points_vector[i](0) << " "
        << sampling_points_vector[i](1) << " "
        << sampling_points_vector[i](2) << " "
        << " .025" << std::endl;
    }
    filestr3.close();
#endif
#undef DEBUG3

    //noisify angles
    if(sampling_noise!=0.0)
    {
        for (size_t n = 0;
             n < sampling_points_vector.size();
             n++)
        {
            FOR_ALL_ELEMENTS_IN_MATRIX1D(sampling_points_vector[n])
            (sampling_points_vector[n])(i) += rnd_gaus(0., sampling_noise);
            sampling_points_vector[n].selfNormalize();
        }
    }

    //#define DEBUG3
#ifdef  DEBUG3
    for (int i = 0;
         i < sampling_points_vector.size();
         i++)
    {
        std::cout  <<  ".color 0 1 0" << std::endl
        <<  ".sphere " << sampling_points_vector[i].transpose()  <<
        " .03" << std::endl;
    }
#endif
#undef DEBUG3

    // store sampling points as angles
    Matrix1D<double> aux(3), aux1(3);
    ZZ(aux) = 0.;
    for (size_t i = 0;
         i < sampling_points_vector.size();
         i++)
    {
        XX(aux) = atan2(YY(sampling_points_vector[i]),
                        XX(sampling_points_vector[i]));
        YY(aux) = acos(ZZ(sampling_points_vector[i]));
        if (YY(aux) < 0.)
            YY(aux) += PI;
        sampling_points_angles.push_back(aux*180. / PI);
    }
#ifdef NEVERDEFINED
    //sort points
    int k;
    int bound = sampling_points_angles.size() - 1;
    int t;
    int last_swap;
    Matrix1D<double> aux_angle(3);
    Matrix1D<double> aux_vector(3);

    while (bound)
    {
        last_swap = 0;
        for (k = 0; k < bound; k++)
        {
            aux = sampling_points_angles[k]; /* aux is a maximum of A[0]..A[k] */
            aux1 = sampling_points_vector[k]; /* aux is a maximum of A[0]..A[k] */
            if (sortFunc(aux, sampling_points_angles[k+1]))
            {
                sampling_points_angles[k] = sampling_points_angles[k+1];
                //                sampling_points_angles[k] = sampling_points_angles[k+1];
                sampling_points_angles[k+1] = aux; /*swap*/
                sampling_points_vector[k] = sampling_points_vector[k+1];
                sampling_points_vector[k] = sampling_points_vector[k+1];
                sampling_points_vector[k+1] = aux1; /*swap*/
                last_swap = k; /* mark the last swap position */
            }//if
        }//for
        bound = last_swap; /*elements after bound already sorted */
    }//while


#endif
    //#define DEBUG3
#ifdef  DEBUG3
    std::ofstream filestr4;
    filestr4.open ("computeSamplingPoints_4.bild");

    filestr4    << ".color yellow"
    << std::endl
    ;
    for (int i = 0;
         i < sampling_points_vector.size();
         i++)
    {
        filestr4 <<  ".sphere "
        << sampling_points_vector[i](0) << " "
        << sampling_points_vector[i](1) << " "
        << sampling_points_vector[i](2) << " "
        << " .025" << std::endl;
    }
    filestr4.close();
#endif
#undef DEBUG3

    numberSamplesAsymmetricUnit=sampling_points_vector.size();
}

createAsymUnitFile()

void Sampling::createAsymUnitFile

(

const FileName &

docfilename

)

save asymmetric unit sampling in a doc file

Definition at line 1440 of file sampling.cpp.

{
    MetaDataVec DF;
    FileName tmp_filename;
    //#define CHIMERA
#ifdef CHIMERA

    std::ofstream filestr;
    filestr.open ("create_asym_unit_file.bild");
    filestr    << ".color white"
    << std::endl
    << ".sphere 0 0 0 .95"
    << std::endl
    ;
    filestr    << ".color green"
    << std::endl
    ;
#endif

    MDRowVec row;
    for (size_t i = 0; i < no_redundant_sampling_points_vector.size(); i++)
    {
#ifdef CHIMERA
        filestr  << ".sphere "
        << no_redundant_sampling_points_vector[i]
        << " 0.018"
        << std::endl
        ;
#endif

        row.setValue(MDL_REF, (int)i);
        row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
        row.setValue(MDL_ANGLE_ROT,XX(no_redundant_sampling_points_angles[i]));
        row.setValue(MDL_ANGLE_TILT,YY(no_redundant_sampling_points_angles[i]));
        row.setValue(MDL_ANGLE_PSI,ZZ(no_redundant_sampling_points_angles[i]));
        row.setValue(MDL_X,XX(no_redundant_sampling_points_vector[i]));
        row.setValue(MDL_Y,YY(no_redundant_sampling_points_vector[i]));
        row.setValue(MDL_Z,ZZ(no_redundant_sampling_points_vector[i]));
        DF.addRow(row);
    }
#ifdef CHIMERA
    filestr.close();
#endif
    #undef CHIMERA

    tmp_filename = docfilename + "_angles.doc";
    DF.setComment("REF refers to the projection directions BEFORE delete those not in a neighborhood, ");
    DF.write(tmp_filename);
}

createSymFile()

void Sampling::createSymFile	(	const FileName &	simFp,
		int	symmetry,
		int	sym_order
	)

create symmetry file from introduced symmetry see SymList class

Definition at line 1319 of file sampling.cpp.

{
    symmetry_file = simFp + ".sym";
    std::ofstream SymFile;
    SymFile.open(symmetry_file.c_str(), std::ios::out);
    if (symmetry == pg_CN)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
    }
    else if (symmetry == pg_CI)
    {
        SymFile << "inversion ";
    }
    else if (symmetry == pg_CS)
    {
        SymFile << "mirror_plane 0 0 1";
    }
    else if (symmetry == pg_CNV)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "mirror_plane 0 1 0";
    }
    else if (symmetry == pg_CNH)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "mirror_plane 0 0 -1";
    }
    else if (symmetry == pg_SN)
    {
        int order = sym_order / 2;
        SymFile << "rot_axis " << order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "inversion ";
    }
    else if (symmetry == pg_DN)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 0 1 0";
    }
    else if (symmetry == pg_DNV)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 0 1 0";
        SymFile << std::endl;
        SymFile << "mirror_plane 0 1 0";
    }
    else if (symmetry == pg_DNH)
    {
        SymFile << "rot_axis " << sym_order << " 0 0 1";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 1 0 0";
        SymFile << std::endl;
        SymFile << "mirror_plane 0 0 1";
    }
    else if (symmetry == pg_T)
    {
        SymFile << "rot_axis " << "3" << "  0. 0. 1.";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 0. 0.816496 0.577350";
    }
    else if (symmetry == pg_TD)
    {
        SymFile << "rot_axis " << "3" << "  0. 0. 1.";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 0. 0.816496 0.577350";
        SymFile << std::endl;
        SymFile << "mirror_plane 1.4142136 2.4494897 0.0000000";
    }
    else if (symmetry == pg_TH)
    {
        SymFile << "rot_axis " << "3" << "  0. 0. 1.";
        SymFile << std::endl;
        SymFile << "rot_axis " << "2" << " 0. -0.816496 -0.577350";
        SymFile << std::endl;
        SymFile << "inversion";
    }
    else if (symmetry == pg_O)
    {
        SymFile << "rot_axis " << "3" << "  .5773502  .5773502 .5773502";
        SymFile << std::endl;
        SymFile << "rot_axis " << "4" << " 0 0 1";
    }
    else if (symmetry == pg_OH)
    {
        SymFile << "rot_axis " << "3" << "  .5773502  .5773502 .5773502";
        SymFile << std::endl;
        SymFile << "rot_axis " << "4" << " 0 0 1";
        SymFile << std::endl;
        SymFile << "mirror_plane 0 1 1";
    }
    else if (symmetry == pg_I)
    {
        SymFile << "rot_axis 2  0             0          1";
        SymFile << std::endl;
        SymFile << "rot_axis 5 -1.618033989  -1           0";
        SymFile << std::endl;
        SymFile << "rot_axis 3 -0.53934467   -1.4120227   0";
    }
    else if (symmetry == pg_IH)
    {
        SymFile << "rot_axis 2  0             0          1";
        SymFile << std::endl;
        SymFile << "rot_axis 5 -1.618033989  -1           0";
        SymFile << std::endl;
        SymFile << "rot_axis 3 -0.53934467   -1.4120227   0";
        SymFile << std::endl;
        SymFile << "mirror_plane 1 0 0";
    }
    else
    {
        std::cerr << "ERROR:: Symmetry " << symmetry  << "is not known" << std::endl;
        exit(0);
    }
    SymFile.close();
    SL.readSymmetryFile(symmetry_file);

}

fillDistance()

void Sampling::fillDistance	(	const Matrix1D< double > &	starting_point,
		const Matrix1D< double > &	ending_point,
		std::vector< Matrix1D< double > > &	edge_vector,
		int	number,
		bool	only_half_spheree,
		double	min_z = -10.,
		double	max_z = +10.
	)

fill distance

Definition at line 631 of file sampling.cpp.

{
    Matrix1D<double> v_aux(3);

    double alpha;
    double beta;
    double gamma;
    // skip first corener, already computed;
    double upsilon = acos(dotProduct(starting_point, ending_point));

    for (int i1 = 1; i1 < my_number_of_samples; i1++)
    {
        gamma  = (double)i1 / my_number_of_samples;
        alpha  = sin((1. - gamma) * upsilon) / (sin(upsilon));
        beta   = sin(gamma * upsilon) / sin(upsilon);
        v_aux = alpha * starting_point + beta * ending_point;
        v_aux.selfNormalize();
        if (   (only_half_sphere && ZZ(v_aux) < 0.0)
               || ZZ(v_aux) < min_z
               || ZZ(v_aux) > max_z
           )
            //if (only_half_sphere && ZZ(v_aux) < 0.0)
            continue;
        else
            sampling_points_vector.push_back(v_aux);
    }
    /*
        //Remove points not in tilt range, check only the edges
        for (int i = 0;
             i < auxCounter;
             i++)
        {
            if (  ZZ(sampling_points_vector[i]) < min_z
               || ZZ(sampling_points_vector[i]) > max_z
               )
            sampling_points_vector[i].remove();
        }
    */
}

fillEdge()

void Sampling::fillEdge	(	const Matrix1D< double > &	starting_point,
		const Matrix1D< double > &	ending_point,
		std::vector< Matrix1D< double > > &	edge_vector,
		bool	FLAG_END
	)

fill edge

Definition at line 606 of file sampling.cpp.

{
    Matrix1D<double> v_aux(3);

    double alpha;
    double beta;
    double gamma;
    // skip first corener, already computed;
    double upsilon = acos(dotProduct(starting_point, ending_point));
    for (size_t i1 = 1; i1 < number_of_samples; i1++)
    {
        gamma  = (double)i1 / (number_of_samples - 1);
        alpha  = sin((1. - gamma) * upsilon) / (sin(upsilon));
        beta   = sin(gamma * upsilon) / sin(upsilon);
        v_aux = alpha * starting_point + beta * ending_point;
        v_aux.selfNormalize();
        if (beta > 0.9999 && END_FLAG)
            continue;
        edge_vector.push_back(v_aux);
    }
}

fillExpDataProjectionDirectionByLR() [1/2]

void Sampling::fillExpDataProjectionDirectionByLR

(

const MetaData &

DFi

)

Precalculate exp_data by symmetry matrices (speeds up calculations)

Definition at line 2256 of file sampling.cpp.

{
    std::vector <Matrix1D<double> > exp_data_projection_direction;
    Matrix1D<double>  direction(3);
    DFi.firstRowId();
    //#define CHIMERA
#ifdef CHIMERA

    std::ofstream filestr;
    filestr.open ("exp_data_projection_direction_by_L_R.bild");
    filestr    << ".color white"
    << std::endl
    << ".sphere 0 0 0 .95"
    << std::endl
    ;
    filestr    << ".color green"
    << std::endl
    ;
#endif

    double img_tilt,img_rot,img_psi;
    FileName imgName;
    exp_data_fileNames.clear();
    for (size_t objId : DFi.ids())
    {
        DFi.getValue(MDL_ANGLE_ROT,img_rot, objId);
        DFi.getValue(MDL_ANGLE_TILT,img_tilt, objId);
        DFi.getValue(MDL_ANGLE_PSI,img_psi, objId);
        Euler_direction(img_rot, img_tilt, img_psi, direction);
        exp_data_projection_direction.push_back(direction);
        DFi.getValue(MDL_IMAGE,imgName, objId);
        exp_data_fileNames.push_back(imgName);
    }

    exp_data_projection_direction_by_L_R.clear();
#ifdef MYPSI

    exp_data_projection_direction_by_L_R_psi.clear();
#endif

    for (size_t i = 0; i < exp_data_projection_direction.size(); i++)
        for (size_t j = 0; j < R_repository.size(); j++)
        {
            direction =  L_repository[j] *
                         (exp_data_projection_direction[i].transpose() *
                          R_repository[j]).transpose();
            exp_data_projection_direction_by_L_R.push_back(direction);
#ifdef MYPSI

            Euler_apply_transf(L_repository[j],
                               R_repository[j], img_rot,
                               img_tilt,
                               img_psi,
                               rotp,
                               tiltp,
                               psip);
            exp_data_projection_direction_by_L_R_psi.push_back(psip);
#endif
            //#define CHIMERA
#ifdef CHIMERA

            filestr << ".sphere " << direction.transpose()
            << " 0.02" << std::endl
            ;
#endif

        }
    //#define CHIMERA
#ifdef CHIMERA
    filestr.close();
#endif
    #undef CHIMERA
}

fillExpDataProjectionDirectionByLR() [2/2]

void Sampling::fillExpDataProjectionDirectionByLR

(

const FileName &

FnexperimentalImages

)

Precalculate exp_data by symmetry matrices (speeds up calculations)

Definition at line 2247 of file sampling.cpp.

{
    //read input files
    MetaDataVec DFi;
    DFi.read(FnexperimentalImages);//experimental points
    fillExpDataProjectionDirectionByLR(DFi);
}

fillLRRepository()

void Sampling::fillLRRepository

(

void

)

fill R and L Repository (vector with symmetry matrices)

Definition at line 2216 of file sampling.cpp.

{
    Matrix2D<double>  L(4, 4), R(4, 4);
    Matrix2D<double>  Identity(3,3);
    Identity.initIdentity();
    //NEXT 2 ROB
    R_repository.clear();
    L_repository.clear();
    R_repository.push_back(Identity);
    L_repository.push_back(Identity);
    for (int isym = 0; isym < SL.symsNo(); isym++)
    {
        SL.getMatrices(isym, L, R);
        R.resize(3, 3);
        L.resize(3, 3);
        R_repository.push_back(R);
        L_repository.push_back(L);
    }
//#define DEBUG3
#ifdef  DEBUG3
    std::cout << "==============================\n" ;
    for (int isym = 0; isym < R_repository.size(); isym++)
    {
        std::cout << R_repository[isym];
        //std::cout << L_repository[isym];
    }
    std::cout << "-------------------------------\n";
#endif
#undef DEBUG3
}

findClosestExperimentalPoint()

void Sampling::findClosestExperimentalPoint

(

)

for each sampling point find the experimental images closer to that point than to any other

Definition at line 2100 of file sampling.cpp.

{
    double my_dotProduct,my_dotProduct_aux;
    Matrix1D<double>  row(3),direction(3);
    int winner_sampling=-1;
    int winner_exp=-1;
    //#define CHIMERA
#ifdef CHIMERA

    std::vector<std::vector<int> >  aux_vec;
    aux_vec.resize(no_redundant_sampling_points_vector.size());
#endif

    std::vector<std::vector<size_t> >  aux_my_exp_img_per_sampling_point;

    //resize vector
    aux_my_exp_img_per_sampling_point.resize(
        no_redundant_sampling_points_vector.size());

    for(size_t i=0,l=0;i< exp_data_projection_direction_by_L_R.size();l++)
    {
        my_dotProduct=-2;
        for (size_t k = 0; k < R_repository.size(); k++,i++)
        {
            for(size_t j=0;j< no_redundant_sampling_points_vector.size();j++)
            {
                my_dotProduct_aux =
                    dotProduct(exp_data_projection_direction_by_L_R[i],
                               no_redundant_sampling_points_vector[j]);

                if ( my_dotProduct_aux > my_dotProduct)
                {
                    my_dotProduct = my_dotProduct_aux;
                    winner_sampling = j;
#ifdef CHIMERA

                    winner_exp_L_R  = i;
#endif

                    winner_exp = l;
                }
            }//for j
        }//for k
        aux_my_exp_img_per_sampling_point[winner_sampling].push_back(winner_exp);
#ifdef CHIMERA

        aux_vec[winner_sampling].push_back(winner_exp_L_R);
#endif

    }//for i aux_my_exp_img_per_sampling_point
    for(size_t i=0;i< aux_my_exp_img_per_sampling_point.size();i++)
        if(aux_my_exp_img_per_sampling_point[i].size()!=0)
            my_exp_img_per_sampling_point.push_back(aux_my_exp_img_per_sampling_point[i]);
#ifdef CHIMERA

    std::ofstream filestr;
    filestr.open ("find_closest_experimental_point.bild");
    filestr    << ".color white"
    << std::endl
    << ".sphere 0 0 0 .95"
    << std::endl
    ;
    filestr    << ".color red"
    << std::endl
    ;
    for (int i = 0;
         i < no_redundant_sampling_points_vector.size();
         i++)
    {
        filestr    <<  ".sphere " << no_redundant_sampling_points_vector[i].transpose()  <<
        " .018" << std::endl;
    }
    int my_sample_point=5;
    filestr    << ".color green"
    << std::endl
    ;
    int ii;
    for (int i = 0;
         i < my_exp_img_per_sampling_point[my_sample_point].size();
         i++)
    {
        ii=aux_vec[my_sample_point][i];
        filestr    << ".sphere "
        << exp_data_projection_direction_by_L_R[ii].transpose()
        << " .017" << std::endl;
    }
    filestr.close();

#endif
    #undef CHIMERA
    //#define DEBUG4
#ifdef DEBUG4

    std::ofstream filestr;
    filestr.open ("find_closest_experimental_point.txt");

    for (int i = 0;
         i < my_exp_img_per_sampling_point.size();
         i++)
    {   //for each sampling point write its experimental images
        filestr << i << std::endl;
        for (int j = 0;
             j < my_exp_img_per_sampling_point[i].size();
             j++)
        {
            filestr    << my_exp_img_per_sampling_point[i][j]
            << " " ;
        }
        filestr << std::endl;
    }
    filestr.close();

#endif
    #undef DEBUG4
}

findClosestSamplingPoint() [1/2]

void Sampling::findClosestSamplingPoint	(	const MetaData &	DFi,
		const FileName &	output_file_root
	)

Find the closest sampling point for a docfile of experimental projections

Definition at line 1991 of file sampling.cpp.

{
    double my_dotProduct,my_dotProduct_aux;
    Matrix1D<double>  row(3),direction(3);
    Matrix1D<double> docline;
    docline.initZeros(7);//three original angles, one winnir, new angles
    Matrix2D<double>  L(4, 4), R(4, 4);
    int winner_sampling=-1;
#if defined(CHIMERA) || defined(MYPSI)

    int winner_exp_L_R=-1;
#endif

    MetaDataVec DFo;
    size_t id;

    DFo.setComment("Original rot, tilt, psi, Xoff, Yoff are stored as comments");

    //#define DEBUG3
#ifdef  DEBUG3

    std::ofstream filestr;
    filestr.open ("find_closest_sampling_point.bild");
    int exp_image=1;
#endif

    auto idIter(DFi.ids().begin());
    for (size_t i = 0; i <  exp_data_projection_direction_by_L_R.size(); )
    {
        my_dotProduct=-2;
        for (size_t k = 0; k < R_repository.size(); k++,i++)
        {
#ifdef  DEBUG3
            //experimental points plus symmetry
            if( i>(exp_image*R_repository.size()-1) && i< ((exp_image+1)*R_repository.size()))
            {
                filestr    <<  ".color red" << std::endl
                <<  ".sphere "   << exp_data_projection_direction_by_L_R[i]
                <<  " .019"      << std::endl;
            }
#endif
            for(size_t j=0;j< no_redundant_sampling_points_vector.size();j++)
            {
                my_dotProduct_aux =
                    dotProduct(exp_data_projection_direction_by_L_R[i],
                               no_redundant_sampling_points_vector[j]);

                if ( my_dotProduct_aux > my_dotProduct)
                {
                    my_dotProduct = my_dotProduct_aux;
                    winner_sampling = j;
#if defined(CHIMERA) || defined(MYPSI)

                    winner_exp_L_R  = i;
#endif

                }
            }//for j
        }//for k
#ifdef  DEBUG3
        if( i==  ((exp_image+1)*R_repository.size()) )
        {
            filestr    <<  ".color yellow" << std::endl
            <<  ".sphere "   << no_redundant_sampling_points_vector[winner_sampling]
            <<  " .020"      << std::endl;
        }
#endif
        //add winner to the DOC fILE
        std::string fnImg, comment;
        double aux;
        size_t objId = *idIter;
        DFi.getValue(MDL_IMAGE, fnImg, objId);
        DFi.getValue(MDL_ANGLE_ROT,aux, objId);
        comment+=floatToString(aux)+" ";
        DFi.getValue(MDL_ANGLE_TILT,aux, objId);
        comment+=floatToString(aux)+" ";
        DFi.getValue(MDL_ANGLE_PSI,aux, objId);
        comment+=floatToString(aux)+" ";
        DFi.getValue(MDL_SHIFT_X,aux, objId);
        comment+=floatToString(aux)+" ";
        DFi.getValue(MDL_SHIFT_Y,aux, objId);
        comment+=floatToString(aux);
        id = DFo.addObject();
        DFo.setValue(MDL_STAR_COMMENT,comment, id);
        DFo.setValue(MDL_IMAGE,fnImg, id);
        DFo.setValue(MDL_REF, winner_sampling, id);
#ifdef MYPSI

        DFo.set(6, exp_data_projection_direction_by_L_R_psi[winner_exp_L_R]);
#endif

        DFo.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[winner_sampling], id);
        DFo.setValue(MDL_ANGLE_ROT,XX(no_redundant_sampling_points_angles[winner_sampling]), id);
        DFo.setValue(MDL_ANGLE_TILT,YY(no_redundant_sampling_points_angles[winner_sampling]), id);
        DFo.setValue(MDL_ANGLE_PSI,ZZ(no_redundant_sampling_points_angles[winner_sampling]), id);

        ++idIter;
    }//for i

    if (output_file_root.size() > 0)
        DFo.write(output_file_root+ "_closest_sampling_points.doc");
#ifdef  DEBUG3

    filestr.close();
#endif
#undef DEBUG3
}

findClosestSamplingPoint() [2/2]

void Sampling::findClosestSamplingPoint	(	const FileName &	FnexperimentalImages,
		const FileName &	output_file_root
	)

Find the closest sampling point for a docfile of experimental projections

Definition at line 1982 of file sampling.cpp.

{
    //read input files
    MetaDataVec DFi;
    DFi.read(FnexperimentalImages);//experimental points
    findClosestSamplingPoint(DFi, output_file_root);

}

operator==()

bool Sampling::operator==

(

const Sampling &

op

)

const

'is equal to' (equality).

Definition at line 107 of file sampling.cpp.

{

    return (XMIPP_EQUAL_REAL(sampling_rate_rad, op.sampling_rate_rad) &&
            XMIPP_EQUAL_REAL(cos_neighborhood_radius, op.cos_neighborhood_radius) &&
            no_redundant_sampling_points_angles == op.no_redundant_sampling_points_angles &&
            no_redundant_sampling_points_vector == op.no_redundant_sampling_points_vector &&
            no_redundant_sampling_points_index == op.no_redundant_sampling_points_index &&
            my_neighbors == op.my_neighbors //&&
            //exp_data_fileNames == op.exp_data_fileNames
           );

}

readSamplingFile()

void Sampling::readSamplingFile	(	const FileName &	infilename,
		bool	read_vectors = true,
		bool	read_sampling_sphere = false
	)

Read neighbors i

This is the object ID in the metadata, usually starts at 1

This is the index of the object, starts at 0

Definition at line 1592 of file sampling.cpp.

{
    //Read extra info
    MetaDataVec md(FN_SAMPLING_EXTRA(fn_base));
    size_t id = md.firstRowId();
    md.getValue(MDL_SAMPLINGRATE, sampling_rate_rad, id);
    md.getValue(MDL_NEIGHBORHOOD_RADIUS, cos_neighborhood_radius, id);
    md.getValue(MDL_POINTSASYMETRICUNIT,numberSamplesAsymmetricUnit,id);

    //Read neighbors
    md.read(FN_SAMPLING_NEI(fn_base));
    my_neighbors.resize(md.size());
    bool readFileName = md.containsLabel(MDL_IMAGE);
    if (readFileName)
    {
        exp_data_fileNames.clear();
        exp_data_fileNames.resize(md.size());
    }

    {
        size_t i = 0, ii = 0;
        for (size_t objId : md.ids())
        {
            if (readFileName)
                md.getValue(MDL_IMAGE,exp_data_fileNames[ii++], objId);
            md.getValue(MDL_NEIGHBORS, my_neighbors[i], objId);
            i++;
        }
    }

    //Read projection directions
    md.read(FN_SAMPLING_PROJ(fn_base));
    size_t size = md.size();
    no_redundant_sampling_points_index.resize(size);
    no_redundant_sampling_points_angles.resize(size);
    if (read_vectors)
        no_redundant_sampling_points_vector.resize(size);

    {
        size_t i = 0;
        for (size_t objId : md.ids())
        {
            //size_t objId;
            //size_t objIndex;

            Matrix1D<double> &angles = no_redundant_sampling_points_angles[i];
            angles.resizeNoCopy(3);
            md.getValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i], objId);
            md.getValue(MDL_ANGLE_ROT, XX(angles), objId);
            md.getValue(MDL_ANGLE_TILT, YY(angles), objId);
            md.getValue(MDL_ANGLE_PSI, ZZ(angles), objId);

            if (read_vectors)
            {
                Matrix1D<double> &vectors = no_redundant_sampling_points_vector[i];
                vectors.resizeNoCopy(3);
                md.getValue(MDL_X, XX(vectors), objId);
                md.getValue(MDL_Y, YY(vectors), objId);
                md.getValue(MDL_Z, ZZ(vectors), objId);
            }

            i++;
        }
    }

//#define DEBUG5
#ifdef  DEBUG5
//This bild file does not make sense since x,y,z are angles
    std::ofstream filestr5;
    filestr5.open ("/tmp/loadsamplingfile_removeRedundantPoints.bild");
    filestr5 << ".color 0 0 1" << std::endl;
    for (int i = 0;
         i < no_redundant_sampling_points_index.size();
         i++)
    {
        filestr5  <<  ".color 1 0 0" << std::endl
        <<  ".sphere "
        << no_redundant_sampling_points_vector[i](0) << " "
        << no_redundant_sampling_points_vector[i](1) << " "
        << no_redundant_sampling_points_vector[i](2) << " "
        << " .03" << std::endl;
    }
    filestr5.close();
#endif
#undef DEBUG5

    if (read_sampling_sphere)
    {
        //Read projection directions
        md.read(FN_SAMPLING_SPHERE(fn_base));
        size_t size = md.size();
        sampling_points_angles.resize(size);
        if (read_vectors)
            sampling_points_vector.resize(size);

        size_t i = 0;
        for (size_t objId : md.ids())
        {
            Matrix1D<double> &angles = sampling_points_angles[i];
            angles.resizeNoCopy(3);
            md.getValue(MDL_ANGLE_ROT, XX(angles), objId);
            md.getValue(MDL_ANGLE_TILT, YY(angles), objId);
            md.getValue(MDL_ANGLE_PSI, ZZ(angles), objId);

            if (read_vectors)
            {
                Matrix1D<double> &vectors = sampling_points_vector[i];
                vectors.resizeNoCopy(3);
                md.getValue(MDL_X, XX(vectors), objId);
                md.getValue(MDL_Y, YY(vectors), objId);
                md.getValue(MDL_Z, ZZ(vectors), objId);
            }

            i++;
        }

    }
}

removePointsFarAwayFromExperimentalData()

void Sampling::removePointsFarAwayFromExperimentalData

(

)

remove all those points that are further away from experimental data than neighborhood_radius_rad

Definition at line 1928 of file sampling.cpp.

{
    double my_dotProduct;
    Matrix1D<double>  row(3),direction(3);
    Matrix2D<double>  L(4, 4), R(4, 4);

    size_t my_end = no_redundant_sampling_points_vector.size() - 1;

    for (size_t i = 0; i <= my_end; i++)
    {
        bool my_delete=true;
        for (size_t j=0; my_delete && j< exp_data_projection_direction_by_L_R.size();j++)
        {
            my_dotProduct = dotProduct(no_redundant_sampling_points_vector[i],
                                       exp_data_projection_direction_by_L_R[j]);
            if (my_dotProduct > cos_neighborhood_radius)
                my_delete=false;
        }//for j
        if(my_delete)
        {
            REMOVE_LAST(no_redundant_sampling_points_vector);
            REMOVE_LAST(no_redundant_sampling_points_angles);
            REMOVE_LAST(no_redundant_sampling_points_index);

            --my_end;
            --i;//since a point has been swaped we should repeat the same index
        }// if(my_delete)
    }//for i end
    //#define CHIMERA
#ifdef CHIMERA
    std::ofstream filestr;
    filestr.open ("remove_points_far_away_from_experimental_data.bild");
    filestr    << ".color white"
    << std::endl
    << ".sphere 0 0 0 .95"
    << std::endl
    ;
    filestr    << ".color green"
    << std::endl
    ;
    //green neighbours
    for (int i = 0;
         i < no_redundant_sampling_points_vector.size();
         i++)
    {
        filestr    <<  ".color green" << std::endl
        <<  ".sphere " << no_redundant_sampling_points_vector[i].transpose()  <<
        " .018" << std::endl;
    }
    filestr.close();

#endif
    #undef CHIMERA
}

removeRedundantPoints()

void Sampling::removeRedundantPoints	(	const int	symmetry,
		int	sym_order
	)

Definition at line 691 of file sampling.cpp.

{
    Matrix2D<double>  L(4, 4), R(4, 4);
    Matrix2D<double>  aux(3, 3);
    Matrix1D<double>  row1(3), row2(3);

    //int j_end=0;
    Matrix1D<double>  row(3);

    CLEAR_VECTORS();

    if (symmetry == pg_CN)
    {//OK
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >= (-180. / sym_order) &&
                XX(sampling_points_angles[i]) <= (180. / sym_order))
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry == pg_CI  ||
             symmetry == pg_CS )
    {//OK
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (YY(sampling_points_angles[i]) <= 90)
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_CNV )
    {//OK
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >=    0. / sym_order &&
                XX(sampling_points_angles[i]) <=  180. / sym_order)
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_CNH )
    {//OK
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >= -180. / sym_order &&
                XX(sampling_points_angles[i]) <=  180. / sym_order &&
                YY(sampling_points_angles[i]) <=    90.
               )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_SN )
    {//OK
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >= -180.*2. / sym_order &&
                XX(sampling_points_angles[i]) <=  180.*2. / sym_order &&
                YY(sampling_points_angles[i]) <=    90.
               )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_DN )
    {
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >= -180. / sym_order  + 90. &&
                XX(sampling_points_angles[i]) <=  180. / sym_order  + 90. &&
                YY(sampling_points_angles[i]) <=    90.
               )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_DNV )
    {
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >=    0. + 90. &&
                XX(sampling_points_angles[i]) <=  180. / sym_order +90. &&
                YY(sampling_points_angles[i]) <=    90.
               )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_DNH )
    {
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >= 0. &&
                XX(sampling_points_angles[i]) <=  180. / sym_order &&
                YY(sampling_points_angles[i]) <=   90.
               )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_T )
    {//OK
        Matrix1D<double>  _3_fold_axis_1_by_3_fold_axis_2(3);
        _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(-0.942809, 0., 0.);
        _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
        Matrix1D<double>  _3_fold_axis_2_by_3_fold_axis_3(3);
        _3_fold_axis_2_by_3_fold_axis_3 = vectorR3(0.471405, 0.272165, 0.7698);
        _3_fold_axis_2_by_3_fold_axis_3.selfNormalize();
        Matrix1D<double>  _3_fold_axis_3_by_3_fold_axis_1(3);
        _3_fold_axis_3_by_3_fold_axis_1 = vectorR3(0.471404, 0.816497, 0.);
        _3_fold_axis_3_by_3_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if ((XX(sampling_points_angles[i]) >=     90. &&
                 XX(sampling_points_angles[i]) <=   150. )||
                 XX(sampling_points_angles[i]) ==     0
                )
                if (
                    dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
                    dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_3_fold_axis_3) >= 0 &&
                    dotProduct(sampling_points_vector[i], _3_fold_axis_3_by_3_fold_axis_1) >= 0
                )
                {
                    no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                    no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
                }
        }// for i
    }
    else if (symmetry  == pg_TD )
    {//OK
        Matrix1D<double>  _2_fold_axis_1_by_3_fold_axis_2(3);
        _2_fold_axis_1_by_3_fold_axis_2 = vectorR3(-0.942809, 0., 0.);
        _2_fold_axis_1_by_3_fold_axis_2.selfNormalize();
        Matrix1D<double>  _3_fold_axis_2_by_3_fold_axis_5(3);
        _3_fold_axis_2_by_3_fold_axis_5 = vectorR3(0.471405, 0.272165, 0.7698);
        _3_fold_axis_2_by_3_fold_axis_5.selfNormalize();
        Matrix1D<double>  _3_fold_axis_5_by_2_fold_axis_1(3);
        _3_fold_axis_5_by_2_fold_axis_1 = vectorR3(0., 0.471405, -0.666667);
        _3_fold_axis_5_by_2_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            //           if ( XX(sampling_points_angles[i])>=     120. &&
            //                 XX(sampling_points_angles[i])<=   150. ||
            //                 XX(sampling_points_angles[i])==     0
            //              )
            if (
                dotProduct(sampling_points_vector[i], _2_fold_axis_1_by_3_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_3_fold_axis_5) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_5_by_2_fold_axis_1) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_TH )
    {//OK
        Matrix1D<double>  _3_fold_axis_1_by_2_fold_axis_1(3);
        _3_fold_axis_1_by_2_fold_axis_1 = vectorR3(-0.816496, 0., 0.);
        _3_fold_axis_1_by_2_fold_axis_1.selfNormalize();
        Matrix1D<double>  _2_fold_axis_1_by_2_fold_axis_2(3);
        _2_fold_axis_1_by_2_fold_axis_2 = vectorR3(0.707107, 0.408248, -0.57735);
        _2_fold_axis_1_by_2_fold_axis_2.selfNormalize();
        Matrix1D<double>  _2_fold_axis_2_by_3_fold_axis_1(3);
        _2_fold_axis_2_by_3_fold_axis_1 = vectorR3(-0.408248, -0.707107, 0.);
        _2_fold_axis_2_by_3_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            //           if ( XX(sampling_points_angles[i])>=     120. &&
            //                 XX(sampling_points_angles[i])<=   150. ||
            //                 XX(sampling_points_angles[i])==     0
            //              )
            if (
                dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_2_fold_axis_1) >= 0 &&
                dotProduct(sampling_points_vector[i], _2_fold_axis_1_by_2_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _2_fold_axis_2_by_3_fold_axis_1) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_O )
    {//OK
        Matrix1D<double>  _3_fold_axis_1_by_3_fold_axis_2(3);
        _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(0., -1., 1.);
        _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
        Matrix1D<double>  _3_fold_axis_2_by_4_fold_axis(3);
        _3_fold_axis_2_by_4_fold_axis = vectorR3(1., 1., 0.);
        _3_fold_axis_2_by_4_fold_axis.selfNormalize();
        Matrix1D<double>  _4_fold_axis_by_3_fold_axis_1(3);
        _4_fold_axis_by_3_fold_axis_1 = vectorR3(-1., 1., 0.);
        _4_fold_axis_by_3_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if ((XX(sampling_points_angles[i]) >=   45. &&
                 XX(sampling_points_angles[i]) <=  135. &&
                 YY(sampling_points_angles[i]) <=  90.) ||
                XX(sampling_points_angles[i]) ==  0.
               )
                if (
                    dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
                    dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_4_fold_axis) >= 0 &&
                    dotProduct(sampling_points_vector[i], _4_fold_axis_by_3_fold_axis_1) >= 0
                )
                {
                    no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                    no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
                }
        }// for i
    }
    else if (symmetry  == pg_OH )
    {//OK
        Matrix1D<double>  _3_fold_axis_1_by_3_fold_axis_2(3);
        _3_fold_axis_1_by_3_fold_axis_2 = vectorR3(0., -1., 1.);
        _3_fold_axis_1_by_3_fold_axis_2.selfNormalize();
        Matrix1D<double>  _3_fold_axis_2_by_4_fold_axis(3);
        _3_fold_axis_2_by_4_fold_axis = vectorR3(1., 1., 0.);
        _3_fold_axis_2_by_4_fold_axis.selfNormalize();
        Matrix1D<double>  _4_fold_axis_by_3_fold_axis_1(3);
        _4_fold_axis_by_3_fold_axis_1 = vectorR3(-1., 1., 0.);
        _4_fold_axis_by_3_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (XX(sampling_points_angles[i]) >=   90. &&
                XX(sampling_points_angles[i]) <=  135. &&
                YY(sampling_points_angles[i]) <=  90.)
                if (
                    dotProduct(sampling_points_vector[i], _3_fold_axis_1_by_3_fold_axis_2) >= 0 &&
                    dotProduct(sampling_points_vector[i], _3_fold_axis_2_by_4_fold_axis) >= 0 &&
                    dotProduct(sampling_points_vector[i], _4_fold_axis_by_3_fold_axis_1) >= 0
                )
                {
                    no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                    no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
                }
        }// for i
    }
    else if (symmetry  == pg_I || symmetry  == pg_I2)
    {//OK
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = vectorR3(0., 1., 0.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = vectorR3(-0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = vectorR3(0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I1)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0., 90., 0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();

        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I3)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0.,31.7174745559,0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();

        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_by_5_fold_axis_1) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I4)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0.,-31.7174745559,0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
                                        -0.303530987314591,
                                        -0.491123477863004);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
                                        0.303530987314591,
                                        -0.491123477863004);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();

        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_2_by_3_fold_axis)   <= 0 &&
                dotProduct(sampling_points_vector[i],
                           _3_fold_axis_by_5_fold_axis_1)   <= 0 &&
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_1_by_5_fold_axis_2) <= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I5)
    {//OK
        std::cerr << "ERROR: Symmetry pg_I5 not implemented" << std::endl;
        exit(0);
    }
    else if (symmetry  == pg_IH || symmetry  == pg_I2H)
    {//OK
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = vectorR3(0., 1., 0.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = vectorR3(-0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = vectorR3(0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_2_fold_axis(3);
        _3_fold_axis_by_2_fold_axis =  vectorR3(1.,0.,0.);
        _3_fold_axis_by_2_fold_axis.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_by_2_fold_axis) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I1H)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0., 90., 0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 1., 0.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(-0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.4999999839058737,
                                        -0.8090170074556163,
                                        0.3090169861701543);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_2_fold_axis(3);
        _3_fold_axis_by_2_fold_axis =  A * vectorR3(1.,0.,0.);
        _3_fold_axis_by_2_fold_axis.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i], _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i], _5_fold_axis_2_by_3_fold_axis) >= 0 &&
                dotProduct(sampling_points_vector[i], _3_fold_axis_by_2_fold_axis) >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I3H)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0.,31.7174745559,0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
                                        -0.303530987314591,
                                        -0.491123477863004);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
                                        0.303530987314591,
                                        -0.491123477863004);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_2_fold_axis(3);
        _3_fold_axis_by_2_fold_axis = vectorR3(0.,1.,0.);
        _3_fold_axis_by_2_fold_axis.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_2_by_3_fold_axis)   >= 0 &&
                dotProduct(sampling_points_vector[i],
                           _3_fold_axis_by_5_fold_axis_1)   >= 0 &&
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_1_by_5_fold_axis_2) >= 0 &&
                dotProduct(sampling_points_vector[i],
                           _3_fold_axis_by_2_fold_axis)     >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I4H)
    {//OK
        Matrix2D<double>  A(3, 3);
        Euler_angles2matrix(0.,-31.7174745559,0., A);
        Matrix1D<double>  _5_fold_axis_1_by_5_fold_axis_2(3);
        _5_fold_axis_1_by_5_fold_axis_2 = A * vectorR3(0., 0., 1.);
        _5_fold_axis_1_by_5_fold_axis_2.selfNormalize();
        Matrix1D<double>  _5_fold_axis_2_by_3_fold_axis(3);
        _5_fold_axis_2_by_3_fold_axis = A * vectorR3(0.187592467856686,
                                        -0.303530987314591,
                                        -0.491123477863004);
        _5_fold_axis_2_by_3_fold_axis.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_5_fold_axis_1(3);
        _3_fold_axis_by_5_fold_axis_1 = A * vectorR3(0.187592467856686,
                                        0.303530987314591,
                                        -0.491123477863004);
        _3_fold_axis_by_5_fold_axis_1.selfNormalize();
        Matrix1D<double>  _3_fold_axis_by_2_fold_axis(3);
        _3_fold_axis_by_2_fold_axis = vectorR3(0.,1.,0.);
        _3_fold_axis_by_2_fold_axis.selfNormalize();
        for (size_t i = 0; i < sampling_points_angles.size(); i++)
        {
            if (
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_2_by_3_fold_axis)   <= 0 &&
                dotProduct(sampling_points_vector[i],
                           _3_fold_axis_by_5_fold_axis_1)   <= 0 &&
                dotProduct(sampling_points_vector[i],
                           _5_fold_axis_1_by_5_fold_axis_2) <= 0 &&
                dotProduct(sampling_points_vector[i],
                           _3_fold_axis_by_2_fold_axis)     >= 0
            )
            {
                no_redundant_sampling_points_angles.push_back(sampling_points_angles[i]);
                no_redundant_sampling_points_vector.push_back(sampling_points_vector[i]);
            }
        }// for i
    }
    else if (symmetry  == pg_I5H)
    {//OK
        std::cerr << "ERROR: pg_I5H Symmetry not implemented" << std::endl;
        exit(0);
    }
    else
    {
        std::cerr << "ERROR: Symmetry " << symmetry  << "is not known" << std::endl;
        exit(0);
    }

    CREATE_INDEXES();
    numberSamplesAsymmetricUnit=no_redundant_sampling_points_vector.size();
    //#define DEBUG5
#ifdef  DEBUG5

    std::ofstream filestr5;
    filestr5.open ("removeRedundantPoints.bild");
    filestr5 << ".color 0 0 1" << std::endl;
    for (int i = 0;
         i < no_redundant_sampling_points_vector.size();
         i++)
    {
        filestr5  <<  ".color 1 0 0" << std::endl
        <<  ".sphere "
        << no_redundant_sampling_points_vector[i](0) << " "
        << no_redundant_sampling_points_vector[i](1) << " "
        << no_redundant_sampling_points_vector[i](2) << " "
        << " .03" << std::endl;
    }
    filestr5.close();
#endif
#undef DEBUG5

}

removeRedundantPointsExhaustive()

void Sampling::removeRedundantPointsExhaustive	(	const int	symmetry,
		int	sym_order,
		bool	only_half_sphere,
		double	max_ang
	)

Definition at line 1253 of file sampling.cpp.

{
    // Maximum distance
    double cos_max_ang = cos(DEG2RAD(max_ang));
    double my_dotProduct;
    //int j_end=0;
    Matrix1D<double>  direction(3), direction1(3);

    // First call to conventional removeRedundantPoints
    removeRedundantPoints(symmetry, sym_order);
    //    for (int isym = 0; isym < no_redundant_sampling_points_vector.size(); isym++)
    //         std::cout << "sampling:" << no_redundant_sampling_points_vector[isym];
    std::vector <Matrix1D<double> > old_vector = no_redundant_sampling_points_vector;
    std::vector <Matrix1D<double> > old_angles = no_redundant_sampling_points_angles;

    CLEAR_VECTORS();

    // Precalculate symmetry matrices
    fillLRRepository();

    // Then check all points versus each other
    for (size_t i = 0; i < old_angles.size(); i++)
    {
        //direction1=(old_vector[i]).transpose();
        direction1=old_vector[i];
        bool uniq = true;
        for (size_t j = 0; j < R_repository.size(); j++)
        {
            for (size_t k = 0; k < no_redundant_sampling_points_vector.size(); k++)
            {
                direction =  L_repository[j] *
                             (no_redundant_sampling_points_vector[k].transpose() *
                              R_repository[j]).transpose();
                //Calculate distance
                my_dotProduct = dotProduct(direction,direction1);
                if (only_half_sphere)
                    my_dotProduct = ABS(my_dotProduct);

                if (my_dotProduct > cos_max_ang)
                {
                    uniq = false;
                    break;
                }
            }// for k
            if (!uniq)
                break;
        } // for j
        if (uniq)
        {
            no_redundant_sampling_points_vector.push_back(old_vector[i]);
            no_redundant_sampling_points_angles.push_back(old_angles[i]);
        }
    } // for i

    CREATE_INDEXES();

}

saveSamplingFile()

void Sampling::saveSamplingFile	(	const FileName &	fn_base,
		bool	write_vectors = true,
		bool	write_sampling_sphere = false
	)

Save neighbors as three metadtada blocks 1) header with sampling rate and angular distance 2) one row for ach experimentald data with neighbours 3) sampling points

Definition at line 1495 of file sampling.cpp.

{
    {
        MetaDataVec md;
        MDRowVec row;
        row.setValue(MDL_SAMPLINGRATE, sampling_rate_rad);
        row.setValue(MDL_NEIGHBORHOOD_RADIUS, cos_neighborhood_radius);
        row.setValue(MDL_POINTSASYMETRICUNIT,numberSamplesAsymmetricUnit);
        md.setComment("data_extra -> sampling description;"\
                      " data_neighbors --> List with order of each"\
                      "experimental images and its neighbors"\
                     );
        md.setColumnFormat(false);
        md.addRow(row);
        md.write(FN_SAMPLING_EXTRA(fn_base), MD_OVERWRITE);
    }

    {
        MetaDataVec md;
        md.setColumnFormat(true);
        size_t size = my_neighbors.size();
        //Write first block with experimental images order and its neighbors
        bool writeFileName = !exp_data_fileNames.empty();
        for(size_t i = 0; i < size; ++i)
        {
            MDRowVec row;
            row.setValue(MDL_ORDER,i+FIRST_IMAGE);
            if (writeFileName)
                row.setValue(MDL_IMAGE,exp_data_fileNames[i]);
            row.setValue(MDL_NEIGHBORS, my_neighbors[i]);
            md.addRow(row);
        }

        md.write(FN_SAMPLING_NEI(fn_base), MD_APPEND);
    }

    {
        //Write projection directions
        MetaDataVec md;
        size_t size = no_redundant_sampling_points_index.size();

        for (size_t i = 0; i < size; ++i)
        {
            MDRowVec row;
            Matrix1D<double> &angles = no_redundant_sampling_points_angles[i];
            row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
            row.setValue(MDL_ANGLE_ROT, XX(angles));
            row.setValue(MDL_ANGLE_TILT, YY(angles));
            row.setValue(MDL_ANGLE_PSI, ZZ(angles));

            if (write_vectors)
            {
                Matrix1D<double> &vectors = no_redundant_sampling_points_vector[i];
                row.setValue(MDL_X, XX(vectors));
                row.setValue(MDL_Y, YY(vectors));
                row.setValue(MDL_Z, ZZ(vectors));
            }
            md.addRow(row);
        }
        md.write(FN_SAMPLING_PROJ(fn_base), MD_APPEND);
    }

    if (write_sampling_sphere)
    {
        MetaDataVec md;
        size_t size = sampling_points_angles.size();

        for (size_t i = 0; i < size; ++i)
        {
            MDRowVec row;
            row.setValue(MDL_NEIGHBOR, no_redundant_sampling_points_index[i]);
            Matrix1D<double> &angles = sampling_points_angles[i];
            row.setValue(MDL_ANGLE_ROT, XX(angles));
            row.setValue(MDL_ANGLE_TILT, YY(angles));
            row.setValue(MDL_ANGLE_PSI, ZZ(angles));


            if (write_vectors)
            {
                Matrix1D<double> &vectors = sampling_points_vector[i];
                row.setValue(MDL_X, XX(vectors));
                row.setValue(MDL_Y, YY(vectors));
                row.setValue(MDL_Z, ZZ(vectors));
            }
            md.addRow(row);
        }
        md.setComment("NEIGHBOR index points to the slice in a stack file in which the  projection projected in the direction defined by ref is stored");
        md.write(FN_SAMPLING_SPHERE(fn_base), MD_APPEND);

    }

}

setNeighborhoodRadius()

void Sampling::setNeighborhoodRadius

(

double

neighborhood

)

set neighborhood distance

Definition at line 140 of file sampling.cpp.

{
    if(neighborhood<0)
        cos_neighborhood_radius=-1.01;
    else if(neighborhood>180.001)
        REPORT_ERROR(ERR_ARG_INCORRECT,"Neighborhood can not be greater than 180. \n \
                     Use any negative value to cover the whole space ");
    else
    {
        neighborhood_radius_rad = DEG2RAD(neighborhood);
        cos_neighborhood_radius = cos(neighborhood_radius_rad);
    }
}

setNoise()

void Sampling::setNoise	(	double	deviation,
		int	my_seed = -1
	)

set sampling noise for projection vectors create in the unit sphere

Definition at line 134 of file sampling.cpp.

{
    sampling_noise = noise_deviation;
    init_random_generator(my_seed);
}

setSampling()

void Sampling::setSampling

(

double

sampling

)

set sampling rate

Definition at line 121 of file sampling.cpp.

{
    sampling_rate_rad = DEG2RAD(sampling);
    number_of_samples = ROUND(cte_w / sampling_rate_rad)+1;
    if (number_of_samples < 3)
    {
        std::cerr << "maximum value of angular sampling rate is "
        << cte_w*0.5*180./PI
        << std::endl;
        exit(1);
    }
}

sortFunc()

int Sampling::sortFunc	(	const Matrix1D< double > &	a,
		const Matrix1D< double > &	b
	)

Definition at line 586 of file sampling.cpp.

{
    if (YY(t) - 0.000001 > YY(a))
    {
        return 1;
    }
    else if (YY(t) + 0.000001 < YY(a))
    {
        return 0;
    };
    if (XX(t)  > XX(a))
    {
        return 1;
    }
    else
    {
        return 0;
    };
}

Member Data Documentation

cos_neighborhood_radius

double Sampling::cos_neighborhood_radius

cosine of neighborhood s

Definition at line 84 of file sampling.h.

exp_data_fileNames

std::vector<FileName > Sampling::exp_data_fileNames

vector with product of experimental images and L and R

Definition at line 110 of file sampling.h.

exp_data_projection_direction_by_L_R

std::vector<Matrix1D<double> > Sampling::exp_data_projection_direction_by_L_R

vector with product of experimental images and L and R

Definition at line 108 of file sampling.h.

L_repository

std::vector<Matrix2D<double> > Sampling::L_repository

Definition at line 106 of file sampling.h.

my_cross_correlation

std::vector<std::vector<double> > Sampling::my_cross_correlation

vector with angular distance between points (dot_product)

Definition at line 98 of file sampling.h.

my_exp_img_per_sampling_point

std::vector<std::vector<size_t> > Sampling::my_exp_img_per_sampling_point

vector with experimental images per sampling point

Definition at line 90 of file sampling.h.

my_neighbors

std::vector<std::vector<size_t> > Sampling::my_neighbors

vector with neighbors

Definition at line 87 of file sampling.h.

neighborhood_radius_rad

double Sampling::neighborhood_radius_rad

neighborhood in radians

Definition at line 81 of file sampling.h.

no_redundant_sampling_points_angles

std::vector<Matrix1D<double> > Sampling::no_redundant_sampling_points_angles

vector with sampling points described by angles, only store the non redundant part

Definition at line 121 of file sampling.h.

no_redundant_sampling_points_index

std::vector<size_t> Sampling::no_redundant_sampling_points_index

vector with the indexes of each sampling point

Definition at line 123 of file sampling.h.

no_redundant_sampling_points_vector

std::vector<Matrix1D<double> > Sampling::no_redundant_sampling_points_vector

vector with sampling points described by vectors, only store the non redundant part

Definition at line 118 of file sampling.h.

number_of_samples

size_t Sampling::number_of_samples

number of samples

Definition at line 75 of file sampling.h.

numberSamplesAsymmetricUnit

size_t Sampling::numberSamplesAsymmetricUnit

number of samples in the asymmetric unit

Definition at line 78 of file sampling.h.

R_repository

std::vector<Matrix2D<double> > Sampling::R_repository

vector with symmetry matrices

Definition at line 105 of file sampling.h.

sampling_noise

double Sampling::sampling_noise

sampling rate for the unit vectors

Definition at line 72 of file sampling.h.

sampling_points_angles

std::vector<Matrix1D<double> > Sampling::sampling_points_angles

vector with sampling points described by angles

Definition at line 103 of file sampling.h.

sampling_points_vector

std::vector<Matrix1D<double> > Sampling::sampling_points_vector

vector with sampling points described by vectors

Definition at line 101 of file sampling.h.

sampling_rate_rad

double Sampling::sampling_rate_rad

sampling rate in radians

Definition at line 69 of file sampling.h.

SL

SymList Sampling::SL

symmetry information

Definition at line 138 of file sampling.h.

symmetry_file

FileName Sampling::symmetry_file

symmetry file

Definition at line 135 of file sampling.h.

verbose

int Sampling::verbose

Verbose

Definition at line 126 of file sampling.h.

vertices_angles

Vect_angles Sampling::vertices_angles

Geographical co-ordinates of the home Vertex of the 10 diamonds as angles

Definition at line 62 of file sampling.h.

vertices_vectors

std::vector<Matrix1D<double> > Sampling::vertices_vectors

Geographical co-ordinates of the home Vertex of the 10 diamonds as vectors

Definition at line 66 of file sampling.h.

---

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

• xmipp/libraries/data/sampling.h
• xmipp/libraries/data/sampling.cpp