Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

EMAN::PointArray Class Reference

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

#include <pointarray.h>

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).
Transform3Dalign_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 Transform3D &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="", Transform3D *transform=0)
void set_from (double *source, int num, const string &sym="", Transform3D *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)
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.

Private Attributes

double * points
size_t n


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

enum EMAN::PointArray::Density2PointsArrayAlgorithm
 

Enumerator:
PEAKS_SUB 
PEAKS_DIV 
KMEANS 

Definition at line 58 of file pointarray.h.

00059                 {
00060                         PEAKS_SUB, PEAKS_DIV, KMEANS
00061                 };


Constructor & Destructor Documentation

PointArray::PointArray  ) 
 

Definition at line 95 of file pointarray.cpp.

References n, and points.

Referenced by copy().

00096 {
00097         points = 0;
00098         n = 0;
00099 }

PointArray::PointArray int  nn  )  [explicit]
 

Definition at line 101 of file pointarray.cpp.

References n, and points.

00102 {
00103         n = nn;
00104         points = (double *) calloc(4 * n, sizeof(double));
00105 }

PointArray::~PointArray  ) 
 

Definition at line 107 of file pointarray.cpp.

References points.

00108 {
00109         if( points )
00110         {
00111                 free(points);
00112                 points = 0;
00113         }
00114 }


Member Function Documentation

Transform3D * 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 Transform3D to map 'this' to 'to'

Definition at line 252 of file pointarray.cpp.

References EMAN::Util::calc_bilinear_least_square(), get_vector_at(), match_points(), EMAN::Transform3D::set_pretrans(), and EMAN::Transform3D::set_rotation().

00252                                                                {
00253 vector<int> match=match_points(to,max_dist);
00254 Transform3D *ret=new Transform3D();
00255 
00256 // we use bilinear least squares to get 3/6 matrix components
00257 unsigned int i,j;
00258 
00259 vector<float> pts;
00260 for (i=0; i<match.size(); i++) {
00261         if (match[i]==-1) continue;
00262 
00263 //      printf("%d -> %d\n",i,match[i]);
00264         pts.push_back(get_vector_at(i)[0]);
00265         pts.push_back(get_vector_at(i)[1]);
00266         pts.push_back(to->get_vector_at(match[i])[0]);
00267 }
00268 
00269 Vec3f vx=Util::calc_bilinear_least_square(pts);
00270 
00271 // then we get the other 3/6
00272 for (i=j=0; i<match.size(); i++) {
00273         if (match[i]==-1) continue;
00274         pts[j*3]  =get_vector_at(i)[0];
00275         pts[j*3+1]=get_vector_at(i)[1];
00276         pts[j*3+2]=to->get_vector_at(match[i])[1];
00277         j++;
00278 }
00279 
00280 Vec3f vy=Util::calc_bilinear_least_square(pts);
00281 
00282 //ret->set_rotation(vx[1],vx[2],0.0,vy[1],vy[2],0.0,0.0,0.0,1.0);
00283 ret->set_rotation(vx[1],vy[1],0.0,vx[2],vy[2],0.0,0.0,0.0,1.0);
00284 ret->set_pretrans(Vec3f(-vx[0],-vy[0],0));
00285 
00286 return ret;
00287 }

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 289 of file pointarray.cpp.

References get_number_points(), get_value_at(), get_vector_at(), and EMAN::Vec3< Type >::length().

00289                                                                                              {
00290 printf("Warning, this is old code. Use align_2d.\n");
00291 
00292 // returns (dx,dy,residual error,n points used)
00293 // dxhint,dyhint should translate this->to
00294 // flags : 1 - use hint values, 2 - center by strongest point (highest 'value')
00295 int na=   get_number_points();
00296 int nb=to->get_number_points();
00297 if (na<=0 || nb<=0) return vector<float>(4,0);
00298 
00299 int *a2b = (int *)malloc(na*sizeof(int));
00300 int *b2a = (int *)malloc(nb*sizeof(int));
00301 
00302 
00303 // find unweighted centers
00304 float cax,cay,cbx,cby;
00305 int i,j;
00306 
00307 if (flags&1) {
00308         cbx=dxhint;
00309         cby=dyhint;
00310         cax=cay=0;
00311 }
00312 else if (flags&2) {
00313         // find the 'a' list peak
00314         float hia=0.0f;
00315         int hina=0;
00316         for (i=0; i<na; i++) {
00317                 if (get_value_at(i)>hia) { hia=static_cast<float>(get_value_at(i)); hina=i; }
00318         }
00319         cax=get_vector_at(hina)[0];
00320         cay=get_vector_at(hina)[1];
00321 
00322         // find the 'b' list peak
00323         float hib=0;
00324         int hinb=0;
00325         for (i=0; i<na; i++) {
00326                 if (to->get_value_at(i)>hib) { hib=static_cast<float>(to->get_value_at(i)); hinb=i; }
00327         }
00328         cbx=to->get_vector_at(hinb)[0];
00329         cby=to->get_vector_at(hinb)[1];
00330 }
00331 else {
00332         cax=cay=cbx=cby=0;
00333 
00334         for (i=0; i<na; i++) { cax+=get_vector_at(i)[0]; cay+=get_vector_at(i)[1]; }
00335         cax/=(float)na;
00336         cay/=(float)na;
00337 
00338         for (i=0; i<nb; i++) { cbx+=to->get_vector_at(i)[0]; cby+=to->get_vector_at(i)[1]; }
00339         cbx/=(float)nb;
00340         cby/=(float)nb;
00341 }
00342 
00343 Vec3f offset(cbx-cax,cby-cay,0);
00344 
00345 // find the nearest point for each x point, taking the estimated centers into account
00346 for (i=0; i<na; i++) {
00347         float rmin=1.0e30f;
00348         for (j=0; j<nb; j++) {
00349                 float r=(get_vector_at(i)+offset-to->get_vector_at(j)).length();
00350                 if (r<rmin) { a2b[i]=j; rmin=r; }
00351         }
00352 }
00353 
00354 // find the nearest point for each y point
00355 for (i=0; i<nb; i++) {
00356         float rmin=1.0e30f;
00357         for (j=0; j<na; j++) {
00358                 float r=(get_vector_at(j)+offset-to->get_vector_at(i)).length();
00359                 if (r<rmin) { b2a[i]=j; rmin=r; }
00360         }
00361 }
00362 
00363 // now keep only points where x->y matches y->x
00364 for (i=0; i<na; i++) {
00365         if (a2b[i]<0) continue;
00366         if (b2a[a2b[i]]!=i) { printf(" #%d!=%d# ",b2a[a2b[i]],i);  b2a[a2b[i]]=-1; a2b[i]=-1; }
00367         printf("%d->%d  ",i,a2b[i]);
00368 }
00369 printf("\n");
00370 
00371 for (i=0; i<nb; i++) {
00372         if (b2a[i]<0) continue;
00373         if (a2b[b2a[i]]!=i) { a2b[b2a[i]]=-1; b2a[i]=-1; }
00374         printf("%d->%d  ",i,b2a[i]);
00375 }
00376 printf("\n");
00377 
00378 // Compute the average translation required to align the points
00379 float dx=0,dy=0,dr=0,nr=0;
00380 for (i=0; i<na; i++) {
00381         if (a2b[i]==-1) continue;
00382         dx+=to->get_vector_at(a2b[i])[0]-get_vector_at(i)[0];
00383         dy+=to->get_vector_at(a2b[i])[1]-get_vector_at(i)[1];
00384         nr+=1.0;
00385 }
00386 //printf("%f %f %f\n",dx,dy,nr);
00387 if (nr<2) return vector<float>(4,0);
00388 dx/=nr;
00389 dy/=nr;
00390 
00391 // Compute the residual error
00392 for (i=0; i<na; i++) {
00393         if (i==-1  || a2b[i]==-1) continue;
00394         dr+=(to->get_vector_at(a2b[i])-get_vector_at(i)-Vec3f(dx,dy,0)).length();
00395 }
00396 dr/=nr;
00397 
00398 free(a2b);
00399 free(b2a);
00400 vector<float> ret(4);
00401 ret[0]=dx;
00402 ret[1]=dy;
00403 ret[2]=dr;
00404 ret[3]=(float)nr;
00405 return ret;
00406 }

void PointArray::center_to_zero  ) 
 

Definition at line 534 of file pointarray.cpp.

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

00535 {
00536         FloatPoint center = get_center();
00537         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00538                 points[i] -= center[0];
00539                 points[i + 1] -= center[1];
00540                 points[i + 2] -= center[2];
00541         }
00542 }

PointArray * PointArray::copy  ) 
 

Definition at line 121 of file pointarray.cpp.

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

Referenced by mask(), and mask_asymmetric_unit().

00122 {
00123         PointArray *pa2 = new PointArray();
00124         pa2->set_number_points(get_number_points());
00125         double *pa2data = pa2->get_points_array();
00126         memcpy(pa2data, get_points_array(), sizeof(double) * 4 * get_number_points());
00127 
00128         return pa2;
00129 }

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 163 of file pointarray.cpp.

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

Referenced by match_points().

00163                                        {
00164 if (n==0) return NULL;
00165 
00166 unsigned int i,j;
00167 
00168 EMData *ret= new EMData(n,n,1);
00169 ret->to_zero();
00170 
00171 for (i=0; i<n; i++) {
00172         for (j=i+1; j<n; j++) {
00173                 float r=(get_vector_at(i)-get_vector_at(j)).length();
00174                 ret->set_value_at(i,j,0,r);
00175                 ret->set_value_at(j,i,0,r);
00176         }
00177 }
00178 
00179 if (sortrows) {
00180         float *data=ret->get_data();
00181         for (i=0; i<n; i++) qsort(&data[i*n],n,sizeof(float),cmp_float);
00182         ret->update();
00183 }
00184 
00185 return ret;
00186 }

Region PointArray::get_bounding_box  ) 
 

Definition at line 544 of file pointarray.cpp.

References get_number_points(), and points.

00545 {
00546         double xmin, ymin, zmin;
00547         double xmax, ymax, zmax;
00548         xmin = xmax = points[0];
00549         ymin = ymax = points[1];
00550         zmin = zmax = points[2];
00551         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00552                 if (points[i] > xmax)
00553                         xmax = points[i];
00554                 if (points[i] < xmin)
00555                         xmin = points[i];
00556                 if (points[i + 1] > ymax)
00557                         ymax = points[i + 1];
00558                 if (points[i + 1] < ymin)
00559                         ymin = points[i + 1];
00560                 if (points[i + 2] > zmax)
00561                         zmax = points[i + 2];
00562                 if (points[i + 2] < zmin)
00563                         zmin = points[i + 2];
00564         }
00565         return Region(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin);
00566 }

FloatPoint PointArray::get_center  ) 
 

Definition at line 515 of file pointarray.cpp.

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

Referenced by center_to_zero().

00516 {
00517         double xc, yc, zc;
00518         xc = yc = zc = 0.0;
00519         double norm = 0.0;
00520         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00521                 xc += points[i] * points[i + 3];
00522                 yc += points[i + 1] * points[i + 3];
00523                 zc += points[i + 2] * points[i + 3];
00524                 norm += points[i + 3];
00525         }
00526         if (norm <= 0) {
00527                 fprintf(stderr, "Abnormal total value (%g) for PointArray, it should be >0\n", norm);
00528                 return FloatPoint(0, 0, 0);
00529         }
00530         else
00531                 return FloatPoint(xc / norm, yc / norm, zc / norm);
00532 }

size_t PointArray::get_number_points  )  const
 

Definition at line 140 of file pointarray.cpp.

References n.

Referenced by align_trans_2d(), center_to_zero(), copy(), get_bounding_box(), get_center(), mask(), mask_asymmetric_unit(), operator=(), opt_from_proj(), pdb2mrc_by_nfft(), pdb2mrc_by_summation(), projection_by_nfft(), projection_by_summation(), read_from_pdb(), save_to_pdb(), set_from(), and set_from_density_map().

00141 {
00142         return n;
00143 }

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 656 of file pointarray.cpp.

References n, and points.

00656                                      {
00657 vector<float> ret;
00658 for (unsigned int i=0; i<n; i++) {
00659         ret.push_back((float)points[i*4]);
00660         ret.push_back((float)points[i*4+1]);
00661         ret.push_back((float)points[i*4+2]);
00662 }
00663 
00664 return ret;
00665 }

double * PointArray::get_points_array  ) 
 

Definition at line 153 of file pointarray.cpp.

References points.

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

00154 {
00155         return points;
00156 }

double PointArray::get_value_at int  i  ) 
 

Definition at line 1653 of file pointarray.cpp.

References points.

Referenced by align_trans_2d().

01654 {
01655 return points[i*4+3];
01656 }

Vec3f PointArray::get_vector_at int  i  ) 
 

Definition at line 1648 of file pointarray.cpp.

References points.

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

01649 {
01650 return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
01651 }

void PointArray::mask double  rmax,
double  rmin = 0.0
 

Definition at line 569 of file pointarray.cpp.

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

00570 {
00571         double rmax2 = rmax * rmax, rmin2 = rmin * rmin;
00572         PointArray *tmp = this->copy();
00573         double *tmp_points = tmp->get_points_array();
00574          int count = 0;
00575         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00576                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v =
00577                         tmp_points[i + 3];
00578                 double r2 = x * x + y * y + z * z;
00579                 if (r2 >= rmin2 && r2 <= rmax2) {
00580                         points[count * 4] = x;
00581                         points[count * 4 + 1] = y;
00582                         points[count * 4 + 2] = z;
00583                         points[count * 4 + 3] = v;
00584                         ++count;
00585                 }
00586         }
00587         set_number_points(count);
00588         if( tmp )
00589         {
00590                 delete tmp;
00591                 tmp = 0;
00592         }
00593 }

void PointArray::mask_asymmetric_unit const string &  sym  ) 
 

Definition at line 596 of file pointarray.cpp.

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

00597 {
00598         if (sym == "c1" || sym == "C1")
00599                 return;                                 // do nothing for C1 symmetry
00600         double alt0 = 0, alt1 = M_PI, alt2 = M_PI;
00601         double az0 = 0, az1 = M_PI;
00602         if (sym[0] == 'c' || sym[0] == 'C') {
00603                 int nsym = atoi(sym.c_str() + 1);
00604                 az1 = 2.0 * M_PI / nsym / 2.0;
00605         }
00606         else if (sym[0] == 'd' || sym[0] == 'D') {
00607                 int nsym = atoi(sym.c_str() + 1);
00608                 alt1 = M_PI / 2.0;
00609                 alt2 = alt1;
00610                 az1 = 2.0 * M_PI / nsym / 2.0;
00611         }
00612         else if (sym == "icos" || sym == "ICOS") {
00613                 alt1 = 0.652358139784368185995; // 5fold to 3fold
00614                 alt2 = 0.55357435889704525151;  // half of edge ie. 5fold to 2fold along the edge
00615                 az1 = 2.0 * M_PI / 5 / 2.0;
00616         }
00617         else {
00618                 LOGERR("PointArray::set_to_asymmetric_unit(): sym = %s is not implemented yet",
00619                            sym.c_str());
00620                 return;
00621         }
00622 #ifdef DEBUG
00623         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);
00624 #endif
00625 
00626         PointArray *tmp = this->copy();
00627         double *tmp_points = tmp->get_points_array();
00628          int count = 0;
00629         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00630                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v = tmp_points[i + 3];
00631                 double az = atan2(y, x);
00632                 double az_abs = fabs(az - az0);
00633                 if (az_abs < (az1 - az0)) {
00634                         double alt_max = alt1 + (alt2 - alt1) * az_abs / (az1 - az0);
00635                         double alt = acos(z / sqrt(x * x + y * y + z * z));
00636                         if (alt < alt_max && alt >= alt0) {
00637 #ifdef DEBUG
00638                                 printf("Point %3d: 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);
00639 #endif
00640                                 points[count * 4] = x;
00641                                 points[count * 4 + 1] = y;
00642                                 points[count * 4 + 2] = z;
00643                                 points[count * 4 + 3] = v;
00644                                 count++;
00645                         }
00646                 }
00647         }
00648         set_number_points(count);
00649         if( tmp )
00650         {
00651                 delete tmp;
00652                 tmp = 0;
00653         }
00654 }

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

00188                                                                   {
00189 EMData *mx0=distmx(1);
00190 EMData *mx1=to->distmx(1);
00191 unsigned int n2=mx1->get_xsize();       // same as get_number_points on to
00192 
00193 if (max_miss<0) max_miss=(float)mx0->get_attr("sigma")/10.0f;
00194 //printf("max error %f\n",max_miss);
00195 
00196 
00197 
00198 vector<int> ret(n,-1);
00199 vector<float> rete(n,0.0);
00200 unsigned int i,j,k,l;
00201 
00202 if (!mx0 || !mx1) {
00203         if (mx0) delete mx0;
00204         if (mx1) delete mx1;
00205         return ret;
00206 }
00207 
00208 // i iterates over elements of 'this', j looks for a match in 'to'
00209 // k and l iterate over the individual distances
00210 for (i=0; i<n; i++) {
00211         int bestn=-1;                   // number of best match in mx1
00212         double bestd=1.0e38;            // residual error distance in best match
00213         for (j=0; j<n2; j++) {
00214                 double d=0;
00215                 int nd=0;
00216                 for (k=l=0; k<n-1 && l<n2-1; k++,l++) {
00217                         float d1=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l,j));
00218                         float d2=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l,j));
00219                         float d3=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l+1,j));
00220                         float d4=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l+1,j));
00221                         if (d2<d1 && d4>d2) { l--; continue; }
00222                         if (d3<d1 && d4>d3) { k--; continue; }
00223                         d+=d1;
00224                         nd++;
00225                 }
00226                 d/=(float)nd;
00227 //              printf("%d -> %d\t%f\t%d\n",i,j,d,nd);
00228                 if (d<bestd) { bestd=d; bestn=j; }
00229         }
00230         ret[i]=bestn;
00231         rete[i]=static_cast<float>(bestd);
00232 }
00233 
00234 // This will remove duplicate assignments, keeping the best one
00235 // also remove any assignments with large errors
00236 for (i=0; i<n; i++) {
00237         for (j=0; j<n; j++) {
00238                 if (rete[i]>max_miss) { ret[i]=-1; break; }
00239                 if (i==j || ret[i]!=ret[j] || ret[i]==-1) continue;
00240                 if (rete[i]>rete[j]) { ret[i]=-1; break; }
00241         }
00242 }
00243 
00244 delete mx0;
00245 delete mx1;
00246 
00247 return ret;
00248 }

PointArray & PointArray::operator= PointArray pa  ) 
 

Definition at line 131 of file pointarray.cpp.

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

00132 {
00133         if (this != &pa) {
00134                 set_number_points(pa.get_number_points());
00135                 memcpy(get_points_array(), pa.get_points_array(), sizeof(double) * 4 * get_number_points());
00136         }
00137         return *this;
00138 }

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:
proj A vector of EMData objects containing projections with orientations
pixres Size of each Gaussian in pixels

Definition at line 1625 of file pointarray.cpp.

References get_number_points(), and LOGWARN.

01625                                                                         {
01626 #ifdef OPTPP
01627         optdata=proj;
01628         optobj=this;
01629         optpixres=pixres;
01630 
01631         FDNLF1 nlf(get_number_points()*4,calc_opt_proj,init_opt_proj);
01632 //      NLF1 nlf(get_number_points()*4,init_opt_proj,calc_opt_projd);
01633         nlf.initFcn();
01634 
01635         OptCG opt(&nlf);
01636 //      OptQNewton opt(&nlf);
01637         opt.setMaxFeval(2000);
01638         opt.optimize();
01639         opt.printStatus("Done");
01640 #else
01641         (void)proj;             //suppress warning message
01642         (void)pixres;   //suppress warning message
01643         LOGWARN("OPT++ support not enabled.\n");
01644         return;
01645 #endif
01646 }

EMData * PointArray::pdb2mrc_by_nfft int  map_size,
float  apix,
float  res
 

Definition at line 1255 of file pointarray.cpp.

References data, EMAN::EMData::do_ift(), EMAN::EMData::get_data(), get_number_points(), LOGWARN, 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().

01256 {
01257 #if defined NFFT
01258         nfft_3D_plan my_plan;           // plan for the nfft
01259 
01261         nfft_3D_init(&my_plan, map_size, get_number_points());
01262 
01264         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01265                 // FFTW and nfft use row major array layout, EMAN uses column major
01266                 my_plan.v[3 * j + 2] = (fftw_real) (points[i] / (apix * map_size));
01267                 my_plan.v[3 * j + 1] = (fftw_real) (points[i + 1] / (apix * map_size));
01268                 my_plan.v[3 * j] = (fftw_real) (points[i + 2] / (apix * map_size));
01269                 my_plan.f[j].re = (fftw_real) (points[i + 3]);
01270                 my_plan.f[j].im = 0.0;
01271         }
01272 
01274         if (my_plan.nfft_flags & PRE_PSI) {
01275                 nfft_3D_precompute_psi(&my_plan);
01276         }
01277 
01278         // compute the uniform Fourier transform
01279         nfft_3D_transpose(&my_plan);
01280 
01281         // copy the Fourier transform to EMData data array
01282         EMData *fft = new EMData();
01283         fft->set_size(map_size + 2, map_size, map_size);
01284         fft->set_complex(true);
01285         fft->set_ri(true);
01286         fft->to_zero();
01287         float *data = fft->get_data();
01288         double norm = 1.0 / (map_size * map_size * map_size);
01289         for (int k = 0; k < map_size; k++) {
01290                 for (int j = 0; j < map_size; j++) {
01291                         for (int i = 0; i < map_size / 2; i++) {
01292                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01293                                         (float) (my_plan.
01294                                                          f_hat[k * map_size * map_size + j * map_size + i +
01295                                                                    map_size / 2].re) * norm;
01296                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01297                                         (float) (my_plan.
01298                                                          f_hat[k * map_size * map_size + j * map_size + i +
01299                                                                    map_size / 2].im) * norm * (-1.0);
01300                         }
01301                 }
01302         }
01304         nfft_3D_finalize(&my_plan);
01305 
01306         // low pass processor
01307         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01308         int index = 0;
01309         for (int k = 0; k < map_size; k++) {
01310                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01311                 for (int j = 0; j < map_size; j++) {
01312                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01313                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01314                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01315                                 data[index] *= val;
01316                                 data[index + 1] *= val;
01317                         }
01318                 }
01319         }
01320         fft->update();
01321         //fft->process_inplace("eman1.filter.lowpass.gaussian",Dict("lowpass", map_size*apix/res));
01322 
01323         fft->process_inplace("xform.phaseorigin.tocorner");     // move phase origin to center of image map_size, instead of at corner
01324         EMData *map = fft->do_ift();
01325         map->set_attr("apix_x", apix);
01326         map->set_attr("apix_y", apix);
01327         map->set_attr("apix_z", apix);
01328         map->set_attr("origin_x", -map_size/2*apix);
01329         map->set_attr("origin_y", -map_size/2*apix);
01330         map->set_attr("origin_z", -map_size/2*apix);
01331         if( fft )
01332         {
01333                 delete fft;
01334                 fft = 0;
01335         }
01336         return map;
01337 #elif defined NFFT2
01338         nfft_plan my_plan;                      // plan for the nfft
01339 
01341         nfft_init_3d(&my_plan, map_size, map_size, map_size, get_number_points());
01342 
01344         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01345                 // FFTW and nfft use row major array layout, EMAN uses column major
01346                 my_plan.x[3 * j + 2] = (double) (points[i] / (apix * map_size));
01347                 my_plan.x[3 * j + 1] = (double) (points[i + 1] / (apix * map_size));
01348                 my_plan.x[3 * j] = (double) (points[i + 2] / (apix * map_size));
01349                 my_plan.f[j][0] = (double) (points[i + 3]);
01350                 my_plan.f[j][1] = 0.0;
01351         }
01352 
01354         if (my_plan.nfft_flags & PRE_PSI) {
01355                 nfft_precompute_psi(&my_plan);
01356                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01357                         nfft_full_psi(&my_plan, pow(10, -10));
01358         }
01359 
01360         // compute the uniform Fourier transform
01361         nfft_transposed(&my_plan);
01362 
01363         // copy the Fourier transform to EMData data array
01364         EMData *fft = new EMData();
01365         fft->set_size(map_size + 2, map_size, map_size);
01366         fft->set_complex(true);
01367         fft->set_ri(true);
01368         fft->to_zero();
01369         float *data = fft->get_data();
01370         double norm = 1.0 / (map_size * map_size * map_size);
01371         for (int k = 0; k < map_size; k++) {
01372                 for (int j = 0; j < map_size; j++) {
01373                         for (int i = 0; i < map_size / 2; i++) {
01374                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01375                                         (float) (my_plan.
01376                                                          f_hat[k * map_size * map_size + j * map_size + i +
01377                                                                    map_size / 2][0]) * norm;
01378                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01379                                         (float) (my_plan.
01380                                                          f_hat[k * map_size * map_size + j * map_size + i +
01381                                                                    map_size / 2][1]) * norm;
01382                         }
01383                 }
01384         }
01386         nfft_finalize(&my_plan);
01387 
01388         // low pass processor
01389         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01390         int index = 0;
01391         for (int k = 0; k < map_size; k++) {
01392                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01393                 for (int j = 0; j < map_size; j++) {
01394                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01395                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01396                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01397                                 data[index] *= val;
01398                                 data[index + 1] *= val;
01399                         }
01400                 }
01401         }
01402         fft->update();
01403         //fft->process_inplace("eman1.filter.lowpass.gaussian",Dict("lowpass", map_size*apix/res));
01404 
01405         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image map_size, instead of at corner
01406         EMData *map = fft->do_ift();
01407         map->set_attr("apix_x", apix);
01408         map->set_attr("apix_y", apix);
01409         map->set_attr("apix_z", apix);
01410         map->set_attr("origin_x", -map_size / 2 * apix);
01411         map->set_attr("origin_y", -map_size / 2 * apix);
01412         map->set_attr("origin_z", -map_size / 2 * apix);
01413         if( fft )
01414         {
01415                 delete fft;
01416                 fft = 0;
01417         }
01418         return map;
01419 #else
01420         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01421         return 0;
01422 #endif
01423 }

EMData * PointArray::pdb2mrc_by_summation int  map_size,
float  apix,
float  res
 

Definition at line 1021 of file pointarray.cpp.

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

01022 {
01023 #ifdef DEBUG
01024         printf("PointArray::pdb2mrc_by_summation(): %d points\tmapsize = %4d\tapix = %g\tres = %g\n",get_number_points(),map_size, apix, res);
01025 #endif
01026         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01027         //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);
01028 
01029         double min_table_val = 1e-7;
01030         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01031 
01032         double table_step_size = 0.001; // number of steps for each pixel
01033         double inv_table_step_size = 1.0 / table_step_size;
01034         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01035         vector<double> table;
01036         table.resize(table_size);
01037         //double *table = (double *) malloc(sizeof(double) * table_size);
01038         for (int i = 0; i < table_size; i++) {
01039                 double x = -i * table_step_size * apix / gauss_real_width;
01040                 table[i] = exp(-x * x);
01041         }
01042 
01043         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01044         if (gbox <= 0)
01045                 gbox = 1;
01046 
01047         sort_by_axis(2);                        // sort by Z-axis
01048 
01049         EMData *map = new EMData();
01050         map->set_size(map_size, map_size, map_size);
01051         map->to_zero();
01052         float *pd = map->get_data();
01053         for ( size_t s = 0; s < get_number_points(); ++s) {
01054                 double xc = points[4 * s] / apix + map_size / 2;
01055                 double yc = points[4 * s + 1] / apix + map_size / 2;
01056                 double zc = points[4 * s + 2] / apix + map_size / 2;
01057                 double fval = points[4 * s + 3];
01058                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01059                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01060                 int kmin = int (zc) - gbox, kmax = int (zc) + gbox;
01061                 if (imin < 0)
01062                         imin = 0;
01063                 if (jmin < 0)
01064                         jmin = 0;
01065                 if (kmin < 0)
01066                         kmin = 0;
01067                 if (imax > map_size)
01068                         imax = map_size;
01069                 if (jmax > map_size)
01070                         jmax = map_size;
01071                 if (kmax > map_size)
01072                         kmax = map_size;
01073 
01074                 for (int k = kmin; k < kmax; k++) {
01075                         size_t table_index_z = size_t (fabs(k - zc) * inv_table_step_size);
01076                         if ( table_index_z >= table.size() ) continue;
01077                         double zval = table[table_index_z];
01078                         size_t pd_index_z = k * map_size * map_size;
01079                         for (int j = jmin; j < jmax; j++) {
01080                                 size_t table_index_y = size_t (fabs(j - yc) * inv_table_step_size);
01081                                 if ( table_index_y >= table.size() ) continue;
01082                                 double yval = table[table_index_y];
01083                                 size_t pd_index = pd_index_z + j * map_size + imin;
01084                                 for (int i = imin; i < imax; i++, pd_index++) {
01085                                         size_t table_index_x = size_t (fabs(i - xc) * inv_table_step_size);
01086                                         if ( table_index_x >= table.size() ) continue;
01087                                         double xval = table[table_index_x];
01088                                         pd[pd_index] += (float) (fval * zval * yval * xval);
01089                                 }
01090                         }
01091                 }
01092         }
01093         //for(int i=0; i<map_size*map_size; i++) pd[i]/=sqrt(M_PI);
01094         map->update();
01095         map->set_attr("apix_x", apix);
01096         map->set_attr("apix_y", apix);
01097         map->set_attr("apix_z", apix);
01098         map->set_attr("origin_x", -map_size/2*apix);
01099         map->set_attr("origin_y", -map_size/2*apix);
01100         map->set_attr("origin_z", -map_size/2*apix);
01101 
01102         return map;
01103 }

EMData * PointArray::projection_by_nfft int  image_size,
float  apix,
float  res = 0
 

Definition at line 1425 of file pointarray.cpp.

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

01426 {
01427 #if defined NFFT
01428         nfft_2D_plan my_plan;           // plan for the nfft
01429         int N[2], n[2];
01430         N[0] = image_size;
01431         n[0] = next_power_of_2(2 * image_size);
01432         N[1] = image_size;
01433         n[1] = next_power_of_2(2 * image_size);
01434 
01436         nfft_2D_init(&my_plan, image_size, get_number_points());
01437         //nfft_2D_init_specific(&my_plan, N, get_number_points(), n, 3,
01438         //                 PRE_PHI_HUT | PRE_PSI,
01439         //                 FFTW_ESTIMATE | FFTW_OUT_OF_PLACE);
01441         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01442                 // FFTW and nfft use row major array layout, EMAN uses column major
01443                 my_plan.v[2 * j + 1] = (fftw_real) (points[i] / (apix * image_size));
01444                 my_plan.v[2 * j] = (fftw_real) (points[i + 1] / (apix * image_size));
01445                 my_plan.f[j].re = (fftw_real) points[i + 3];
01446                 my_plan.f[j].im = 0.0;
01447         }
01448 
01450         if (my_plan.nfft_flags & PRE_PSI) {
01451                 nfft_2D_precompute_psi(&my_plan);
01452         }
01453 
01454         // compute the uniform Fourier transform
01455         nfft_2D_transpose(&my_plan);
01456 
01457         // copy the Fourier transform to EMData data array
01458         EMData *fft = new EMData();
01459         fft->set_size(image_size + 2, image_size, 1);
01460         fft->set_complex(true);
01461         fft->set_ri(true);
01462         fft->to_zero();
01463         float *data = fft->get_data();
01464         double norm = 1.0 / (image_size * image_size);
01465         for (int j = 0; j < image_size; j++) {
01466                 for (int i = 0; i < image_size / 2; i++) {
01467                         data[j * (image_size + 2) + 2 * i] =
01468                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].re) * norm;
01469                         data[j * (image_size + 2) + 2 * i + 1] =
01470                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].im) * norm * (-1.0);
01471                 }
01472         }
01474         nfft_2D_finalize(&my_plan);
01475 
01476         if (res > 0) {
01477                 // Gaussian low pass processor
01478                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01479                 int index = 0;
01480                 for (int j = 0; j < image_size; j++) {
01481                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01482                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01483                                 double val = exp(-(i * i + RY2) / sigma2);
01484                                 data[index] *= val;
01485                                 data[index + 1] *= val;
01486                         }
01487                 }
01488         }
01489         fft->update();
01490         //fft->process_inplace("eman1.filter.lowpass.gaussian",Dict("lowpass", box*apix/res));
01491 
01492         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01493 
01494         return fft;
01495 #elif defined NFFT2
01496         nfft_plan my_plan;                      // plan for the nfft
01497         int N[2], n[2];
01498         N[0] = image_size;
01499         n[0] = next_power_of_2(2 * image_size);
01500         N[1] = image_size;
01501         n[1] = next_power_of_2(2 * image_size);
01502 
01504         //nfft_init_2d(&my_plan,image_size,image_size,get_number_points());
01505         nfft_init_specific(&my_plan, 2, N, get_number_points(), n, 3,
01506                                            PRE_PHI_HUT | PRE_PSI |
01507                                            MALLOC_X | MALLOC_F_HAT | MALLOC_F, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01509         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01510                 // FFTW and nfft use row major array layout, EMAN uses column major
01511                 my_plan.x[2 * j + 1] = (double) (points[i] / (apix * image_size));
01512                 my_plan.x[2 * j] = (double) (points[i + 1] / (apix * image_size));
01513                 my_plan.f[j][0] = (double) points[i + 3];
01514                 my_plan.f[j][1] = 0.0;
01515         }
01516 
01518         if (my_plan.nfft_flags & PRE_PSI) {
01519                 nfft_precompute_psi(&my_plan);
01520                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01521                         nfft_full_psi(&my_plan, pow(10, -6));
01522         }
01523 
01524         // compute the uniform Fourier transform
01525         nfft_transposed(&my_plan);
01526 
01527         // copy the Fourier transform to EMData data array
01528         EMData *fft = new EMData();
01529         fft->set_size(image_size + 2, image_size, 1);
01530         fft->set_complex(true);
01531         fft->set_ri(true);
01532         fft->to_zero();
01533         float *data = fft->get_data();
01534         double norm = 1.0 / (image_size * image_size);
01535         for (int j = 0; j < image_size; j++) {
01536                 for (int i = 0; i < image_size / 2; i++) {
01537                         data[j * (image_size + 2) + 2 * i] =
01538                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][0]) * norm;
01539                         data[j * (image_size + 2) + 2 * i + 1] =
01540                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][1]) * norm;
01541                 }
01542         }
01544         nfft_finalize(&my_plan);
01545 
01546         if (res > 0) {
01547                 // Gaussian low pass process
01548                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01549                 int index = 0;
01550                 for (int j = 0; j < image_size; j++) {
01551                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01552                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01553                                 double val = exp(-(i * i + RY2) / sigma2);
01554                                 data[index] *= val;
01555                                 data[index + 1] *= val;
01556                         }
01557                 }
01558         }
01559         fft->update();
01560         //fft->process_inplace("eman1.filter.lowpass.gaussian",Dict("lowpass", box*apix/res));
01561 
01562         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01563 
01564         return fft;
01565 #else
01566         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01567         return 0;
01568 #endif
01569 }

EMData * PointArray::projection_by_summation int  image_size,
float  apix,
float  res
 

Definition at line 1106 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.

01107 {
01108         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01109         //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);
01110 
01111         double min_table_val = 1e-7;
01112         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01113 
01114         //double table_step_size = 0.001;    // number of steps for x=[0,1] in exp(-x*x)
01115         //int table_size = int(max_table_x / table_step_size *1.25);
01116         //double* table = (double*)malloc(sizeof(double) * table_size);
01117         //for(int i=0; i<table_size; i++) table[i]=exp(-i*i*table_step_size*table_step_size);
01118 
01119         double table_step_size = 0.001; // number of steps for each pixel
01120         double inv_table_step_size = 1.0 / table_step_size;
01121         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01122         double *table = (double *) malloc(sizeof(double) * table_size);
01123         for (int i = 0; i < table_size; i++) {
01124                 double x = -i * table_step_size * apix / gauss_real_width;
01125                 table[i] = exp(-x * x);
01126         }
01127 
01128         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01129         if (gbox <= 0)
01130                 gbox = 1;
01131         EMData *proj = new EMData();
01132         proj->set_size(image_size, image_size, 1);
01133         proj->to_zero();
01134         float *pd = proj->get_data();
01135         for ( size_t s = 0; s < get_number_points(); ++s) {
01136                 double xc = points[4 * s] / apix + image_size / 2;
01137                 double yc = points[4 * s + 1] / apix + image_size / 2;
01138                 double fval = points[4 * s + 3];
01139                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01140                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01141                 if (imin < 0)
01142                         imin = 0;
01143                 if (jmin < 0)
01144                         jmin = 0;
01145                 if (imax > image_size)
01146                         imax = image_size;
01147                 if (jmax > image_size)
01148                         jmax = image_size;
01149 
01150                 for (int j = jmin; j < jmax; j++) {
01151                         //int table_index_y = int(fabs(j-yc)*apix/gauss_real_width/table_step_size);
01152                         int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01153                         double yval = table[table_index_y];
01154 #ifdef DEBUG
01155                         //double yval2 = exp( - (j-yc)*(j-yc)*apix*apix/(gauss_real_width*gauss_real_width));
01156                         //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);
01157 #endif
01158                         int pd_index = j * image_size + imin;
01159                         for (int i = imin; i < imax; i++, pd_index++) {
01160                                 //int table_index_x = int(fabs(i-xc)*apix/gauss_real_width/table_step_size);
01161                                 int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01162                                 double xval = table[table_index_x];
01163 #ifdef DEBUG
01164                                 //double xval2 = exp( - (i-xc)*(i-xc)*apix*apix/(gauss_real_width*gauss_real_width));
01165                                 //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);
01166 #endif
01167                                 pd[pd_index] += (float)(fval * yval * xval);
01168                         }
01169                 }
01170         }
01171         for (int i = 0; i < image_size * image_size; i++)
01172                 pd[i] /= sqrt(M_PI);
01173         proj->update();
01174         return proj;
01175 }

bool PointArray::read_from_pdb const char *  file  ) 
 

Definition at line 409 of file pointarray.cpp.

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

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

00410 {
00411         struct stat filestat;
00412         stat(file, &filestat);
00413         set_number_points(( int)(filestat.st_size / 80 + 1));
00414 #ifdef DEBUG
00415         printf("PointArray::read_from_pdb(): try %4d atoms first\n", get_number_points());
00416 #endif
00417 
00418         FILE *fp = fopen(file, "r");
00419         if(!fp) {
00420                 fprintf(stderr,"ERROR in PointArray::read_from_pdb(): cannot open file %s\n",file);
00421                 throw;
00422         }
00423         char s[200];
00424         size_t count = 0;
00425         while ((fgets(s, 200, fp) != NULL)) {
00426                 if (strncmp(s, "ENDMDL", 6) == 0)
00427                         break;
00428                 if (strncmp(s, "ATOM", 4) != 0)
00429                         continue;
00430 
00431                 if (s[13] == ' ')
00432                         s[13] = s[14];
00433                 if (s[13] == ' ')
00434                         s[13] = s[15];
00435 
00436                 float e = 0;
00437                 char ctt, ctt2 = ' ';
00438                 if (s[13] == ' ')
00439                         ctt = s[14];
00440                 else if (s[12] == ' ') {
00441                         ctt = s[13];
00442                         ctt2 = s[14];
00443                 }
00444                 else {
00445                         ctt = s[12];
00446                         ctt2 = s[13];
00447                 }
00448 
00449                 switch (ctt) {
00450                 case 'H':
00451                         e = 1.0;
00452                         break;
00453                 case 'C':
00454                         e = 6.0;
00455                         break;
00456                 case 'A':
00457                         if (ctt2 == 'U') {
00458                                 e = 79.0;
00459                                 break;
00460                         }
00461                         // treat 'A'mbiguous atoms as N, not perfect, but good enough
00462                 case 'N':
00463                         e = 7.0;
00464                         break;
00465                 case 'O':
00466                         e = 8.0;
00467                         break;
00468                 case 'P':
00469                         e = 15.0;
00470                         break;
00471                 case 'S':
00472                         e = 16.0;
00473                         break;
00474                 case 'W':
00475                         e = 18.0;
00476                         break;                          // ficticious water 'atom'
00477                 default:
00478                         fprintf(stderr, "Unknown atom %c%c\n", ctt, ctt2);
00479                         e = 0;
00480                 }
00481                 if (e == 0)
00482                         continue;
00483 
00484                 float x, y, z;
00485                 sscanf(&s[28], " %f %f %f", &x, &y, &z);
00486 
00487                 if (count + 1 > get_number_points())
00488                         set_number_points(2 * (count + 1));    //makes sure point array is big enough
00489 #ifdef DEBUG
00490                 printf("Atom %4d: x,y,z = %8g,%8g,%8g\te = %g\n", count, x, y, z, e);
00491 #endif
00492                 points[4 * count] = x;
00493                 points[4 * count + 1] = y;
00494                 points[4 * count + 2] = z;
00495                 points[4 * count + 3] = e;
00496                 count++;
00497         }
00498         fclose(fp);
00499         set_number_points(count);
00500         return true;
00501 }

void PointArray::replace_by_summation EMData image,
int  i,
Vec3f  vec,
float  amp,
float  apix,
float  res
 

Definition at line 1177 of file pointarray.cpp.

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

01178 {
01179         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01180 
01181         double min_table_val = 1e-7;
01182         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01183 
01184         double table_step_size = 0.001; // number of steps for each pixel
01185         double inv_table_step_size = 1.0 / table_step_size;
01186         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01187         double *table = (double *) malloc(sizeof(double) * table_size);
01188         for (int i = 0; i < table_size; i++) {
01189                 double x = -i * table_step_size * apix / gauss_real_width;
01190                 table[i] = exp(-x * x)/pow((float)M_PI,.25f);
01191         }
01192         int image_size=proj->get_xsize();
01193 
01194         // subtract the old point
01195         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01196         if (gbox <= 0)
01197                 gbox = 1;
01198         float *pd = proj->get_data();
01199          int s = ind;
01200         double xc = points[4 * s] / apix + image_size / 2;
01201         double yc = points[4 * s + 1] / apix + image_size / 2;
01202         double fval = points[4 * s + 3];
01203         int imin = int (xc) - gbox, imax = int (xc) + gbox;
01204         int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01205 
01206         if (imin < 0) imin = 0;
01207         if (jmin < 0) jmin = 0;
01208         if (imax > image_size) imax = image_size;
01209         if (jmax > image_size) jmax = image_size;
01210 
01211         for (int j = jmin; j < jmax; j++) {
01212                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01213                 double yval = table[table_index_y];
01214                 int pd_index = j * image_size + imin;
01215                 for (int i = imin; i < imax; i++, pd_index++) {
01216                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01217                         double xval = table[table_index_x];
01218                         pd[pd_index] -= (float)(fval * yval * xval);
01219                 }
01220         }
01221 
01222         // add the new point
01223         gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
01224         if (gbox <= 0)
01225                 gbox = 1;
01226         pd = proj->get_data();
01227         s = ind;
01228         xc = vec[0] / apix + image_size / 2;
01229         yc = vec[1] / apix + image_size / 2;
01230         fval = amp;
01231         imin = int (xc) - gbox, imax = int (xc) + gbox;
01232         jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01233 
01234         if (imin < 0) imin = 0;
01235         if (jmin < 0) jmin = 0;
01236         if (imax > image_size) imax = image_size;
01237         if (jmax > image_size) jmax = image_size;
01238 
01239         for (int j = jmin; j < jmax; j++) {
01240                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01241                 double yval = table[table_index_y];
01242                 int pd_index = j * image_size + imin;
01243                 for (int i = imin; i < imax; i++, pd_index++) {
01244                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01245                         double xval = table[table_index_x];
01246                         pd[pd_index] -= (float)(fval * yval * xval);
01247                 }
01248         }
01249 
01250         proj->update();
01251         return;
01252 }

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 679 of file pointarray.cpp.

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

00679                                                            {
00680         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00681                 Transform s = transform.transpose();
00682                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00683                 v= s*v;
00684                 points[i]  =v[0];
00685                 points[i+1]=v[1];
00686                 points[i+2]=v[2];
00687         }
00688 
00689 }

void PointArray::save_to_pdb const char *  file  ) 
 

Definition at line 504 of file pointarray.cpp.

References get_number_points(), and points.

00505 {
00506         FILE *fp = fopen(file, "w");
00507         for ( size_t i = 0; i < get_number_points(); i++) {
00508                 fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
00509                                 points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
00510         }
00511         fclose(fp);
00512 }

void PointArray::set_from double *  source,
int  num,
const string &  sym = "",
Transform3D transform = 0
 

Definition at line 696 of file pointarray.cpp.

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

00697 {
00698          int nsym = xform->get_nsym(sym);
00699 
00700         if (get_number_points() != (size_t)nsym * num)
00701                 set_number_points((size_t)nsym * num);
00702 
00703         double *target = get_points_array();
00704 
00705         for ( int s = 0; s < nsym; s++) {
00706                 int index = s * 4 * num;
00707                 for ( int i = 0; i < 4 * num; i += 4, index += 4) {
00708                         Vec3f v((float)src[i],(float)src[i+1],(float)src[i+2]);
00709                         v=v*xform->get_sym(sym,s);
00710                         target[index]  =v[0];
00711                         target[index+1]=v[1];
00712                         target[index+2]=v[2];
00713                         target[index+3]=src[i+3];
00714                 }
00715         }
00716 }

void PointArray::set_from PointArray source,
const string &  sym = "",
Transform3D transform = 0
 

Definition at line 690 of file pointarray.cpp.

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

00691 {
00692         set_from(source->get_points_array(), source->get_number_points(), sym, transform);
00693 
00694 }

void PointArray::set_from vector< float >   ) 
 

Definition at line 718 of file pointarray.cpp.

References points, and set_number_points().

Referenced by set_from().

00718                                            {
00719         set_number_points(pts.size()/4);
00720         for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];
00721 
00722 }

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

Definition at line 724 of file pointarray.cpp.

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

00726 {
00727         if (mode == PEAKS_SUB || mode == PEAKS_DIV) {
00728                 // find out how many voxels are useful voxels
00729                 int num_voxels = 0;
00730                 int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00731                 size_t size = nx * ny * nz;
00732                 EMData *tmp_map = map->copy();
00733                 float *pd = tmp_map->get_data();
00734                 for (size_t i = 0; i < size; ++i) {
00735                         if (pd[i] > thresh)
00736                                 ++num_voxels;
00737                 }
00738 
00739                 double pointvol = double (num_voxels) / double (num);
00740                 double gauss_real_width = pow(pointvol, 1. / 3.);       // in pixels
00741 #ifdef DEBUG
00742                 printf("Average point range radius = %g pixels for %d points from %d used voxels\n",
00743                            gauss_real_width, num, num_voxels);
00744 #endif
00745 
00746                 double min_table_val = 1e-4;
00747                 double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
00748 
00749                 double table_step_size = 1.;    // number of steps for each pixel
00750                 double inv_table_step_size = 1.0 / table_step_size;
00751                 int table_size = int (max_table_x * gauss_real_width / (table_step_size) * 1.25) + 1;
00752                 double *table = (double *) malloc(sizeof(double) * table_size);
00753                 for (int i = 0; i < table_size; i++) {
00754                         double x = i * table_step_size / gauss_real_width;
00755                         table[i] = exp(-x * x);
00756                 }
00757 
00758                 int gbox = int (max_table_x * gauss_real_width);        // local box half size in pixels to consider for each point
00759                 if (gbox <= 0)
00760                         gbox = 1;
00761 
00762                 set_number_points(num);
00763                 for (int count = 0; count < num; count++) {
00764                         float cmax = pd[0];
00765                         int cmaxpos = 0;
00766                         for (size_t i = 0; i < size; ++i) {
00767                                 if (pd[i] > cmax) {
00768                                         cmax = pd[i];
00769                                         cmaxpos = i;
00770                                 }
00771                         }
00772                         int iz = cmaxpos / (nx * ny), iy = (cmaxpos - iz * nx * ny) / nx, ix =
00773                                 cmaxpos - iz * nx * ny - iy * nx;
00774 
00775                         // update coordinates in pixels
00776                         points[4*count  ] = ix;
00777                         points[4*count+1] = iy;
00778                         points[4*count+2] = iz;
00779                         points[4*count+3] = cmax;
00780 #ifdef DEBUG
00781                         printf("Point %d: val = %g\tat  %d, %d, %d\n", count, cmax, ix, iy, iz);
00782 #endif
00783 
00784                         int imin = ix - gbox, imax = ix + gbox;
00785                         int jmin = iy - gbox, jmax = iy + gbox;
00786                         int kmin = iz - gbox, kmax = iz + gbox;
00787                         if (imin < 0)
00788                                 imin = 0;
00789                         if (jmin < 0)
00790                                 jmin = 0;
00791                         if (kmin < 0)
00792                                 kmin = 0;
00793                         if (imax > nx)
00794                                 imax = nx;
00795                         if (jmax > ny)
00796                                 jmax = ny;
00797                         if (kmax > nz)
00798                                 kmax = nz;
00799 
00800                         for (int k = kmin; k < kmax; ++k) {
00801                                 int table_index_z = int (fabs(double (k - iz)) * inv_table_step_size);
00802                                 double zval = table[table_index_z];
00803                                 size_t pd_index_z = k * nx * ny;
00804                                 //printf("k = %8d\tx = %8g\tval = %8g\n", k, float(k-iz), zval);
00805                                 for (int j = jmin; j < jmax; ++j) {
00806                                         int table_index_y = int (fabs(double (j - iy)) * inv_table_step_size);
00807                                         double yval = table[table_index_y];
00808                                         size_t pd_index = pd_index_z + j * nx + imin;
00809                                         for (int i = imin; i < imax; ++i, ++pd_index) {
00810                                                 int table_index_x = int (fabs(double (i - ix)) * inv_table_step_size);
00811                                                 double xval = table[table_index_x];
00812                                                 if (mode == PEAKS_SUB)
00813                                                         pd[pd_index] -= (float)(cmax * zval * yval * xval);
00814                                                 else
00815                                                         pd[pd_index] *= (float)(1.0 - zval * yval * xval);      // mode == PEAKS_DIV
00816                                         }
00817                                 }
00818                         }
00819                 }
00820                 set_number_points(num);
00821                 tmp_map->update();
00822                 if( tmp_map )
00823                 {
00824                         delete tmp_map;
00825                         tmp_map = 0;
00826                 }
00827         }
00828         else if (mode == KMEANS) {
00829                 set_number_points(num);
00830                 zero();
00831 
00832                 PointArray tmp_pa;
00833                 tmp_pa.set_number_points(num);
00834                 tmp_pa.zero();
00835 
00836                  int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00837                 float *pd = map->get_data();
00838 
00839                 // initialize segments with random centers at pixels with values > thresh
00840 #ifdef DEBUG
00841                 printf("Start initial random seeding\n");
00842 #endif
00843                 for ( size_t i = 0; i < get_number_points(); i++) {
00844                          int x, y, z;
00845                         double v;
00846                         do {
00847                                 x = ( int) Util::get_frand(0, nx - 1);
00848                                 y = ( int) Util::get_frand(0, ny - 1);
00849                                 z = ( int) Util::get_frand(0, nz - 1);
00850                                 v = pd[z * nx * ny + y * nx + x];
00851 #ifdef DEBUG
00852                                 printf("Trying Point %d: val = %g\tat  %d, %d, %d\tfrom map (%d,%d,%d)\n", i, v, x,
00853                                            y, z, nx, ny, nz);
00854 #endif
00855                         } while (v <= thresh);
00856                         points[4 * i] = (double) x;
00857                         points[4 * i + 1] = (double) y;
00858                         points[4 * i + 2] = (double) z;
00859                         points[4 * i + 3] = (double) v;
00860 #ifdef DEBUG
00861                         printf("Point %d: val = %g\tat  %g, %g, %g\n", i, points[4 * i + 3], points[4 * i],
00862                                    points[4 * i + 1], points[4 * i + 2]);
00863 #endif
00864                 }
00865 
00866                 double min_dcen = 1e0;  // minimal mean segment center shift as convergence criterion
00867                 double dcen = 0.0;
00868                 int iter = 0, max_iter = 100;
00869                 do {
00870 #ifdef DEBUG
00871                         printf("Iteration %3d, start\n", iter);
00872 #endif
00873                         double *tmp_points = tmp_pa.get_points_array();
00874 
00875                         // reassign each pixel to the best segment
00876                         for ( int k = 0; k < nz; k++) {
00877                                 for ( int j = 0; j < ny; j++) {
00878                                         for ( int i = 0; i < nx; i++) {
00879                                                 size_t idx = k * nx * ny + j * nx + i;
00880                                                 if (pd[idx] > thresh) {
00881                                                         double min_dist = 1e60; // just a large distance
00882                                                          int min_s = 0;
00883                                                         for ( size_t s = 0; s < get_number_points(); ++s) {
00884                                                                 double x = points[4 * s];
00885                                                                 double y = points[4 * s + 1];
00886                                                                 double z = points[4 * s + 2];
00887                                                                 double dist =
00888                                                                         (k - z) * (k - z) + (j - y) * (j - y) + (i - x) * (i - x);
00889                                                                 if (dist < min_dist) {
00890                                                                         min_dist = dist;
00891                                                                         min_s = s;
00892                                                                 }
00893                                                         }
00894                                                         tmp_points[4 * min_s] += i;
00895                                                         tmp_points[4 * min_s + 1] += j;
00896                                                         tmp_points[4 * min_s + 2] += k;
00897                                                         tmp_points[4 * min_s + 3] += 1.0;
00898                                                 }
00899                                         }
00900                                 }
00901                         }
00902 #ifdef DEBUG
00903                         printf("Iteration %3d, finished reassigning segments\n", iter);
00904 #endif
00905                         // update each segment's center
00906                         dcen = 0.0;
00907                         for ( size_t s = 0; s < get_number_points(); ++s) {
00908                                 if (tmp_points[4 * s + 3]) {
00909                                         tmp_points[4 * s] /= tmp_points[4 * s + 3];
00910                                         tmp_points[4 * s + 1] /= tmp_points[4 * s + 3];
00911                                         tmp_points[4 * s + 2] /= tmp_points[4 * s + 3];
00912 #ifdef DEBUG
00913                                         printf("Iteration %3d, Point %3d at %8g, %8g, %8g -> %8g, %8g, %8g\n", iter, s,
00914                                                    points[4 * s], points[4 * s + 1], points[4 * s + 2], tmp_points[4 * s],
00915                                                    tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00916 #endif
00917                                 }
00918                                 else {                  // empty segments are reseeded
00919                                          int x, y, z;
00920                                         double v;
00921                                         do {
00922                                                 x = ( int) Util::get_frand(0, nx - 1);
00923                                                 y = ( int) Util::get_frand(0, ny - 1);
00924                                                 z = ( int) Util::get_frand(0, nz - 1);
00925                                                 v = pd[z * nx * ny + y * nx + x];
00926                                         } while (v <= thresh);
00927                                         tmp_points[4 * s] = (double) x;
00928                                         tmp_points[4 * s + 1] = (double) y;
00929                                         tmp_points[4 * s + 2] = (double) z;
00930                                         tmp_points[4 * s + 3] = (double) v;
00931 #ifdef DEBUG
00932                                         printf
00933                                                 ("Iteration %3d, Point %3d reseeded from %8g, %8g, %8g -> %8g, %8g, %8g\n",
00934                                                  iter, s, points[4 * s], points[4 * s + 1], points[4 * s + 2],
00935                                                  tmp_points[4 * s], tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00936 #endif
00937                                 }
00938                                 double dx = tmp_points[4 * s] - points[4 * s];
00939                                 double dy = tmp_points[4 * s + 1] - points[4 * s + 1];
00940                                 double dz = tmp_points[4 * s + 2] - points[4 * s + 2];
00941                                 dcen += dx * dx + dy * dy + dz * dz;
00942                         }
00943                         dcen = sqrt(dcen / get_number_points());
00944                         //swap pointter, faster but risky
00945 #ifdef DEBUG
00946                         printf("before swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
00947                                    (long)tmp_pa.get_points_array());
00948 #endif
00949                         double *tp = get_points_array();
00950                         set_points_array(tmp_points);
00951                         tmp_pa.set_points_array(tp);
00952                         tmp_pa.zero();
00953 #ifdef DEBUG
00954                         printf("after  swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
00955                                    (long)tmp_pa.get_points_array());
00956                         printf("Iteration %3d, finished updating segment centers with dcen = %g pixels\n", iter,
00957                                    dcen);
00958 #endif
00959 
00960                         iter++;
00961                 } while (dcen > min_dcen && iter <= max_iter);
00962                 map->update();
00963 
00964                 sort_by_axis(2);        // x,y,z axes = 0, 1, 2
00965         }
00966         else {
00967                 LOGERR("PointArray::set_from_density_map(): mode = %d is not implemented yet", mode);
00968         }
00969         //update to use apix and origin
00970         int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00971         float origx, origy, origz;
00972         try {
00973                 origx = map->get_attr("origin_x");
00974                 origy = map->get_attr("origin_y");
00975                 origz = map->get_attr("origin_z");
00976         }
00977         catch(...) {
00978                 origx = -nx / 2 * apix;
00979                 origy = -ny / 2 * apix;
00980                 origz = -nz / 2 * apix;
00981         }
00982 
00983 #ifdef DEBUG
00984         printf("Apix = %g\torigin x,y,z = %8g,%8g,%8g\n",apix, origx, origy, origz);
00985 #endif
00986 
00987         float *pd = map->get_data();
00988         for ( size_t i = 0; i < get_number_points(); ++i) {
00989 #ifdef DEBUG
00990                 printf("Point %4d: x,y,z,v = %8g,%8g,%8g,%8g",i, points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
00991 #endif
00992                 points[4 * i + 3] =
00993                         pd[(int) points[4 * i + 2] * nx * ny + (int) points[4 * i + 1] * nx +
00994                            (int) points[4 * i]];
00995                 points[4 * i] = points[4 * i] * apix + origx;
00996                 points[4 * i + 1] = points[4 * i + 1] * apix + origy;
00997                 points[4 * i + 2] = points[4 * i + 2] * apix + origz;
00998 #ifdef DEBUG
00999                 printf("\t->\t%8g,%8g,%8g,%8g\n",points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
01000 #endif
01001         }
01002         map->update();
01003 }

void PointArray::set_number_points size_t  nn  ) 
 

Definition at line 145 of file pointarray.cpp.

References n, and points.

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

00146 {
00147         if (n != nn) {
00148                 n = nn;
00149                 points = (double *) realloc(points, 4 * n * sizeof(double));
00150         }
00151 }

void PointArray::set_points_array double *  p  ) 
 

Definition at line 158 of file pointarray.cpp.

References points.

Referenced by set_from_density_map().

00159 {
00160         points = p;
00161 }

void PointArray::set_vector_at int  i,
vector< double > 
 

Definition at line 1666 of file pointarray.cpp.

References points.

01667 {
01668 points[i*4]  =v[0];
01669 points[i*4+1]=v[1];
01670 points[i*4+2]=v[2];
01671 if (v.size()>=4) points[i*4+3]=v[3];
01672 }

void PointArray::set_vector_at int  i,
Vec3f  vec,
double  value
 

Definition at line 1658 of file pointarray.cpp.

References points.

01659 {
01660 points[i*4]=vec[0];
01661 points[i*4+1]=vec[1];
01662 points[i*4+2]=vec[2];
01663 points[i*4+3]=value;
01664 }

void PointArray::sort_by_axis int  axis = 1  ) 
 

Definition at line 1008 of file pointarray.cpp.

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

Referenced by pdb2mrc_by_summation(), and set_from_density_map().

01009 {
01010         if (axis == 0)
01011                 qsort(points, n, sizeof(double) * 4, cmp_axis_x);
01012         else if (axis == 1)
01013                 qsort(points, n, sizeof(double) * 4, cmp_axis_y);
01014         else if (axis == 2)
01015                 qsort(points, n, sizeof(double) * 4, cmp_axis_z);
01016         else
01017                 qsort(points, n, sizeof(double) * 4, cmp_val);
01018 }

void PointArray::transform const Transform3D transform  ) 
 

Definition at line 667 of file pointarray.cpp.

References n, points, and v.

00667                                                 {
00668 
00669         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00670                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00671                 v=v*xf;
00672                 points[i]  =v[0];
00673                 points[i+1]=v[1];
00674                 points[i+2]=v[2];
00675         }
00676 
00677 }

void PointArray::zero  ) 
 

Definition at line 116 of file pointarray.cpp.

References n, and points.

Referenced by set_from_density_map().

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


Member Data Documentation

size_t EMAN::PointArray::n [private]
 

Definition at line 157 of file pointarray.h.

Referenced by distmx(), get_number_points(), get_points(), match_points(), PointArray(), projection_by_nfft(), right_transform(), set_number_points(), sort_by_axis(), transform(), and zero().

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

Definition at line 156 of file pointarray.h.

Referenced by center_to_zero(), get_bounding_box(), get_center(), get_points(), get_points_array(), get_value_at(), get_vector_at(), mask(), mask_asymmetric_unit(), pdb2mrc_by_nfft(), pdb2mrc_by_summation(), PointArray(), projection_by_nfft(), projection_by_summation(), read_from_pdb(), replace_by_summation(), right_transform(), save_to_pdb(), set_from(), set_from_density_map(), set_number_points(), set_points_array(), set_vector_at(), sort_by_axis(), transform(), zero(), and ~PointArray().


The documentation for this class was generated from the following files:
Generated on Tue May 25 17:36:42 2010 for EMAN2 by  doxygen 1.4.4