EMAN2
Public Types | Public Member Functions | Private Attributes
EMAN::PointArray Class Reference

PointArray defines a double array of points with values in a 3D space. More...

#include <pointarray.h>

Collaboration diagram for EMAN::PointArray:
Collaboration graph
[legend]

List of all members.

Public Types

enum  Density2PointsArrayAlgorithm { PEAKS_SUB, PEAKS_DIV, KMEANS }

Public Member Functions

 PointArray ()
 PointArray (int nn)
 ~PointArray ()
void zero ()
PointArraycopy ()
PointArrayoperator= (PointArray &pa)
size_t get_number_points () const
void set_number_points (size_t nn)
bool read_from_pdb (const char *file)
void save_to_pdb (const char *file)
FloatPoint get_center ()
void center_to_zero ()
Region get_bounding_box ()
double * get_points_array ()
Vec3f get_vector_at (int i)
double get_value_at (int i)
void set_vector_at (int i, Vec3f vec, double value)
void set_vector_at (int i, vector< double >)
void set_points_array (double *p)
vector< float > get_points ()
 Returns all x,y,z triplets packed into a vector<float>
EMDatadistmx (int sortrows=0)
 Calculates a (symmetrized) distance matrix for the current PointArray.
vector< int > match_points (PointArray *to, float max_miss=-1.0)
 Will try to establish a 1-1 correspondence between points in two different PointArray objects (this and to).
Transformalign_2d (PointArray *to, float max_dist)
 Aligns one PointArray to another in 2 dimensions.
vector< float > align_trans_2d (PointArray *to, int flags=0, float dxhint=0, float dyhint=0)
 Translationally aligns one PointArray to another in 2 dimensions.
void mask (double rmax, double rmin=0.0)
void mask_asymmetric_unit (const string &sym)
void transform (const Transform &transform)
void right_transform (const Transform &transform)
 Does Transform*v as opposed to v*Transform (as in the transform function)
void set_from (vector< float >)
void set_from (PointArray *source, const string &sym="", Transform *transform=0)
void set_from (double *source, int num, const string &sym="", Transform *transform=0)
void set_from_density_map (EMData *map, int num, float thresh, float apix, Density2PointsArrayAlgorithm mode=PEAKS_DIV)
void sort_by_axis (int axis=1)
EMDatapdb2mrc_by_nfft (int map_size, float apix, float res)
EMDatapdb2mrc_by_summation (int map_size, float apix, float res, int addpdbbfactor)
EMDataprojection_by_nfft (int image_size, float apix, float res=0)
EMDataprojection_by_summation (int image_size, float apix, float res)
void replace_by_summation (EMData *image, int i, Vec3f vec, float amp, float apix, float res)
void opt_from_proj (const vector< EMData * > &proj, float pixres)
 Optimizes a pointarray based on a set of projection images (EMData objects) This is effectively a 3D reconstruction algorithm.
double sim_pointpotential (double dist, double ang, double dihed)
 Computes a potential value for a single point from a set of angles/distances using current energy settings.
double sim_potential ()
 Computes overall potential for the configuration.
double sim_potentiald (int i)
 Compute a single point potential value.
double sim_potentialdxyz (int i, double dx, double dy, double dz)
 Compute a potential value for a perturbed point, including +-2 nearest neighbors which will also be impacted.
void sim_updategeom ()
 Updates the dist,ang,dihed parameters.
Vec3f sim_descent (int i)
 returns a vector pointing downhill for a single point
void sim_minstep (double maxshift)
 Takes a step to minimize the potential.
void sim_minstep_seq (double meanshift)
 Takes a step to minimize the potential.
void sim_rescale ()
 rescale the entire set so the mean bond length matches dist0
void sim_printstat ()
 prints some statistics to the screen
void sim_set_pot_parms (double pdist0, double pdistc, double pangc, double pdihed0, double pdihedc, double pmapc, EMData *pmap, double pmindistc, double pdistpenc)
 Sets the parameters for the energy function.
void sim_add_point_double ()
 Add more points during the simulation.
void sim_add_point_one ()
double calc_total_length ()
 Calculate total length.
vector< double > fit_helix (EMData *pmap, int minlength, float mindensity, vector< int > edge, int twodir)
 Fit helix from peptide chains.
vector< float > do_pca (int start, int end)
 Do principal component analysis to (a subset of) the point array.
vector< float > do_filter (vector< float > pts, float *ft, int num)
vector< double > construct_helix (int start, int end, float phs, float &score, bool &dir)
void reverse_chain ()
void save_pdb_with_helix (const char *file, vector< float > hlxid)
void remove_helix_from_map (EMData *m, vector< float > hlxid)

Private Attributes

double * points
size_t n
double * bfactor
double * adist
double * aang
double * adihed
double dist0
 Used for simplistic loop dynamics simulation Assumes all points are connected sequetially in a closed loop, potential includes distance, angle and dihedral terms, with optional denisty based terms.
double distc
double angc
double dihed0
double dihedc
double mapc
double apix
double mindistc
double distpenc
EMDatamap
EMDatagradx
EMDatagrady
EMDatagradz
vector< Vec3foldshifts

Detailed Description

PointArray defines a double array of points with values in a 3D space.

Definition at line 55 of file pointarray.h.


Member Enumeration Documentation

Enumerator:
PEAKS_SUB 
PEAKS_DIV 
KMEANS 

Definition at line 58 of file pointarray.h.


Constructor & Destructor Documentation

PointArray::PointArray ( )

Definition at line 108 of file pointarray.cpp.

References aang, adihed, adist, bfactor, gradx, grady, gradz, map, n, and points.

Referenced by copy().

{
        points = 0;
        bfactor = 0;
        n = 0;
        
        adist=0;
        aang=0;
        adihed=0;
        
        map=gradx=grady=gradz=0;
}
PointArray::PointArray ( int  nn) [explicit]

Definition at line 121 of file pointarray.cpp.

References aang, adihed, adist, gradx, grady, gradz, map, n, nn(), and points.

{
        n = nn;
        points = (double *) calloc(4 * n, sizeof(double));
        
        adist=0;
        aang=0;
        adihed=0;
        map=gradx=grady=gradz=0;
}
PointArray::~PointArray ( )

Definition at line 132 of file pointarray.cpp.

References aang, adihed, adist, gradx, grady, gradz, and points.

{
        if( points )
        {
                free(points);
                points = 0;
        }
        
        if (adist) free(adist);
        if (aang) free(aang);
        if (adihed) free(adihed);
//      if (map!=0) delete map;
        if (gradx!=0) delete gradx;
        if (grady!=0) delete grady;
        if (gradz!=0) delete gradz;
}

Member Function Documentation

Transform * PointArray::align_2d ( PointArray to,
float  max_dist 
)

Aligns one PointArray to another in 2 dimensions.

Parameters:
toAnother PointArray to align to
max_dist
Returns:
a Transform to map 'this' to 'to'

Definition at line 286 of file pointarray.cpp.

References EMAN::Util::calc_bilinear_least_square(), get_vector_at(), match_points(), EMAN::Transform::set(), and EMAN::Transform::set_pre_trans().

                                                             {
vector<int> match=match_points(to,max_dist);
Transform *ret=new Transform();

// we use bilinear least squares to get 3/6 matrix components
unsigned int i,j;

vector<float> pts;
for (i=0; i<match.size(); i++) {
        if (match[i]==-1) continue;

//      printf("%d -> %d\n",i,match[i]);
        pts.push_back(get_vector_at(i)[0]);
        pts.push_back(get_vector_at(i)[1]);
        pts.push_back(to->get_vector_at(match[i])[0]);
}

Vec3f vx=Util::calc_bilinear_least_square(pts);

// then we get the other 3/6
for (i=j=0; i<match.size(); i++) {
        if (match[i]==-1) continue;
        pts[j*3]  =get_vector_at(i)[0];
        pts[j*3+1]=get_vector_at(i)[1];
        pts[j*3+2]=to->get_vector_at(match[i])[1];
        j++;
}

Vec3f vy=Util::calc_bilinear_least_square(pts);

//ret->set_rotation(vx[1],vx[2],0.0,vy[1],vy[2],0.0,0.0,0.0,1.0);
//ret->set_rotation(vx[1],vy[1],0.0,vx[2],vy[2],0.0,0.0,0.0,1.0);
//ret->set_pretrans(Vec3f(-vx[0],-vy[0],0));

ret->set(0, 0, vx[1]);
ret->set(0, 1, vy[1]);
ret->set(0, 2, 0.0f);
ret->set(1, 0, vx[2]);
ret->set(1, 1, vy[2]);
ret->set(1, 2, 0.0f);
ret->set(2, 0, 0.0f);
ret->set(2, 1, 0.0f);
ret->set(2, 2, 1.0f);
ret->set_pre_trans(Vec3f(-vx[0],-vy[0],0));

return ret;
}
vector< float > PointArray::align_trans_2d ( PointArray to,
int  flags = 0,
float  dxhint = 0,
float  dyhint = 0 
)

Translationally aligns one PointArray to another in 2 dimensions.

Parameters:
toAnother PointArray to align to
flags
dxhint
dyhint
Returns:
A Pixel containing dx,dy and a quality factor (smaller better)

Definition at line 334 of file pointarray.cpp.

References get_number_points(), get_value_at(), get_vector_at(), and EMAN::length().

                                                                                             {
printf("Warning, this is old code. Use align_2d.\n");

// returns (dx,dy,residual error,n points used)
// dxhint,dyhint should translate this->to
// flags : 1 - use hint values, 2 - center by strongest point (highest 'value')
int na=   get_number_points();
int nb=to->get_number_points();
if (na<=0 || nb<=0) return vector<float>(4,0);

int *a2b = (int *)malloc(na*sizeof(int));
int *b2a = (int *)malloc(nb*sizeof(int));


// find unweighted centers
float cax,cay,cbx,cby;
int i,j;

if (flags&1) {
        cbx=dxhint;
        cby=dyhint;
        cax=cay=0;
}
else if (flags&2) {
        // find the 'a' list peak
        float hia=0.0f;
        int hina=0;
        for (i=0; i<na; i++) {
                if (get_value_at(i)>hia) { hia=static_cast<float>(get_value_at(i)); hina=i; }
        }
        cax=get_vector_at(hina)[0];
        cay=get_vector_at(hina)[1];

        // find the 'b' list peak
        float hib=0;
        int hinb=0;
        for (i=0; i<na; i++) {
                if (to->get_value_at(i)>hib) { hib=static_cast<float>(to->get_value_at(i)); hinb=i; }
        }
        cbx=to->get_vector_at(hinb)[0];
        cby=to->get_vector_at(hinb)[1];
}
else {
        cax=cay=cbx=cby=0;

        for (i=0; i<na; i++) { cax+=get_vector_at(i)[0]; cay+=get_vector_at(i)[1]; }
        cax/=(float)na;
        cay/=(float)na;

        for (i=0; i<nb; i++) { cbx+=to->get_vector_at(i)[0]; cby+=to->get_vector_at(i)[1]; }
        cbx/=(float)nb;
        cby/=(float)nb;
}

Vec3f offset(cbx-cax,cby-cay,0);

// find the nearest point for each x point, taking the estimated centers into account
for (i=0; i<na; i++) {
        float rmin=1.0e30f;
        for (j=0; j<nb; j++) {
                float r=(get_vector_at(i)+offset-to->get_vector_at(j)).length();
                if (r<rmin) { a2b[i]=j; rmin=r; }
        }
}

// find the nearest point for each y point
for (i=0; i<nb; i++) {
        float rmin=1.0e30f;
        for (j=0; j<na; j++) {
                float r=(get_vector_at(j)+offset-to->get_vector_at(i)).length();
                if (r<rmin) { b2a[i]=j; rmin=r; }
        }
}

// now keep only points where x->y matches y->x
for (i=0; i<na; i++) {
        if (a2b[i]<0) continue;
        if (b2a[a2b[i]]!=i) { printf(" #%d!=%d# ",b2a[a2b[i]],i);  b2a[a2b[i]]=-1; a2b[i]=-1; }
        printf("%d->%d  ",i,a2b[i]);
}
printf("\n");

for (i=0; i<nb; i++) {
        if (b2a[i]<0) continue;
        if (a2b[b2a[i]]!=i) { a2b[b2a[i]]=-1; b2a[i]=-1; }
        printf("%d->%d  ",i,b2a[i]);
}
printf("\n");

// Compute the average translation required to align the points
float dx=0,dy=0,dr=0,nr=0;
for (i=0; i<na; i++) {
        if (a2b[i]==-1) continue;
        dx+=to->get_vector_at(a2b[i])[0]-get_vector_at(i)[0];
        dy+=to->get_vector_at(a2b[i])[1]-get_vector_at(i)[1];
        nr+=1.0;
}
//printf("%f %f %f\n",dx,dy,nr);
if (nr<2) return vector<float>(4,0);
dx/=nr;
dy/=nr;

// Compute the residual error
for (i=0; i<na; i++) {
        if (i==-1  || a2b[i]==-1) continue;
        dr+=(to->get_vector_at(a2b[i])-get_vector_at(i)-Vec3f(dx,dy,0)).length();
}
dr/=nr;

free(a2b);
free(b2a);
vector<float> ret(4);
ret[0]=dx;
ret[1]=dy;
ret[2]=dr;
ret[3]=(float)nr;
return ret;
}
double PointArray::calc_total_length ( )

Calculate total length.

Definition at line 1210 of file pointarray.cpp.

References dist(), n, points, and sqrt().

                                    {
        double dist=0;
        for(int i=0; i<n; i++){
                int k=(i+1)%n;
                double d=(points[i*4]-points[k*4])*(points[i*4]-points[k*4])+(points[i*4+1]-points[k*4+1])*(points[i*4+1]-points[k*4+1])+(points[i*4+2]-points[k*4+2])*(points[i*4+2]-points[k*4+2]);
                d=sqrt(d);
                dist+=d;
        }
        return dist;
}
void PointArray::center_to_zero ( )

Definition at line 584 of file pointarray.cpp.

References get_center(), get_number_points(), and points.

{
        FloatPoint center = get_center();
        for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
                points[i] -= center[0];
                points[i + 1] -= center[1];
                points[i + 2] -= center[2];
        }
}
vector< double > PointArray::construct_helix ( int  start,
int  end,
float  phs,
float &  score,
bool &  dir 
)

Definition at line 1980 of file pointarray.cpp.

References abs, coveig(), do_pca(), eigval, eigvec, EMAN::EMData::get_attr(), EMAN::EMData::get_value_at(), EMAN::Vec3< Type >::length(), map, mean(), and points.

Referenced by fit_helix().

                                                                                               {
        // calculate length
        Vec3f d(points[end*4]-points[start*4],points[end*4+1]-points[start*4+1],points[end*4+2]-points[start*4+2]);
        double len=d.length();
        int nh=int(len/1.54)+2;
        vector<double> helix(nh*3);
        vector<float> eigvec=do_pca(start,end); 
        float eigval[3],vec[9];

        Util::coveig(3,&eigvec[0],eigval,vec);
        float maxeigv=0;
        int maxvi=-1;
        for(int i=0; i<3; i++){
                if(abs(eigval[i])>maxeigv){
                        maxeigv=abs(eigval[i]);
                        maxvi=i;
                }
        }
        dir=eigval[maxvi]>0;
//      printf("%f\t",eigval[maxvi]);
//      vector<float> eigv(eigvec,eigvec+sizeof(eigvec)/sizeof(float));
        // create helix
        helix[0]=0;helix[1]=0;helix[2]=0;
        helix[nh*3-3]=.0;helix[nh*3-2]=0;helix[nh*3-1]=len+.83;
        
        for (int i=0; i<nh-2; i++){
                if(dir){
                        helix[(i+1)*3+0]=cos(((phs+(100*i))*M_PI)/180)*2.3;
                        helix[(i+1)*3+1]=sin(((phs+(100*i))*M_PI)/180)*2.3;
                }
                else{
                        helix[(i+1)*3+1]=cos(((phs+(100*i))*M_PI)/180)*2.3;
                        helix[(i+1)*3+0]=sin(((phs+(100*i))*M_PI)/180)*2.3;
                }       
                helix[(i+1)*3+2]=i*1.54+.83;
        }
        // transform to correct position
        vector<double> pts(nh*3);
        float mean[3];
        for (int k=0; k<3; k++) mean[k]=0;
        for (int k=0; k<nh; k++){
                for (int l=0; l<3; l++){
                        pts[k*3+l]=helix[k*3+0]*eigvec[0*3+l]+helix[k*3+1]*eigvec[1*3+l]+helix[k*3+2]*eigvec[2*3+l];
                        mean[l]+=pts[k*3+l];
                }
        }
        for (int k=0; k<3; k++) mean[k]/=nh;
        for (int k=0; k<nh; k++){
                for (int l=0; l<3; l++){
                        pts[k*3+l]-=mean[l];
                }
        }
        for (int k=0; k<3; k++) mean[k]=0;
        for (int k=start; k<end; k++){
                for (int l=0; l<3; l++){
                        mean[l]+=points[k*4+l];
                }
        }
        for (int k=0; k<3; k++) mean[k]/=(end-start);   
        for (int k=0; k<nh; k++){
                for (int l=0; l<3; l++){
                        pts[k*3+l]+=mean[l];
                }
        }
        
        // correct direction
        Vec3f d1(pts[0]-points[start*4],pts[1]-points[start*4+1],pts[2]-points[start*4+2]);
        Vec3f d2(pts[0]-points[end*4],pts[1]-points[end*4+1],pts[2]-points[end*4+2]);
                
        if (d1.length()>d2.length()) { //do reverse
                double tmp;
                for (int i=0; i<nh/2; i++){
                        for(int j=0; j<3; j++){
                                tmp=pts[i*3+j];
                                pts[i*3+j]=pts[(nh-i-1)*3+j];
                                pts[(nh-i-1)*3+j]=tmp;
                        }
                }
        }
        
        // calculate score
//      int sx=map->get_xsize(),sy=map->get_ysize(),sz=map->get_zsize();
        int sx=0,sy=0,sz=0;
        float ax=map->get_attr("apix_x"),ay=map->get_attr("apix_y"),az=map->get_attr("apix_z");
        score=0;
        for (int i=1; i<nh-1; i++){
                score+=map->get_value_at(int(pts[i*3]/ax+sx/2),int(pts[i*3+1]/ay+sy/2),int(pts[i*3+2]/az+sz/2));
        }
        score/=(nh-2);
                
        return pts;

        
}
PointArray * PointArray::copy ( )

Definition at line 154 of file pointarray.cpp.

References get_number_points(), get_points_array(), PointArray(), and set_number_points().

Referenced by mask(), and mask_asymmetric_unit().

{
        PointArray *pa2 = new PointArray();
        pa2->set_number_points(get_number_points());
        double *pa2data = pa2->get_points_array();
        memcpy(pa2data, get_points_array(), sizeof(double) * 4 * get_number_points());

        return pa2;
}
EMData * PointArray::distmx ( int  sortrows = 0)

Calculates a (symmetrized) distance matrix for the current PointArray.

Parameters:
sortrowsif set, will sort the values in each row. The return will no longer be a true similarity matrix.
Returns:
An EMData object containing the similarity matrix

Definition at line 197 of file pointarray.cpp.

References cmp_float(), data, EMAN::EMData::get_data(), get_vector_at(), EMAN::length(), n, EMAN::EMData::set_value_at(), EMAN::EMData::to_zero(), and EMAN::EMData::update().

Referenced by match_points().

                                       {
if (n==0) return NULL;

unsigned int i,j;

EMData *ret= new EMData(n,n,1);
ret->to_zero();

for (i=0; i<n; i++) {
        for (j=i+1; j<n; j++) {
                float r=(get_vector_at(i)-get_vector_at(j)).length();
                ret->set_value_at(i,j,0,r);
                ret->set_value_at(j,i,0,r);
        }
}

if (sortrows) {
        float *data=ret->get_data();
        for (i=0; i<n; i++) qsort(&data[i*n],n,sizeof(float),cmp_float);
        ret->update();
}

return ret;
}
vector< float > PointArray::do_filter ( vector< float >  pts,
float *  ft,
int  num 
)

Definition at line 1646 of file pointarray.cpp.

Referenced by fit_helix().

                                                                        {
        // filter a 1D array
        vector<float> result(pts);
        for (uint i=0; i<pts.size(); i++)
                result[i]=0;
        for (int i=(num-1)/2; i<pts.size()-(num-1)/2; i++){
                for (int j=0; j<num; j++){
                        int k=i+j-(num-1)/2;
                        result[i]+=pts[k]*ft[j];
                }
        }
        return result;
        
}
vector< float > PointArray::do_pca ( int  start = 0,
int  end = -1 
)

Do principal component analysis to (a subset of) the point array.

Definition at line 1609 of file pointarray.cpp.

References coveig(), covmat, eigval, eigvec, mean(), n, and points.

Referenced by construct_helix(), and fit_helix().

                                                       {
        
        if (end==-1) end=n;
        float covmat[9],mean[3];
        for (int i=0; i<3; i++) mean[i]=0;
        for (int i=start; i<end; i++){
                for (int j=0; j<3; j++){
                        mean[j]+=points[i*4+j];
                }
        }
        for (int i=0; i<3; i++) mean[i]/=end-start;
        
        for (int i=0; i<3; i++){
                for (int j=0; j<3; j++){
                        if (j<i){
                                covmat[i*3+j]=covmat[j*3+i];
                        }
                        else{
                                covmat[i*3+j]=0;
                                for (int k=start; k<end; k++)
                                {
                                        covmat[i*3+j]+=(points[k*4+i]-mean[i])*(points[k*4+j]-mean[j]);
                                }
                        }
                        
//                      printf("%f\t",covmat[i*3+j]);
                }
//              printf("\n");
        }
        
        float eigval[3],eigvec[9];
        Util::coveig(3,covmat,eigval,eigvec);
        vector<float> eigv(eigvec,eigvec+sizeof(eigvec)/sizeof(float));
//      printf(" %f,%f,%f\n %f,%f,%f\n %f,%f,%f\n",eigvec[0],eigvec[1],eigvec[2],eigvec[3],eigvec[4],eigvec[5],eigvec[6],eigvec[7],eigvec[8]);
        return eigv;
}
vector< double > PointArray::fit_helix ( EMData pmap,
int  minlength = 13,
float  mindensity = 4,
vector< int >  edge = vector<int>(),
int  twodir = 0 
)

Fit helix from peptide chains.

Definition at line 1661 of file pointarray.cpp.

References abs, construct_helix(), dist(), do_filter(), do_pca(), eigvec, map, mean(), n, points, reverse_chain(), set_number_points(), sqrt(), and v.

{
        vector<float> hlxlen(n);
        vector<int> helix;
        map=pmap;
        float ft[7]={0.0044,0.0540,0.2420,0.3989,0.2420,0.0540,0.0044};
//      float ft[7]={0,0,0,1,0,0,0};

        // search for long rods in the point array globally
        
        for (int dir=twodir; dir<2; dir++){ 
                // search in both directions and combine the result
                if( twodir==0)
                  reverse_chain();
                for (uint i=0; i<n; i++){
                        vector<float> dist(50);
                        // for each point, search the following 50 points, find the longest rod
                        for (int len=5; len<50; len++){
                                int pos=i+len;
                                if (pos>=n || pos<0)    break;
                                vector<float> eigvec=do_pca(i,pos); // align the points 
                                vector<float> pts((len+1)*3);
                                float mean[3];
                                for (int k=0; k<3; k++) mean[k]=0;
                                for (int k=i; k<pos; k++){
                                        for (int l=0; l<3; l++){
                                                pts[(k-i)*3+l]=points[k*4+0]*eigvec[l*3+0]+points[k*4+1]*eigvec[l*3+1]+points[k*4+2]*eigvec[l*3+2];
                                                mean[l]+=pts[(k-i)*3+l];
                                        }
                                }
                                for (int k=0; k<3; k++) mean[k]/=len;
                                float dst=0;
                                // distance to the center axis
                                for (int k=0; k<len; k++){
                                        dst+=abs((pts[k*3]-mean[0])*(pts[k*3]-mean[0])+(pts[k*3+1]-mean[1])*(pts[k*3+1]-mean[1]));
                                }
                                dist[len]=1-dst/len/len;
                        }
                        
                        vector<float> nd=do_filter(dist,ft,7);
                        nd=do_filter(nd,ft,7);
                        // length of the first rod
                        for (int j=7; j<49; j++){
                                if(nd[j]>nd[j-1] && nd[j]>nd[j+1]){
                                        hlxlen[i]=j-6;
                                        break;
                                }
                        }
                        
                        if(hlxlen[i]>25) hlxlen[i]=0;
                }
                // filter the array before further process
//              hlxlen[50]=100;
//              for (int i=0; i<n; i++) printf("%d\t%f\n",i,hlxlen[i]);
//              for (int i=0; i<3; i++) hlxlen=do_filter(hlxlen,ft,7);
//              for (int i=0; i<n; i++) printf("%d\t%f\n",i,hlxlen[i]);
                vector<float> ishlx(n);
                int hlx=0;
                float up=minlength; // rod length threshold
                // record position of possible helixes
                for (uint i=0; i<n; i++){
                        if(hlx<=0){
                                if(hlxlen[i]>up){
                                        hlx=hlxlen[i];
                                        helix.push_back(i);
                                        helix.push_back(i+hlxlen[i]-5);
                                }
                                
                        }
                        else{
                                hlx--;
                                ishlx[i]=hlxlen[i];
                        }
                }
                // while counting reversely
                if(dir==0){
                        for (uint i=0; i<helix.size(); i++) helix[i]=n-1-helix[i];
                        for (uint i=0; i<helix.size()/2; i++){
                                int tmp=helix[i*2+1];
                                helix[i*2+1]=helix[i*2];
                                helix[i*2]=tmp;
                        }
                }
                        

        }

#ifdef DEBUG
        printf("potential helix counting from both sides: \n");
        for (uint i=0; i<helix.size()/2; i++){
                printf("%d\t%d\n",helix[i*2],helix[i*2+1]);
        }       
        printf("\n\n");
#endif


        // Combine the result from both side
        for (uint i=0; i<helix.size()/2; i++){
                int change=1;
                while(change==1){
                        change=0;
                        for (uint j=i+1; j<helix.size()/2; j++){
                                if(helix[j*2]==0) continue;
                                if(helix[j*2]-2<helix[i*2+1] && helix[j*2+1]+2>helix[i*2]){
                                        helix[i*2]=(helix[i*2]<helix[j*2])?helix[i*2]:helix[j*2];
                                        helix[i*2+1]=(helix[i*2+1]>helix[j*2+1])?helix[i*2+1]:helix[j*2+1];
                                        helix[j*2]=0;
                                        helix[j*2+1]=0;
                                        change=1;
                                }
                        }       
                }
        }
        
        vector<int> allhlx;
        int minid=1;
        while (minid>=0){
                int mins=10000;
                minid=-1;
                for (uint i=0;i<helix.size()/2; i++){
                        if(helix[i*2]<.1) continue;
                        if(helix[i*2]<mins){
                                mins=helix[i*2];
                                minid=i;
                        }
                }
                if(minid>=0){
                        allhlx.push_back(helix[minid*2]);
                        allhlx.push_back(helix[minid*2+1]);
                        helix[minid*2]=-1;
                }               
        }
        
#ifdef DEBUG
        printf("combined result: \n");  
        for (uint i=0; i<allhlx.size()/2; i++){
                printf("%d\t%d\n",allhlx[i*2],allhlx[i*2+1]);
        }       
        printf("\n\n");
#endif
        
        // local search to decide the start and end point of each helix
//      vector<float> allscore(allhlx.size()/2);
        for (uint i=0; i<allhlx.size()/2; i++){
                int sz=10;
                int start=allhlx[i*2]-sz,end=allhlx[i*2+1]+sz;
                start=start>0?start:0;
                end=end<n?end:n;
                float minscr=100000;
                int mj=0,mk=0;
                
                for (int j=start; j<end; j++){
                        for (int k=j+6; k<end; k++){
                                vector<float> eigvec=do_pca(j,k);
                                vector<float> pts((k-j)*3);
                                float mean[3];
                                for (int u=0; u<3; u++) mean[u]=0;
                                for (int u=j; u<k; u++){
                                        for (int v=0; v<3; v++){
                                                pts[(u-j)*3+v]=points[u*4+0]*eigvec[v*3+0]+points[u*4+1]*eigvec[v*3+1]+points[u*4+2]*eigvec[v*3+2];
                                                mean[v]+=pts[(u-j)*3+v];
                                        }
                                }
                                for (int u=0; u<3; u++) mean[u]/=(k-j);
                                float dst=0;
                                // distance to the center axis
                                for (int u=0; u<k-j; u++){
                                        dst+=sqrt((pts[u*3]-mean[0])*(pts[u*3]-mean[0])+(pts[u*3+1]-mean[1])*(pts[u*3+1]-mean[1]));
                                }
                                float len=k-j;
                                float scr=dst/len/len;
                                if (scr<minscr){
//                                      printf("%f\t%d\t%d\n",scr,j,k);
                                        minscr=scr;
                                        mj=j;
                                        mk=k;
                                }
                        }
                }

//              printf("%d\t%d\n",mj,mk);
                
                allhlx[i*2]=mj;
                allhlx[i*2+1]=mk;        
//              allscore[i]=minscr;
//              if (mk-mj>60)
//                      allscore[i]=100;
        }
        
        for (uint i=0; i<edge.size()/2; i++){
                allhlx.push_back(edge[i*2]);
                allhlx.push_back(edge[i*2+1]);
        }
        
        
        vector<int> allhlx2;
        minid=1;
        while (minid>=0){
                int mins=10000;
                minid=-1;
                for (uint i=0;i<allhlx.size()/2; i++){
                        if(allhlx[i*2]<.1) continue;
                        if(allhlx[i*2]<mins){
                                mins=allhlx[i*2];
                                minid=i;
                        }
                }
                if(minid>=0){
                        allhlx2.push_back(allhlx[minid*2]<allhlx[minid*2+1]?allhlx[minid*2]:allhlx[minid*2+1]);
                        allhlx2.push_back(allhlx[minid*2]>allhlx[minid*2+1]?allhlx[minid*2]:allhlx[minid*2+1]);
                        allhlx[minid*2]=-1;
                }               
        }
        allhlx=allhlx2;
        
#ifdef DEBUG
        printf("Fitted helixes: \n");
        for (uint i=0; i<allhlx.size()/2; i++){
                printf("%d\t%d\n",allhlx[i*2],allhlx[i*2+1]);
        }       
        printf("\n\n");
#endif
        // create ideal helix
        uint ia=0,ka=0;
        bool dir;
        vector<double> finalhlx;
        vector<double> hlxid;
        printf("Confirming helix... \n");
        while(ia<n){
                if (ia==allhlx[ka*2]){
//                      int sz=(allhlx[ka*2+1]-allhlx[ka*2])>10?5:(allhlx[ka*2+1]-allhlx[ka*2])/2;
                        int sz=3;
                        float score=0,maxscr=0;
                        float bestphs=0,phsscore=0,pscore=0;
                        int mi,mj;
                        for (int i=0; i<sz; i++){
                                for (int j=0; j<sz; j++){
                                        int start=allhlx[ka*2]+i,end=allhlx[ka*2+1]-j;
                                        phsscore=0;bestphs=-1;
                                        for (float phs=-180; phs<180; phs+=10){ //search for phase
                                                construct_helix(start,end,phs,pscore,dir);
                                                if (pscore>phsscore){
                                                        phsscore=pscore;
                                                        bestphs=phs;
                                                }
                                        }
//                                      printf("%f\t",bestphs);
                                        construct_helix(start,end,bestphs,score,dir);
                                        if (score>maxscr){
                                                maxscr=score;
                                                mi=i;
                                                mj=j;
                                        }
                                }
                        }
                        for (int i=0; i<mi; i++){                       
                                finalhlx.push_back(points[(i+ia)*4]);
                                finalhlx.push_back(points[(i+ia)*4+1]);
                                finalhlx.push_back(points[(i+ia)*4+2]);
                        }
                        int start=allhlx[ka*2]+mi,end=allhlx[ka*2+1]-mj;
                        printf("%d\t%d\t%f\t%d\t",start,end,maxscr,maxscr>mindensity);
                        if (maxscr>mindensity){
                                phsscore=0;
                                for (float phs=-180; phs<180; phs+=10){ //search for phase
                                                construct_helix(start,end,phs,pscore,dir);
                                                if (pscore>phsscore){
                                                        phsscore=pscore;
                                                        bestphs=phs;
                                                }
                                }
                                vector<double> pts=construct_helix(start,end,bestphs,score,dir);
                                int lendiff=end-start-pts.size()/3-2;
                                printf("%d\t",dir);
                                if (pts.size()/3-2>9 && lendiff>-5){
                                        
                                        hlxid.push_back(finalhlx.size()/3+1);
                                        printf("%d\t",finalhlx.size()/3+1);
                                        for (uint j=3; j<pts.size()-3; j++)
                                                finalhlx.push_back(pts[j]);
                                        hlxid.push_back(finalhlx.size()/3-2);
                                        printf("%d\t",finalhlx.size()/3-2);
                                        for (uint j=0; j<3; j++)
                                                hlxid.push_back(pts[j]);
                                        for (uint j=pts.size()-3; j<pts.size(); j++)
                                                hlxid.push_back(pts[j]);
                                        ia=end;
                                }
                                else{
                                        printf("%d\t",pts.size()/3-2);
                                }
                        }
                        printf("\n");
                        ka++;
                        while(allhlx[ka*2]<ia)
                                ka++;
                }
                else{
                        finalhlx.push_back(points[ia*4]);
                        finalhlx.push_back(points[ia*4+1]);
                        finalhlx.push_back(points[ia*4+2]);
                        ia++;                   
                }
        }
        
        set_number_points(finalhlx.size()/3);
        
        for (uint i=0; i<n; i++){
                for (uint j=0; j<3; j++)
                        points[i*4+j]=finalhlx[i*3+j];
                points[i*4+3]=0;
        }
                        

        
        printf("\n\n");
        return hlxid;
}
Region PointArray::get_bounding_box ( )

Definition at line 594 of file pointarray.cpp.

References get_number_points(), and points.

{
        double xmin, ymin, zmin;
        double xmax, ymax, zmax;
        xmin = xmax = points[0];
        ymin = ymax = points[1];
        zmin = zmax = points[2];
        for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
                if (points[i] > xmax)
                        xmax = points[i];
                if (points[i] < xmin)
                        xmin = points[i];
                if (points[i + 1] > ymax)
                        ymax = points[i + 1];
                if (points[i + 1] < ymin)
                        ymin = points[i + 1];
                if (points[i + 2] > zmax)
                        zmax = points[i + 2];
                if (points[i + 2] < zmin)
                        zmin = points[i + 2];
        }
        return Region(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin);
}
FloatPoint PointArray::get_center ( )

Definition at line 565 of file pointarray.cpp.

References get_number_points(), norm(), and points.

Referenced by center_to_zero().

{
        double xc, yc, zc;
        xc = yc = zc = 0.0;
        double norm = 0.0;
        for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
                xc += points[i] * points[i + 3];
                yc += points[i + 1] * points[i + 3];
                zc += points[i + 2] * points[i + 3];
                norm += points[i + 3];
        }
        if (norm <= 0) {
                fprintf(stderr, "Abnormal total value (%g) for PointArray, it should be >0\n", norm);
                return FloatPoint(0, 0, 0);
        }
        else
                return FloatPoint(xc / norm, yc / norm, zc / norm);
}
size_t PointArray::get_number_points ( ) const
vector< float > PointArray::get_points ( )

Returns all x,y,z triplets packed into a vector<float>

Returns:
All points packed into a vector<float>

Definition at line 706 of file pointarray.cpp.

References n, and points.

                                     {
vector<float> ret;
for (unsigned int i=0; i<n; i++) {
        ret.push_back((float)points[i*4]);
        ret.push_back((float)points[i*4+1]);
        ret.push_back((float)points[i*4+2]);
}

return ret;
}
double * PointArray::get_points_array ( )

Definition at line 187 of file pointarray.cpp.

References points.

Referenced by copy(), mask(), operator=(), set_from(), and set_from_density_map().

{
        return points;
}
double PointArray::get_value_at ( int  i)

Definition at line 2809 of file pointarray.cpp.

References points.

Referenced by align_trans_2d().

{
return points[i*4+3];
}
Vec3f PointArray::get_vector_at ( int  i)

Definition at line 2804 of file pointarray.cpp.

References points.

Referenced by align_2d(), align_trans_2d(), and distmx().

{
return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
}
void PointArray::mask ( double  rmax,
double  rmin = 0.0 
)

Definition at line 619 of file pointarray.cpp.

References copy(), get_number_points(), get_points_array(), points, set_number_points(), v, x, and y.

{
        double rmax2 = rmax * rmax, rmin2 = rmin * rmin;
        PointArray *tmp = this->copy();
        double *tmp_points = tmp->get_points_array();
         int count = 0;
        for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
                double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v =
                        tmp_points[i + 3];
                double r2 = x * x + y * y + z * z;
                if (r2 >= rmin2 && r2 <= rmax2) {
                        points[count * 4] = x;
                        points[count * 4 + 1] = y;
                        points[count * 4 + 2] = z;
                        points[count * 4 + 3] = v;
                        ++count;
                }
        }
        set_number_points(count);
        if( tmp )
        {
                delete tmp;
                tmp = 0;
        }
}
void PointArray::mask_asymmetric_unit ( const string &  sym)

Definition at line 646 of file pointarray.cpp.

References copy(), LOGERR, points, set_number_points(), sqrt(), v, x, and y.

{
        if (sym == "c1" || sym == "C1")
                return;                                 // do nothing for C1 symmetry
        double alt0 = 0, alt1 = M_PI, alt2 = M_PI;
        double az0 = 0, az1 = M_PI;
        if (sym[0] == 'c' || sym[0] == 'C') {
                int nsym = atoi(sym.c_str() + 1);
                az1 = 2.0 * M_PI / nsym / 2.0;
        }
        else if (sym[0] == 'd' || sym[0] == 'D') {
                int nsym = atoi(sym.c_str() + 1);
                alt1 = M_PI / 2.0;
                alt2 = alt1;
                az1 = 2.0 * M_PI / nsym / 2.0;
        }
        else if (sym == "icos" || sym == "ICOS") {
                alt1 = 0.652358139784368185995; // 5fold to 3fold
                alt2 = 0.55357435889704525151;  // half of edge ie. 5fold to 2fold along the edge
                az1 = 2.0 * M_PI / 5 / 2.0;
        }
        else {
                LOGERR("PointArray::set_to_asymmetric_unit(): sym = %s is not implemented yet",
                           sym.c_str());
                return;
        }
#ifdef DEBUG
        printf("Sym %s: alt0 = %8g\talt1 = %8g\talt2 = %8g\taz0 = %8g\taz1 = %8g\n", sym.c_str(), alt0*180.0/M_PI, alt1*180.0/M_PI, alt2*180.0/M_PI, az0*180.0/M_PI, az1*180.0/M_PI);
#endif

        PointArray *tmp = this->copy();
        double *tmp_points = tmp->get_points_array();
         int count = 0;
        for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
                double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v = tmp_points[i + 3];
                double az = atan2(y, x);
                double az_abs = fabs(az - az0);
                if (az_abs < (az1 - az0)) {
                        double alt_max = alt1 + (alt2 - alt1) * az_abs / (az1 - az0);
                        double alt = acos(z / sqrt(x * x + y * y + z * z));
                        if (alt < alt_max && alt >= alt0) {
#ifdef DEBUG
                                printf("Point %3lu: x,y,z = %8g,%8g,%8g\taz = %8g\talt = %8g\n",i/4,x,y,z,az*180.0/M_PI, alt*180.0/M_PI);
#endif
                                points[count * 4] = x;
                                points[count * 4 + 1] = y;
                                points[count * 4 + 2] = z;
                                points[count * 4 + 3] = v;
                                count++;
                        }
                }
        }
        set_number_points(count);
        if( tmp )
        {
                delete tmp;
                tmp = 0;
        }
}
vector< int > PointArray::match_points ( PointArray to,
float  max_miss = -1.0 
)

Will try to establish a 1-1 correspondence between points in two different PointArray objects (this and to).

Returns a vector<int> where the index is addresses the points in 'this' and the value addresses points in 'to'. A value of -1 means there was no match for that point.

Returns:
A vector<int> with the mapping of points from 'this' to 'to'. e.g. - ret[2]=5 means point 2 in 'this' matches point 5 in 'to'

Definition at line 222 of file pointarray.cpp.

References distmx(), EMAN::EMObject::f, EMAN::EMData::get_attr(), EMAN::EMData::get_value_at(), EMAN::EMData::get_xsize(), and n.

Referenced by align_2d().

                                                                  {
EMData *mx0=distmx(1);
EMData *mx1=to->distmx(1);
unsigned int n2=mx1->get_xsize();       // same as get_number_points on to

if (max_miss<0) max_miss=(float)mx0->get_attr("sigma")/10.0f;
//printf("max error %f\n",max_miss);



vector<int> ret(n,-1);
vector<float> rete(n,0.0);
unsigned int i,j,k,l;

if (!mx0 || !mx1) {
        if (mx0) delete mx0;
        if (mx1) delete mx1;
        return ret;
}

// i iterates over elements of 'this', j looks for a match in 'to'
// k and l iterate over the individual distances
for (i=0; i<n; i++) {
        int bestn=-1;                   // number of best match in mx1
        double bestd=1.0e38;            // residual error distance in best match
        for (j=0; j<n2; j++) {
                double d=0;
                int nd=0;
                for (k=l=0; k<n-1 && l<n2-1; k++,l++) {
                        float d1=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l,j));
                        float d2=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l,j));
                        float d3=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l+1,j));
                        float d4=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l+1,j));
                        if (d2<d1 && d4>d2) { l--; continue; }
                        if (d3<d1 && d4>d3) { k--; continue; }
                        d+=d1;
                        nd++;
                }
                d/=(float)nd;
//              printf("%d -> %d\t%f\t%d\n",i,j,d,nd);
                if (d<bestd) { bestd=d; bestn=j; }
        }
        ret[i]=bestn;
        rete[i]=static_cast<float>(bestd);
}

// This will remove duplicate assignments, keeping the best one
// also remove any assignments with large errors
for (i=0; i<n; i++) {
        for (j=0; j<n; j++) {
                if (rete[i]>max_miss) { ret[i]=-1; break; }
                if (i==j || ret[i]!=ret[j] || ret[i]==-1) continue;
                if (rete[i]>rete[j]) { ret[i]=-1; break; }
        }
}

delete mx0;
delete mx1;

return ret;
}
PointArray & PointArray::operator= ( PointArray pa)

Definition at line 164 of file pointarray.cpp.

References get_number_points(), get_points_array(), and set_number_points().

{
        if (this != &pa) {
                set_number_points(pa.get_number_points());
                memcpy(get_points_array(), pa.get_points_array(), sizeof(double) * 4 * get_number_points());
        }
        return *this;
}
void PointArray::opt_from_proj ( const vector< EMData * > &  proj,
float  pixres 
)

Optimizes a pointarray based on a set of projection images (EMData objects) This is effectively a 3D reconstruction algorithm.

Author:
Steve Ludtke 11/27/2004
Parameters:
projA vector of EMData objects containing projections with orientations
pixresSize of each Gaussian in pixels

Definition at line 2781 of file pointarray.cpp.

References get_number_points(), LOGWARN, and proj.

                                                                        {
#ifdef OPTPP
        optdata=proj;
        optobj=this;
        optpixres=pixres;

        FDNLF1 nlf(get_number_points()*4,calc_opt_proj,init_opt_proj);
//      NLF1 nlf(get_number_points()*4,init_opt_proj,calc_opt_projd);
        nlf.initFcn();

        OptCG opt(&nlf);
//      OptQNewton opt(&nlf);
        opt.setMaxFeval(2000);
        opt.optimize();
        opt.printStatus("Done");
#else
        (void)proj;             //suppress warning message
        (void)pixres;   //suppress warning message
        LOGWARN("OPT++ support not enabled.\n");
        return;
#endif
}
EMData * PointArray::pdb2mrc_by_nfft ( int  map_size,
float  apix,
float  res 
)

Definition at line 2411 of file pointarray.cpp.

References apix, data, EMAN::EMData::do_ift(), EMAN::EMData::get_data(), get_number_points(), LOGWARN, map, norm(), points, EMAN::EMData::process_inplace(), EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), EMAN::EMData::to_zero(), and EMAN::EMData::update().

{
#if defined NFFT
        nfft_3D_plan my_plan;           // plan for the nfft

        nfft_3D_init(&my_plan, map_size, get_number_points());

        for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
                // FFTW and nfft use row major array layout, EMAN uses column major
                my_plan.v[3 * j + 2] = (fftw_real) (points[i] / (apix * map_size));
                my_plan.v[3 * j + 1] = (fftw_real) (points[i + 1] / (apix * map_size));
                my_plan.v[3 * j] = (fftw_real) (points[i + 2] / (apix * map_size));
                my_plan.f[j].re = (fftw_real) (points[i + 3]);
                my_plan.f[j].im = 0.0;
        }

        if (my_plan.nfft_flags & PRE_PSI) {
                nfft_3D_precompute_psi(&my_plan);
        }

        // compute the uniform Fourier transform
        nfft_3D_transpose(&my_plan);

        // copy the Fourier transform to EMData data array
        EMData *fft = new EMData();
        fft->set_size(map_size + 2, map_size, map_size);
        fft->set_complex(true);
        fft->set_ri(true);
        fft->to_zero();
        float *data = fft->get_data();
        double norm = 1.0 / (map_size * map_size * map_size);
        for (int k = 0; k < map_size; k++) {
                for (int j = 0; j < map_size; j++) {
                        for (int i = 0; i < map_size / 2; i++) {
                                data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
                                        (float) (my_plan.
                                                         f_hat[k * map_size * map_size + j * map_size + i +
                                                                   map_size / 2].re) * norm;
                                data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
                                        (float) (my_plan.
                                                         f_hat[k * map_size * map_size + j * map_size + i +
                                                                   map_size / 2].im) * norm * (-1.0);
                        }
                }
        }
        nfft_3D_finalize(&my_plan);

        // low pass processor
        double sigma2 = (map_size * apix / res) * (map_size * apix / res);
        int index = 0;
        for (int k = 0; k < map_size; k++) {
                double RZ2 = (k - map_size / 2) * (k - map_size / 2);
                for (int j = 0; j < map_size; j++) {
                        double RY2 = (j - map_size / 2) * (j - map_size / 2);
                        for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
                                float val = exp(-(i * i + RY2 + RZ2) / sigma2);
                                data[index] *= val;
                                data[index + 1] *= val;
                        }
                }
        }
        fft->update();
        //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));

        fft->process_inplace("xform.phaseorigin.tocorner");     // move phase origin to center of image map_size, instead of at corner
        EMData *map = fft->do_ift();
        map->set_attr("apix_x", apix);
        map->set_attr("apix_y", apix);
        map->set_attr("apix_z", apix);
        map->set_attr("origin_x", -map_size/2*apix);
        map->set_attr("origin_y", -map_size/2*apix);
        map->set_attr("origin_z", -map_size/2*apix);
        if( fft )
        {
                delete fft;
                fft = 0;
        }
        return map;
#elif defined NFFT2
        nfft_plan my_plan;                      // plan for the nfft

        nfft_init_3d(&my_plan, map_size, map_size, map_size, get_number_points());

        for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
                // FFTW and nfft use row major array layout, EMAN uses column major
                my_plan.x[3 * j + 2] = (double) (points[i] / (apix * map_size));
                my_plan.x[3 * j + 1] = (double) (points[i + 1] / (apix * map_size));
                my_plan.x[3 * j] = (double) (points[i + 2] / (apix * map_size));
                my_plan.f[j][0] = (double) (points[i + 3]);
                my_plan.f[j][1] = 0.0;
        }

        if (my_plan.nfft_flags & PRE_PSI) {
                nfft_precompute_psi(&my_plan);
                if (my_plan.nfft_flags & PRE_FULL_PSI)
                        nfft_full_psi(&my_plan, pow(10, -10));
        }

        // compute the uniform Fourier transform
        nfft_transposed(&my_plan);

        // copy the Fourier transform to EMData data array
        EMData *fft = new EMData();
        fft->set_size(map_size + 2, map_size, map_size);
        fft->set_complex(true);
        fft->set_ri(true);
        fft->to_zero();
        float *data = fft->get_data();
        double norm = 1.0 / (map_size * map_size * map_size);
        for (int k = 0; k < map_size; k++) {
                for (int j = 0; j < map_size; j++) {
                        for (int i = 0; i < map_size / 2; i++) {
                                data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
                                        (float) (my_plan.
                                                         f_hat[k * map_size * map_size + j * map_size + i +
                                                                   map_size / 2][0]) * norm;
                                data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
                                        (float) (my_plan.
                                                         f_hat[k * map_size * map_size + j * map_size + i +
                                                                   map_size / 2][1]) * norm;
                        }
                }
        }
        nfft_finalize(&my_plan);

        // low pass processor
        double sigma2 = (map_size * apix / res) * (map_size * apix / res);
        int index = 0;
        for (int k = 0; k < map_size; k++) {
                double RZ2 = (k - map_size / 2) * (k - map_size / 2);
                for (int j = 0; j < map_size; j++) {
                        double RY2 = (j - map_size / 2) * (j - map_size / 2);
                        for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
                                float val = exp(-(i * i + RY2 + RZ2) / sigma2);
                                data[index] *= val;
                                data[index + 1] *= val;
                        }
                }
        }
        fft->update();
        //fft->process_inplace(".filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));

        fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image map_size, instead of at corner
        EMData *map = fft->do_ift();
        map->set_attr("apix_x", apix);
        map->set_attr("apix_y", apix);
        map->set_attr("apix_z", apix);
        map->set_attr("origin_x", -map_size / 2 * apix);
        map->set_attr("origin_y", -map_size / 2 * apix);
        map->set_attr("origin_z", -map_size / 2 * apix);
        if( fft )
        {
                delete fft;
                fft = 0;
        }
        return map;
#else
        LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
        return 0;
#endif
}
EMData * PointArray::pdb2mrc_by_summation ( int  map_size,
float  apix,
float  res,
int  addpdbbfactor 
)

Definition at line 2148 of file pointarray.cpp.

References bfactor, EMAN::EMData::get_data(), get_number_points(), log(), map, points, EMAN::EMData::set_attr(), EMAN::EMData::set_size(), sqrt(), EMAN::EMData::to_zero(), EMAN::EMData::update(), and x.

{
#ifdef DEBUG
        printf("PointArray::pdb2mrc_by_summation(): %lu points\tmapsize = %4d\tapix = %g\tres = %g\n",get_number_points(),map_size, apix, res);
#endif
        //if ( gauss_real_width < apix) LOGERR("PointArray::projection_by_summation(): apix(%g) is too large for resolution (%g Angstrom in Fourier space) with %g pixels of 1/e half width", apix, res, gauss_real_width);

        double min_table_val = 1e-7;
        double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)

        double table_step_size = 0.001; // number of steps for each pixel
        double inv_table_step_size = 1.0 / table_step_size;
        
//      sort_by_axis(2);                        // sort by Z-axis

        EMData *map = new EMData();
        map->set_size(map_size, map_size, map_size);
        map->to_zero();
        float *pd = map->get_data();
        
        vector<double> table;
        double gauss_real_width;
        int table_size;
        int gbox;
        
        
        for ( size_t s = 0; s < get_number_points(); ++s) {
                double xc = points[4 * s] / apix + map_size / 2;
                double yc = points[4 * s + 1] / apix + map_size / 2;
                double zc = points[4 * s + 2] / apix + map_size / 2;
                double fval = points[4 * s + 3];
                
                //printf("\n bfactor=%f",bfactor[s]);
                
                
                
                
                
                if(addpdbbfactor==-1){
                        gauss_real_width = res/M_PI;    // in Angstrom, res is in Angstrom
                }else{
                        gauss_real_width = (bfactor[s])/(4*sqrt(2.0)*M_PI);     // in Angstrom, res is in Angstrom
                }
                
                
                table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
                table.resize(table_size);
                for (int i = 0; i < table_size; i++){
                        table[i] = 0;
                }
                
                for (int i = 0; i < table_size; i++) {                                          
                        double x = -i * table_step_size * apix / gauss_real_width;
                        if(addpdbbfactor==-1){
                                table[i] =  exp(-x * x);        
                        }
                        else{
                        table[i] =  exp(-x * x)/sqrt(gauss_real_width * gauss_real_width * 2 * M_PI);
                        }
                }

                gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
                if (gbox <= 0)
                        gbox = 1;
                
                int imin = int (xc) - gbox, imax = int (xc) + gbox;
                int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
                int kmin = int (zc) - gbox, kmax = int (zc) + gbox;
                if (imin < 0)
                        imin = 0;
                if (jmin < 0)
                        jmin = 0;
                if (kmin < 0)
                        kmin = 0;
                if (imax > map_size)
                        imax = map_size;
                if (jmax > map_size)
                        jmax = map_size;
                if (kmax > map_size)
                        kmax = map_size;                
                
                for (int k = kmin; k < kmax; k++) {
                        size_t table_index_z = size_t (fabs(k - zc) * inv_table_step_size);
                        if ( table_index_z >= table.size() ) continue;
                        double zval = table[table_index_z];
                        size_t pd_index_z = k * map_size * map_size;
                        
                        for (int j = jmin; j < jmax; j++) {
                                size_t table_index_y = size_t (fabs(j - yc) * inv_table_step_size);
                                if ( table_index_y >= table.size() ) continue;
                                double yval = table[table_index_y];
                                size_t pd_index = pd_index_z + j * map_size + imin;
                                for (int i = imin; i < imax; i++, pd_index++) {
                                        size_t table_index_x = size_t (fabs(i - xc) * inv_table_step_size);
                                        if ( table_index_x >= table.size() ) continue;
                                        double xval = table[table_index_x];
                                        pd[pd_index] += (float) (fval * zval * yval * xval);
                                }
                        }
                }
        }

        map->update();
        map->set_attr("apix_x", apix);
        map->set_attr("apix_y", apix);
        map->set_attr("apix_z", apix);
        map->set_attr("origin_x", -map_size/2*apix);
        map->set_attr("origin_y", -map_size/2*apix);
        map->set_attr("origin_z", -map_size/2*apix);

        return map;
}
EMData * PointArray::projection_by_nfft ( int  image_size,
float  apix,
float  res = 0 
)

Definition at line 2581 of file pointarray.cpp.

References apix, data, EMAN::EMData::get_data(), get_number_points(), LOGWARN, n, norm(), points, EMAN::EMData::process_inplace(), EMAN::EMData::set_complex(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), EMAN::EMData::to_zero(), and EMAN::EMData::update().

{
#if defined NFFT
        nfft_2D_plan my_plan;           // plan for the nfft
        int N[2], n[2];
        N[0] = image_size;
        n[0] = next_power_of_2(2 * image_size);
        N[1] = image_size;
        n[1] = next_power_of_2(2 * image_size);

        nfft_2D_init(&my_plan, image_size, get_number_points());
        //nfft_2D_init_specific(&my_plan, N, get_number_points(), n, 3,
        //                 PRE_PHI_HUT | PRE_PSI,
        //                 FFTW_ESTIMATE | FFTW_OUT_OF_PLACE);
        for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
                // FFTW and nfft use row major array layout, EMAN uses column major
                my_plan.v[2 * j + 1] = (fftw_real) (points[i] / (apix * image_size));
                my_plan.v[2 * j] = (fftw_real) (points[i + 1] / (apix * image_size));
                my_plan.f[j].re = (fftw_real) points[i + 3];
                my_plan.f[j].im = 0.0;
        }

        if (my_plan.nfft_flags & PRE_PSI) {
                nfft_2D_precompute_psi(&my_plan);
        }

        // compute the uniform Fourier transform
        nfft_2D_transpose(&my_plan);

        // copy the Fourier transform to EMData data array
        EMData *fft = new EMData();
        fft->set_size(image_size + 2, image_size, 1);
        fft->set_complex(true);
        fft->set_ri(true);
        fft->to_zero();
        float *data = fft->get_data();
        double norm = 1.0 / (image_size * image_size);
        for (int j = 0; j < image_size; j++) {
                for (int i = 0; i < image_size / 2; i++) {
                        data[j * (image_size + 2) + 2 * i] =
                                (float) (my_plan.f_hat[j * image_size + i + image_size / 2].re) * norm;
                        data[j * (image_size + 2) + 2 * i + 1] =
                                (float) (my_plan.f_hat[j * image_size + i + image_size / 2].im) * norm * (-1.0);
                }
        }
        nfft_2D_finalize(&my_plan);

        if (res > 0) {
                // Gaussian low pass processor
                double sigma2 = (image_size * apix / res) * (image_size * apix / res);
                int index = 0;
                for (int j = 0; j < image_size; j++) {
                        double RY2 = (j - image_size / 2) * (j - image_size / 2);
                        for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
                                double val = exp(-(i * i + RY2) / sigma2);
                                data[index] *= val;
                                data[index + 1] *= val;
                        }
                }
        }
        fft->update();
        //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));

        fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner

        return fft;
#elif defined NFFT2
        nfft_plan my_plan;                      // plan for the nfft
        int N[2], n[2];
        N[0] = image_size;
        n[0] = next_power_of_2(2 * image_size);
        N[1] = image_size;
        n[1] = next_power_of_2(2 * image_size);

        //nfft_init_2d(&my_plan,image_size,image_size,get_number_points());
        nfft_init_specific(&my_plan, 2, N, get_number_points(), n, 3,
                                           PRE_PHI_HUT | PRE_PSI |
                                           MALLOC_X | MALLOC_F_HAT | MALLOC_F, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
                // FFTW and nfft use row major array layout, EMAN uses column major
                my_plan.x[2 * j + 1] = (double) (points[i] / (apix * image_size));
                my_plan.x[2 * j] = (double) (points[i + 1] / (apix * image_size));
                my_plan.f[j][0] = (double) points[i + 3];
                my_plan.f[j][1] = 0.0;
        }

        if (my_plan.nfft_flags & PRE_PSI) {
                nfft_precompute_psi(&my_plan);
                if (my_plan.nfft_flags & PRE_FULL_PSI)
                        nfft_full_psi(&my_plan, pow(10, -6));
        }

        // compute the uniform Fourier transform
        nfft_transposed(&my_plan);

        // copy the Fourier transform to EMData data array
        EMData *fft = new EMData();
        fft->set_size(image_size + 2, image_size, 1);
        fft->set_complex(true);
        fft->set_ri(true);
        fft->to_zero();
        float *data = fft->get_data();
        double norm = 1.0 / (image_size * image_size);
        for (int j = 0; j < image_size; j++) {
                for (int i = 0; i < image_size / 2; i++) {
                        data[j * (image_size + 2) + 2 * i] =
                                (float) (my_plan.f_hat[j * image_size + i + image_size / 2][0]) * norm;
                        data[j * (image_size + 2) + 2 * i + 1] =
                                (float) (my_plan.f_hat[j * image_size + i + image_size / 2][1]) * norm;
                }
        }
        nfft_finalize(&my_plan);

        if (res > 0) {
                // Gaussian low pass process
                double sigma2 = (image_size * apix / res) * (image_size * apix / res);
                int index = 0;
                for (int j = 0; j < image_size; j++) {
                        double RY2 = (j - image_size / 2) * (j - image_size / 2);
                        for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
                                double val = exp(-(i * i + RY2) / sigma2);
                                data[index] *= val;
                                data[index + 1] *= val;
                        }
                }
        }
        fft->update();
        //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));

        fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner

        return fft;
#else
        LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
        return 0;
#endif
}
EMData * PointArray::projection_by_summation ( int  image_size,
float  apix,
float  res 
)

Definition at line 2262 of file pointarray.cpp.

References EMAN::EMData::get_data(), get_number_points(), log(), points, proj, EMAN::EMData::set_size(), sqrt(), EMAN::EMData::to_zero(), EMAN::EMData::update(), and x.

{
        double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
        //if ( gauss_real_width < apix) LOGERR("PointArray::projection_by_summation(): apix(%g) is too large for resolution (%g Angstrom in Fourier space) with %g pixels of 1/e half width", apix, res, gauss_real_width);

        double min_table_val = 1e-7;
        double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)

        //double table_step_size = 0.001;    // number of steps for x=[0,1] in exp(-x*x)
        //int table_size = int(max_table_x / table_step_size *1.25);
        //double* table = (double*)malloc(sizeof(double) * table_size);
        //for(int i=0; i<table_size; i++) table[i]=exp(-i*i*table_step_size*table_step_size);

        double table_step_size = 0.001; // number of steps for each pixel
        double inv_table_step_size = 1.0 / table_step_size;
        int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
        double *table = (double *) malloc(sizeof(double) * table_size);
        for (int i = 0; i < table_size; i++) {
                double x = -i * table_step_size * apix / gauss_real_width;
                table[i] = exp(-x * x);
        }

        int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
        if (gbox <= 0)
                gbox = 1;
        EMData *proj = new EMData();
        proj->set_size(image_size, image_size, 1);
        proj->to_zero();
        float *pd = proj->get_data();
        for ( size_t s = 0; s < get_number_points(); ++s) {
                double xc = points[4 * s] / apix + image_size / 2;
                double yc = points[4 * s + 1] / apix + image_size / 2;
                double fval = points[4 * s + 3];
                int imin = int (xc) - gbox, imax = int (xc) + gbox;
                int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
                if (imin < 0)
                        imin = 0;
                if (jmin < 0)
                        jmin = 0;
                if (imax > image_size)
                        imax = image_size;
                if (jmax > image_size)
                        jmax = image_size;

                for (int j = jmin; j < jmax; j++) {
                        //int table_index_y = int(fabs(j-yc)*apix/gauss_real_width/table_step_size);
                        int table_index_y = int (fabs(j - yc) * inv_table_step_size);
                        double yval = table[table_index_y];
#ifdef DEBUG
                        //double yval2 = exp( - (j-yc)*(j-yc)*apix*apix/(gauss_real_width*gauss_real_width));
                        //if(fabs(yval2-yval)/yval2>1e-2) printf("s=%d\txc,yc=%g,%g\tyval,yval2=%g,%g\tdiff=%g\n",s,xc,yc,yval,yval2,fabs(yval2-yval)/yval2);
#endif
                        int pd_index = j * image_size + imin;
                        for (int i = imin; i < imax; i++, pd_index++) {
                                //int table_index_x = int(fabs(i-xc)*apix/gauss_real_width/table_step_size);
                                int table_index_x = int (fabs(i - xc) * inv_table_step_size);
                                double xval = table[table_index_x];
#ifdef DEBUG
                                //double xval2 = exp( - (i-xc)*(i-xc)*apix*apix/(gauss_real_width*gauss_real_width));
                                //if(fabs(xval2-xval)/xval2>1e-2) printf("\ts=%d\txc,yc=%g,%g\txval,xval2=%g,%g\tdiff=%g\n",s,xc,yc,xval,xval2,fabs(xval2-xval)/xval2);
#endif
                                pd[pd_index] += (float)(fval * yval * xval);
                        }
                }
        }
        for (int i = 0; i < image_size * image_size; i++)
                pd[i] /= sqrt(M_PI);
        proj->update();
        return proj;
}
bool PointArray::read_from_pdb ( const char *  file)

Definition at line 454 of file pointarray.cpp.

References bfactor, get_number_points(), points, q, set_number_points(), x, and y.

Referenced by EMAN::PDBReader::makePointArray().

{
        struct stat filestat;
        stat(file, &filestat);
        set_number_points(( int)(filestat.st_size / 80 + 1));

        #ifdef DEBUG
        printf("PointArray::read_from_pdb(): try %4lu atoms first\n", get_number_points());
        #endif

        FILE *fp = fopen(file, "r");
        if(!fp) {
                fprintf(stderr,"ERROR in PointArray::read_from_pdb(): cannot open file %s\n",file);
                throw;
        }
        char s[200];
        size_t count = 0;
        
        while ((fgets(s, 200, fp) != NULL)) {
                if (strncmp(s, "ENDMDL", 6) == 0)
                        break;
                if (strncmp(s, "ATOM", 4) != 0)
                        continue;

                if (s[13] == ' ')
                        s[13] = s[14];
                if (s[13] == ' ')
                        s[13] = s[15];

                float e = 0;
                char ctt, ctt2 = ' ';
                if (s[13] == ' ')
                        ctt = s[14];
                else if (s[12] == ' ') {
                        ctt = s[13];
                        ctt2 = s[14];
                }
                else {
                        ctt = s[12];
                        ctt2 = s[13];
                }

                switch (ctt) {
                case 'H':
                        e = 1.0;
                        break;
                case 'C':
                        e = 6.0;
                        break;
                case 'A':
                        if (ctt2 == 'U') {
                                e = 79.0;
                                break;
                        }
                        // treat 'A'mbiguous atoms as N, not perfect, but good enough
                case 'N':
                        e = 7.0;
                        break;
                case 'O':
                        e = 8.0;
                        break;
                case 'P':
                        e = 15.0;
                        break;
                case 'S':
                        e = 16.0;
                        break;
                case 'W':
                        e = 18.0;
                        break;                          // ficticious water 'atom'
                default:
                        fprintf(stderr, "Unknown atom %c%c\n", ctt, ctt2);
                        e = 0;
                }
                if (e == 0)
                        continue;

                float x, y, z, q;
                sscanf(&s[28], " %f %f %f", &x, &y, &z);
                sscanf(&s[60], " %f", &q);

                if (count + 1 > get_number_points())
                        set_number_points(2 * (count + 1));    //makes sure point array is big enough
                
#ifdef DEBUG
                printf("Atom %4lu: x,y,z = %8g,%8g,%8g\te = %g\n", count, x, y, z, e);
#endif
                points[4 * count] = x;
                points[4 * count + 1] = y;
                points[4 * count + 2] = z;
                points[4 * count + 3] = e;
                bfactor[count] = q;
                count++;
        }
        fclose(fp);
        set_number_points(count);
        return true;
}
void PointArray::remove_helix_from_map ( EMData m,
vector< float >  hlxid 
)

Definition at line 2091 of file pointarray.cpp.

References EMAN::cross(), EMAN::dot(), EMAN::EMData::get_attr(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::length(), EMAN::Vec3< Type >::length(), EMAN::EMData::set_value_at(), t, x, and y.

                                                                    {
        
        int sx=m->get_xsize(),sy=m->get_ysize(),sz=m->get_zsize();
        float ax=m->get_attr("apix_x"),ay=m->get_attr("apix_y"),az=m->get_attr("apix_z");
        for (int x=0; x<sx; x++){
                for (int y=0; y<sy; y++){
                        for (int z=0; z<sz; z++){
                                Vec3f p0((x)*ax,(y)*ay,(z)*az);
                                bool inhlx=false;
                                for (uint i=0; i<hlxid.size()/8; i++){
                                        Vec3f p1(hlxid[i*8+2],hlxid[i*8+3],hlxid[i*8+4]),p2(hlxid[i*8+5],hlxid[i*8+6],hlxid[i*8+7]);
                                        Vec3f dp=p2-p1;
                                        float l=dp.length();
                                        float d=((p0-p1).cross(p0-p2)).length()/l;
                                        float t=-(p1-p0).dot(p2-p1)/(l*l);
                                        if (d<5 && t>0 && t<1){
                                                inhlx=true;
                                                break;
                                        }
                                }
                                if(inhlx){
                                        m->set_value_at(x,y,z,0);
                                }
                                        
                        }
                }
        }
}
void PointArray::replace_by_summation ( EMData image,
int  i,
Vec3f  vec,
float  amp,
float  apix,
float  res 
)

Definition at line 2333 of file pointarray.cpp.

References EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), log(), points, sqrt(), EMAN::EMData::update(), and x.

{
        double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom

        double min_table_val = 1e-7;
        double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)

        double table_step_size = 0.001; // number of steps for each pixel
        double inv_table_step_size = 1.0 / table_step_size;
        int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
        double *table = (double *) malloc(sizeof(double) * table_size);
        for (int i = 0; i < table_size; i++) {
                double x = -i * table_step_size * apix / gauss_real_width;
                table[i] = exp(-x * x)/pow((float)M_PI,.25f);
        }
        int image_size=proj->get_xsize();

        // subtract the old point
        int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
        if (gbox <= 0)
                gbox = 1;
        float *pd = proj->get_data();
         int s = ind;
        double xc = points[4 * s] / apix + image_size / 2;
        double yc = points[4 * s + 1] / apix + image_size / 2;
        double fval = points[4 * s + 3];
        int imin = int (xc) - gbox, imax = int (xc) + gbox;
        int jmin = int (yc) - gbox, jmax = int (yc) + gbox;

        if (imin < 0) imin = 0;
        if (jmin < 0) jmin = 0;
        if (imax > image_size) imax = image_size;
        if (jmax > image_size) jmax = image_size;

        for (int j = jmin; j < jmax; j++) {
                int table_index_y = int (fabs(j - yc) * inv_table_step_size);
                double yval = table[table_index_y];
                int pd_index = j * image_size + imin;
                for (int i = imin; i < imax; i++, pd_index++) {
                        int table_index_x = int (fabs(i - xc) * inv_table_step_size);
                        double xval = table[table_index_x];
                        pd[pd_index] -= (float)(fval * yval * xval);
                }
        }

        // add the new point
        gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
        if (gbox <= 0)
                gbox = 1;
        pd = proj->get_data();
        s = ind;
        xc = vec[0] / apix + image_size / 2;
        yc = vec[1] / apix + image_size / 2;
        fval = amp;
        imin = int (xc) - gbox, imax = int (xc) + gbox;
        jmin = int (yc) - gbox, jmax = int (yc) + gbox;

        if (imin < 0) imin = 0;
        if (jmin < 0) jmin = 0;
        if (imax > image_size) imax = image_size;
        if (jmax > image_size) jmax = image_size;

        for (int j = jmin; j < jmax; j++) {
                int table_index_y = int (fabs(j - yc) * inv_table_step_size);
                double yval = table[table_index_y];
                int pd_index = j * image_size + imin;
                for (int i = imin; i < imax; i++, pd_index++) {
                        int table_index_x = int (fabs(i - xc) * inv_table_step_size);
                        double xval = table[table_index_x];
                        pd[pd_index] -= (float)(fval * yval * xval);
                }
        }

        proj->update();
        return;
}
void PointArray::reverse_chain ( )

Definition at line 2120 of file pointarray.cpp.

References n, and points.

Referenced by fit_helix().

                              {
        // reverse the point array chain, from the last to the first point
        double tmp;
//      for(int i=0; i<n/2; i++){
//              for (int j=0; j<4; j++){
//                      printf("%f\t",
        for(int i=0; i<n/2; i++){
                for (int j=0; j<4; j++){
                        tmp=points[(n-1-i)*4+j];
                        points[(n-1-i)*4+j]=points[i*4+j];
                        points[i*4+j]=tmp;
                }
        }
}
void PointArray::right_transform ( const Transform transform)

Does Transform*v as opposed to v*Transform (as in the transform function)

Parameters:
transforman EMAN2 Transform object

Definition at line 729 of file pointarray.cpp.

References n, points, EMAN::Transform::transpose(), and v.

                                                           {
        for ( unsigned int i = 0; i < 4 * n; i += 4) {
                Transform s = transform.transpose();
                Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
                v= s*v;
                points[i]  =v[0];
                points[i+1]=v[1];
                points[i+2]=v[2];
        }

}
void PointArray::save_pdb_with_helix ( const char *  file,
vector< float >  hlxid 
)

Definition at line 2075 of file pointarray.cpp.

References get_number_points(), and points.

{

        FILE *fp = fopen(file, "w");
        
        for (uint i=0; i<hlxid.size()/8; i++){
                fprintf(fp, "HELIX%5lu   A ALA A%5lu  ALA A%5lu  1                        %5lu\n", 
                                i, (int)hlxid[i*8], (int)hlxid[i*8+1], int(hlxid[i*8+1]-hlxid[i*8]+4));
        }
        for ( size_t i = 0; i < get_number_points(); i++) {
                fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
                                points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
        }
        fclose(fp);
}
void PointArray::save_to_pdb ( const char *  file)

Definition at line 554 of file pointarray.cpp.

References get_number_points(), and points.

{
        FILE *fp = fopen(file, "w");
        for ( size_t i = 0; i < get_number_points(); i++) {
                fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
                                points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
        }
        fclose(fp);
}
void PointArray::set_from ( PointArray source,
const string &  sym = "",
Transform transform = 0 
)

Definition at line 740 of file pointarray.cpp.

References get_number_points(), get_points_array(), set_from(), and transform().

{
        set_from(source->get_points_array(), source->get_number_points(), sym, transform);

}
void PointArray::set_from ( double *  source,
int  num,
const string &  sym = "",
Transform transform = 0 
)

Definition at line 746 of file pointarray.cpp.

References EMAN::Transform::get_nsym(), get_number_points(), get_points_array(), EMAN::Transform::get_sym(), set_number_points(), and v.

{
         int nsym = xform->get_nsym(sym);

        if (get_number_points() != (size_t)nsym * num)
                set_number_points((size_t)nsym * num);

        double *target = get_points_array();

        for ( int s = 0; s < nsym; s++) {
                int index = s * 4 * num;
                for ( int i = 0; i < 4 * num; i += 4, index += 4) {
                        Vec3f v((float)src[i],(float)src[i+1],(float)src[i+2]);
                        v=v*xform->get_sym(sym,s);
                        target[index]  =v[0];
                        target[index+1]=v[1];
                        target[index+2]=v[2];
                        target[index+3]=src[i+3];
                }
        }
}
void PointArray::set_from ( vector< float >  pts)

Definition at line 768 of file pointarray.cpp.

References points, and set_number_points().

Referenced by set_from().

                                           {
        set_number_points(pts.size()/4);
        for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];

}
void PointArray::set_from_density_map ( EMData map,
int  num,
float  thresh,
float  apix,
Density2PointsArrayAlgorithm  mode = PEAKS_DIV 
)

Definition at line 774 of file pointarray.cpp.

References apix, EMAN::EMData::copy(), dist(), EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::Util::get_frand(), get_number_points(), get_points_array(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), KMEANS, log(), LOGERR, nx, ny, PEAKS_DIV, PEAKS_SUB, points, set_number_points(), set_points_array(), sort_by_axis(), sqrt(), EMAN::EMData::update(), v, x, y, and zero().

{
        if (mode == PEAKS_SUB || mode == PEAKS_DIV) {
                // find out how many voxels are useful voxels
                int num_voxels = 0;
                int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
                size_t size = (size_t)nx * ny * nz;
                EMData *tmp_map = map->copy();
                float *pd = tmp_map->get_data();
                for (size_t i = 0; i < size; ++i) {
                        if (pd[i] > thresh)
                                ++num_voxels;
                }

                double pointvol = double (num_voxels) / double (num);
                double gauss_real_width = pow(pointvol, 1. / 3.);       // in pixels
#ifdef DEBUG
                printf("Average point range radius = %g pixels for %d points from %d used voxels\n",
                           gauss_real_width, num, num_voxels);
#endif

                double min_table_val = 1e-4;
                double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)

                double table_step_size = 1.;    // number of steps for each pixel
                double inv_table_step_size = 1.0 / table_step_size;
                int table_size = int (max_table_x * gauss_real_width / (table_step_size) * 1.25) + 1;
                double *table = (double *) malloc(sizeof(double) * table_size);
                for (int i = 0; i < table_size; i++) {
                        double x = i * table_step_size / gauss_real_width;
                        table[i] = exp(-x * x);
                }

                int gbox = int (max_table_x * gauss_real_width);        // local box half size in pixels to consider for each point
                if (gbox <= 0)
                        gbox = 1;

                set_number_points(num);
                for (int count = 0; count < num; count++) {
                        float cmax = pd[0];
                        int cmaxpos = 0;
                        for (size_t i = 0; i < size; ++i) {
                                if (pd[i] > cmax) {
                                        cmax = pd[i];
                                        cmaxpos = i;
                                }
                        }
                        int iz = cmaxpos / (nx * ny), iy = (cmaxpos - iz * nx * ny) / nx, ix =
                                cmaxpos - iz * nx * ny - iy * nx;

                        // update coordinates in pixels
                        points[4*count  ] = ix;
                        points[4*count+1] = iy;
                        points[4*count+2] = iz;
                        points[4*count+3] = cmax;
#ifdef DEBUG
                        printf("Point %d: val = %g\tat  %d, %d, %d\n", count, cmax, ix, iy, iz);
#endif

                        int imin = ix - gbox, imax = ix + gbox;
                        int jmin = iy - gbox, jmax = iy + gbox;
                        int kmin = iz - gbox, kmax = iz + gbox;
                        if (imin < 0)
                                imin = 0;
                        if (jmin < 0)
                                jmin = 0;
                        if (kmin < 0)
                                kmin = 0;
                        if (imax > nx)
                                imax = nx;
                        if (jmax > ny)
                                jmax = ny;
                        if (kmax > nz)
                                kmax = nz;

                        for (int k = kmin; k < kmax; ++k) {
                                int table_index_z = int (fabs(double (k - iz)) * inv_table_step_size);
                                double zval = table[table_index_z];
                                size_t pd_index_z = k * nx * ny;
                                //printf("k = %8d\tx = %8g\tval = %8g\n", k, float(k-iz), zval);
                                for (int j = jmin; j < jmax; ++j) {
                                        int table_index_y = int (fabs(double (j - iy)) * inv_table_step_size);
                                        double yval = table[table_index_y];
                                        size_t pd_index = pd_index_z + j * nx + imin;
                                        for (int i = imin; i < imax; ++i, ++pd_index) {
                                                int table_index_x = int (fabs(double (i - ix)) * inv_table_step_size);
                                                double xval = table[table_index_x];
                                                if (mode == PEAKS_SUB)
                                                        pd[pd_index] -= (float)(cmax * zval * yval * xval);
                                                else
                                                        pd[pd_index] *= (float)(1.0 - zval * yval * xval);      // mode == PEAKS_DIV
                                        }
                                }
                        }
                }
                set_number_points(num);
                tmp_map->update();
                if( tmp_map )
                {
                        delete tmp_map;
                        tmp_map = 0;
                }
        }
        else if (mode == KMEANS) {
                set_number_points(num);
                zero();

                PointArray tmp_pa;
                tmp_pa.set_number_points(num);
                tmp_pa.zero();

                 int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
                float *pd = map->get_data();

                // initialize segments with random centers at pixels with values > thresh
#ifdef DEBUG
                printf("Start initial random seeding\n");
#endif
                for ( size_t i = 0; i < get_number_points(); i++) {
                         int x, y, z;
                        double v;
                        do {
                                x = ( int) Util::get_frand(0, nx - 1);
                                y = ( int) Util::get_frand(0, ny - 1);
                                z = ( int) Util::get_frand(0, nz - 1);
                                v = pd[z * nx * ny + y * nx + x];
#ifdef DEBUG
                                printf("Trying Point %lu: val = %g\tat  %d, %d, %d\tfrom map (%d,%d,%d)\n", i, v, x,
                                           y, z, nx, ny, nz);
#endif
                        } while (v <= thresh);
                        points[4 * i] = (double) x;
                        points[4 * i + 1] = (double) y;
                        points[4 * i + 2] = (double) z;
                        points[4 * i + 3] = (double) v;
#ifdef DEBUG
                        printf("Point %lu: val = %g\tat  %g, %g, %g\n", i, points[4 * i + 3], points[4 * i],
                                   points[4 * i + 1], points[4 * i + 2]);
#endif
                }

                double min_dcen = 1e0;  // minimal mean segment center shift as convergence criterion
                double dcen = 0.0;
                int iter = 0, max_iter = 100;
                do {
#ifdef DEBUG
                        printf("Iteration %3d, start\n", iter);
#endif
                        double *tmp_points = tmp_pa.get_points_array();

                        // reassign each pixel to the best segment
                        for ( int k = 0; k < nz; k++) {
                                for ( int j = 0; j < ny; j++) {
                                        for ( int i = 0; i < nx; i++) {
                                                size_t idx = k * nx * ny + j * nx + i;
                                                if (pd[idx] > thresh) {
                                                        double min_dist = 1e60; // just a large distance
                                                         int min_s = 0;
                                                        for ( size_t s = 0; s < get_number_points(); ++s) {
                                                                double x = points[4 * s];
                                                                double y = points[4 * s + 1];
                                                                double z = points[4 * s + 2];
                                                                double dist =
                                                                        (k - z) * (k - z) + (j - y) * (j - y) + (i - x) * (i - x);
                                                                if (dist < min_dist) {
                                                                        min_dist = dist;
                                                                        min_s = s;
                                                                }
                                                        }
                                                        tmp_points[4 * min_s] += i;
                                                        tmp_points[4 * min_s + 1] += j;
                                                        tmp_points[4 * min_s + 2] += k;
                                                        tmp_points[4 * min_s + 3] += 1.0;
                                                }
                                        }
                                }
                        }
#ifdef DEBUG
                        printf("Iteration %3d, finished reassigning segments\n", iter);
#endif
                        // update each segment's center
                        dcen = 0.0;
                        for ( size_t s = 0; s < get_number_points(); ++s) {
                                if (tmp_points[4 * s + 3]) {
                                        tmp_points[4 * s] /= tmp_points[4 * s + 3];
                                        tmp_points[4 * s + 1] /= tmp_points[4 * s + 3];
                                        tmp_points[4 * s + 2] /= tmp_points[4 * s + 3];
#ifdef DEBUG
                                        printf("Iteration %3d, Point %3lu at %8g, %8g, %8g -> %8g, %8g, %8g\n", iter, s,
                                                   points[4 * s], points[4 * s + 1], points[4 * s + 2], tmp_points[4 * s],
                                                   tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
#endif
                                }
                                else {                  // empty segments are reseeded
                                         int x, y, z;
                                        double v;
                                        do {
                                                x = ( int) Util::get_frand(0, nx - 1);
                                                y = ( int) Util::get_frand(0, ny - 1);
                                                z = ( int) Util::get_frand(0, nz - 1);
                                                v = pd[z * nx * ny + y * nx + x];
                                        } while (v <= thresh);
                                        tmp_points[4 * s] = (double) x;
                                        tmp_points[4 * s + 1] = (double) y;
                                        tmp_points[4 * s + 2] = (double) z;
                                        tmp_points[4 * s + 3] = (double) v;
#ifdef DEBUG
                                        printf
                                                ("Iteration %3d, Point %3lu reseeded from %8g, %8g, %8g -> %8g, %8g, %8g\n",
                                                 iter, s, points[4 * s], points[4 * s + 1], points[4 * s + 2],
                                                 tmp_points[4 * s], tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
#endif
                                }
                                double dx = tmp_points[4 * s] - points[4 * s];
                                double dy = tmp_points[4 * s + 1] - points[4 * s + 1];
                                double dz = tmp_points[4 * s + 2] - points[4 * s + 2];
                                dcen += dx * dx + dy * dy + dz * dz;
                        }
                        dcen = sqrt(dcen / get_number_points());
                        //swap pointter, faster but risky
#ifdef DEBUG
                        printf("before swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
                                   (long)tmp_pa.get_points_array());
#endif
                        double *tp = get_points_array();
                        set_points_array(tmp_points);
                        tmp_pa.set_points_array(tp);
                        tmp_pa.zero();
#ifdef DEBUG
                        printf("after  swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
                                   (long)tmp_pa.get_points_array());
                        printf("Iteration %3d, finished updating segment centers with dcen = %g pixels\n", iter,
                                   dcen);
#endif

                        iter++;
                } while (dcen > min_dcen && iter <= max_iter);
                map->update();

                sort_by_axis(2);        // x,y,z axes = 0, 1, 2
        }
        else {
                LOGERR("PointArray::set_from_density_map(): mode = %d is not implemented yet", mode);
        }
        //update to use apix and origin
        int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
        float origx, origy, origz;
        try {
                origx = map->get_attr("origin_x");
                origy = map->get_attr("origin_y");
                origz = map->get_attr("origin_z");
        }
        catch(...) {
                origx = -nx / 2 * apix;
                origy = -ny / 2 * apix;
                origz = -nz / 2 * apix;
        }

#ifdef DEBUG
        printf("Apix = %g\torigin x,y,z = %8g,%8g,%8g\n",apix, origx, origy, origz);
#endif

        float *pd = map->get_data();
        for ( size_t i = 0; i < get_number_points(); ++i) {
#ifdef DEBUG
                printf("Point %4lu: x,y,z,v = %8g,%8g,%8g,%8g",i, points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
#endif
                points[4 * i + 3] =
                        pd[(int) points[4 * i + 2] * nx * ny + (int) points[4 * i + 1] * nx +
                           (int) points[4 * i]];
                points[4 * i] = points[4 * i] * apix + origx;
                points[4 * i + 1] = points[4 * i + 1] * apix + origy;
                points[4 * i + 2] = points[4 * i + 2] * apix + origz;
#ifdef DEBUG
                printf("\t->\t%8g,%8g,%8g,%8g\n",points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
#endif
        }
        map->update();
}
void PointArray::set_number_points ( size_t  nn)

Definition at line 178 of file pointarray.cpp.

References bfactor, n, nn(), and points.

Referenced by copy(), fit_helix(), mask(), mask_asymmetric_unit(), operator=(), read_from_pdb(), set_from(), set_from_density_map(), and sim_add_point_double().

{       
        if (n != nn) {
                n = nn;         
                points = (double *) realloc(points, 4 * n * sizeof(double));
                bfactor = (double *) realloc(bfactor, n * sizeof(double));
        }
}
void PointArray::set_points_array ( double *  p)

Definition at line 192 of file pointarray.cpp.

References points.

Referenced by set_from_density_map(), sim_add_point_double(), and sim_add_point_one().

{
        points = p;
}
void PointArray::set_vector_at ( int  i,
Vec3f  vec,
double  value 
)

Definition at line 2814 of file pointarray.cpp.

References points.

{
points[i*4]=vec[0];
points[i*4+1]=vec[1];
points[i*4+2]=vec[2];
points[i*4+3]=value;
}
void PointArray::set_vector_at ( int  i,
vector< double >  v 
)

Definition at line 2822 of file pointarray.cpp.

References points.

{
points[i*4]  =v[0];
points[i*4+1]=v[1];
points[i*4+2]=v[2];
if (v.size()>=4) points[i*4+3]=v[3];
}
void PointArray::sim_add_point_double ( )

Add more points during the simulation.

Definition at line 1440 of file pointarray.cpp.

References aang, adihed, adist, n, nn(), points, set_number_points(), set_points_array(), and sim_updategeom().

                                      {
        
        int nn=n*2;
        int tmpn=n;
        set_number_points(nn);
        double* pa2data=(double *) calloc(4 * nn, sizeof(double));
        bool *newpt=new bool[nn];
        for (int i=0;i<nn;i++){
                if (i%2==0)
                        newpt[i]=1;
                else
                        newpt[i]=0;
        }
        int i=0;
        for (int ii=0;ii<nn;ii++){
                if (newpt[ii]) {
                        
                        pa2data[ii*4]=points[i*4];
                        pa2data[ii*4+1]=points[i*4+1];
                        pa2data[ii*4+2]=points[i*4+2];
                        pa2data[ii*4+3]=1;
                        i++;
                }
                else{
                        int k;
                        if (i<tmpn)
                                k=i;
                        else
                                k=0;
                        pa2data[ii*4]=(points[k*4]+points[(i-1)*4])/2;
                        pa2data[ii*4+1]=(points[k*4+1]+points[(i-1)*4+1])/2;
                        pa2data[ii*4+2]=(points[k*4+2]+points[(i-1)*4+2])/2;
                        pa2data[ii*4+3]=1;
                }
                        
        }
                
        delete []newpt;
        free(points);
        set_points_array(pa2data);
        
        if (adist) free(adist);
        if (aang) free(aang);
        if (adihed) free(adihed);
        adist=aang=adihed=0;
        sim_updategeom();
}
void PointArray::sim_add_point_one ( )

Definition at line 1490 of file pointarray.cpp.

References aang, adihed, adist, apix, EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), map, mapc, n, points, set_points_array(), EMAN::EMData::sget_value_at_interp(), sim_pointpotential(), sim_potentiald(), and sim_updategeom().

                                   {

        
        double maxpot=-1000000,pot,meanpot=0;
        int ipt=-1;
        bool onedge=0;
        // Find the highest potential point
        for (int i=0; i<n; i++) {
                meanpot+=sim_pointpotential(adist[i],aang[i],adihed[i]);
                pot=/*sim_pointpotential(adist[i],aang[i],adihed[i])*/-mapc*map->sget_value_at_interp(points[i*4]/apix+map->get_xsize()/2,points[i*4+1]/apix+map->get_ysize()/2,points[i*4+2]/apix+map->get_zsize()/2);
                if (pot>maxpot){
                        maxpot=pot;
                        ipt=i;
                }
        }
        meanpot/=n;
        
        for (int i=0; i<n; i++) {
                int k=(i+1)%n;
                double pt0,pt1,pt2;
                pt0=(points[k*4]+points[i*4])/2;
                pt1=(points[k*4+1]+points[i*4+1])/2;
                pt2=(points[k*4+2]+points[i*4+2])/2;
                pot=/*meanpot*/-mapc*map->sget_value_at_interp(pt0/apix+map->get_xsize()/2,pt1/apix+map->get_ysize()/2,pt2/apix+map->get_zsize()/2);
                if (pot>maxpot){
                        maxpot=pot;
                        ipt=i;
                        onedge=1;
                }
        
        }
        
        // The rest points remain the same
        int i;
        double* pa2data=(double *) calloc(4 * (n+1), sizeof(double));
        for (int ii=0; ii<n+1; ii++) {
                if(ii!=ipt && ii!=ipt+1){
                        if(ii<ipt)
                                i=ii;
                        else    // shift the points after the adding position
                                i=ii-1;
                        
                        pa2data[ii*4]=points[i*4];
                        pa2data[ii*4+1]=points[i*4+1];
                        pa2data[ii*4+2]=points[i*4+2];
                        pa2data[ii*4+3]=1;
                }
        }
        // Adding points
        if( onedge ) {
                int k0,k1;
                k0=((ipt+n-1)%n);
                k1=((ipt+1)%n);
                pa2data[ipt*4]=points[ipt*4];
                pa2data[ipt*4+1]=points[ipt*4+1];
                pa2data[ipt*4+2]=points[ipt*4+2];
                pa2data[ipt*4+3]=1;
                
                pa2data[(ipt+1)*4]=(points[ipt*4]+points[k1*4])/2;
                pa2data[(ipt+1)*4+1]=(points[ipt*4+1]+points[k1*4+1])/2;
                pa2data[(ipt+1)*4+2]=(points[ipt*4+2]+points[k1*4+2])/2;
                pa2data[(ipt+1)*4+3]=1; 
                
        }
        else {
                int k0,k1;
                k0=((ipt+n-1)%n);
                k1=((ipt+1)%n);
                pa2data[ipt*4]=(points[ipt*4]+points[k0*4])/2;
                pa2data[ipt*4+1]=(points[ipt*4+1]+points[k0*4+1])/2;
                pa2data[ipt*4+2]=(points[ipt*4+2]+points[k0*4+2])/2;
                pa2data[ipt*4+3]=1;
                
                pa2data[(ipt+1)*4]=(points[ipt*4]+points[k1*4])/2;
                pa2data[(ipt+1)*4+1]=(points[ipt*4+1]+points[k1*4+1])/2;
                pa2data[(ipt+1)*4+2]=(points[ipt*4+2]+points[k1*4+2])/2;
                pa2data[(ipt+1)*4+3]=1; 
        
        }
        free(points);
        n++;
        set_points_array(pa2data);
        
        if (adist) free(adist);
        if (aang) free(aang);
        if (adihed) free(adihed);
        adist=aang=adihed=0;
        sim_updategeom();
        
        // search for the best position for the new points
        if (onedge){
                i=ipt+1;
                double bestpot=10000,nowpot;
                Vec3f old(points[i*4],points[(i+1)*4],points[(i+2)*4]);
                Vec3f newpt(points[i*4],points[(i+1)*4],points[(i+2)*4]);
                for (int ii=0;ii<5000;ii++){
                        // Try to minimize potential globally
                        boost::mt19937 rng; 
                        boost::normal_distribution<> nd(0.0, 0.0);
                        boost::variate_generator<boost::mt19937&,boost::normal_distribution<> > var_nor(rng, nd);
                        points[i*4]=old[0]+var_nor();
                        points[i*4+1]=old[1]+var_nor();
                        points[i*4+2]=old[2]+var_nor();
                        nowpot=sim_potentiald(i);
                        if (nowpot<bestpot) {
                                bestpot=nowpot;
                                newpt[0]=points[i*4];
                                newpt[1]=points[(i+1)*4];
                                newpt[2]=points[(i+2)*4];
                        }
                                
                }
                points[i*4]=newpt[0];
                points[i*4+1]=newpt[1];
                points[i*4+2]=newpt[2];
        }
        
}
Vec3f PointArray::sim_descent ( int  i)

returns a vector pointing downhill for a single point

Definition at line 1222 of file pointarray.cpp.

References points, and sim_potentiald().

Referenced by sim_minstep(), and sim_minstep_seq().

                                   {
        Vec3f old(points[i*4],points[i*4+1],points[i*4+2]);
        double pot=0.,potx=0.,poty=0.,potz=0.;
        double stepsz=0.01;
        
        for (int ii=i-2; ii<=i+2; ii++) pot+=sim_potentiald(ii);
//      pot=potential();
        
        points[i*4]=old[0]+stepsz;
        for (int ii=i-2; ii<=i+2; ii++) potx+=sim_potentiald(ii);
//      potx=potential();
        points[i*4]=old[0];

        points[i*4+1]=old[1]+stepsz;
        for (int ii=i-2; ii<=i+2; ii++) poty+=sim_potentiald(ii);
//      poty=potential();
        points[i*4+1]=old[1];

        points[i*4+2]=old[2]+stepsz;
        for (int ii=i-2; ii<=i+2; ii++) potz+=sim_potentiald(ii);
//      potz=potential();
        points[i*4+2]=old[2];
        
        //      printf("%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\n",pot,potx,poty,potz,(pot-potx),(pot-poty),(pot-potz));
        
//      if (pot==potx) potx=pot+1000000.0;
//      if (pot==potx) potx=pot+1000000.0;
//      if (pot==potx) potx=pot+1000000.0;
        return Vec3f((pot-potx),(pot-poty),(pot-potz));
}
void PointArray::sim_minstep ( double  maxshift)

Takes a step to minimize the potential.

Definition at line 1254 of file pointarray.cpp.

References max, mean(), n, oldshifts, points, shifts, and sim_descent().

                                            { 
        vector<Vec3f> shifts;
        
        double max=0.0;
        double mean=0.0;
        for (uint i=0; i<n; i++) {
                if (oldshifts.size()==n) shifts.push_back((sim_descent(i)+oldshifts[i])/2.0);
                else shifts.push_back(sim_descent(i));
                float len=shifts[i].length();
                if (len>max) max=len;
                mean+=len;
        }
        oldshifts=shifts;
        
//      printf("max vec %1.2f\tmean %1.3f\n",max,mean/n);
        
        for (uint i=0; i<n; i++) {
//              if (std::isnan(shifts[i][0]) ||std::isnan(shifts[i][1]) ||std::isnan(shifts[i][2])) { printf("Nan: %d\n",i); shifts[i]=Vec3f(max,max,max); }
                points[i*4]+=shifts[i][0]*maxshift/max;
                points[i*4+1]+=shifts[i][1]*maxshift/max;
                points[i*4+2]+=shifts[i][2]*maxshift/max;
//              printf("%d. %1.2f\t%1.2f\t%1.2f\n",i,shifts[i][0]*maxshift/max,shifts[i][1]*maxshift/max,shifts[i][2]*maxshift/max);
        }
}
void PointArray::sim_minstep_seq ( double  meanshift)

Takes a step to minimize the potential.

Definition at line 1280 of file pointarray.cpp.

References EMAN::Util::get_irand(), EMAN::Vec3< Type >::length(), mean(), n, points, and sim_descent().

                                                 {
        /*
        // Try to minimize potential globally
        boost::mt19937 rng; 
        boost::normal_distribution<> nd(0.0, 20.0);
        boost::variate_generator<boost::mt19937&,boost::normal_distribution<> > var_nor(rng, nd);
        double *oldpts=new double[4*get_number_points()];
        double *bestpts=new double[4*get_number_points()];
        double best_pot,new_pot;
        memcpy(oldpts, get_points_array(), sizeof(double) * 4 * get_number_points());
        memcpy(bestpts, get_points_array(), sizeof(double) * 4 * get_number_points());
        best_pot=sim_potential();
        double disttmp=0;
        for (int k=0; k<n; k++) disttmp+=adist[k];
        best_pot+=distc*pow((disttmp-336*3.3),2.0);
        for (int i=0; i<1000; i++){
                for (int j=0; j<n; j++){
                        points[4*j]=oldpts[4*j]+var_nor();
                        points[4*j+1]=oldpts[4*j+1]+var_nor();
                        points[4*j+2]=oldpts[4*j+2]+var_nor();
                }
                new_pot=sim_potential();
                disttmp=0;
                for (int k=0; k<n; k++) disttmp+=adist[k];
                new_pot+=distc*pow((disttmp-336*3.3),2.0);
                if (new_pot<best_pot){
                        memcpy(bestpts, get_points_array(), sizeof(double) * 4 * get_number_points());
                        best_pot=new_pot;
                        printf("%f\t",best_pot);
                }
        }
        memcpy(get_points_array(),bestpts, sizeof(double) * 4 * get_number_points());
        
        delete []oldpts;
        delete []bestpts;
        */
        // we compute 10 random gradients and use these to adjust stepsize
        double mean=0.0;
        for (int i=0; i<10; i++) {
//              Vec3f shift=sim_descent(random()%n);
                Vec3f shift=sim_descent(Util::get_irand(0,n-1));
                mean+=shift.length();
        }
        mean/=10.0;
        double stepadj=meanshift/mean;
//      printf("\t%1.4g\n",stepadj);

        // Now we go through all points sequentially and move each downhill
        // The trick here is the sequential part, as each point is impacted by the point already moved before it.
        // This may create a "seam" at the first point which won't be adjusted to compensate for the last point (wraparound)
        // until the next cycle
        Vec3f oshifts;
        for (uint ii=0; ii<n; ii++) {
                uint i=2*(ii%(n/2))+2*ii/n;     // this maps a linear sequence to an all-even -> all odd sequence
                Vec3f shift,d;
                if (ii==n) {
                        d=sim_descent(i);
                        shift=(d+oshifts)/2.0;
                        oshifts=d;
                }
                else {
                        shift=sim_descent(i);
                        oshifts=shift;
                }
                
//              double p2=potential();
//              double pot=sim_potentialdxyz(i,0.0,0.0,0.0);
//              double pots=sim_potentialdxyz(i,shift[0]*stepadj,shift[1]*stepadj,shift[2]*stepadj);

                // only step if it actually improves the potential for this particle (note that it does not need to improve the overall potential)
//              if (pots<pot) {
                        points[i*4]+=shift[0]*stepadj;
                        points[i*4+1]+=shift[1]*stepadj;
                        points[i*4+2]+=shift[2]*stepadj;
//                      printf("%d. %1.4g -> %1.4g  %1.3g %1.3g %1.3g %1.3g\n",i,pot,pots,shift[0],shift[1],shift[2],stepadj);
//                      if (potential()>p2) printf("%d. %1.4g %1.4g\t%1.4g %1.4g\n",i,pot,pots,p2,potential());
//              }
                
        }       
}
double EMAN::PointArray::sim_pointpotential ( double  dist,
double  ang,
double  dihed 
) [inline]

Computes a potential value for a single point from a set of angles/distances using current energy settings.

Definition at line 156 of file pointarray.h.

References angc, dihed0, dihedc, dist0, and distc.

Referenced by sim_add_point_one(), sim_potential(), and sim_potentiald().

                                                                                        {
                        return pow(dist-dist0,2.0)*distc+ang*ang*angc+pow(dihed-dihed0,2.0)*dihedc;
                }
double PointArray::sim_potential ( )

Computes overall potential for the configuration.

Definition at line 1099 of file pointarray.cpp.

References aang, adihed, adist, apix, EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), map, mapc, n, points, EMAN::EMData::sget_value_at_interp(), sim_pointpotential(), and sim_updategeom().

Referenced by sim_printstat().

                                 {
        double ret=0;
        sim_updategeom();
        
        if (map &&mapc) {
                for (uint i=0; i<n; i++) ret+=sim_pointpotential(adist[i],aang[i],adihed[i])-mapc*map->sget_value_at_interp(points[i*4]/apix+map->get_xsize()/2,points[i*4+1]/apix+map->get_ysize()/2,points[i*4+2]/apix+map->get_zsize()/2);
        }
        else {
                for (uint i=0; i<n; i++) ret+=sim_pointpotential(adist[i],aang[i],adihed[i]);
        }

#ifdef _WIN32
//      if (_isnan(ret/n))
#else
//      if (std::isnan(ret/n))
#endif
//              printf("%f             %f\n",ret,n);
        
        return ret/n;
}
double PointArray::sim_potentiald ( int  i)

Compute a single point potential value.

Definition at line 1121 of file pointarray.cpp.

References aang, adihed, adist, apix, b, EMAN::Vec3< Type >::cross(), dist(), distpenc, EMAN::Vec3< Type >::dot(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Vec3< Type >::length(), map, mapc, mindistc, n, points, EMAN::EMData::sget_value_at_interp(), sim_pointpotential(), and sim_updategeom().

Referenced by sim_add_point_one(), sim_descent(), and sim_potentialdxyz().

                                       {
        if (!adist) sim_updategeom();           // wasteful, but only once
        
//      if (i<0 || i>=n) throw InvalidParameterException("Point number out of range");
        if (i<0) i-=n*(i/n-1);
        if (i>=n) i=i%n;
        
        // how expensive is % ?  Worth replacing ?
        int ib=4*((i+n-1)%n);           // point before i with wraparound
        int ibb=4*((i+n-2)%n);  // 2 points before i with wraparound
        int ia=4*((i+1)%n);             // 1 point after
        int ii=i*4;
        
        Vec3f a(points[ib]-points[ibb],points[ib+1]-points[ibb+1],points[ib+2]-points[ibb+2]);                  // -2 -> -1
        Vec3f b(points[ii]-points[ib],points[ii+1]-points[ib+1],points[ii+2]-points[ib+2]);             // -1 -> 0
        Vec3f c(points[ia]-points[ii],points[ia+1]-points[ii+1],points[ia+2]-points[ii+2]);             // 0 -> 1
        double dist=b.length();
        adist[i]=dist;
        // Angle, tests should avoid isnan being necessary
        double ang=b.dot(c);
        if (ang!=0.0) {                                 // if b.dot(c) is 0, we set it to the last determined value...
                ang/=(dist*c.length());
                if (ang>1.0) ang=1.0;           // should never happen, but just in case of roundoff error
                if (ang<-1.0) ang=-1.0;
                ang=acos(ang);
        }
        else ang=aang[i];
        aang[i]=ang;
        
        // Dihedral
        Vec3f cr1=a.cross(b);
        Vec3f cr2=b.cross(c);
        double dihed;
        double denom=cr1.length()*cr2.length();
        if (denom==0) dihed=adihed[i];                                          // set the dihedral to the last determined value if indeterminate
        else {
                dihed=cr1.dot(cr2)/denom;
                if (dihed>1.0) dihed=1.0;                               // should never happen, but just in case of roundoff error
                if (dihed<-1.0) dihed=-1.0;
                dihed=acos(dihed);
        }
        adihed[i]=dihed;
//*     Do not need for small amount of points
        // Distance to the closest neighbor
        double mindist=10000;
        for (int j=0;j<n;j++){
                if(j==i)
                        continue;
                int ja=4*j;
                Vec3f d(points[ii]-points[ja],points[ii+1]-points[ja+1],points[ii+2]-points[ja+2]);
                double jdst=d.length();
                if(jdst<mindist)
                        mindist=jdst;
        }
        double distpen=0;
        if (mindist<mindistc)
                distpen=distpenc/mindist;
//*/    
        
        
//      if (std::isnan(dist) || std::isnan(ang) || std::isnan(dihed)) printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\n",i,dist,ang,dihed,b.length(),c.length(),b.dot(c)/(dist*c.length()));
//      if (std::isnan(dihed)) dihed=dihed0;
//      if (std::isnan(ang)) ang=0;
//      if (std::isnan(dist)) dist=3.3;
//      if (isnan(dihed)) dihed=dihed0;
//      if (isnan(ang)) ang=0;
        if (map && mapc) {
                return distpen+sim_pointpotential(dist,ang,dihed)-mapc*map->sget_value_at_interp(points[ii]/apix+map->get_xsize()/2,points[ii+1]/apix+map->get_ysize()/2,points[ii+2]/apix+map->get_zsize()/2);
        }
        return sim_pointpotential(dist,ang,dihed);
}
double PointArray::sim_potentialdxyz ( int  i,
double  dx,
double  dy,
double  dz 
)

Compute a potential value for a perturbed point, including +-2 nearest neighbors which will also be impacted.

Definition at line 1194 of file pointarray.cpp.

References points, and sim_potentiald().

                                                                           {
        Vec3f old(points[i*4],points[i*4+1],points[i*4+2]);
        double potd=0.;
                
        points[i*4]=old[0]+dx;
        points[i*4+1]=old[1]+dy;
        points[i*4+2]=old[2]+dz;
        for (int ii=i-2; ii<=i+2; ii++) potd+=sim_potentiald(ii);
//      potd=potential();
        points[i*4]=old[0];
        points[i*4+1]=old[1];
        points[i*4+2]=old[2];
        
        return potd;
}
void PointArray::sim_printstat ( )

prints some statistics to the screen

Definition at line 1388 of file pointarray.cpp.

References aang, adihed, adist, n, sim_potential(), and sim_updategeom().

                               {
        sim_updategeom();       
        
        double mdist=0.0,mang=0.0,mdihed=0.0;
        double midist=1000.0,miang=M_PI*2,midihed=M_PI*2;
        double madist=0.0,maang=0.0,madihed=0.0;
        
        for (int i=0; i<n; i++) {
                mdist+=adist[i];
                mang+=aang[i];
                mdihed+=adihed[i];

                midist=adist[i]<midist?adist[i]:midist;
                madist=adist[i]>madist?adist[i]:madist;
                
                miang=aang[i]<miang?aang[i]:miang;
                maang=aang[i]>maang?aang[i]:maang;
                
                midihed=adihed[i]<midihed?adihed[i]:midihed;
                madihed=adihed[i]>madihed?adihed[i]:madihed;
        }
        double p=sim_potential();
        double anorm = 180.0/M_PI;
        printf(" potential: %1.1f\t dist: %1.2f / %1.2f / %1.2f\tang: %1.2f / %1.2f / %1.2f\tdihed: %1.2f / %1.2f / %1.2f  ln=%1.1f\n",p,midist,mdist/n,madist,miang*anorm,mang/n*anorm,maang*anorm,midihed*anorm,mdihed/n*anorm,madihed*anorm,mdihed/(M_PI*2.0)-n/10.0);
                
}
void PointArray::sim_rescale ( )

rescale the entire set so the mean bond length matches dist0

Definition at line 1362 of file pointarray.cpp.

References b, dist0, EMAN::Vec3< Type >::length(), max, min, n, and points.

                             {
        double meandist=0.0;
        double max=0.0,min=dist0*1000.0;
        
        for (uint ii=0; ii<n; ii++) {
                int ib=4*((ii+n-1)%n);          // point before i with wraparound
                int i=4*ii;

                Vec3f b(points[i]-points[ib],points[i+1]-points[ib+1],points[i+2]-points[ib+2]);                // -1 -> 0
                double len=b.length();
                meandist+=len;
                max=len>max?len:max;
                min=len<min?len:min;
        }
        meandist/=n;
        double scale=dist0/meandist;
        
        printf("mean = %1.3f rescaled: %1.3f - %1.3f\n",meandist,min*scale,max*scale);
        
        for (uint i=0; i<n; i++) {
                points[i*4]*=scale;
                points[i*4+1]*=scale;
                points[i*4+2]*=scale;
        }

}       
void PointArray::sim_set_pot_parms ( double  pdist0,
double  pdistc,
double  pangc,
double  pdihed0,
double  pdihedc,
double  pmapc,
EMData pmap,
double  pmindistc,
double  pdistpenc 
)

Sets the parameters for the energy function.

Definition at line 1415 of file pointarray.cpp.

References angc, apix, dihed0, dihedc, dist0, distc, distpenc, EMAN::EMData::get_attr(), gradx, grady, gradz, map, mapc, and mindistc.

                                                                                                                                                                          {
        dist0=pdist0;
        distc=pdistc;
        angc=pangc;
        dihed0=pdihed0;
        dihedc=pdihedc;
        mapc=pmapc;
        mindistc=pmindistc;
        distpenc=pdistpenc;
        if (pmap!=0 && pmap!=map) {
//              if (map!=0) delete map;
                if (gradx!=0) delete gradx;
                if (grady!=0) delete grady;
                if (gradz!=0) delete gradz;
                
                map=pmap;
                apix=map->get_attr("apix_x");
//              gradx=map->process("math.edge.xgradient");  // we compute the gradient to make the minimization easier
//              grady=map->process("math.edge.ygradient");
//              gradz=map->process("math.edge.zgradient");
        }
                
}
void PointArray::sim_updategeom ( )

Updates the dist,ang,dihed parameters.

Definition at line 1056 of file pointarray.cpp.

References aang, adihed, adist, b, EMAN::Vec3< Type >::cross(), EMAN::Vec3< Type >::dot(), EMAN::Vec3< Type >::length(), n, and points.

Referenced by sim_add_point_double(), sim_add_point_one(), sim_potential(), sim_potentiald(), and sim_printstat().

                                {
        if (!adist) adist=(double *)malloc(sizeof(double)*n);
        if (!aang) aang=(double *)malloc(sizeof(double)*n);
        if (!adihed) adihed=(double *)malloc(sizeof(double)*n);
        
        for (uint ii=0; ii<n; ii++) {
                // how expensive is % ?  Worth replacing ?
                int ib=4*((ii+n-1)%n);          // point before i with wraparound
                int ibb=4*((ii+n-2)%n); // 2 points before i with wraparound
                int ia=4*((ii+1)%n);            // 1 point after
                int i=4*ii;

                Vec3f a(points[ib]-points[ibb],points[ib+1]-points[ibb+1],points[ib+2]-points[ibb+2]);  // -2 -> -1
                Vec3f b(points[i]-points[ib],points[i+1]-points[ib+1],points[i+2]-points[ib+2]);                // -1 -> 0
                Vec3f c(points[ia]-points[i],points[ia+1]-points[i+1],points[ia+2]-points[i+2]);                // 0 -> 1
                adist[ii]=b.length();
                
                // angle
                aang[ii]=b.dot(c);
                if (aang[ii]!=0) {
                        aang[ii]/=(adist[ii]*c.length());
                        if (aang[ii]>=1.0) aang[ii]=0.0;
                        else if (aang[ii]<=-1.0) aang[ii]=M_PI;
                        else aang[ii]=acos(aang[ii]);
                }

                // dihedral
                Vec3f cr1=a.cross(b);
                Vec3f cr2=b.cross(c);
                double denom=cr1.length()*cr2.length();
                if (denom==0) adihed[ii]=0;
                else {
                        double tmp=cr1.dot(cr2)/(denom);
                        if (tmp>1) tmp=1;
                        if (tmp<-1) tmp=-1;
                        adihed[ii]=acos(tmp);
                }
                
//              if (std::isnan(ang[ii])) ang[ii]=0;
                
        }
}
void PointArray::sort_by_axis ( int  axis = 1)

Definition at line 2135 of file pointarray.cpp.

References cmp_axis_x(), cmp_axis_y(), cmp_axis_z(), cmp_val(), n, and points.

Referenced by set_from_density_map().

{
        if (axis == 0)
                qsort(points, n, sizeof(double) * 4, cmp_axis_x);
        else if (axis == 1)
                qsort(points, n, sizeof(double) * 4, cmp_axis_y);
        else if (axis == 2)
                qsort(points, n, sizeof(double) * 4, cmp_axis_z);
        else
                qsort(points, n, sizeof(double) * 4, cmp_val);
}
void PointArray::transform ( const Transform transform)

Definition at line 717 of file pointarray.cpp.

References n, points, and v.

Referenced by set_from().

                                              {

        for ( unsigned int i = 0; i < 4 * n; i += 4) {
                Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
                v=v*xf;
                points[i]  =v[0];
                points[i+1]=v[1];
                points[i+2]=v[2];
        }

}
void PointArray::zero ( )

Definition at line 149 of file pointarray.cpp.

References n, and points.

Referenced by set_from_density_map().

{
        memset((void *) points, 0, 4 * n * sizeof(double));
}

Member Data Documentation

double* EMAN::PointArray::aang [private]
double* EMAN::PointArray::adihed [private]
double* EMAN::PointArray::adist [private]
double EMAN::PointArray::angc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_pointpotential(), and sim_set_pot_parms().

double EMAN::PointArray::apix [private]
double* EMAN::PointArray::bfactor [private]

Definition at line 214 of file pointarray.h.

Referenced by pdb2mrc_by_summation(), PointArray(), read_from_pdb(), and set_number_points().

double EMAN::PointArray::dihed0 [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_pointpotential(), and sim_set_pot_parms().

double EMAN::PointArray::dihedc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_pointpotential(), and sim_set_pot_parms().

double EMAN::PointArray::dist0 [private]

Used for simplistic loop dynamics simulation Assumes all points are connected sequetially in a closed loop, potential includes distance, angle and dihedral terms, with optional denisty based terms.

Parameters:
dist0- center distance for quadratic energy term
distc- quadratic distance coefficient c(dist-dist0)^2
angc- quadratic angle coefficient c*(180-ang)^2
dihed0- dihedral center angle
dihedc- dihedral angle coefficient c*(dihed-dihed0)^2
mapc- coefficient for map energy
map- EMData representing map to match/fit
mindistc- mininum distance between two points
distpenc- penalty for close points

Definition at line 235 of file pointarray.h.

Referenced by sim_pointpotential(), sim_rescale(), and sim_set_pot_parms().

double EMAN::PointArray::distc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_pointpotential(), and sim_set_pot_parms().

double EMAN::PointArray::distpenc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_potentiald(), and sim_set_pot_parms().

Definition at line 237 of file pointarray.h.

Referenced by PointArray(), sim_set_pot_parms(), and ~PointArray().

Definition at line 237 of file pointarray.h.

Referenced by PointArray(), sim_set_pot_parms(), and ~PointArray().

Definition at line 237 of file pointarray.h.

Referenced by PointArray(), sim_set_pot_parms(), and ~PointArray().

double EMAN::PointArray::mapc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_add_point_one(), sim_potential(), sim_potentiald(), and sim_set_pot_parms().

double EMAN::PointArray::mindistc [private]

Definition at line 235 of file pointarray.h.

Referenced by sim_potentiald(), and sim_set_pot_parms().

size_t EMAN::PointArray::n [private]

Definition at line 238 of file pointarray.h.

Referenced by sim_minstep().

double* EMAN::PointArray::points [private]

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