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).
Transformalign_2d (PointArray *to, float max_dist)
 Aligns one PointArray to another in 2 dimensions.
vector< float > align_trans_2d (PointArray *to, int flags=0, float dxhint=0, float dyhint=0)
 Translationally aligns one PointArray to another in 2 dimensions.
void mask (double rmax, double rmin=0.0)
void mask_asymmetric_unit (const string &sym)
void transform (const Transform &transform)
void right_transform (const Transform &transform)
 Does Transform*v as opposed to v*Transform (as in the transform function).
void set_from (vector< float >)
void set_from (PointArray *source, const string &sym="", Transform *transform=0)
void set_from (double *source, int num, const string &sym="", Transform *transform=0)
void set_from_density_map (EMData *map, int num, float thresh, float apix, Density2PointsArrayAlgorithm mode=PEAKS_DIV)
void sort_by_axis (int axis=1)
EMDatapdb2mrc_by_nfft (int map_size, float apix, float res)
EMDatapdb2mrc_by_summation (int map_size, float apix, float res, int addpdbbfactor)
EMDataprojection_by_nfft (int image_size, float apix, float res=0)
EMDataprojection_by_summation (int image_size, float apix, float res)
void replace_by_summation (EMData *image, int i, Vec3f vec, float amp, float apix, float res)
void opt_from_proj (const vector< EMData * > &proj, float pixres)
 Optimizes a pointarray based on a set of projection images (EMData objects) This is effectively a 3D reconstruction algorithm.

Private Attributes

double * points
size_t n
double * bfactor


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
 

Enumeration values:
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 bfactor, n, and points.

Referenced by copy().

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

PointArray::PointArray int  nn  )  [explicit]
 

Definition at line 102 of file pointarray.cpp.

References n, and points.

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

PointArray::~PointArray  ) 
 

Definition at line 108 of file pointarray.cpp.

References points.

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


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

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

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

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

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

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

void PointArray::center_to_zero  ) 
 

Definition at line 552 of file pointarray.cpp.

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

00553 {
00554         FloatPoint center = get_center();
00555         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00556                 points[i] -= center[0];
00557                 points[i + 1] -= center[1];
00558                 points[i + 2] -= center[2];
00559         }
00560 }

PointArray * PointArray::copy  ) 
 

Definition at line 122 of file pointarray.cpp.

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

Referenced by mask(), and mask_asymmetric_unit().

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

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

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

Region PointArray::get_bounding_box  ) 
 

Definition at line 562 of file pointarray.cpp.

References get_number_points(), and points.

00563 {
00564         double xmin, ymin, zmin;
00565         double xmax, ymax, zmax;
00566         xmin = xmax = points[0];
00567         ymin = ymax = points[1];
00568         zmin = zmax = points[2];
00569         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00570                 if (points[i] > xmax)
00571                         xmax = points[i];
00572                 if (points[i] < xmin)
00573                         xmin = points[i];
00574                 if (points[i + 1] > ymax)
00575                         ymax = points[i + 1];
00576                 if (points[i + 1] < ymin)
00577                         ymin = points[i + 1];
00578                 if (points[i + 2] > zmax)
00579                         zmax = points[i + 2];
00580                 if (points[i + 2] < zmin)
00581                         zmin = points[i + 2];
00582         }
00583         return Region(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin);
00584 }

FloatPoint PointArray::get_center  ) 
 

Definition at line 533 of file pointarray.cpp.

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

Referenced by center_to_zero().

00534 {
00535         double xc, yc, zc;
00536         xc = yc = zc = 0.0;
00537         double norm = 0.0;
00538         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00539                 xc += points[i] * points[i + 3];
00540                 yc += points[i + 1] * points[i + 3];
00541                 zc += points[i + 2] * points[i + 3];
00542                 norm += points[i + 3];
00543         }
00544         if (norm <= 0) {
00545                 fprintf(stderr, "Abnormal total value (%g) for PointArray, it should be >0\n", norm);
00546                 return FloatPoint(0, 0, 0);
00547         }
00548         else
00549                 return FloatPoint(xc / norm, yc / norm, zc / norm);
00550 }

size_t PointArray::get_number_points  )  const
 

Definition at line 141 of file pointarray.cpp.

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

00142 {
00143         return n;
00144 }

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

References points.

00674                                      {
00675 vector<float> ret;
00676 for (unsigned int i=0; i<n; i++) {
00677         ret.push_back((float)points[i*4]);
00678         ret.push_back((float)points[i*4+1]);
00679         ret.push_back((float)points[i*4+2]);
00680 }
00681 
00682 return ret;
00683 }

double * PointArray::get_points_array  ) 
 

Definition at line 155 of file pointarray.cpp.

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

00156 {
00157         return points;
00158 }

double PointArray::get_value_at int  i  ) 
 

Definition at line 1700 of file pointarray.cpp.

References points.

Referenced by align_trans_2d().

01701 {
01702 return points[i*4+3];
01703 }

Vec3f PointArray::get_vector_at int  i  ) 
 

Definition at line 1695 of file pointarray.cpp.

References points, and EMAN::Vec3f.

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

01696 {
01697 return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
01698 }

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

Definition at line 587 of file pointarray.cpp.

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

00588 {
00589         double rmax2 = rmax * rmax, rmin2 = rmin * rmin;
00590         PointArray *tmp = this->copy();
00591         double *tmp_points = tmp->get_points_array();
00592          int count = 0;
00593         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00594                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v =
00595                         tmp_points[i + 3];
00596                 double r2 = x * x + y * y + z * z;
00597                 if (r2 >= rmin2 && r2 <= rmax2) {
00598                         points[count * 4] = x;
00599                         points[count * 4 + 1] = y;
00600                         points[count * 4 + 2] = z;
00601                         points[count * 4 + 3] = v;
00602                         ++count;
00603                 }
00604         }
00605         set_number_points(count);
00606         if( tmp )
00607         {
00608                 delete tmp;
00609                 tmp = 0;
00610         }
00611 }

void PointArray::mask_asymmetric_unit const string &  sym  ) 
 

Definition at line 614 of file pointarray.cpp.

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

00615 {
00616         if (sym == "c1" || sym == "C1")
00617                 return;                                 // do nothing for C1 symmetry
00618         double alt0 = 0, alt1 = M_PI, alt2 = M_PI;
00619         double az0 = 0, az1 = M_PI;
00620         if (sym[0] == 'c' || sym[0] == 'C') {
00621                 int nsym = atoi(sym.c_str() + 1);
00622                 az1 = 2.0 * M_PI / nsym / 2.0;
00623         }
00624         else if (sym[0] == 'd' || sym[0] == 'D') {
00625                 int nsym = atoi(sym.c_str() + 1);
00626                 alt1 = M_PI / 2.0;
00627                 alt2 = alt1;
00628                 az1 = 2.0 * M_PI / nsym / 2.0;
00629         }
00630         else if (sym == "icos" || sym == "ICOS") {
00631                 alt1 = 0.652358139784368185995; // 5fold to 3fold
00632                 alt2 = 0.55357435889704525151;  // half of edge ie. 5fold to 2fold along the edge
00633                 az1 = 2.0 * M_PI / 5 / 2.0;
00634         }
00635         else {
00636                 LOGERR("PointArray::set_to_asymmetric_unit(): sym = %s is not implemented yet",
00637                            sym.c_str());
00638                 return;
00639         }
00640 #ifdef DEBUG
00641         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);
00642 #endif
00643 
00644         PointArray *tmp = this->copy();
00645         double *tmp_points = tmp->get_points_array();
00646          int count = 0;
00647         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00648                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v = tmp_points[i + 3];
00649                 double az = atan2(y, x);
00650                 double az_abs = fabs(az - az0);
00651                 if (az_abs < (az1 - az0)) {
00652                         double alt_max = alt1 + (alt2 - alt1) * az_abs / (az1 - az0);
00653                         double alt = acos(z / sqrt(x * x + y * y + z * z));
00654                         if (alt < alt_max && alt >= alt0) {
00655 #ifdef DEBUG
00656                                 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);
00657 #endif
00658                                 points[count * 4] = x;
00659                                 points[count * 4 + 1] = y;
00660                                 points[count * 4 + 2] = z;
00661                                 points[count * 4 + 3] = v;
00662                                 count++;
00663                         }
00664                 }
00665         }
00666         set_number_points(count);
00667         if( tmp )
00668         {
00669                 delete tmp;
00670                 tmp = 0;
00671         }
00672 }

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

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

PointArray & PointArray::operator= PointArray pa  ) 
 

Definition at line 132 of file pointarray.cpp.

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

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

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

References get_number_points(), and LOGWARN.

01672                                                                         {
01673 #ifdef OPTPP
01674         optdata=proj;
01675         optobj=this;
01676         optpixres=pixres;
01677 
01678         FDNLF1 nlf(get_number_points()*4,calc_opt_proj,init_opt_proj);
01679 //      NLF1 nlf(get_number_points()*4,init_opt_proj,calc_opt_projd);
01680         nlf.initFcn();
01681 
01682         OptCG opt(&nlf);
01683 //      OptQNewton opt(&nlf);
01684         opt.setMaxFeval(2000);
01685         opt.optimize();
01686         opt.printStatus("Done");
01687 #else
01688         (void)proj;             //suppress warning message
01689         (void)pixres;   //suppress warning message
01690         LOGWARN("OPT++ support not enabled.\n");
01691         return;
01692 #endif
01693 }

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

Definition at line 1302 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().

01303 {
01304 #if defined NFFT
01305         nfft_3D_plan my_plan;           // plan for the nfft
01306 
01308         nfft_3D_init(&my_plan, map_size, get_number_points());
01309 
01311         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01312                 // FFTW and nfft use row major array layout, EMAN uses column major
01313                 my_plan.v[3 * j + 2] = (fftw_real) (points[i] / (apix * map_size));
01314                 my_plan.v[3 * j + 1] = (fftw_real) (points[i + 1] / (apix * map_size));
01315                 my_plan.v[3 * j] = (fftw_real) (points[i + 2] / (apix * map_size));
01316                 my_plan.f[j].re = (fftw_real) (points[i + 3]);
01317                 my_plan.f[j].im = 0.0;
01318         }
01319 
01321         if (my_plan.nfft_flags & PRE_PSI) {
01322                 nfft_3D_precompute_psi(&my_plan);
01323         }
01324 
01325         // compute the uniform Fourier transform
01326         nfft_3D_transpose(&my_plan);
01327 
01328         // copy the Fourier transform to EMData data array
01329         EMData *fft = new EMData();
01330         fft->set_size(map_size + 2, map_size, map_size);
01331         fft->set_complex(true);
01332         fft->set_ri(true);
01333         fft->to_zero();
01334         float *data = fft->get_data();
01335         double norm = 1.0 / (map_size * map_size * map_size);
01336         for (int k = 0; k < map_size; k++) {
01337                 for (int j = 0; j < map_size; j++) {
01338                         for (int i = 0; i < map_size / 2; i++) {
01339                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01340                                         (float) (my_plan.
01341                                                          f_hat[k * map_size * map_size + j * map_size + i +
01342                                                                    map_size / 2].re) * norm;
01343                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01344                                         (float) (my_plan.
01345                                                          f_hat[k * map_size * map_size + j * map_size + i +
01346                                                                    map_size / 2].im) * norm * (-1.0);
01347                         }
01348                 }
01349         }
01351         nfft_3D_finalize(&my_plan);
01352 
01353         // low pass processor
01354         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01355         int index = 0;
01356         for (int k = 0; k < map_size; k++) {
01357                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01358                 for (int j = 0; j < map_size; j++) {
01359                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01360                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01361                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01362                                 data[index] *= val;
01363                                 data[index + 1] *= val;
01364                         }
01365                 }
01366         }
01367         fft->update();
01368         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
01369 
01370         fft->process_inplace("xform.phaseorigin.tocorner");     // move phase origin to center of image map_size, instead of at corner
01371         EMData *map = fft->do_ift();
01372         map->set_attr("apix_x", apix);
01373         map->set_attr("apix_y", apix);
01374         map->set_attr("apix_z", apix);
01375         map->set_attr("origin_x", -map_size/2*apix);
01376         map->set_attr("origin_y", -map_size/2*apix);
01377         map->set_attr("origin_z", -map_size/2*apix);
01378         if( fft )
01379         {
01380                 delete fft;
01381                 fft = 0;
01382         }
01383         return map;
01384 #elif defined NFFT2
01385         nfft_plan my_plan;                      // plan for the nfft
01386 
01388         nfft_init_3d(&my_plan, map_size, map_size, map_size, get_number_points());
01389 
01391         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01392                 // FFTW and nfft use row major array layout, EMAN uses column major
01393                 my_plan.x[3 * j + 2] = (double) (points[i] / (apix * map_size));
01394                 my_plan.x[3 * j + 1] = (double) (points[i + 1] / (apix * map_size));
01395                 my_plan.x[3 * j] = (double) (points[i + 2] / (apix * map_size));
01396                 my_plan.f[j][0] = (double) (points[i + 3]);
01397                 my_plan.f[j][1] = 0.0;
01398         }
01399 
01401         if (my_plan.nfft_flags & PRE_PSI) {
01402                 nfft_precompute_psi(&my_plan);
01403                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01404                         nfft_full_psi(&my_plan, pow(10, -10));
01405         }
01406 
01407         // compute the uniform Fourier transform
01408         nfft_transposed(&my_plan);
01409 
01410         // copy the Fourier transform to EMData data array
01411         EMData *fft = new EMData();
01412         fft->set_size(map_size + 2, map_size, map_size);
01413         fft->set_complex(true);
01414         fft->set_ri(true);
01415         fft->to_zero();
01416         float *data = fft->get_data();
01417         double norm = 1.0 / (map_size * map_size * map_size);
01418         for (int k = 0; k < map_size; k++) {
01419                 for (int j = 0; j < map_size; j++) {
01420                         for (int i = 0; i < map_size / 2; i++) {
01421                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01422                                         (float) (my_plan.
01423                                                          f_hat[k * map_size * map_size + j * map_size + i +
01424                                                                    map_size / 2][0]) * norm;
01425                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01426                                         (float) (my_plan.
01427                                                          f_hat[k * map_size * map_size + j * map_size + i +
01428                                                                    map_size / 2][1]) * norm;
01429                         }
01430                 }
01431         }
01433         nfft_finalize(&my_plan);
01434 
01435         // low pass processor
01436         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01437         int index = 0;
01438         for (int k = 0; k < map_size; k++) {
01439                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01440                 for (int j = 0; j < map_size; j++) {
01441                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01442                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01443                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01444                                 data[index] *= val;
01445                                 data[index + 1] *= val;
01446                         }
01447                 }
01448         }
01449         fft->update();
01450         //fft->process_inplace(".filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
01451 
01452         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image map_size, instead of at corner
01453         EMData *map = fft->do_ift();
01454         map->set_attr("apix_x", apix);
01455         map->set_attr("apix_y", apix);
01456         map->set_attr("apix_z", apix);
01457         map->set_attr("origin_x", -map_size / 2 * apix);
01458         map->set_attr("origin_y", -map_size / 2 * apix);
01459         map->set_attr("origin_z", -map_size / 2 * apix);
01460         if( fft )
01461         {
01462                 delete fft;
01463                 fft = 0;
01464         }
01465         return map;
01466 #else
01467         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01468         return 0;
01469 #endif
01470 }

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

Definition at line 1039 of file pointarray.cpp.

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

01040 {
01041 #ifdef DEBUG
01042         printf("PointArray::pdb2mrc_by_summation(): %lu points\tmapsize = %4d\tapix = %g\tres = %g\n",get_number_points(),map_size, apix, res);
01043 #endif
01044         //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);
01045 
01046         double min_table_val = 1e-7;
01047         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01048 
01049         double table_step_size = 0.001; // number of steps for each pixel
01050         double inv_table_step_size = 1.0 / table_step_size;
01051         
01052 //      sort_by_axis(2);                        // sort by Z-axis
01053 
01054         EMData *map = new EMData();
01055         map->set_size(map_size, map_size, map_size);
01056         map->to_zero();
01057         float *pd = map->get_data();
01058         
01059         vector<double> table;
01060         double gauss_real_width;
01061         int table_size;
01062         int gbox;
01063         
01064         
01065         for ( size_t s = 0; s < get_number_points(); ++s) {
01066                 double xc = points[4 * s] / apix + map_size / 2;
01067                 double yc = points[4 * s + 1] / apix + map_size / 2;
01068                 double zc = points[4 * s + 2] / apix + map_size / 2;
01069                 double fval = points[4 * s + 3];
01070                 
01071                 //printf("\n bfactor=%f",bfactor[s]);
01072                 
01073                 
01074                 
01075                 
01076                 
01077                 if(addpdbbfactor==-1){
01078                         gauss_real_width = res/M_PI;    // in Angstrom, res is in Angstrom
01079                 }else{
01080                         gauss_real_width = (bfactor[s])/(4*sqrt(2.0)*M_PI);     // in Angstrom, res is in Angstrom
01081                 }
01082                 
01083                 
01084                 table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01085                 table.resize(table_size);
01086                 for (int i = 0; i < table_size; i++){
01087                         table[i] = 0;
01088                 }
01089                 
01090                 for (int i = 0; i < table_size; i++) {                                          
01091                         double x = -i * table_step_size * apix / gauss_real_width;
01092                         if(addpdbbfactor==-1){
01093                                 table[i] =  exp(-x * x);        
01094                         }
01095                         else{
01096                         table[i] =  exp(-x * x)/sqrt(gauss_real_width * gauss_real_width * 2 * M_PI);
01097                         }
01098                 }
01099 
01100                 gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
01101                 if (gbox <= 0)
01102                         gbox = 1;
01103                 
01104                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01105                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01106                 int kmin = int (zc) - gbox, kmax = int (zc) + gbox;
01107                 if (imin < 0)
01108                         imin = 0;
01109                 if (jmin < 0)
01110                         jmin = 0;
01111                 if (kmin < 0)
01112                         kmin = 0;
01113                 if (imax > map_size)
01114                         imax = map_size;
01115                 if (jmax > map_size)
01116                         jmax = map_size;
01117                 if (kmax > map_size)
01118                         kmax = map_size;                
01119                 
01120                 for (int k = kmin; k < kmax; k++) {
01121                         size_t table_index_z = size_t (fabs(k - zc) * inv_table_step_size);
01122                         if ( table_index_z >= table.size() ) continue;
01123                         double zval = table[table_index_z];
01124                         size_t pd_index_z = k * map_size * map_size;
01125                         
01126                         for (int j = jmin; j < jmax; j++) {
01127                                 size_t table_index_y = size_t (fabs(j - yc) * inv_table_step_size);
01128                                 if ( table_index_y >= table.size() ) continue;
01129                                 double yval = table[table_index_y];
01130                                 size_t pd_index = pd_index_z + j * map_size + imin;
01131                                 for (int i = imin; i < imax; i++, pd_index++) {
01132                                         size_t table_index_x = size_t (fabs(i - xc) * inv_table_step_size);
01133                                         if ( table_index_x >= table.size() ) continue;
01134                                         double xval = table[table_index_x];
01135                                         pd[pd_index] += (float) (fval * zval * yval * xval);
01136                                 }
01137                         }
01138                 }
01139         }
01140 
01141         map->update();
01142         map->set_attr("apix_x", apix);
01143         map->set_attr("apix_y", apix);
01144         map->set_attr("apix_z", apix);
01145         map->set_attr("origin_x", -map_size/2*apix);
01146         map->set_attr("origin_y", -map_size/2*apix);
01147         map->set_attr("origin_z", -map_size/2*apix);
01148 
01149         return map;
01150 }

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

Definition at line 1472 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().

01473 {
01474 #if defined NFFT
01475         nfft_2D_plan my_plan;           // plan for the nfft
01476         int N[2], n[2];
01477         N[0] = image_size;
01478         n[0] = next_power_of_2(2 * image_size);
01479         N[1] = image_size;
01480         n[1] = next_power_of_2(2 * image_size);
01481 
01483         nfft_2D_init(&my_plan, image_size, get_number_points());
01484         //nfft_2D_init_specific(&my_plan, N, get_number_points(), n, 3,
01485         //                 PRE_PHI_HUT | PRE_PSI,
01486         //                 FFTW_ESTIMATE | FFTW_OUT_OF_PLACE);
01488         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01489                 // FFTW and nfft use row major array layout, EMAN uses column major
01490                 my_plan.v[2 * j + 1] = (fftw_real) (points[i] / (apix * image_size));
01491                 my_plan.v[2 * j] = (fftw_real) (points[i + 1] / (apix * image_size));
01492                 my_plan.f[j].re = (fftw_real) points[i + 3];
01493                 my_plan.f[j].im = 0.0;
01494         }
01495 
01497         if (my_plan.nfft_flags & PRE_PSI) {
01498                 nfft_2D_precompute_psi(&my_plan);
01499         }
01500 
01501         // compute the uniform Fourier transform
01502         nfft_2D_transpose(&my_plan);
01503 
01504         // copy the Fourier transform to EMData data array
01505         EMData *fft = new EMData();
01506         fft->set_size(image_size + 2, image_size, 1);
01507         fft->set_complex(true);
01508         fft->set_ri(true);
01509         fft->to_zero();
01510         float *data = fft->get_data();
01511         double norm = 1.0 / (image_size * image_size);
01512         for (int j = 0; j < image_size; j++) {
01513                 for (int i = 0; i < image_size / 2; i++) {
01514                         data[j * (image_size + 2) + 2 * i] =
01515                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].re) * norm;
01516                         data[j * (image_size + 2) + 2 * i + 1] =
01517                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].im) * norm * (-1.0);
01518                 }
01519         }
01521         nfft_2D_finalize(&my_plan);
01522 
01523         if (res > 0) {
01524                 // Gaussian low pass processor
01525                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01526                 int index = 0;
01527                 for (int j = 0; j < image_size; j++) {
01528                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01529                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01530                                 double val = exp(-(i * i + RY2) / sigma2);
01531                                 data[index] *= val;
01532                                 data[index + 1] *= val;
01533                         }
01534                 }
01535         }
01536         fft->update();
01537         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
01538 
01539         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01540 
01541         return fft;
01542 #elif defined NFFT2
01543         nfft_plan my_plan;                      // plan for the nfft
01544         int N[2], n[2];
01545         N[0] = image_size;
01546         n[0] = next_power_of_2(2 * image_size);
01547         N[1] = image_size;
01548         n[1] = next_power_of_2(2 * image_size);
01549 
01551         //nfft_init_2d(&my_plan,image_size,image_size,get_number_points());
01552         nfft_init_specific(&my_plan, 2, N, get_number_points(), n, 3,
01553                                            PRE_PHI_HUT | PRE_PSI |
01554                                            MALLOC_X | MALLOC_F_HAT | MALLOC_F, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01556         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01557                 // FFTW and nfft use row major array layout, EMAN uses column major
01558                 my_plan.x[2 * j + 1] = (double) (points[i] / (apix * image_size));
01559                 my_plan.x[2 * j] = (double) (points[i + 1] / (apix * image_size));
01560                 my_plan.f[j][0] = (double) points[i + 3];
01561                 my_plan.f[j][1] = 0.0;
01562         }
01563 
01565         if (my_plan.nfft_flags & PRE_PSI) {
01566                 nfft_precompute_psi(&my_plan);
01567                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01568                         nfft_full_psi(&my_plan, pow(10, -6));
01569         }
01570 
01571         // compute the uniform Fourier transform
01572         nfft_transposed(&my_plan);
01573 
01574         // copy the Fourier transform to EMData data array
01575         EMData *fft = new EMData();
01576         fft->set_size(image_size + 2, image_size, 1);
01577         fft->set_complex(true);
01578         fft->set_ri(true);
01579         fft->to_zero();
01580         float *data = fft->get_data();
01581         double norm = 1.0 / (image_size * image_size);
01582         for (int j = 0; j < image_size; j++) {
01583                 for (int i = 0; i < image_size / 2; i++) {
01584                         data[j * (image_size + 2) + 2 * i] =
01585                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][0]) * norm;
01586                         data[j * (image_size + 2) + 2 * i + 1] =
01587                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][1]) * norm;
01588                 }
01589         }
01591         nfft_finalize(&my_plan);
01592 
01593         if (res > 0) {
01594                 // Gaussian low pass process
01595                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01596                 int index = 0;
01597                 for (int j = 0; j < image_size; j++) {
01598                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01599                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01600                                 double val = exp(-(i * i + RY2) / sigma2);
01601                                 data[index] *= val;
01602                                 data[index + 1] *= val;
01603                         }
01604                 }
01605         }
01606         fft->update();
01607         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
01608 
01609         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01610 
01611         return fft;
01612 #else
01613         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01614         return 0;
01615 #endif
01616 }

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

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

01154 {
01155         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01156         //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);
01157 
01158         double min_table_val = 1e-7;
01159         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01160 
01161         //double table_step_size = 0.001;    // number of steps for x=[0,1] in exp(-x*x)
01162         //int table_size = int(max_table_x / table_step_size *1.25);
01163         //double* table = (double*)malloc(sizeof(double) * table_size);
01164         //for(int i=0; i<table_size; i++) table[i]=exp(-i*i*table_step_size*table_step_size);
01165 
01166         double table_step_size = 0.001; // number of steps for each pixel
01167         double inv_table_step_size = 1.0 / table_step_size;
01168         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01169         double *table = (double *) malloc(sizeof(double) * table_size);
01170         for (int i = 0; i < table_size; i++) {
01171                 double x = -i * table_step_size * apix / gauss_real_width;
01172                 table[i] = exp(-x * x);
01173         }
01174 
01175         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01176         if (gbox <= 0)
01177                 gbox = 1;
01178         EMData *proj = new EMData();
01179         proj->set_size(image_size, image_size, 1);
01180         proj->to_zero();
01181         float *pd = proj->get_data();
01182         for ( size_t s = 0; s < get_number_points(); ++s) {
01183                 double xc = points[4 * s] / apix + image_size / 2;
01184                 double yc = points[4 * s + 1] / apix + image_size / 2;
01185                 double fval = points[4 * s + 3];
01186                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01187                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01188                 if (imin < 0)
01189                         imin = 0;
01190                 if (jmin < 0)
01191                         jmin = 0;
01192                 if (imax > image_size)
01193                         imax = image_size;
01194                 if (jmax > image_size)
01195                         jmax = image_size;
01196 
01197                 for (int j = jmin; j < jmax; j++) {
01198                         //int table_index_y = int(fabs(j-yc)*apix/gauss_real_width/table_step_size);
01199                         int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01200                         double yval = table[table_index_y];
01201 #ifdef DEBUG
01202                         //double yval2 = exp( - (j-yc)*(j-yc)*apix*apix/(gauss_real_width*gauss_real_width));
01203                         //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);
01204 #endif
01205                         int pd_index = j * image_size + imin;
01206                         for (int i = imin; i < imax; i++, pd_index++) {
01207                                 //int table_index_x = int(fabs(i-xc)*apix/gauss_real_width/table_step_size);
01208                                 int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01209                                 double xval = table[table_index_x];
01210 #ifdef DEBUG
01211                                 //double xval2 = exp( - (i-xc)*(i-xc)*apix*apix/(gauss_real_width*gauss_real_width));
01212                                 //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);
01213 #endif
01214                                 pd[pd_index] += (float)(fval * yval * xval);
01215                         }
01216                 }
01217         }
01218         for (int i = 0; i < image_size * image_size; i++)
01219                 pd[i] /= sqrt(M_PI);
01220         proj->update();
01221         return proj;
01222 }

bool PointArray::read_from_pdb const char *  file  ) 
 

Definition at line 422 of file pointarray.cpp.

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

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

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

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

Definition at line 1224 of file pointarray.cpp.

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

01225 {
01226         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01227 
01228         double min_table_val = 1e-7;
01229         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01230 
01231         double table_step_size = 0.001; // number of steps for each pixel
01232         double inv_table_step_size = 1.0 / table_step_size;
01233         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01234         double *table = (double *) malloc(sizeof(double) * table_size);
01235         for (int i = 0; i < table_size; i++) {
01236                 double x = -i * table_step_size * apix / gauss_real_width;
01237                 table[i] = exp(-x * x)/pow((float)M_PI,.25f);
01238         }
01239         int image_size=proj->get_xsize();
01240 
01241         // subtract the old point
01242         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01243         if (gbox <= 0)
01244                 gbox = 1;
01245         float *pd = proj->get_data();
01246          int s = ind;
01247         double xc = points[4 * s] / apix + image_size / 2;
01248         double yc = points[4 * s + 1] / apix + image_size / 2;
01249         double fval = points[4 * s + 3];
01250         int imin = int (xc) - gbox, imax = int (xc) + gbox;
01251         int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01252 
01253         if (imin < 0) imin = 0;
01254         if (jmin < 0) jmin = 0;
01255         if (imax > image_size) imax = image_size;
01256         if (jmax > image_size) jmax = image_size;
01257 
01258         for (int j = jmin; j < jmax; j++) {
01259                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01260                 double yval = table[table_index_y];
01261                 int pd_index = j * image_size + imin;
01262                 for (int i = imin; i < imax; i++, pd_index++) {
01263                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01264                         double xval = table[table_index_x];
01265                         pd[pd_index] -= (float)(fval * yval * xval);
01266                 }
01267         }
01268 
01269         // add the new point
01270         gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
01271         if (gbox <= 0)
01272                 gbox = 1;
01273         pd = proj->get_data();
01274         s = ind;
01275         xc = vec[0] / apix + image_size / 2;
01276         yc = vec[1] / apix + image_size / 2;
01277         fval = amp;
01278         imin = int (xc) - gbox, imax = int (xc) + gbox;
01279         jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01280 
01281         if (imin < 0) imin = 0;
01282         if (jmin < 0) jmin = 0;
01283         if (imax > image_size) imax = image_size;
01284         if (jmax > image_size) jmax = image_size;
01285 
01286         for (int j = jmin; j < jmax; j++) {
01287                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01288                 double yval = table[table_index_y];
01289                 int pd_index = j * image_size + imin;
01290                 for (int i = imin; i < imax; i++, pd_index++) {
01291                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01292                         double xval = table[table_index_x];
01293                         pd[pd_index] -= (float)(fval * yval * xval);
01294                 }
01295         }
01296 
01297         proj->update();
01298         return;
01299 }

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

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

00697                                                            {
00698         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00699                 Transform s = transform.transpose();
00700                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00701                 v= s*v;
00702                 points[i]  =v[0];
00703                 points[i+1]=v[1];
00704                 points[i+2]=v[2];
00705         }
00706 
00707 }

void PointArray::save_to_pdb const char *  file  ) 
 

Definition at line 522 of file pointarray.cpp.

References get_number_points(), and points.

00523 {
00524         FILE *fp = fopen(file, "w");
00525         for ( size_t i = 0; i < get_number_points(); i++) {
00526                 fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
00527                                 points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
00528         }
00529         fclose(fp);
00530 }

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

Definition at line 714 of file pointarray.cpp.

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

00715 {
00716          int nsym = xform->get_nsym(sym);
00717 
00718         if (get_number_points() != (size_t)nsym * num)
00719                 set_number_points((size_t)nsym * num);
00720 
00721         double *target = get_points_array();
00722 
00723         for ( int s = 0; s < nsym; s++) {
00724                 int index = s * 4 * num;
00725                 for ( int i = 0; i < 4 * num; i += 4, index += 4) {
00726                         Vec3f v((float)src[i],(float)src[i+1],(float)src[i+2]);
00727                         v=v*xform->get_sym(sym,s);
00728                         target[index]  =v[0];
00729                         target[index+1]=v[1];
00730                         target[index+2]=v[2];
00731                         target[index+3]=src[i+3];
00732                 }
00733         }
00734 }

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

Definition at line 708 of file pointarray.cpp.

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

00709 {
00710         set_from(source->get_points_array(), source->get_number_points(), sym, transform);
00711 
00712 }

void PointArray::set_from vector< float >   ) 
 

Definition at line 736 of file pointarray.cpp.

References points, and set_number_points().

Referenced by set_from().

00736                                            {
00737         set_number_points(pts.size()/4);
00738         for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];
00739 
00740 }

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

Definition at line 742 of file pointarray.cpp.

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

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

void PointArray::set_number_points size_t  nn  ) 
 

Definition at line 146 of file pointarray.cpp.

References bfactor, n, and points.

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

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

void PointArray::set_points_array double *  p  ) 
 

Definition at line 160 of file pointarray.cpp.

References points.

Referenced by set_from_density_map().

00161 {
00162         points = p;
00163 }

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

Definition at line 1713 of file pointarray.cpp.

References points, and v.

01714 {
01715 points[i*4]  =v[0];
01716 points[i*4+1]=v[1];
01717 points[i*4+2]=v[2];
01718 if (v.size()>=4) points[i*4+3]=v[3];
01719 }

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

Definition at line 1705 of file pointarray.cpp.

References points, and EMAN::Vec3f.

01706 {
01707 points[i*4]=vec[0];
01708 points[i*4+1]=vec[1];
01709 points[i*4+2]=vec[2];
01710 points[i*4+3]=value;
01711 }

void PointArray::sort_by_axis int  axis = 1  ) 
 

Definition at line 1026 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().

01027 {
01028         if (axis == 0)
01029                 qsort(points, n, sizeof(double) * 4, cmp_axis_x);
01030         else if (axis == 1)
01031                 qsort(points, n, sizeof(double) * 4, cmp_axis_y);
01032         else if (axis == 2)
01033                 qsort(points, n, sizeof(double) * 4, cmp_axis_z);
01034         else
01035                 qsort(points, n, sizeof(double) * 4, cmp_val);
01036 }

void PointArray::transform const Transform transform  ) 
 

Definition at line 685 of file pointarray.cpp.

References points, v, and EMAN::Vec3f.

00685                                               {
00686 
00687         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00688                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00689                 v=v*xf;
00690                 points[i]  =v[0];
00691                 points[i+1]=v[1];
00692                 points[i+2]=v[2];
00693         }
00694 
00695 }

void PointArray::zero  ) 
 

Definition at line 117 of file pointarray.cpp.

References n, and points.

Referenced by set_from_density_map().

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


Member Data Documentation

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

Definition at line 158 of file pointarray.h.

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

size_t EMAN::PointArray::n [private]
 

Definition at line 157 of file pointarray.h.

Referenced by distmx(), match_points(), PointArray(), projection_by_nfft(), set_number_points(), sort_by_axis(), 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_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 Jun 11 13:48:14 2013 for EMAN2 by  doxygen 1.3.9.1