pointarray.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Wen Jiang, 08/11/2004 (jiang12@purdue.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "pointarray.h"
00037 #include <vector>
00038 #include <cstring>
00039 
00040 using namespace EMAN;
00041 
00042 int cmp_axis_x(const void *a, const void *b)
00043 {
00044         double diff = ((double *) a)[0] - ((double *) b)[0];
00045         if (diff < 0.0)
00046                 return -1;
00047         else if (diff > 0.0)
00048                 return 1;
00049         else
00050                 return 0;
00051 }
00052 int cmp_axis_y(const void *a, const void *b)
00053 {
00054         double diff = ((double *) a)[1] - ((double *) b)[1];
00055         if (diff < 0.0)
00056                 return -1;
00057         else if (diff > 0.0)
00058                 return 1;
00059         else
00060                 return 0;
00061 }
00062 int cmp_axis_z(const void *a, const void *b)
00063 {
00064         double diff = ((double *) a)[2] - ((double *) b)[2];
00065         if (diff < 0.0)
00066                 return -1;
00067         else if (diff > 0.0)
00068                 return 1;
00069         else
00070                 return 0;
00071 }
00072 int cmp_val(const void *a, const void *b)
00073 {
00074         double diff = ((double *) a)[3] - ((double *) b)[3];
00075         if (diff < 0.0)
00076                 return -1;
00077         else if (diff > 0.0)
00078                 return 1;
00079         else
00080                 return 0;
00081 }
00082 // This will do a sort in descending order
00083 int cmp_float(const void *a, const void *b)
00084 {
00085         double diff = *((float *) a) - *((float *) b);
00086         if (diff < 0.0)
00087                 return 1;
00088         else if (diff > 0.0)
00089                 return -1;
00090         else
00091                 return 0;
00092 }
00093 
00094 
00095 PointArray::PointArray()
00096 {
00097         points = 0;
00098         n = 0;
00099 }
00100 
00101 PointArray::PointArray( int nn)
00102 {
00103         n = nn;
00104         points = (double *) calloc(4 * n, sizeof(double));
00105 }
00106 
00107 PointArray::~PointArray()
00108 {
00109         if( points )
00110         {
00111                 free(points);
00112                 points = 0;
00113         }
00114 }
00115 
00116 void PointArray::zero()
00117 {
00118         memset((void *) points, 0, 4 * n * sizeof(double));
00119 }
00120 
00121 PointArray *PointArray::copy()
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 }
00130 
00131 PointArray & PointArray::operator=(PointArray & pa)
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 }
00139 
00140  size_t PointArray::get_number_points() const
00141 {
00142         return n;
00143 }
00144 
00145 void PointArray::set_number_points(size_t nn)
00146 {
00147         if (n != nn) {
00148                 n = nn;
00149                 points = (double *) realloc(points, 4 * n * sizeof(double));
00150         }
00151 }
00152 
00153 double *PointArray::get_points_array()
00154 {
00155         return points;
00156 }
00157 
00158 void PointArray::set_points_array(double *p)
00159 {
00160         points = p;
00161 }
00162 
00163 EMData *PointArray::distmx(int sortrows) {
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 }
00187 
00188 vector<int> PointArray::match_points(PointArray *to,float max_miss) {
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 }
00249 
00250 // uses bilinear least-squares to generate a transformation
00251 // matrix for pairs of points
00252 Transform3D *PointArray::align_2d(PointArray *to,float max_dist) {
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 }
00288 
00289 vector<float> PointArray::align_trans_2d(PointArray *to, int flags, float dxhint,float dyhint) {
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 }
00407 
00408 
00409 bool PointArray::read_from_pdb(const char *file)
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 }
00502 
00503 
00504 void PointArray::save_to_pdb(const char *file)
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 }
00513 
00514 
00515 FloatPoint PointArray::get_center()
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 }
00533 
00534 void PointArray::center_to_zero()
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 }
00543 
00544 Region PointArray::get_bounding_box()
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 }
00567 
00568 
00569 void PointArray::mask(double rmax, double rmin)
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 }
00594 
00595 
00596 void PointArray::mask_asymmetric_unit(const string & sym)
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 }
00655 
00656 vector<float> PointArray::get_points() {
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 }
00666 
00667 void PointArray::transform(const Transform3D& xf) {
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 }
00678 
00679 void PointArray::right_transform(const Transform& transform) {
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 }
00690 void PointArray::set_from(PointArray * source, const string & sym, Transform3D *transform)
00691 {
00692         set_from(source->get_points_array(), source->get_number_points(), sym, transform);
00693 
00694 }
00695 
00696 void PointArray::set_from(double *src,  int num, const string & sym, Transform3D *xform)
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 }
00717 
00718 void PointArray::set_from(vector<float> pts) {
00719         set_number_points(pts.size()/4);
00720         for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];
00721 
00722 }
00723 
00724 void PointArray::set_from_density_map(EMData * map, int num, float thresh, float apix,
00725                                                                           Density2PointsArrayAlgorithm mode)
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 }
01004 
01005 
01006 
01007 
01008 void PointArray::sort_by_axis(int axis)
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 }
01019 
01020 
01021 EMData *PointArray::pdb2mrc_by_summation(int map_size, float apix, float res)
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 }
01104 
01105 
01106 EMData *PointArray::projection_by_summation(int image_size, float apix, float res)
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 }
01176 
01177 void PointArray::replace_by_summation(EMData *proj, int ind, Vec3f vec, float amp, float apix, float res)
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 }
01253 
01254 
01255 EMData *PointArray::pdb2mrc_by_nfft(int , float , float )
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 }
01424 
01425 EMData *PointArray::projection_by_nfft(int , float , float )
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 }
01570 
01571 #ifdef OPTPP
01572 #include "NLF.h"
01573 #include "BoundConstraint.h"
01574 #include "OptCG.h"
01575 //#include "OptNewton.h"
01576 #include "newmatap.h"
01577 
01578 vector<EMData*> optdata;
01579 PointArray *optobj;
01580 float optpixres;
01581 
01582 void init_opt_proj(int ndim, ColumnVector& x)
01583 {
01584 int i;
01585 double *data=optobj->get_points_array();
01586 
01587 for (i=0; i<ndim; i++) x(i+1)=data[i];
01588 }
01589 
01590 void calc_opt_proj(int n, const ColumnVector& x, double& fx, int& result)
01591 {
01592         int i;
01593         PointArray pa;
01594         Transform3D xform;
01595         int size=optdata[0]->get_xsize();
01596         fx=0;
01597 
01598         for (i=0; i<optdata.size(); i++) {
01599                 xform=(optdata[i]->get_transform());
01600                 pa.set_from((double *)x.nric()+1,n/4,std::string("c1"),&xform);
01601                 EMData *p=pa.projection_by_summation(size,1.0,optpixres);
01602                 p->process_inplace("normalize.unitlen");
01603                 fx-=sqrt(p->cmp("dot",EMObject(optdata[i]),Dict()));
01604         }
01605 
01606         result=NLPFunction;
01607 
01608         printf("%g\t%1.1f %1.1f %1.1f %g\t%1.1f %1.1f %1.1f %g\t%1.1f %1.1f %1.1f %g\n",fx,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12));
01609 }
01610 
01611 void calc_opt_projd(int mode,int n, const ColumnVector& x, double& fx, ColumnVector& gx, int& result)
01612 {
01613 
01614 if (mode & NLPFunction) {
01615 
01616 }
01617 
01618 if (mode & NLPGradient) {
01619 
01620 }
01621 
01622 }
01623 #endif
01624 
01625 void PointArray::opt_from_proj(const vector<EMData*> & proj,float pixres) {
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 }
01647 
01648 Vec3f PointArray::get_vector_at(int i)
01649 {
01650 return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
01651 }
01652 
01653 double PointArray::get_value_at(int i)
01654 {
01655 return points[i*4+3];
01656 }
01657 
01658 void PointArray::set_vector_at(int i,Vec3f vec,double value)
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 }
01665 
01666 void PointArray::set_vector_at(int i,vector<double> v)
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 }

Generated on Tue May 25 17:13:33 2010 for EMAN2 by  doxygen 1.4.7