EMAN2
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:
[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)
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>
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 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 ()
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.
Add more points during the simulation.
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 * aang
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
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.

```                {
PEAKS_SUB, PEAKS_DIV, KMEANS
};
```

## Constructor & Destructor Documentation

 PointArray::PointArray ( )

Definition at line 108 of file pointarray.cpp.

Referenced by copy().

```{
points = 0;
bfactor = 0;
n = 0;

aang=0;

}
```
 PointArray::PointArray ( int nn ) ` [explicit]`

Definition at line 121 of file pointarray.cpp.

```{
n = nn;
points = (double *) calloc(4 * n, sizeof(double));

aang=0;
}
```
 PointArray::~PointArray ( )

Definition at line 132 of file pointarray.cpp.

```{
if( points )
{
free(points);
points = 0;
}

if (aang) free(aang);
//      if (map!=0) delete map;
}
```

## Member Function Documentation

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

Aligns one PointArray to another in 2 dimensions.

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

Definition at line 286 of file pointarray.cpp.

```                                                             {
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:
 to Another 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.

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().

```{
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:
 sortrows if 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.

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 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

Definition at line 173 of file pointarray.cpp.

References n.

```{
return n;
}
```
 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.

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.

Parameters:
 proj A vector of EMData objects containing projections with orientations pixres Size 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.

```{
#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.

```{
#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]);

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;
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.

```{
#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.

```{
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.

```{
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.

```                                                                    {

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.

```{
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);
}
}

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:
 transform an 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.

```{
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.

```{
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.

```{
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.

```{
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];
}
```

Add more points during the simulation.

Definition at line 1440 of file pointarray.cpp.

```                                      {

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 (aang) free(aang);
sim_updategeom();
}
```

Definition at line 1490 of file pointarray.cpp.

```                                   {

double maxpot=-1000000,pot,meanpot=0;
int ipt=-1;
bool onedge=0;
// Find the highest potential point
for (int i=0; i<n; i++) {
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;
}
}
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 (aang) free(aang);
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;

// 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);

// 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) {
//                      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.

Referenced by sim_printstat().

```                                 {
double ret=0;
sim_updategeom();

if (map &&mapc) {
}
else {
}

#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.

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();
// 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);
}
//*     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.

```                               {
sim_updategeom();

double mdist=0.0,mang=0.0,mdihed=0.0;
double midist=1000.0,miang=M_PI*2,midihed=M_PI*2;

for (int i=0; i<n; i++) {
mang+=aang[i];

miang=aang[i]<miang?aang[i]:miang;
maang=aang[i]>maang?aang[i]:maang;

}
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.

```                                                                                                                                                                          {
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;

map=pmap;
apix=map->get_attr("apix_x");
}

}
```
 void PointArray::sim_updategeom ( )

Definition at line 1056 of file pointarray.cpp.

```                                {
if (!aang) aang=(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

// angle
aang[ii]=b.dot(c);
if (aang[ii]!=0) {
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();
else {
double tmp=cr1.dot(cr2)/(denom);
if (tmp>1) tmp=1;
if (tmp<-1) tmp=-1;
}

//              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]`

Definition at line 219 of file pointarray.h.

 double* EMAN::PointArray::adihed` [private]`

Definition at line 220 of file pointarray.h.

 double* EMAN::PointArray::adist` [private]`

Definition at line 218 of file pointarray.h.

 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]`

Definition at line 235 of file pointarray.h.

 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().

 EMData* EMAN::PointArray::gradx` [private]`

Definition at line 237 of file pointarray.h.

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

 EMData * EMAN::PointArray::grady` [private]`

Definition at line 237 of file pointarray.h.

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

 EMData * EMAN::PointArray::gradz` [private]`

Definition at line 237 of file pointarray.h.

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

 EMData* EMAN::PointArray::map` [private]`

Definition at line 236 of file pointarray.h.

 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 213 of file pointarray.h.

 vector EMAN::PointArray::oldshifts` [private]`

Definition at line 238 of file pointarray.h.

Referenced by sim_minstep().

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