EMAN2
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 
00037 #include "pointarray.h"
00038 #include "util.h"
00039 #include "vec3.h"
00040 #include <vector>
00041 #include <cstring>
00042 #include <boost/random.hpp>
00043 #include <boost/random/normal_distribution.hpp>
00044 
00045 #ifdef __APPLE__
00046         typedef unsigned int uint;
00047 #endif  //__APPLE__
00048 
00049 #ifdef _WIN32
00050         typedef unsigned int uint;
00051 #endif  //_WIN32
00052 
00053 using namespace EMAN;
00054 
00055 int cmp_axis_x(const void *a, const void *b)
00056 {
00057         double diff = ((double *) a)[0] - ((double *) b)[0];
00058         if (diff < 0.0)
00059                 return -1;
00060         else if (diff > 0.0)
00061                 return 1;
00062         else
00063                 return 0;
00064 }
00065 int cmp_axis_y(const void *a, const void *b)
00066 {
00067         double diff = ((double *) a)[1] - ((double *) b)[1];
00068         if (diff < 0.0)
00069                 return -1;
00070         else if (diff > 0.0)
00071                 return 1;
00072         else
00073                 return 0;
00074 }
00075 int cmp_axis_z(const void *a, const void *b)
00076 {
00077         double diff = ((double *) a)[2] - ((double *) b)[2];
00078         if (diff < 0.0)
00079                 return -1;
00080         else if (diff > 0.0)
00081                 return 1;
00082         else
00083                 return 0;
00084 }
00085 int cmp_val(const void *a, const void *b)
00086 {
00087         double diff = ((double *) a)[3] - ((double *) b)[3];
00088         if (diff < 0.0)
00089                 return -1;
00090         else if (diff > 0.0)
00091                 return 1;
00092         else
00093                 return 0;
00094 }
00095 // This will do a sort in descending order
00096 int cmp_float(const void *a, const void *b)
00097 {
00098         double diff = *((float *) a) - *((float *) b);
00099         if (diff < 0.0)
00100                 return 1;
00101         else if (diff > 0.0)
00102                 return -1;
00103         else
00104                 return 0;
00105 }
00106 
00107 
00108 PointArray::PointArray()
00109 {
00110         points = 0;
00111         bfactor = 0;
00112         n = 0;
00113         
00114         adist=0;
00115         aang=0;
00116         adihed=0;
00117         
00118         map=gradx=grady=gradz=0;
00119 }
00120 
00121 PointArray::PointArray( int nn)
00122 {
00123         n = nn;
00124         points = (double *) calloc(4 * n, sizeof(double));
00125         
00126         adist=0;
00127         aang=0;
00128         adihed=0;
00129         map=gradx=grady=gradz=0;
00130 }
00131 
00132 PointArray::~PointArray()
00133 {
00134         if( points )
00135         {
00136                 free(points);
00137                 points = 0;
00138         }
00139         
00140         if (adist) free(adist);
00141         if (aang) free(aang);
00142         if (adihed) free(adihed);
00143 //      if (map!=0) delete map;
00144         if (gradx!=0) delete gradx;
00145         if (grady!=0) delete grady;
00146         if (gradz!=0) delete gradz;
00147 }
00148 
00149 void PointArray::zero()
00150 {
00151         memset((void *) points, 0, 4 * n * sizeof(double));
00152 }
00153 
00154 PointArray *PointArray::copy()
00155 {
00156         PointArray *pa2 = new PointArray();
00157         pa2->set_number_points(get_number_points());
00158         double *pa2data = pa2->get_points_array();
00159         memcpy(pa2data, get_points_array(), sizeof(double) * 4 * get_number_points());
00160 
00161         return pa2;
00162 }
00163 
00164 PointArray & PointArray::operator=(PointArray & pa)
00165 {
00166         if (this != &pa) {
00167                 set_number_points(pa.get_number_points());
00168                 memcpy(get_points_array(), pa.get_points_array(), sizeof(double) * 4 * get_number_points());
00169         }
00170         return *this;
00171 }
00172 
00173  size_t PointArray::get_number_points() const
00174 {
00175         return n;
00176 }
00177 
00178 void PointArray::set_number_points(size_t nn)
00179 {       
00180         if (n != nn) {
00181                 n = nn;         
00182                 points = (double *) realloc(points, 4 * n * sizeof(double));
00183                 bfactor = (double *) realloc(bfactor, n * sizeof(double));
00184         }
00185 }
00186 
00187 double *PointArray::get_points_array()
00188 {
00189         return points;
00190 }
00191 
00192 void PointArray::set_points_array(double *p)
00193 {
00194         points = p;
00195 }
00196 
00197 EMData *PointArray::distmx(int sortrows) {
00198 if (n==0) return NULL;
00199 
00200 unsigned int i,j;
00201 
00202 EMData *ret= new EMData(n,n,1);
00203 ret->to_zero();
00204 
00205 for (i=0; i<n; i++) {
00206         for (j=i+1; j<n; j++) {
00207                 float r=(get_vector_at(i)-get_vector_at(j)).length();
00208                 ret->set_value_at(i,j,0,r);
00209                 ret->set_value_at(j,i,0,r);
00210         }
00211 }
00212 
00213 if (sortrows) {
00214         float *data=ret->get_data();
00215         for (i=0; i<n; i++) qsort(&data[i*n],n,sizeof(float),cmp_float);
00216         ret->update();
00217 }
00218 
00219 return ret;
00220 }
00221 
00222 vector<int> PointArray::match_points(PointArray *to,float max_miss) {
00223 EMData *mx0=distmx(1);
00224 EMData *mx1=to->distmx(1);
00225 unsigned int n2=mx1->get_xsize();       // same as get_number_points on to
00226 
00227 if (max_miss<0) max_miss=(float)mx0->get_attr("sigma")/10.0f;
00228 //printf("max error %f\n",max_miss);
00229 
00230 
00231 
00232 vector<int> ret(n,-1);
00233 vector<float> rete(n,0.0);
00234 unsigned int i,j,k,l;
00235 
00236 if (!mx0 || !mx1) {
00237         if (mx0) delete mx0;
00238         if (mx1) delete mx1;
00239         return ret;
00240 }
00241 
00242 // i iterates over elements of 'this', j looks for a match in 'to'
00243 // k and l iterate over the individual distances
00244 for (i=0; i<n; i++) {
00245         int bestn=-1;                   // number of best match in mx1
00246         double bestd=1.0e38;            // residual error distance in best match
00247         for (j=0; j<n2; j++) {
00248                 double d=0;
00249                 int nd=0;
00250                 for (k=l=0; k<n-1 && l<n2-1; k++,l++) {
00251                         float d1=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l,j));
00252                         float d2=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l,j));
00253                         float d3=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l+1,j));
00254                         float d4=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l+1,j));
00255                         if (d2<d1 && d4>d2) { l--; continue; }
00256                         if (d3<d1 && d4>d3) { k--; continue; }
00257                         d+=d1;
00258                         nd++;
00259                 }
00260                 d/=(float)nd;
00261 //              printf("%d -> %d\t%f\t%d\n",i,j,d,nd);
00262                 if (d<bestd) { bestd=d; bestn=j; }
00263         }
00264         ret[i]=bestn;
00265         rete[i]=static_cast<float>(bestd);
00266 }
00267 
00268 // This will remove duplicate assignments, keeping the best one
00269 // also remove any assignments with large errors
00270 for (i=0; i<n; i++) {
00271         for (j=0; j<n; j++) {
00272                 if (rete[i]>max_miss) { ret[i]=-1; break; }
00273                 if (i==j || ret[i]!=ret[j] || ret[i]==-1) continue;
00274                 if (rete[i]>rete[j]) { ret[i]=-1; break; }
00275         }
00276 }
00277 
00278 delete mx0;
00279 delete mx1;
00280 
00281 return ret;
00282 }
00283 
00284 // uses bilinear least-squares to generate a transformation
00285 // matrix for pairs of points
00286 Transform *PointArray::align_2d(PointArray *to,float max_dist) {
00287 vector<int> match=match_points(to,max_dist);
00288 Transform *ret=new Transform();
00289 
00290 // we use bilinear least squares to get 3/6 matrix components
00291 unsigned int i,j;
00292 
00293 vector<float> pts;
00294 for (i=0; i<match.size(); i++) {
00295         if (match[i]==-1) continue;
00296 
00297 //      printf("%d -> %d\n",i,match[i]);
00298         pts.push_back(get_vector_at(i)[0]);
00299         pts.push_back(get_vector_at(i)[1]);
00300         pts.push_back(to->get_vector_at(match[i])[0]);
00301 }
00302 
00303 Vec3f vx=Util::calc_bilinear_least_square(pts);
00304 
00305 // then we get the other 3/6
00306 for (i=j=0; i<match.size(); i++) {
00307         if (match[i]==-1) continue;
00308         pts[j*3]  =get_vector_at(i)[0];
00309         pts[j*3+1]=get_vector_at(i)[1];
00310         pts[j*3+2]=to->get_vector_at(match[i])[1];
00311         j++;
00312 }
00313 
00314 Vec3f vy=Util::calc_bilinear_least_square(pts);
00315 
00316 //ret->set_rotation(vx[1],vx[2],0.0,vy[1],vy[2],0.0,0.0,0.0,1.0);
00317 //ret->set_rotation(vx[1],vy[1],0.0,vx[2],vy[2],0.0,0.0,0.0,1.0);
00318 //ret->set_pretrans(Vec3f(-vx[0],-vy[0],0));
00319 
00320 ret->set(0, 0, vx[1]);
00321 ret->set(0, 1, vy[1]);
00322 ret->set(0, 2, 0.0f);
00323 ret->set(1, 0, vx[2]);
00324 ret->set(1, 1, vy[2]);
00325 ret->set(1, 2, 0.0f);
00326 ret->set(2, 0, 0.0f);
00327 ret->set(2, 1, 0.0f);
00328 ret->set(2, 2, 1.0f);
00329 ret->set_pre_trans(Vec3f(-vx[0],-vy[0],0));
00330 
00331 return ret;
00332 }
00333 
00334 vector<float> PointArray::align_trans_2d(PointArray *to, int flags, float dxhint,float dyhint) {
00335 printf("Warning, this is old code. Use align_2d.\n");
00336 
00337 // returns (dx,dy,residual error,n points used)
00338 // dxhint,dyhint should translate this->to
00339 // flags : 1 - use hint values, 2 - center by strongest point (highest 'value')
00340 int na=   get_number_points();
00341 int nb=to->get_number_points();
00342 if (na<=0 || nb<=0) return vector<float>(4,0);
00343 
00344 int *a2b = (int *)malloc(na*sizeof(int));
00345 int *b2a = (int *)malloc(nb*sizeof(int));
00346 
00347 
00348 // find unweighted centers
00349 float cax,cay,cbx,cby;
00350 int i,j;
00351 
00352 if (flags&1) {
00353         cbx=dxhint;
00354         cby=dyhint;
00355         cax=cay=0;
00356 }
00357 else if (flags&2) {
00358         // find the 'a' list peak
00359         float hia=0.0f;
00360         int hina=0;
00361         for (i=0; i<na; i++) {
00362                 if (get_value_at(i)>hia) { hia=static_cast<float>(get_value_at(i)); hina=i; }
00363         }
00364         cax=get_vector_at(hina)[0];
00365         cay=get_vector_at(hina)[1];
00366 
00367         // find the 'b' list peak
00368         float hib=0;
00369         int hinb=0;
00370         for (i=0; i<na; i++) {
00371                 if (to->get_value_at(i)>hib) { hib=static_cast<float>(to->get_value_at(i)); hinb=i; }
00372         }
00373         cbx=to->get_vector_at(hinb)[0];
00374         cby=to->get_vector_at(hinb)[1];
00375 }
00376 else {
00377         cax=cay=cbx=cby=0;
00378 
00379         for (i=0; i<na; i++) { cax+=get_vector_at(i)[0]; cay+=get_vector_at(i)[1]; }
00380         cax/=(float)na;
00381         cay/=(float)na;
00382 
00383         for (i=0; i<nb; i++) { cbx+=to->get_vector_at(i)[0]; cby+=to->get_vector_at(i)[1]; }
00384         cbx/=(float)nb;
00385         cby/=(float)nb;
00386 }
00387 
00388 Vec3f offset(cbx-cax,cby-cay,0);
00389 
00390 // find the nearest point for each x point, taking the estimated centers into account
00391 for (i=0; i<na; i++) {
00392         float rmin=1.0e30f;
00393         for (j=0; j<nb; j++) {
00394                 float r=(get_vector_at(i)+offset-to->get_vector_at(j)).length();
00395                 if (r<rmin) { a2b[i]=j; rmin=r; }
00396         }
00397 }
00398 
00399 // find the nearest point for each y point
00400 for (i=0; i<nb; i++) {
00401         float rmin=1.0e30f;
00402         for (j=0; j<na; j++) {
00403                 float r=(get_vector_at(j)+offset-to->get_vector_at(i)).length();
00404                 if (r<rmin) { b2a[i]=j; rmin=r; }
00405         }
00406 }
00407 
00408 // now keep only points where x->y matches y->x
00409 for (i=0; i<na; i++) {
00410         if (a2b[i]<0) continue;
00411         if (b2a[a2b[i]]!=i) { printf(" #%d!=%d# ",b2a[a2b[i]],i);  b2a[a2b[i]]=-1; a2b[i]=-1; }
00412         printf("%d->%d  ",i,a2b[i]);
00413 }
00414 printf("\n");
00415 
00416 for (i=0; i<nb; i++) {
00417         if (b2a[i]<0) continue;
00418         if (a2b[b2a[i]]!=i) { a2b[b2a[i]]=-1; b2a[i]=-1; }
00419         printf("%d->%d  ",i,b2a[i]);
00420 }
00421 printf("\n");
00422 
00423 // Compute the average translation required to align the points
00424 float dx=0,dy=0,dr=0,nr=0;
00425 for (i=0; i<na; i++) {
00426         if (a2b[i]==-1) continue;
00427         dx+=to->get_vector_at(a2b[i])[0]-get_vector_at(i)[0];
00428         dy+=to->get_vector_at(a2b[i])[1]-get_vector_at(i)[1];
00429         nr+=1.0;
00430 }
00431 //printf("%f %f %f\n",dx,dy,nr);
00432 if (nr<2) return vector<float>(4,0);
00433 dx/=nr;
00434 dy/=nr;
00435 
00436 // Compute the residual error
00437 for (i=0; i<na; i++) {
00438         if (i==-1  || a2b[i]==-1) continue;
00439         dr+=(to->get_vector_at(a2b[i])-get_vector_at(i)-Vec3f(dx,dy,0)).length();
00440 }
00441 dr/=nr;
00442 
00443 free(a2b);
00444 free(b2a);
00445 vector<float> ret(4);
00446 ret[0]=dx;
00447 ret[1]=dy;
00448 ret[2]=dr;
00449 ret[3]=(float)nr;
00450 return ret;
00451 }
00452 
00453 
00454 bool PointArray::read_from_pdb(const char *file)
00455 {
00456         struct stat filestat;
00457         stat(file, &filestat);
00458         set_number_points(( int)(filestat.st_size / 80 + 1));
00459 
00460         #ifdef DEBUG
00461         printf("PointArray::read_from_pdb(): try %4lu atoms first\n", get_number_points());
00462         #endif
00463 
00464         FILE *fp = fopen(file, "r");
00465         if(!fp) {
00466                 fprintf(stderr,"ERROR in PointArray::read_from_pdb(): cannot open file %s\n",file);
00467                 throw;
00468         }
00469         char s[200];
00470         size_t count = 0;
00471         
00472         while ((fgets(s, 200, fp) != NULL)) {
00473                 if (strncmp(s, "ENDMDL", 6) == 0)
00474                         break;
00475                 if (strncmp(s, "ATOM", 4) != 0)
00476                         continue;
00477 
00478                 if (s[13] == ' ')
00479                         s[13] = s[14];
00480                 if (s[13] == ' ')
00481                         s[13] = s[15];
00482 
00483                 float e = 0;
00484                 char ctt, ctt2 = ' ';
00485                 if (s[13] == ' ')
00486                         ctt = s[14];
00487                 else if (s[12] == ' ') {
00488                         ctt = s[13];
00489                         ctt2 = s[14];
00490                 }
00491                 else {
00492                         ctt = s[12];
00493                         ctt2 = s[13];
00494                 }
00495 
00496                 switch (ctt) {
00497                 case 'H':
00498                         e = 1.0;
00499                         break;
00500                 case 'C':
00501                         e = 6.0;
00502                         break;
00503                 case 'A':
00504                         if (ctt2 == 'U') {
00505                                 e = 79.0;
00506                                 break;
00507                         }
00508                         // treat 'A'mbiguous atoms as N, not perfect, but good enough
00509                 case 'N':
00510                         e = 7.0;
00511                         break;
00512                 case 'O':
00513                         e = 8.0;
00514                         break;
00515                 case 'P':
00516                         e = 15.0;
00517                         break;
00518                 case 'S':
00519                         e = 16.0;
00520                         break;
00521                 case 'W':
00522                         e = 18.0;
00523                         break;                          // ficticious water 'atom'
00524                 default:
00525                         fprintf(stderr, "Unknown atom %c%c\n", ctt, ctt2);
00526                         e = 0;
00527                 }
00528                 if (e == 0)
00529                         continue;
00530 
00531                 float x, y, z, q;
00532                 sscanf(&s[28], " %f %f %f", &x, &y, &z);
00533                 sscanf(&s[60], " %f", &q);
00534 
00535                 if (count + 1 > get_number_points())
00536                         set_number_points(2 * (count + 1));    //makes sure point array is big enough
00537                 
00538 #ifdef DEBUG
00539                 printf("Atom %4lu: x,y,z = %8g,%8g,%8g\te = %g\n", count, x, y, z, e);
00540 #endif
00541                 points[4 * count] = x;
00542                 points[4 * count + 1] = y;
00543                 points[4 * count + 2] = z;
00544                 points[4 * count + 3] = e;
00545                 bfactor[count] = q;
00546                 count++;
00547         }
00548         fclose(fp);
00549         set_number_points(count);
00550         return true;
00551 }
00552 
00553 
00554 void PointArray::save_to_pdb(const char *file)
00555 {
00556         FILE *fp = fopen(file, "w");
00557         for ( size_t i = 0; i < get_number_points(); i++) {
00558                 fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
00559                                 points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
00560         }
00561         fclose(fp);
00562 }
00563 
00564 
00565 FloatPoint PointArray::get_center()
00566 {
00567         double xc, yc, zc;
00568         xc = yc = zc = 0.0;
00569         double norm = 0.0;
00570         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00571                 xc += points[i] * points[i + 3];
00572                 yc += points[i + 1] * points[i + 3];
00573                 zc += points[i + 2] * points[i + 3];
00574                 norm += points[i + 3];
00575         }
00576         if (norm <= 0) {
00577                 fprintf(stderr, "Abnormal total value (%g) for PointArray, it should be >0\n", norm);
00578                 return FloatPoint(0, 0, 0);
00579         }
00580         else
00581                 return FloatPoint(xc / norm, yc / norm, zc / norm);
00582 }
00583 
00584 void PointArray::center_to_zero()
00585 {
00586         FloatPoint center = get_center();
00587         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00588                 points[i] -= center[0];
00589                 points[i + 1] -= center[1];
00590                 points[i + 2] -= center[2];
00591         }
00592 }
00593 
00594 Region PointArray::get_bounding_box()
00595 {
00596         double xmin, ymin, zmin;
00597         double xmax, ymax, zmax;
00598         xmin = xmax = points[0];
00599         ymin = ymax = points[1];
00600         zmin = zmax = points[2];
00601         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00602                 if (points[i] > xmax)
00603                         xmax = points[i];
00604                 if (points[i] < xmin)
00605                         xmin = points[i];
00606                 if (points[i + 1] > ymax)
00607                         ymax = points[i + 1];
00608                 if (points[i + 1] < ymin)
00609                         ymin = points[i + 1];
00610                 if (points[i + 2] > zmax)
00611                         zmax = points[i + 2];
00612                 if (points[i + 2] < zmin)
00613                         zmin = points[i + 2];
00614         }
00615         return Region(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin);
00616 }
00617 
00618 
00619 void PointArray::mask(double rmax, double rmin)
00620 {
00621         double rmax2 = rmax * rmax, rmin2 = rmin * rmin;
00622         PointArray *tmp = this->copy();
00623         double *tmp_points = tmp->get_points_array();
00624          int count = 0;
00625         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00626                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v =
00627                         tmp_points[i + 3];
00628                 double r2 = x * x + y * y + z * z;
00629                 if (r2 >= rmin2 && r2 <= rmax2) {
00630                         points[count * 4] = x;
00631                         points[count * 4 + 1] = y;
00632                         points[count * 4 + 2] = z;
00633                         points[count * 4 + 3] = v;
00634                         ++count;
00635                 }
00636         }
00637         set_number_points(count);
00638         if( tmp )
00639         {
00640                 delete tmp;
00641                 tmp = 0;
00642         }
00643 }
00644 
00645 
00646 void PointArray::mask_asymmetric_unit(const string & sym)
00647 {
00648         if (sym == "c1" || sym == "C1")
00649                 return;                                 // do nothing for C1 symmetry
00650         double alt0 = 0, alt1 = M_PI, alt2 = M_PI;
00651         double az0 = 0, az1 = M_PI;
00652         if (sym[0] == 'c' || sym[0] == 'C') {
00653                 int nsym = atoi(sym.c_str() + 1);
00654                 az1 = 2.0 * M_PI / nsym / 2.0;
00655         }
00656         else if (sym[0] == 'd' || sym[0] == 'D') {
00657                 int nsym = atoi(sym.c_str() + 1);
00658                 alt1 = M_PI / 2.0;
00659                 alt2 = alt1;
00660                 az1 = 2.0 * M_PI / nsym / 2.0;
00661         }
00662         else if (sym == "icos" || sym == "ICOS") {
00663                 alt1 = 0.652358139784368185995; // 5fold to 3fold
00664                 alt2 = 0.55357435889704525151;  // half of edge ie. 5fold to 2fold along the edge
00665                 az1 = 2.0 * M_PI / 5 / 2.0;
00666         }
00667         else {
00668                 LOGERR("PointArray::set_to_asymmetric_unit(): sym = %s is not implemented yet",
00669                            sym.c_str());
00670                 return;
00671         }
00672 #ifdef DEBUG
00673         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);
00674 #endif
00675 
00676         PointArray *tmp = this->copy();
00677         double *tmp_points = tmp->get_points_array();
00678          int count = 0;
00679         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00680                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v = tmp_points[i + 3];
00681                 double az = atan2(y, x);
00682                 double az_abs = fabs(az - az0);
00683                 if (az_abs < (az1 - az0)) {
00684                         double alt_max = alt1 + (alt2 - alt1) * az_abs / (az1 - az0);
00685                         double alt = acos(z / sqrt(x * x + y * y + z * z));
00686                         if (alt < alt_max && alt >= alt0) {
00687 #ifdef DEBUG
00688                                 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);
00689 #endif
00690                                 points[count * 4] = x;
00691                                 points[count * 4 + 1] = y;
00692                                 points[count * 4 + 2] = z;
00693                                 points[count * 4 + 3] = v;
00694                                 count++;
00695                         }
00696                 }
00697         }
00698         set_number_points(count);
00699         if( tmp )
00700         {
00701                 delete tmp;
00702                 tmp = 0;
00703         }
00704 }
00705 
00706 vector<float> PointArray::get_points() {
00707 vector<float> ret;
00708 for (unsigned int i=0; i<n; i++) {
00709         ret.push_back((float)points[i*4]);
00710         ret.push_back((float)points[i*4+1]);
00711         ret.push_back((float)points[i*4+2]);
00712 }
00713 
00714 return ret;
00715 }
00716 
00717 void PointArray::transform(const Transform& xf) {
00718 
00719         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00720                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00721                 v=v*xf;
00722                 points[i]  =v[0];
00723                 points[i+1]=v[1];
00724                 points[i+2]=v[2];
00725         }
00726 
00727 }
00728 
00729 void PointArray::right_transform(const Transform& transform) {
00730         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00731                 Transform s = transform.transpose();
00732                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00733                 v= s*v;
00734                 points[i]  =v[0];
00735                 points[i+1]=v[1];
00736                 points[i+2]=v[2];
00737         }
00738 
00739 }
00740 void PointArray::set_from(PointArray * source, const string & sym, Transform *transform)
00741 {
00742         set_from(source->get_points_array(), source->get_number_points(), sym, transform);
00743 
00744 }
00745 
00746 void PointArray::set_from(double *src,  int num, const string & sym, Transform *xform)
00747 {
00748          int nsym = xform->get_nsym(sym);
00749 
00750         if (get_number_points() != (size_t)nsym * num)
00751                 set_number_points((size_t)nsym * num);
00752 
00753         double *target = get_points_array();
00754 
00755         for ( int s = 0; s < nsym; s++) {
00756                 int index = s * 4 * num;
00757                 for ( int i = 0; i < 4 * num; i += 4, index += 4) {
00758                         Vec3f v((float)src[i],(float)src[i+1],(float)src[i+2]);
00759                         v=v*xform->get_sym(sym,s);
00760                         target[index]  =v[0];
00761                         target[index+1]=v[1];
00762                         target[index+2]=v[2];
00763                         target[index+3]=src[i+3];
00764                 }
00765         }
00766 }
00767 
00768 void PointArray::set_from(vector<float> pts) {
00769         set_number_points(pts.size()/4);
00770         for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];
00771 
00772 }
00773 
00774 void PointArray::set_from_density_map(EMData * map, int num, float thresh, float apix,
00775                                                                           Density2PointsArrayAlgorithm mode)
00776 {
00777         if (mode == PEAKS_SUB || mode == PEAKS_DIV) {
00778                 // find out how many voxels are useful voxels
00779                 int num_voxels = 0;
00780                 int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00781                 size_t size = (size_t)nx * ny * nz;
00782                 EMData *tmp_map = map->copy();
00783                 float *pd = tmp_map->get_data();
00784                 for (size_t i = 0; i < size; ++i) {
00785                         if (pd[i] > thresh)
00786                                 ++num_voxels;
00787                 }
00788 
00789                 double pointvol = double (num_voxels) / double (num);
00790                 double gauss_real_width = pow(pointvol, 1. / 3.);       // in pixels
00791 #ifdef DEBUG
00792                 printf("Average point range radius = %g pixels for %d points from %d used voxels\n",
00793                            gauss_real_width, num, num_voxels);
00794 #endif
00795 
00796                 double min_table_val = 1e-4;
00797                 double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
00798 
00799                 double table_step_size = 1.;    // number of steps for each pixel
00800                 double inv_table_step_size = 1.0 / table_step_size;
00801                 int table_size = int (max_table_x * gauss_real_width / (table_step_size) * 1.25) + 1;
00802                 double *table = (double *) malloc(sizeof(double) * table_size);
00803                 for (int i = 0; i < table_size; i++) {
00804                         double x = i * table_step_size / gauss_real_width;
00805                         table[i] = exp(-x * x);
00806                 }
00807 
00808                 int gbox = int (max_table_x * gauss_real_width);        // local box half size in pixels to consider for each point
00809                 if (gbox <= 0)
00810                         gbox = 1;
00811 
00812                 set_number_points(num);
00813                 for (int count = 0; count < num; count++) {
00814                         float cmax = pd[0];
00815                         int cmaxpos = 0;
00816                         for (size_t i = 0; i < size; ++i) {
00817                                 if (pd[i] > cmax) {
00818                                         cmax = pd[i];
00819                                         cmaxpos = i;
00820                                 }
00821                         }
00822                         int iz = cmaxpos / (nx * ny), iy = (cmaxpos - iz * nx * ny) / nx, ix =
00823                                 cmaxpos - iz * nx * ny - iy * nx;
00824 
00825                         // update coordinates in pixels
00826                         points[4*count  ] = ix;
00827                         points[4*count+1] = iy;
00828                         points[4*count+2] = iz;
00829                         points[4*count+3] = cmax;
00830 #ifdef DEBUG
00831                         printf("Point %d: val = %g\tat  %d, %d, %d\n", count, cmax, ix, iy, iz);
00832 #endif
00833 
00834                         int imin = ix - gbox, imax = ix + gbox;
00835                         int jmin = iy - gbox, jmax = iy + gbox;
00836                         int kmin = iz - gbox, kmax = iz + gbox;
00837                         if (imin < 0)
00838                                 imin = 0;
00839                         if (jmin < 0)
00840                                 jmin = 0;
00841                         if (kmin < 0)
00842                                 kmin = 0;
00843                         if (imax > nx)
00844                                 imax = nx;
00845                         if (jmax > ny)
00846                                 jmax = ny;
00847                         if (kmax > nz)
00848                                 kmax = nz;
00849 
00850                         for (int k = kmin; k < kmax; ++k) {
00851                                 int table_index_z = int (fabs(double (k - iz)) * inv_table_step_size);
00852                                 double zval = table[table_index_z];
00853                                 size_t pd_index_z = k * nx * ny;
00854                                 //printf("k = %8d\tx = %8g\tval = %8g\n", k, float(k-iz), zval);
00855                                 for (int j = jmin; j < jmax; ++j) {
00856                                         int table_index_y = int (fabs(double (j - iy)) * inv_table_step_size);
00857                                         double yval = table[table_index_y];
00858                                         size_t pd_index = pd_index_z + j * nx + imin;
00859                                         for (int i = imin; i < imax; ++i, ++pd_index) {
00860                                                 int table_index_x = int (fabs(double (i - ix)) * inv_table_step_size);
00861                                                 double xval = table[table_index_x];
00862                                                 if (mode == PEAKS_SUB)
00863                                                         pd[pd_index] -= (float)(cmax * zval * yval * xval);
00864                                                 else
00865                                                         pd[pd_index] *= (float)(1.0 - zval * yval * xval);      // mode == PEAKS_DIV
00866                                         }
00867                                 }
00868                         }
00869                 }
00870                 set_number_points(num);
00871                 tmp_map->update();
00872                 if( tmp_map )
00873                 {
00874                         delete tmp_map;
00875                         tmp_map = 0;
00876                 }
00877         }
00878         else if (mode == KMEANS) {
00879                 set_number_points(num);
00880                 zero();
00881 
00882                 PointArray tmp_pa;
00883                 tmp_pa.set_number_points(num);
00884                 tmp_pa.zero();
00885 
00886                  int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00887                 float *pd = map->get_data();
00888 
00889                 // initialize segments with random centers at pixels with values > thresh
00890 #ifdef DEBUG
00891                 printf("Start initial random seeding\n");
00892 #endif
00893                 for ( size_t i = 0; i < get_number_points(); i++) {
00894                          int x, y, z;
00895                         double v;
00896                         do {
00897                                 x = ( int) Util::get_frand(0, nx - 1);
00898                                 y = ( int) Util::get_frand(0, ny - 1);
00899                                 z = ( int) Util::get_frand(0, nz - 1);
00900                                 v = pd[z * nx * ny + y * nx + x];
00901 #ifdef DEBUG
00902                                 printf("Trying Point %lu: val = %g\tat  %d, %d, %d\tfrom map (%d,%d,%d)\n", i, v, x,
00903                                            y, z, nx, ny, nz);
00904 #endif
00905                         } while (v <= thresh);
00906                         points[4 * i] = (double) x;
00907                         points[4 * i + 1] = (double) y;
00908                         points[4 * i + 2] = (double) z;
00909                         points[4 * i + 3] = (double) v;
00910 #ifdef DEBUG
00911                         printf("Point %lu: val = %g\tat  %g, %g, %g\n", i, points[4 * i + 3], points[4 * i],
00912                                    points[4 * i + 1], points[4 * i + 2]);
00913 #endif
00914                 }
00915 
00916                 double min_dcen = 1e0;  // minimal mean segment center shift as convergence criterion
00917                 double dcen = 0.0;
00918                 int iter = 0, max_iter = 100;
00919                 do {
00920 #ifdef DEBUG
00921                         printf("Iteration %3d, start\n", iter);
00922 #endif
00923                         double *tmp_points = tmp_pa.get_points_array();
00924 
00925                         // reassign each pixel to the best segment
00926                         for ( int k = 0; k < nz; k++) {
00927                                 for ( int j = 0; j < ny; j++) {
00928                                         for ( int i = 0; i < nx; i++) {
00929                                                 size_t idx = k * nx * ny + j * nx + i;
00930                                                 if (pd[idx] > thresh) {
00931                                                         double min_dist = 1e60; // just a large distance
00932                                                          int min_s = 0;
00933                                                         for ( size_t s = 0; s < get_number_points(); ++s) {
00934                                                                 double x = points[4 * s];
00935                                                                 double y = points[4 * s + 1];
00936                                                                 double z = points[4 * s + 2];
00937                                                                 double dist =
00938                                                                         (k - z) * (k - z) + (j - y) * (j - y) + (i - x) * (i - x);
00939                                                                 if (dist < min_dist) {
00940                                                                         min_dist = dist;
00941                                                                         min_s = s;
00942                                                                 }
00943                                                         }
00944                                                         tmp_points[4 * min_s] += i;
00945                                                         tmp_points[4 * min_s + 1] += j;
00946                                                         tmp_points[4 * min_s + 2] += k;
00947                                                         tmp_points[4 * min_s + 3] += 1.0;
00948                                                 }
00949                                         }
00950                                 }
00951                         }
00952 #ifdef DEBUG
00953                         printf("Iteration %3d, finished reassigning segments\n", iter);
00954 #endif
00955                         // update each segment's center
00956                         dcen = 0.0;
00957                         for ( size_t s = 0; s < get_number_points(); ++s) {
00958                                 if (tmp_points[4 * s + 3]) {
00959                                         tmp_points[4 * s] /= tmp_points[4 * s + 3];
00960                                         tmp_points[4 * s + 1] /= tmp_points[4 * s + 3];
00961                                         tmp_points[4 * s + 2] /= tmp_points[4 * s + 3];
00962 #ifdef DEBUG
00963                                         printf("Iteration %3d, Point %3lu at %8g, %8g, %8g -> %8g, %8g, %8g\n", iter, s,
00964                                                    points[4 * s], points[4 * s + 1], points[4 * s + 2], tmp_points[4 * s],
00965                                                    tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00966 #endif
00967                                 }
00968                                 else {                  // empty segments are reseeded
00969                                          int x, y, z;
00970                                         double v;
00971                                         do {
00972                                                 x = ( int) Util::get_frand(0, nx - 1);
00973                                                 y = ( int) Util::get_frand(0, ny - 1);
00974                                                 z = ( int) Util::get_frand(0, nz - 1);
00975                                                 v = pd[z * nx * ny + y * nx + x];
00976                                         } while (v <= thresh);
00977                                         tmp_points[4 * s] = (double) x;
00978                                         tmp_points[4 * s + 1] = (double) y;
00979                                         tmp_points[4 * s + 2] = (double) z;
00980                                         tmp_points[4 * s + 3] = (double) v;
00981 #ifdef DEBUG
00982                                         printf
00983                                                 ("Iteration %3d, Point %3lu reseeded from %8g, %8g, %8g -> %8g, %8g, %8g\n",
00984                                                  iter, s, points[4 * s], points[4 * s + 1], points[4 * s + 2],
00985                                                  tmp_points[4 * s], tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00986 #endif
00987                                 }
00988                                 double dx = tmp_points[4 * s] - points[4 * s];
00989                                 double dy = tmp_points[4 * s + 1] - points[4 * s + 1];
00990                                 double dz = tmp_points[4 * s + 2] - points[4 * s + 2];
00991                                 dcen += dx * dx + dy * dy + dz * dz;
00992                         }
00993                         dcen = sqrt(dcen / get_number_points());
00994                         //swap pointter, faster but risky
00995 #ifdef DEBUG
00996                         printf("before swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
00997                                    (long)tmp_pa.get_points_array());
00998 #endif
00999                         double *tp = get_points_array();
01000                         set_points_array(tmp_points);
01001                         tmp_pa.set_points_array(tp);
01002                         tmp_pa.zero();
01003 #ifdef DEBUG
01004                         printf("after  swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
01005                                    (long)tmp_pa.get_points_array());
01006                         printf("Iteration %3d, finished updating segment centers with dcen = %g pixels\n", iter,
01007                                    dcen);
01008 #endif
01009 
01010                         iter++;
01011                 } while (dcen > min_dcen && iter <= max_iter);
01012                 map->update();
01013 
01014                 sort_by_axis(2);        // x,y,z axes = 0, 1, 2
01015         }
01016         else {
01017                 LOGERR("PointArray::set_from_density_map(): mode = %d is not implemented yet", mode);
01018         }
01019         //update to use apix and origin
01020         int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
01021         float origx, origy, origz;
01022         try {
01023                 origx = map->get_attr("origin_x");
01024                 origy = map->get_attr("origin_y");
01025                 origz = map->get_attr("origin_z");
01026         }
01027         catch(...) {
01028                 origx = -nx / 2 * apix;
01029                 origy = -ny / 2 * apix;
01030                 origz = -nz / 2 * apix;
01031         }
01032 
01033 #ifdef DEBUG
01034         printf("Apix = %g\torigin x,y,z = %8g,%8g,%8g\n",apix, origx, origy, origz);
01035 #endif
01036 
01037         float *pd = map->get_data();
01038         for ( size_t i = 0; i < get_number_points(); ++i) {
01039 #ifdef DEBUG
01040                 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]);
01041 #endif
01042                 points[4 * i + 3] =
01043                         pd[(int) points[4 * i + 2] * nx * ny + (int) points[4 * i + 1] * nx +
01044                            (int) points[4 * i]];
01045                 points[4 * i] = points[4 * i] * apix + origx;
01046                 points[4 * i + 1] = points[4 * i + 1] * apix + origy;
01047                 points[4 * i + 2] = points[4 * i + 2] * apix + origz;
01048 #ifdef DEBUG
01049                 printf("\t->\t%8g,%8g,%8g,%8g\n",points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
01050 #endif
01051         }
01052         map->update();
01053 }
01054 
01056 void PointArray::sim_updategeom() {
01057         if (!adist) adist=(double *)malloc(sizeof(double)*n);
01058         if (!aang) aang=(double *)malloc(sizeof(double)*n);
01059         if (!adihed) adihed=(double *)malloc(sizeof(double)*n);
01060         
01061         for (uint ii=0; ii<n; ii++) {
01062                 // how expensive is % ?  Worth replacing ?
01063                 int ib=4*((ii+n-1)%n);          // point before i with wraparound
01064                 int ibb=4*((ii+n-2)%n); // 2 points before i with wraparound
01065                 int ia=4*((ii+1)%n);            // 1 point after
01066                 int i=4*ii;
01067 
01068                 Vec3f a(points[ib]-points[ibb],points[ib+1]-points[ibb+1],points[ib+2]-points[ibb+2]);  // -2 -> -1
01069                 Vec3f b(points[i]-points[ib],points[i+1]-points[ib+1],points[i+2]-points[ib+2]);                // -1 -> 0
01070                 Vec3f c(points[ia]-points[i],points[ia+1]-points[i+1],points[ia+2]-points[i+2]);                // 0 -> 1
01071                 adist[ii]=b.length();
01072                 
01073                 // angle
01074                 aang[ii]=b.dot(c);
01075                 if (aang[ii]!=0) {
01076                         aang[ii]/=(adist[ii]*c.length());
01077                         if (aang[ii]>=1.0) aang[ii]=0.0;
01078                         else if (aang[ii]<=-1.0) aang[ii]=M_PI;
01079                         else aang[ii]=acos(aang[ii]);
01080                 }
01081 
01082                 // dihedral
01083                 Vec3f cr1=a.cross(b);
01084                 Vec3f cr2=b.cross(c);
01085                 double denom=cr1.length()*cr2.length();
01086                 if (denom==0) adihed[ii]=0;
01087                 else {
01088                         double tmp=cr1.dot(cr2)/(denom);
01089                         if (tmp>1) tmp=1;
01090                         if (tmp<-1) tmp=-1;
01091                         adihed[ii]=acos(tmp);
01092                 }
01093                 
01094 //              if (std::isnan(ang[ii])) ang[ii]=0;
01095                 
01096         }
01097 }
01098 
01099 double PointArray::sim_potential() {
01100         double ret=0;
01101         sim_updategeom();
01102         
01103         if (map &&mapc) {
01104                 for (uint i=0; i<n; i++) ret+=sim_pointpotential(adist[i],aang[i],adihed[i])-mapc*map->sget_value_at_interp(points[i*4]/apix+map->get_xsize()/2,points[i*4+1]/apix+map->get_ysize()/2,points[i*4+2]/apix+map->get_zsize()/2);
01105         }
01106         else {
01107                 for (uint i=0; i<n; i++) ret+=sim_pointpotential(adist[i],aang[i],adihed[i]);
01108         }
01109 
01110 #ifdef _WIN32
01111 //      if (_isnan(ret/n))
01112 #else
01113 //      if (std::isnan(ret/n))
01114 #endif
01115 //              printf("%f             %f\n",ret,n);
01116         
01117         return ret/n;
01118 }
01119 
01120 // potential for a single point. Note that if a point moves, it will impact energies +-2 from its position. This function computes only for the point i
01121 double PointArray::sim_potentiald(int i) {
01122         if (!adist) sim_updategeom();           // wasteful, but only once
01123         
01124 //      if (i<0 || i>=n) throw InvalidParameterException("Point number out of range");
01125         if (i<0) i-=n*(i/n-1);
01126         if (i>=n) i=i%n;
01127         
01128         // how expensive is % ?  Worth replacing ?
01129         int ib=4*((i+n-1)%n);           // point before i with wraparound
01130         int ibb=4*((i+n-2)%n);  // 2 points before i with wraparound
01131         int ia=4*((i+1)%n);             // 1 point after
01132         int ii=i*4;
01133         
01134         Vec3f a(points[ib]-points[ibb],points[ib+1]-points[ibb+1],points[ib+2]-points[ibb+2]);                  // -2 -> -1
01135         Vec3f b(points[ii]-points[ib],points[ii+1]-points[ib+1],points[ii+2]-points[ib+2]);             // -1 -> 0
01136         Vec3f c(points[ia]-points[ii],points[ia+1]-points[ii+1],points[ia+2]-points[ii+2]);             // 0 -> 1
01137         double dist=b.length();
01138         adist[i]=dist;
01139         // Angle, tests should avoid isnan being necessary
01140         double ang=b.dot(c);
01141         if (ang!=0.0) {                                 // if b.dot(c) is 0, we set it to the last determined value...
01142                 ang/=(dist*c.length());
01143                 if (ang>1.0) ang=1.0;           // should never happen, but just in case of roundoff error
01144                 if (ang<-1.0) ang=-1.0;
01145                 ang=acos(ang);
01146         }
01147         else ang=aang[i];
01148         aang[i]=ang;
01149         
01150         // Dihedral
01151         Vec3f cr1=a.cross(b);
01152         Vec3f cr2=b.cross(c);
01153         double dihed;
01154         double denom=cr1.length()*cr2.length();
01155         if (denom==0) dihed=adihed[i];                                          // set the dihedral to the last determined value if indeterminate
01156         else {
01157                 dihed=cr1.dot(cr2)/denom;
01158                 if (dihed>1.0) dihed=1.0;                               // should never happen, but just in case of roundoff error
01159                 if (dihed<-1.0) dihed=-1.0;
01160                 dihed=acos(dihed);
01161         }
01162         adihed[i]=dihed;
01163 //*     Do not need for small amount of points
01164         // Distance to the closest neighbor
01165         double mindist=10000;
01166         for (int j=0;j<n;j++){
01167                 if(j==i)
01168                         continue;
01169                 int ja=4*j;
01170                 Vec3f d(points[ii]-points[ja],points[ii+1]-points[ja+1],points[ii+2]-points[ja+2]);
01171                 double jdst=d.length();
01172                 if(jdst<mindist)
01173                         mindist=jdst;
01174         }
01175         double distpen=0;
01176         if (mindist<mindistc)
01177                 distpen=distpenc/mindist;
01178 //*/    
01179         
01180         
01181 //      if (std::isnan(dist) || std::isnan(ang) || std::isnan(dihed)) printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\n",i,dist,ang,dihed,b.length(),c.length(),b.dot(c)/(dist*c.length()));
01182 //      if (std::isnan(dihed)) dihed=dihed0;
01183 //      if (std::isnan(ang)) ang=0;
01184 //      if (std::isnan(dist)) dist=3.3;
01185 //      if (isnan(dihed)) dihed=dihed0;
01186 //      if (isnan(ang)) ang=0;
01187         if (map && mapc) {
01188                 return distpen+sim_pointpotential(dist,ang,dihed)-mapc*map->sget_value_at_interp(points[ii]/apix+map->get_xsize()/2,points[ii+1]/apix+map->get_ysize()/2,points[ii+2]/apix+map->get_zsize()/2);
01189         }
01190         return sim_pointpotential(dist,ang,dihed);
01191 }
01192 
01193 // Computes a potential for a point and +-2 nearest neighbors under perturbation of the central point location
01194 double PointArray::sim_potentialdxyz(int i, double dx, double dy, double dz) {
01195         Vec3f old(points[i*4],points[i*4+1],points[i*4+2]);
01196         double potd=0.;
01197                 
01198         points[i*4]=old[0]+dx;
01199         points[i*4+1]=old[1]+dy;
01200         points[i*4+2]=old[2]+dz;
01201         for (int ii=i-2; ii<=i+2; ii++) potd+=sim_potentiald(ii);
01202 //      potd=potential();
01203         points[i*4]=old[0];
01204         points[i*4+1]=old[1];
01205         points[i*4+2]=old[2];
01206         
01207         return potd;
01208 }
01209 
01210 double PointArray::calc_total_length(){
01211         double dist=0;
01212         for(int i=0; i<n; i++){
01213                 int k=(i+1)%n;
01214                 double d=(points[i*4]-points[k*4])*(points[i*4]-points[k*4])+(points[i*4+1]-points[k*4+1])*(points[i*4+1]-points[k*4+1])+(points[i*4+2]-points[k*4+2])*(points[i*4+2]-points[k*4+2]);
01215                 d=sqrt(d);
01216                 dist+=d;
01217         }
01218         return dist;
01219 }
01220 
01221 // Computes a gradient of the potential for a single point, including impact on +-2 nearest neighbors
01222 Vec3f PointArray::sim_descent(int i) {
01223         Vec3f old(points[i*4],points[i*4+1],points[i*4+2]);
01224         double pot=0.,potx=0.,poty=0.,potz=0.;
01225         double stepsz=0.01;
01226         
01227         for (int ii=i-2; ii<=i+2; ii++) pot+=sim_potentiald(ii);
01228 //      pot=potential();
01229         
01230         points[i*4]=old[0]+stepsz;
01231         for (int ii=i-2; ii<=i+2; ii++) potx+=sim_potentiald(ii);
01232 //      potx=potential();
01233         points[i*4]=old[0];
01234 
01235         points[i*4+1]=old[1]+stepsz;
01236         for (int ii=i-2; ii<=i+2; ii++) poty+=sim_potentiald(ii);
01237 //      poty=potential();
01238         points[i*4+1]=old[1];
01239 
01240         points[i*4+2]=old[2]+stepsz;
01241         for (int ii=i-2; ii<=i+2; ii++) potz+=sim_potentiald(ii);
01242 //      potz=potential();
01243         points[i*4+2]=old[2];
01244         
01245         //      printf("%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\n",pot,potx,poty,potz,(pot-potx),(pot-poty),(pot-potz));
01246         
01247 //      if (pot==potx) potx=pot+1000000.0;
01248 //      if (pot==potx) potx=pot+1000000.0;
01249 //      if (pot==potx) potx=pot+1000000.0;
01250         return Vec3f((pot-potx),(pot-poty),(pot-potz));
01251 }
01252 
01254 void PointArray::sim_minstep(double maxshift) { 
01255         vector<Vec3f> shifts;
01256         
01257         double max=0.0;
01258         double mean=0.0;
01259         for (uint i=0; i<n; i++) {
01260                 if (oldshifts.size()==n) shifts.push_back((sim_descent(i)+oldshifts[i])/2.0);
01261                 else shifts.push_back(sim_descent(i));
01262                 float len=shifts[i].length();
01263                 if (len>max) max=len;
01264                 mean+=len;
01265         }
01266         oldshifts=shifts;
01267         
01268 //      printf("max vec %1.2f\tmean %1.3f\n",max,mean/n);
01269         
01270         for (uint i=0; i<n; i++) {
01271 //              if (std::isnan(shifts[i][0]) ||std::isnan(shifts[i][1]) ||std::isnan(shifts[i][2])) { printf("Nan: %d\n",i); shifts[i]=Vec3f(max,max,max); }
01272                 points[i*4]+=shifts[i][0]*maxshift/max;
01273                 points[i*4+1]+=shifts[i][1]*maxshift/max;
01274                 points[i*4+2]+=shifts[i][2]*maxshift/max;
01275 //              printf("%d. %1.2f\t%1.2f\t%1.2f\n",i,shifts[i][0]*maxshift/max,shifts[i][1]*maxshift/max,shifts[i][2]*maxshift/max);
01276         }
01277 }
01278 
01280 void PointArray::sim_minstep_seq(double meanshift) {
01281         /*
01282         // Try to minimize potential globally
01283         boost::mt19937 rng; 
01284         boost::normal_distribution<> nd(0.0, 20.0);
01285         boost::variate_generator<boost::mt19937&,boost::normal_distribution<> > var_nor(rng, nd);
01286         double *oldpts=new double[4*get_number_points()];
01287         double *bestpts=new double[4*get_number_points()];
01288         double best_pot,new_pot;
01289         memcpy(oldpts, get_points_array(), sizeof(double) * 4 * get_number_points());
01290         memcpy(bestpts, get_points_array(), sizeof(double) * 4 * get_number_points());
01291         best_pot=sim_potential();
01292         double disttmp=0;
01293         for (int k=0; k<n; k++) disttmp+=adist[k];
01294         best_pot+=distc*pow((disttmp-336*3.3),2.0);
01295         for (int i=0; i<1000; i++){
01296                 for (int j=0; j<n; j++){
01297                         points[4*j]=oldpts[4*j]+var_nor();
01298                         points[4*j+1]=oldpts[4*j+1]+var_nor();
01299                         points[4*j+2]=oldpts[4*j+2]+var_nor();
01300                 }
01301                 new_pot=sim_potential();
01302                 disttmp=0;
01303                 for (int k=0; k<n; k++) disttmp+=adist[k];
01304                 new_pot+=distc*pow((disttmp-336*3.3),2.0);
01305                 if (new_pot<best_pot){
01306                         memcpy(bestpts, get_points_array(), sizeof(double) * 4 * get_number_points());
01307                         best_pot=new_pot;
01308                         printf("%f\t",best_pot);
01309                 }
01310         }
01311         memcpy(get_points_array(),bestpts, sizeof(double) * 4 * get_number_points());
01312         
01313         delete []oldpts;
01314         delete []bestpts;
01315         */
01316         // we compute 10 random gradients and use these to adjust stepsize
01317         double mean=0.0;
01318         for (int i=0; i<10; i++) {
01319 //              Vec3f shift=sim_descent(random()%n);
01320                 Vec3f shift=sim_descent(Util::get_irand(0,n-1));
01321                 mean+=shift.length();
01322         }
01323         mean/=10.0;
01324         double stepadj=meanshift/mean;
01325 //      printf("\t%1.4g\n",stepadj);
01326 
01327         // Now we go through all points sequentially and move each downhill
01328         // The trick here is the sequential part, as each point is impacted by the point already moved before it.
01329         // This may create a "seam" at the first point which won't be adjusted to compensate for the last point (wraparound)
01330         // until the next cycle
01331         Vec3f oshifts;
01332         for (uint ii=0; ii<n; ii++) {
01333                 uint i=2*(ii%(n/2))+2*ii/n;     // this maps a linear sequence to an all-even -> all odd sequence
01334                 Vec3f shift,d;
01335                 if (ii==n) {
01336                         d=sim_descent(i);
01337                         shift=(d+oshifts)/2.0;
01338                         oshifts=d;
01339                 }
01340                 else {
01341                         shift=sim_descent(i);
01342                         oshifts=shift;
01343                 }
01344                 
01345 //              double p2=potential();
01346 //              double pot=sim_potentialdxyz(i,0.0,0.0,0.0);
01347 //              double pots=sim_potentialdxyz(i,shift[0]*stepadj,shift[1]*stepadj,shift[2]*stepadj);
01348 
01349                 // only step if it actually improves the potential for this particle (note that it does not need to improve the overall potential)
01350 //              if (pots<pot) {
01351                         points[i*4]+=shift[0]*stepadj;
01352                         points[i*4+1]+=shift[1]*stepadj;
01353                         points[i*4+2]+=shift[2]*stepadj;
01354 //                      printf("%d. %1.4g -> %1.4g  %1.3g %1.3g %1.3g %1.3g\n",i,pot,pots,shift[0],shift[1],shift[2],stepadj);
01355 //                      if (potential()>p2) printf("%d. %1.4g %1.4g\t%1.4g %1.4g\n",i,pot,pots,p2,potential());
01356 //              }
01357                 
01358         }       
01359 }
01360 
01361 
01362 void PointArray::sim_rescale() {
01363         double meandist=0.0;
01364         double max=0.0,min=dist0*1000.0;
01365         
01366         for (uint ii=0; ii<n; ii++) {
01367                 int ib=4*((ii+n-1)%n);          // point before i with wraparound
01368                 int i=4*ii;
01369 
01370                 Vec3f b(points[i]-points[ib],points[i+1]-points[ib+1],points[i+2]-points[ib+2]);                // -1 -> 0
01371                 double len=b.length();
01372                 meandist+=len;
01373                 max=len>max?len:max;
01374                 min=len<min?len:min;
01375         }
01376         meandist/=n;
01377         double scale=dist0/meandist;
01378         
01379         printf("mean = %1.3f rescaled: %1.3f - %1.3f\n",meandist,min*scale,max*scale);
01380         
01381         for (uint i=0; i<n; i++) {
01382                 points[i*4]*=scale;
01383                 points[i*4+1]*=scale;
01384                 points[i*4+2]*=scale;
01385         }
01386 
01387 }       
01388 void PointArray::sim_printstat() {
01389         sim_updategeom();       
01390         
01391         double mdist=0.0,mang=0.0,mdihed=0.0;
01392         double midist=1000.0,miang=M_PI*2,midihed=M_PI*2;
01393         double madist=0.0,maang=0.0,madihed=0.0;
01394         
01395         for (int i=0; i<n; i++) {
01396                 mdist+=adist[i];
01397                 mang+=aang[i];
01398                 mdihed+=adihed[i];
01399 
01400                 midist=adist[i]<midist?adist[i]:midist;
01401                 madist=adist[i]>madist?adist[i]:madist;
01402                 
01403                 miang=aang[i]<miang?aang[i]:miang;
01404                 maang=aang[i]>maang?aang[i]:maang;
01405                 
01406                 midihed=adihed[i]<midihed?adihed[i]:midihed;
01407                 madihed=adihed[i]>madihed?adihed[i]:madihed;
01408         }
01409         double p=sim_potential();
01410         double anorm = 180.0/M_PI;
01411         printf(" potential: %1.1f\t dist: %1.2f / %1.2f / %1.2f\tang: %1.2f / %1.2f / %1.2f\tdihed: %1.2f / %1.2f / %1.2f  ln=%1.1f\n",p,midist,mdist/n,madist,miang*anorm,mang/n*anorm,maang*anorm,midihed*anorm,mdihed/n*anorm,madihed*anorm,mdihed/(M_PI*2.0)-n/10.0);
01412                 
01413 }
01414 
01415 void PointArray::sim_set_pot_parms(double pdist0,double pdistc,double pangc, double pdihed0, double pdihedc, double pmapc, EMData *pmap, double pmindistc,double pdistpenc) {
01416         dist0=pdist0;
01417         distc=pdistc;
01418         angc=pangc;
01419         dihed0=pdihed0;
01420         dihedc=pdihedc;
01421         mapc=pmapc;
01422         mindistc=pmindistc;
01423         distpenc=pdistpenc;
01424         if (pmap!=0 && pmap!=map) {
01425 //              if (map!=0) delete map;
01426                 if (gradx!=0) delete gradx;
01427                 if (grady!=0) delete grady;
01428                 if (gradz!=0) delete gradz;
01429                 
01430                 map=pmap;
01431                 apix=map->get_attr("apix_x");
01432 //              gradx=map->process("math.edge.xgradient");  // we compute the gradient to make the minimization easier
01433 //              grady=map->process("math.edge.ygradient");
01434 //              gradz=map->process("math.edge.zgradient");
01435         }
01436                 
01437 }
01438 
01439 // Double the number of points by adding new points on the center of edges
01440 void PointArray::sim_add_point_double() {
01441         
01442         int nn=n*2;
01443         int tmpn=n;
01444         set_number_points(nn);
01445         double* pa2data=(double *) calloc(4 * nn, sizeof(double));
01446         bool *newpt=new bool[nn];
01447         for (int i=0;i<nn;i++){
01448                 if (i%2==0)
01449                         newpt[i]=1;
01450                 else
01451                         newpt[i]=0;
01452         }
01453         int i=0;
01454         for (int ii=0;ii<nn;ii++){
01455                 if (newpt[ii]) {
01456                         
01457                         pa2data[ii*4]=points[i*4];
01458                         pa2data[ii*4+1]=points[i*4+1];
01459                         pa2data[ii*4+2]=points[i*4+2];
01460                         pa2data[ii*4+3]=1;
01461                         i++;
01462                 }
01463                 else{
01464                         int k;
01465                         if (i<tmpn)
01466                                 k=i;
01467                         else
01468                                 k=0;
01469                         pa2data[ii*4]=(points[k*4]+points[(i-1)*4])/2;
01470                         pa2data[ii*4+1]=(points[k*4+1]+points[(i-1)*4+1])/2;
01471                         pa2data[ii*4+2]=(points[k*4+2]+points[(i-1)*4+2])/2;
01472                         pa2data[ii*4+3]=1;
01473                 }
01474                         
01475         }
01476                 
01477         delete []newpt;
01478         free(points);
01479         set_points_array(pa2data);
01480         
01481         if (adist) free(adist);
01482         if (aang) free(aang);
01483         if (adihed) free(adihed);
01484         adist=aang=adihed=0;
01485         sim_updategeom();
01486 }
01487 
01488 // Delete one point with lowest density, and add two points on the edges to that one.
01489 // Or add one point on the edge with lowest density
01490 void PointArray::sim_add_point_one() {
01491 
01492         
01493         double maxpot=-1000000,pot,meanpot=0;
01494         int ipt=-1;
01495         bool onedge=0;
01496         // Find the highest potential point
01497         for (int i=0; i<n; i++) {
01498                 meanpot+=sim_pointpotential(adist[i],aang[i],adihed[i]);
01499                 pot=/*sim_pointpotential(adist[i],aang[i],adihed[i])*/-mapc*map->sget_value_at_interp(points[i*4]/apix+map->get_xsize()/2,points[i*4+1]/apix+map->get_ysize()/2,points[i*4+2]/apix+map->get_zsize()/2);
01500                 if (pot>maxpot){
01501                         maxpot=pot;
01502                         ipt=i;
01503                 }
01504         }
01505         meanpot/=n;
01506         
01507         for (int i=0; i<n; i++) {
01508                 int k=(i+1)%n;
01509                 double pt0,pt1,pt2;
01510                 pt0=(points[k*4]+points[i*4])/2;
01511                 pt1=(points[k*4+1]+points[i*4+1])/2;
01512                 pt2=(points[k*4+2]+points[i*4+2])/2;
01513                 pot=/*meanpot*/-mapc*map->sget_value_at_interp(pt0/apix+map->get_xsize()/2,pt1/apix+map->get_ysize()/2,pt2/apix+map->get_zsize()/2);
01514                 if (pot>maxpot){
01515                         maxpot=pot;
01516                         ipt=i;
01517                         onedge=1;
01518                 }
01519         
01520         }
01521         
01522         // The rest points remain the same
01523         int i;
01524         double* pa2data=(double *) calloc(4 * (n+1), sizeof(double));
01525         for (int ii=0; ii<n+1; ii++) {
01526                 if(ii!=ipt && ii!=ipt+1){
01527                         if(ii<ipt)
01528                                 i=ii;
01529                         else    // shift the points after the adding position
01530                                 i=ii-1;
01531                         
01532                         pa2data[ii*4]=points[i*4];
01533                         pa2data[ii*4+1]=points[i*4+1];
01534                         pa2data[ii*4+2]=points[i*4+2];
01535                         pa2data[ii*4+3]=1;
01536                 }
01537         }
01538         // Adding points
01539         if( onedge ) {
01540                 int k0,k1;
01541                 k0=((ipt+n-1)%n);
01542                 k1=((ipt+1)%n);
01543                 pa2data[ipt*4]=points[ipt*4];
01544                 pa2data[ipt*4+1]=points[ipt*4+1];
01545                 pa2data[ipt*4+2]=points[ipt*4+2];
01546                 pa2data[ipt*4+3]=1;
01547                 
01548                 pa2data[(ipt+1)*4]=(points[ipt*4]+points[k1*4])/2;
01549                 pa2data[(ipt+1)*4+1]=(points[ipt*4+1]+points[k1*4+1])/2;
01550                 pa2data[(ipt+1)*4+2]=(points[ipt*4+2]+points[k1*4+2])/2;
01551                 pa2data[(ipt+1)*4+3]=1; 
01552                 
01553         }
01554         else {
01555                 int k0,k1;
01556                 k0=((ipt+n-1)%n);
01557                 k1=((ipt+1)%n);
01558                 pa2data[ipt*4]=(points[ipt*4]+points[k0*4])/2;
01559                 pa2data[ipt*4+1]=(points[ipt*4+1]+points[k0*4+1])/2;
01560                 pa2data[ipt*4+2]=(points[ipt*4+2]+points[k0*4+2])/2;
01561                 pa2data[ipt*4+3]=1;
01562                 
01563                 pa2data[(ipt+1)*4]=(points[ipt*4]+points[k1*4])/2;
01564                 pa2data[(ipt+1)*4+1]=(points[ipt*4+1]+points[k1*4+1])/2;
01565                 pa2data[(ipt+1)*4+2]=(points[ipt*4+2]+points[k1*4+2])/2;
01566                 pa2data[(ipt+1)*4+3]=1; 
01567         
01568         }
01569         free(points);
01570         n++;
01571         set_points_array(pa2data);
01572         
01573         if (adist) free(adist);
01574         if (aang) free(aang);
01575         if (adihed) free(adihed);
01576         adist=aang=adihed=0;
01577         sim_updategeom();
01578         
01579         // search for the best position for the new points
01580         if (onedge){
01581                 i=ipt+1;
01582                 double bestpot=10000,nowpot;
01583                 Vec3f old(points[i*4],points[(i+1)*4],points[(i+2)*4]);
01584                 Vec3f newpt(points[i*4],points[(i+1)*4],points[(i+2)*4]);
01585                 for (int ii=0;ii<5000;ii++){
01586                         // Try to minimize potential globally
01587                         boost::mt19937 rng; 
01588                         boost::normal_distribution<> nd(0.0, 0.0);
01589                         boost::variate_generator<boost::mt19937&,boost::normal_distribution<> > var_nor(rng, nd);
01590                         points[i*4]=old[0]+var_nor();
01591                         points[i*4+1]=old[1]+var_nor();
01592                         points[i*4+2]=old[2]+var_nor();
01593                         nowpot=sim_potentiald(i);
01594                         if (nowpot<bestpot) {
01595                                 bestpot=nowpot;
01596                                 newpt[0]=points[i*4];
01597                                 newpt[1]=points[(i+1)*4];
01598                                 newpt[2]=points[(i+2)*4];
01599                         }
01600                                 
01601                 }
01602                 points[i*4]=newpt[0];
01603                 points[i*4+1]=newpt[1];
01604                 points[i*4+2]=newpt[2];
01605         }
01606         
01607 }
01608 
01609 vector<float> PointArray::do_pca(int start=0, int end=-1){
01610         
01611         if (end==-1) end=n;
01612         float covmat[9],mean[3];
01613         for (int i=0; i<3; i++) mean[i]=0;
01614         for (int i=start; i<end; i++){
01615                 for (int j=0; j<3; j++){
01616                         mean[j]+=points[i*4+j];
01617                 }
01618         }
01619         for (int i=0; i<3; i++) mean[i]/=end-start;
01620         
01621         for (int i=0; i<3; i++){
01622                 for (int j=0; j<3; j++){
01623                         if (j<i){
01624                                 covmat[i*3+j]=covmat[j*3+i];
01625                         }
01626                         else{
01627                                 covmat[i*3+j]=0;
01628                                 for (int k=start; k<end; k++)
01629                                 {
01630                                         covmat[i*3+j]+=(points[k*4+i]-mean[i])*(points[k*4+j]-mean[j]);
01631                                 }
01632                         }
01633                         
01634 //                      printf("%f\t",covmat[i*3+j]);
01635                 }
01636 //              printf("\n");
01637         }
01638         
01639         float eigval[3],eigvec[9];
01640         Util::coveig(3,covmat,eigval,eigvec);
01641         vector<float> eigv(eigvec,eigvec+sizeof(eigvec)/sizeof(float));
01642 //      printf(" %f,%f,%f\n %f,%f,%f\n %f,%f,%f\n",eigvec[0],eigvec[1],eigvec[2],eigvec[3],eigvec[4],eigvec[5],eigvec[6],eigvec[7],eigvec[8]);
01643         return eigv;
01644 }
01645 
01646 vector<float> PointArray::do_filter(vector<float> pts, float *ft, int num){
01647         // filter a 1D array
01648         vector<float> result(pts);
01649         for (uint i=0; i<pts.size(); i++)
01650                 result[i]=0;
01651         for (int i=(num-1)/2; i<pts.size()-(num-1)/2; i++){
01652                 for (int j=0; j<num; j++){
01653                         int k=i+j-(num-1)/2;
01654                         result[i]+=pts[k]*ft[j];
01655                 }
01656         }
01657         return result;
01658         
01659 }
01660 
01661 vector<double> PointArray::fit_helix(EMData* pmap,int minlength=13,float mindensity=4, vector<int> edge=vector<int>(),int twodir=0)
01662 {
01663         vector<float> hlxlen(n);
01664         vector<int> helix;
01665         map=pmap;
01666         float ft[7]={0.0044,0.0540,0.2420,0.3989,0.2420,0.0540,0.0044};
01667 //      float ft[7]={0,0,0,1,0,0,0};
01668 
01669         // search for long rods in the point array globally
01670         
01671         for (int dir=twodir; dir<2; dir++){ 
01672                 // search in both directions and combine the result
01673                 if( twodir==0)
01674                   reverse_chain();
01675                 for (uint i=0; i<n; i++){
01676                         vector<float> dist(50);
01677                         // for each point, search the following 50 points, find the longest rod
01678                         for (int len=5; len<50; len++){
01679                                 int pos=i+len;
01680                                 if (pos>=n || pos<0)    break;
01681                                 vector<float> eigvec=do_pca(i,pos); // align the points 
01682                                 vector<float> pts((len+1)*3);
01683                                 float mean[3];
01684                                 for (int k=0; k<3; k++) mean[k]=0;
01685                                 for (int k=i; k<pos; k++){
01686                                         for (int l=0; l<3; l++){
01687                                                 pts[(k-i)*3+l]=points[k*4+0]*eigvec[l*3+0]+points[k*4+1]*eigvec[l*3+1]+points[k*4+2]*eigvec[l*3+2];
01688                                                 mean[l]+=pts[(k-i)*3+l];
01689                                         }
01690                                 }
01691                                 for (int k=0; k<3; k++) mean[k]/=len;
01692                                 float dst=0;
01693                                 // distance to the center axis
01694                                 for (int k=0; k<len; k++){
01695                                         dst+=abs((pts[k*3]-mean[0])*(pts[k*3]-mean[0])+(pts[k*3+1]-mean[1])*(pts[k*3+1]-mean[1]));
01696                                 }
01697                                 dist[len]=1-dst/len/len;
01698                         }
01699                         
01700                         vector<float> nd=do_filter(dist,ft,7);
01701                         nd=do_filter(nd,ft,7);
01702                         // length of the first rod
01703                         for (int j=7; j<49; j++){
01704                                 if(nd[j]>nd[j-1] && nd[j]>nd[j+1]){
01705                                         hlxlen[i]=j-6;
01706                                         break;
01707                                 }
01708                         }
01709                         
01710                         if(hlxlen[i]>25) hlxlen[i]=0;
01711                 }
01712                 // filter the array before further process
01713 //              hlxlen[50]=100;
01714 //              for (int i=0; i<n; i++) printf("%d\t%f\n",i,hlxlen[i]);
01715 //              for (int i=0; i<3; i++) hlxlen=do_filter(hlxlen,ft,7);
01716 //              for (int i=0; i<n; i++) printf("%d\t%f\n",i,hlxlen[i]);
01717                 vector<float> ishlx(n);
01718                 int hlx=0;
01719                 float up=minlength; // rod length threshold
01720                 // record position of possible helixes
01721                 for (uint i=0; i<n; i++){
01722                         if(hlx<=0){
01723                                 if(hlxlen[i]>up){
01724                                         hlx=hlxlen[i];
01725                                         helix.push_back(i);
01726                                         helix.push_back(i+hlxlen[i]-5);
01727                                 }
01728                                 
01729                         }
01730                         else{
01731                                 hlx--;
01732                                 ishlx[i]=hlxlen[i];
01733                         }
01734                 }
01735                 // while counting reversely
01736                 if(dir==0){
01737                         for (uint i=0; i<helix.size(); i++) helix[i]=n-1-helix[i];
01738                         for (uint i=0; i<helix.size()/2; i++){
01739                                 int tmp=helix[i*2+1];
01740                                 helix[i*2+1]=helix[i*2];
01741                                 helix[i*2]=tmp;
01742                         }
01743                 }
01744                         
01745 
01746         }
01747 
01748 #ifdef DEBUG
01749         printf("potential helix counting from both sides: \n");
01750         for (uint i=0; i<helix.size()/2; i++){
01751                 printf("%d\t%d\n",helix[i*2],helix[i*2+1]);
01752         }       
01753         printf("\n\n");
01754 #endif
01755 
01756 
01757         // Combine the result from both side
01758         for (uint i=0; i<helix.size()/2; i++){
01759                 int change=1;
01760                 while(change==1){
01761                         change=0;
01762                         for (uint j=i+1; j<helix.size()/2; j++){
01763                                 if(helix[j*2]==0) continue;
01764                                 if(helix[j*2]-2<helix[i*2+1] && helix[j*2+1]+2>helix[i*2]){
01765                                         helix[i*2]=(helix[i*2]<helix[j*2])?helix[i*2]:helix[j*2];
01766                                         helix[i*2+1]=(helix[i*2+1]>helix[j*2+1])?helix[i*2+1]:helix[j*2+1];
01767                                         helix[j*2]=0;
01768                                         helix[j*2+1]=0;
01769                                         change=1;
01770                                 }
01771                         }       
01772                 }
01773         }
01774         
01775         vector<int> allhlx;
01776         int minid=1;
01777         while (minid>=0){
01778                 int mins=10000;
01779                 minid=-1;
01780                 for (uint i=0;i<helix.size()/2; i++){
01781                         if(helix[i*2]<.1) continue;
01782                         if(helix[i*2]<mins){
01783                                 mins=helix[i*2];
01784                                 minid=i;
01785                         }
01786                 }
01787                 if(minid>=0){
01788                         allhlx.push_back(helix[minid*2]);
01789                         allhlx.push_back(helix[minid*2+1]);
01790                         helix[minid*2]=-1;
01791                 }               
01792         }
01793         
01794 #ifdef DEBUG
01795         printf("combined result: \n");  
01796         for (uint i=0; i<allhlx.size()/2; i++){
01797                 printf("%d\t%d\n",allhlx[i*2],allhlx[i*2+1]);
01798         }       
01799         printf("\n\n");
01800 #endif
01801         
01802         // local search to decide the start and end point of each helix
01803 //      vector<float> allscore(allhlx.size()/2);
01804         for (uint i=0; i<allhlx.size()/2; i++){
01805                 int sz=10;
01806                 int start=allhlx[i*2]-sz,end=allhlx[i*2+1]+sz;
01807                 start=start>0?start:0;
01808                 end=end<n?end:n;
01809                 float minscr=100000;
01810                 int mj=0,mk=0;
01811                 
01812                 for (int j=start; j<end; j++){
01813                         for (int k=j+6; k<end; k++){
01814                                 vector<float> eigvec=do_pca(j,k);
01815                                 vector<float> pts((k-j)*3);
01816                                 float mean[3];
01817                                 for (int u=0; u<3; u++) mean[u]=0;
01818                                 for (int u=j; u<k; u++){
01819                                         for (int v=0; v<3; v++){
01820                                                 pts[(u-j)*3+v]=points[u*4+0]*eigvec[v*3+0]+points[u*4+1]*eigvec[v*3+1]+points[u*4+2]*eigvec[v*3+2];
01821                                                 mean[v]+=pts[(u-j)*3+v];
01822                                         }
01823                                 }
01824                                 for (int u=0; u<3; u++) mean[u]/=(k-j);
01825                                 float dst=0;
01826                                 // distance to the center axis
01827                                 for (int u=0; u<k-j; u++){
01828                                         dst+=sqrt((pts[u*3]-mean[0])*(pts[u*3]-mean[0])+(pts[u*3+1]-mean[1])*(pts[u*3+1]-mean[1]));
01829                                 }
01830                                 float len=k-j;
01831                                 float scr=dst/len/len;
01832                                 if (scr<minscr){
01833 //                                      printf("%f\t%d\t%d\n",scr,j,k);
01834                                         minscr=scr;
01835                                         mj=j;
01836                                         mk=k;
01837                                 }
01838                         }
01839                 }
01840 
01841 //              printf("%d\t%d\n",mj,mk);
01842                 
01843                 allhlx[i*2]=mj;
01844                 allhlx[i*2+1]=mk;        
01845 //              allscore[i]=minscr;
01846 //              if (mk-mj>60)
01847 //                      allscore[i]=100;
01848         }
01849         
01850         for (uint i=0; i<edge.size()/2; i++){
01851                 allhlx.push_back(edge[i*2]);
01852                 allhlx.push_back(edge[i*2+1]);
01853         }
01854         
01855         
01856         vector<int> allhlx2;
01857         minid=1;
01858         while (minid>=0){
01859                 int mins=10000;
01860                 minid=-1;
01861                 for (uint i=0;i<allhlx.size()/2; i++){
01862                         if(allhlx[i*2]<.1) continue;
01863                         if(allhlx[i*2]<mins){
01864                                 mins=allhlx[i*2];
01865                                 minid=i;
01866                         }
01867                 }
01868                 if(minid>=0){
01869                         allhlx2.push_back(allhlx[minid*2]<allhlx[minid*2+1]?allhlx[minid*2]:allhlx[minid*2+1]);
01870                         allhlx2.push_back(allhlx[minid*2]>allhlx[minid*2+1]?allhlx[minid*2]:allhlx[minid*2+1]);
01871                         allhlx[minid*2]=-1;
01872                 }               
01873         }
01874         allhlx=allhlx2;
01875         
01876 #ifdef DEBUG
01877         printf("Fitted helixes: \n");
01878         for (uint i=0; i<allhlx.size()/2; i++){
01879                 printf("%d\t%d\n",allhlx[i*2],allhlx[i*2+1]);
01880         }       
01881         printf("\n\n");
01882 #endif
01883         // create ideal helix
01884         uint ia=0,ka=0;
01885         bool dir;
01886         vector<double> finalhlx;
01887         vector<double> hlxid;
01888         printf("Confirming helix... \n");
01889         while(ia<n){
01890                 if (ia==allhlx[ka*2]){
01891 //                      int sz=(allhlx[ka*2+1]-allhlx[ka*2])>10?5:(allhlx[ka*2+1]-allhlx[ka*2])/2;
01892                         int sz=3;
01893                         float score=0,maxscr=0;
01894                         float bestphs=0,phsscore=0,pscore=0;
01895                         int mi,mj;
01896                         for (int i=0; i<sz; i++){
01897                                 for (int j=0; j<sz; j++){
01898                                         int start=allhlx[ka*2]+i,end=allhlx[ka*2+1]-j;
01899                                         phsscore=0;bestphs=-1;
01900                                         for (float phs=-180; phs<180; phs+=10){ //search for phase
01901                                                 construct_helix(start,end,phs,pscore,dir);
01902                                                 if (pscore>phsscore){
01903                                                         phsscore=pscore;
01904                                                         bestphs=phs;
01905                                                 }
01906                                         }
01907 //                                      printf("%f\t",bestphs);
01908                                         construct_helix(start,end,bestphs,score,dir);
01909                                         if (score>maxscr){
01910                                                 maxscr=score;
01911                                                 mi=i;
01912                                                 mj=j;
01913                                         }
01914                                 }
01915                         }
01916                         for (int i=0; i<mi; i++){                       
01917                                 finalhlx.push_back(points[(i+ia)*4]);
01918                                 finalhlx.push_back(points[(i+ia)*4+1]);
01919                                 finalhlx.push_back(points[(i+ia)*4+2]);
01920                         }
01921                         int start=allhlx[ka*2]+mi,end=allhlx[ka*2+1]-mj;
01922                         printf("%d\t%d\t%f\t%d\t",start,end,maxscr,maxscr>mindensity);
01923                         if (maxscr>mindensity){
01924                                 phsscore=0;
01925                                 for (float phs=-180; phs<180; phs+=10){ //search for phase
01926                                                 construct_helix(start,end,phs,pscore,dir);
01927                                                 if (pscore>phsscore){
01928                                                         phsscore=pscore;
01929                                                         bestphs=phs;
01930                                                 }
01931                                 }
01932                                 vector<double> pts=construct_helix(start,end,bestphs,score,dir);
01933                                 int lendiff=end-start-pts.size()/3-2;
01934                                 printf("%d\t",dir);
01935                                 if (pts.size()/3-2>9 && lendiff>-5){
01936                                         
01937                                         hlxid.push_back(finalhlx.size()/3+1);
01938                                         printf("%d\t",finalhlx.size()/3+1);
01939                                         for (uint j=3; j<pts.size()-3; j++)
01940                                                 finalhlx.push_back(pts[j]);
01941                                         hlxid.push_back(finalhlx.size()/3-2);
01942                                         printf("%d\t",finalhlx.size()/3-2);
01943                                         for (uint j=0; j<3; j++)
01944                                                 hlxid.push_back(pts[j]);
01945                                         for (uint j=pts.size()-3; j<pts.size(); j++)
01946                                                 hlxid.push_back(pts[j]);
01947                                         ia=end;
01948                                 }
01949                                 else{
01950                                         printf("%d\t",pts.size()/3-2);
01951                                 }
01952                         }
01953                         printf("\n");
01954                         ka++;
01955                         while(allhlx[ka*2]<ia)
01956                                 ka++;
01957                 }
01958                 else{
01959                         finalhlx.push_back(points[ia*4]);
01960                         finalhlx.push_back(points[ia*4+1]);
01961                         finalhlx.push_back(points[ia*4+2]);
01962                         ia++;                   
01963                 }
01964         }
01965         
01966         set_number_points(finalhlx.size()/3);
01967         
01968         for (uint i=0; i<n; i++){
01969                 for (uint j=0; j<3; j++)
01970                         points[i*4+j]=finalhlx[i*3+j];
01971                 points[i*4+3]=0;
01972         }
01973                         
01974 
01975         
01976         printf("\n\n");
01977         return hlxid;
01978 }
01979 
01980 vector<double> PointArray::construct_helix(int start,int end, float phs, float &score, bool &dir){
01981         // calculate length
01982         Vec3f d(points[end*4]-points[start*4],points[end*4+1]-points[start*4+1],points[end*4+2]-points[start*4+2]);
01983         double len=d.length();
01984         int nh=int(len/1.54)+2;
01985         vector<double> helix(nh*3);
01986         vector<float> eigvec=do_pca(start,end); 
01987         float eigval[3],vec[9];
01988 
01989         Util::coveig(3,&eigvec[0],eigval,vec);
01990         float maxeigv=0;
01991         int maxvi=-1;
01992         for(int i=0; i<3; i++){
01993                 if(abs(eigval[i])>maxeigv){
01994                         maxeigv=abs(eigval[i]);
01995                         maxvi=i;
01996                 }
01997         }
01998         dir=eigval[maxvi]>0;
01999 //      printf("%f\t",eigval[maxvi]);
02000 //      vector<float> eigv(eigvec,eigvec+sizeof(eigvec)/sizeof(float));
02001         // create helix
02002         helix[0]=0;helix[1]=0;helix[2]=0;
02003         helix[nh*3-3]=.0;helix[nh*3-2]=0;helix[nh*3-1]=len+.83;
02004         
02005         for (int i=0; i<nh-2; i++){
02006                 if(dir){
02007                         helix[(i+1)*3+0]=cos(((phs+(100*i))*M_PI)/180)*2.3;
02008                         helix[(i+1)*3+1]=sin(((phs+(100*i))*M_PI)/180)*2.3;
02009                 }
02010                 else{
02011                         helix[(i+1)*3+1]=cos(((phs+(100*i))*M_PI)/180)*2.3;
02012                         helix[(i+1)*3+0]=sin(((phs+(100*i))*M_PI)/180)*2.3;
02013                 }       
02014                 helix[(i+1)*3+2]=i*1.54+.83;
02015         }
02016         // transform to correct position
02017         vector<double> pts(nh*3);
02018         float mean[3];
02019         for (int k=0; k<3; k++) mean[k]=0;
02020         for (int k=0; k<nh; k++){
02021                 for (int l=0; l<3; l++){
02022                         pts[k*3+l]=helix[k*3+0]*eigvec[0*3+l]+helix[k*3+1]*eigvec[1*3+l]+helix[k*3+2]*eigvec[2*3+l];
02023                         mean[l]+=pts[k*3+l];
02024                 }
02025         }
02026         for (int k=0; k<3; k++) mean[k]/=nh;
02027         for (int k=0; k<nh; k++){
02028                 for (int l=0; l<3; l++){
02029                         pts[k*3+l]-=mean[l];
02030                 }
02031         }
02032         for (int k=0; k<3; k++) mean[k]=0;
02033         for (int k=start; k<end; k++){
02034                 for (int l=0; l<3; l++){
02035                         mean[l]+=points[k*4+l];
02036                 }
02037         }
02038         for (int k=0; k<3; k++) mean[k]/=(end-start);   
02039         for (int k=0; k<nh; k++){
02040                 for (int l=0; l<3; l++){
02041                         pts[k*3+l]+=mean[l];
02042                 }
02043         }
02044         
02045         // correct direction
02046         Vec3f d1(pts[0]-points[start*4],pts[1]-points[start*4+1],pts[2]-points[start*4+2]);
02047         Vec3f d2(pts[0]-points[end*4],pts[1]-points[end*4+1],pts[2]-points[end*4+2]);
02048                 
02049         if (d1.length()>d2.length()) { //do reverse
02050                 double tmp;
02051                 for (int i=0; i<nh/2; i++){
02052                         for(int j=0; j<3; j++){
02053                                 tmp=pts[i*3+j];
02054                                 pts[i*3+j]=pts[(nh-i-1)*3+j];
02055                                 pts[(nh-i-1)*3+j]=tmp;
02056                         }
02057                 }
02058         }
02059         
02060         // calculate score
02061 //      int sx=map->get_xsize(),sy=map->get_ysize(),sz=map->get_zsize();
02062         int sx=0,sy=0,sz=0;
02063         float ax=map->get_attr("apix_x"),ay=map->get_attr("apix_y"),az=map->get_attr("apix_z");
02064         score=0;
02065         for (int i=1; i<nh-1; i++){
02066                 score+=map->get_value_at(int(pts[i*3]/ax+sx/2),int(pts[i*3+1]/ay+sy/2),int(pts[i*3+2]/az+sz/2));
02067         }
02068         score/=(nh-2);
02069                 
02070         return pts;
02071 
02072         
02073 }
02074 
02075 void PointArray::save_pdb_with_helix(const char *file, vector<float> hlxid)
02076 {
02077 
02078         FILE *fp = fopen(file, "w");
02079         
02080         for (uint i=0; i<hlxid.size()/8; i++){
02081                 fprintf(fp, "HELIX%5lu   A ALA A%5lu  ALA A%5lu  1                        %5lu\n", 
02082                                 i, (int)hlxid[i*8], (int)hlxid[i*8+1], int(hlxid[i*8+1]-hlxid[i*8]+4));
02083         }
02084         for ( size_t i = 0; i < get_number_points(); i++) {
02085                 fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
02086                                 points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
02087         }
02088         fclose(fp);
02089 }
02090 
02091 void PointArray::remove_helix_from_map(EMData *m, vector<float> hlxid){
02092         
02093         int sx=m->get_xsize(),sy=m->get_ysize(),sz=m->get_zsize();
02094         float ax=m->get_attr("apix_x"),ay=m->get_attr("apix_y"),az=m->get_attr("apix_z");
02095         for (int x=0; x<sx; x++){
02096                 for (int y=0; y<sy; y++){
02097                         for (int z=0; z<sz; z++){
02098                                 Vec3f p0((x)*ax,(y)*ay,(z)*az);
02099                                 bool inhlx=false;
02100                                 for (uint i=0; i<hlxid.size()/8; i++){
02101                                         Vec3f p1(hlxid[i*8+2],hlxid[i*8+3],hlxid[i*8+4]),p2(hlxid[i*8+5],hlxid[i*8+6],hlxid[i*8+7]);
02102                                         Vec3f dp=p2-p1;
02103                                         float l=dp.length();
02104                                         float d=((p0-p1).cross(p0-p2)).length()/l;
02105                                         float t=-(p1-p0).dot(p2-p1)/(l*l);
02106                                         if (d<5 && t>0 && t<1){
02107                                                 inhlx=true;
02108                                                 break;
02109                                         }
02110                                 }
02111                                 if(inhlx){
02112                                         m->set_value_at(x,y,z,0);
02113                                 }
02114                                         
02115                         }
02116                 }
02117         }
02118 }
02119 
02120 void PointArray::reverse_chain(){
02121         // reverse the point array chain, from the last to the first point
02122         double tmp;
02123 //      for(int i=0; i<n/2; i++){
02124 //              for (int j=0; j<4; j++){
02125 //                      printf("%f\t",
02126         for(int i=0; i<n/2; i++){
02127                 for (int j=0; j<4; j++){
02128                         tmp=points[(n-1-i)*4+j];
02129                         points[(n-1-i)*4+j]=points[i*4+j];
02130                         points[i*4+j]=tmp;
02131                 }
02132         }
02133 }
02134 
02135 void PointArray::sort_by_axis(int axis)
02136 {
02137         if (axis == 0)
02138                 qsort(points, n, sizeof(double) * 4, cmp_axis_x);
02139         else if (axis == 1)
02140                 qsort(points, n, sizeof(double) * 4, cmp_axis_y);
02141         else if (axis == 2)
02142                 qsort(points, n, sizeof(double) * 4, cmp_axis_z);
02143         else
02144                 qsort(points, n, sizeof(double) * 4, cmp_val);
02145 }
02146 
02147 
02148 EMData *PointArray::pdb2mrc_by_summation(int map_size, float apix, float res, int addpdbbfactor)
02149 {
02150 #ifdef DEBUG
02151         printf("PointArray::pdb2mrc_by_summation(): %lu points\tmapsize = %4d\tapix = %g\tres = %g\n",get_number_points(),map_size, apix, res);
02152 #endif
02153         //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);
02154 
02155         double min_table_val = 1e-7;
02156         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
02157 
02158         double table_step_size = 0.001; // number of steps for each pixel
02159         double inv_table_step_size = 1.0 / table_step_size;
02160         
02161 //      sort_by_axis(2);                        // sort by Z-axis
02162 
02163         EMData *map = new EMData();
02164         map->set_size(map_size, map_size, map_size);
02165         map->to_zero();
02166         float *pd = map->get_data();
02167         
02168         vector<double> table;
02169         double gauss_real_width;
02170         int table_size;
02171         int gbox;
02172         
02173         
02174         for ( size_t s = 0; s < get_number_points(); ++s) {
02175                 double xc = points[4 * s] / apix + map_size / 2;
02176                 double yc = points[4 * s + 1] / apix + map_size / 2;
02177                 double zc = points[4 * s + 2] / apix + map_size / 2;
02178                 double fval = points[4 * s + 3];
02179                 
02180                 //printf("\n bfactor=%f",bfactor[s]);
02181                 
02182                 
02183                 
02184                 
02185                 
02186                 if(addpdbbfactor==-1){
02187                         gauss_real_width = res/M_PI;    // in Angstrom, res is in Angstrom
02188                 }else{
02189                         gauss_real_width = (bfactor[s])/(4*sqrt(2.0)*M_PI);     // in Angstrom, res is in Angstrom
02190                 }
02191                 
02192                 
02193                 table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
02194                 table.resize(table_size);
02195                 for (int i = 0; i < table_size; i++){
02196                         table[i] = 0;
02197                 }
02198                 
02199                 for (int i = 0; i < table_size; i++) {                                          
02200                         double x = -i * table_step_size * apix / gauss_real_width;
02201                         if(addpdbbfactor==-1){
02202                                 table[i] =  exp(-x * x);        
02203                         }
02204                         else{
02205                         table[i] =  exp(-x * x)/sqrt(gauss_real_width * gauss_real_width * 2 * M_PI);
02206                         }
02207                 }
02208 
02209                 gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
02210                 if (gbox <= 0)
02211                         gbox = 1;
02212                 
02213                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
02214                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
02215                 int kmin = int (zc) - gbox, kmax = int (zc) + gbox;
02216                 if (imin < 0)
02217                         imin = 0;
02218                 if (jmin < 0)
02219                         jmin = 0;
02220                 if (kmin < 0)
02221                         kmin = 0;
02222                 if (imax > map_size)
02223                         imax = map_size;
02224                 if (jmax > map_size)
02225                         jmax = map_size;
02226                 if (kmax > map_size)
02227                         kmax = map_size;                
02228                 
02229                 for (int k = kmin; k < kmax; k++) {
02230                         size_t table_index_z = size_t (fabs(k - zc) * inv_table_step_size);
02231                         if ( table_index_z >= table.size() ) continue;
02232                         double zval = table[table_index_z];
02233                         size_t pd_index_z = k * map_size * map_size;
02234                         
02235                         for (int j = jmin; j < jmax; j++) {
02236                                 size_t table_index_y = size_t (fabs(j - yc) * inv_table_step_size);
02237                                 if ( table_index_y >= table.size() ) continue;
02238                                 double yval = table[table_index_y];
02239                                 size_t pd_index = pd_index_z + j * map_size + imin;
02240                                 for (int i = imin; i < imax; i++, pd_index++) {
02241                                         size_t table_index_x = size_t (fabs(i - xc) * inv_table_step_size);
02242                                         if ( table_index_x >= table.size() ) continue;
02243                                         double xval = table[table_index_x];
02244                                         pd[pd_index] += (float) (fval * zval * yval * xval);
02245                                 }
02246                         }
02247                 }
02248         }
02249 
02250         map->update();
02251         map->set_attr("apix_x", apix);
02252         map->set_attr("apix_y", apix);
02253         map->set_attr("apix_z", apix);
02254         map->set_attr("origin_x", -map_size/2*apix);
02255         map->set_attr("origin_y", -map_size/2*apix);
02256         map->set_attr("origin_z", -map_size/2*apix);
02257 
02258         return map;
02259 }
02260 
02261 
02262 EMData *PointArray::projection_by_summation(int image_size, float apix, float res)
02263 {
02264         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
02265         //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);
02266 
02267         double min_table_val = 1e-7;
02268         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
02269 
02270         //double table_step_size = 0.001;    // number of steps for x=[0,1] in exp(-x*x)
02271         //int table_size = int(max_table_x / table_step_size *1.25);
02272         //double* table = (double*)malloc(sizeof(double) * table_size);
02273         //for(int i=0; i<table_size; i++) table[i]=exp(-i*i*table_step_size*table_step_size);
02274 
02275         double table_step_size = 0.001; // number of steps for each pixel
02276         double inv_table_step_size = 1.0 / table_step_size;
02277         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
02278         double *table = (double *) malloc(sizeof(double) * table_size);
02279         for (int i = 0; i < table_size; i++) {
02280                 double x = -i * table_step_size * apix / gauss_real_width;
02281                 table[i] = exp(-x * x);
02282         }
02283 
02284         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
02285         if (gbox <= 0)
02286                 gbox = 1;
02287         EMData *proj = new EMData();
02288         proj->set_size(image_size, image_size, 1);
02289         proj->to_zero();
02290         float *pd = proj->get_data();
02291         for ( size_t s = 0; s < get_number_points(); ++s) {
02292                 double xc = points[4 * s] / apix + image_size / 2;
02293                 double yc = points[4 * s + 1] / apix + image_size / 2;
02294                 double fval = points[4 * s + 3];
02295                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
02296                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
02297                 if (imin < 0)
02298                         imin = 0;
02299                 if (jmin < 0)
02300                         jmin = 0;
02301                 if (imax > image_size)
02302                         imax = image_size;
02303                 if (jmax > image_size)
02304                         jmax = image_size;
02305 
02306                 for (int j = jmin; j < jmax; j++) {
02307                         //int table_index_y = int(fabs(j-yc)*apix/gauss_real_width/table_step_size);
02308                         int table_index_y = int (fabs(j - yc) * inv_table_step_size);
02309                         double yval = table[table_index_y];
02310 #ifdef DEBUG
02311                         //double yval2 = exp( - (j-yc)*(j-yc)*apix*apix/(gauss_real_width*gauss_real_width));
02312                         //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);
02313 #endif
02314                         int pd_index = j * image_size + imin;
02315                         for (int i = imin; i < imax; i++, pd_index++) {
02316                                 //int table_index_x = int(fabs(i-xc)*apix/gauss_real_width/table_step_size);
02317                                 int table_index_x = int (fabs(i - xc) * inv_table_step_size);
02318                                 double xval = table[table_index_x];
02319 #ifdef DEBUG
02320                                 //double xval2 = exp( - (i-xc)*(i-xc)*apix*apix/(gauss_real_width*gauss_real_width));
02321                                 //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);
02322 #endif
02323                                 pd[pd_index] += (float)(fval * yval * xval);
02324                         }
02325                 }
02326         }
02327         for (int i = 0; i < image_size * image_size; i++)
02328                 pd[i] /= sqrt(M_PI);
02329         proj->update();
02330         return proj;
02331 }
02332 
02333 void PointArray::replace_by_summation(EMData *proj, int ind, Vec3f vec, float amp, float apix, float res)
02334 {
02335         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
02336 
02337         double min_table_val = 1e-7;
02338         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
02339 
02340         double table_step_size = 0.001; // number of steps for each pixel
02341         double inv_table_step_size = 1.0 / table_step_size;
02342         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
02343         double *table = (double *) malloc(sizeof(double) * table_size);
02344         for (int i = 0; i < table_size; i++) {
02345                 double x = -i * table_step_size * apix / gauss_real_width;
02346                 table[i] = exp(-x * x)/pow((float)M_PI,.25f);
02347         }
02348         int image_size=proj->get_xsize();
02349 
02350         // subtract the old point
02351         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
02352         if (gbox <= 0)
02353                 gbox = 1;
02354         float *pd = proj->get_data();
02355          int s = ind;
02356         double xc = points[4 * s] / apix + image_size / 2;
02357         double yc = points[4 * s + 1] / apix + image_size / 2;
02358         double fval = points[4 * s + 3];
02359         int imin = int (xc) - gbox, imax = int (xc) + gbox;
02360         int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
02361 
02362         if (imin < 0) imin = 0;
02363         if (jmin < 0) jmin = 0;
02364         if (imax > image_size) imax = image_size;
02365         if (jmax > image_size) jmax = image_size;
02366 
02367         for (int j = jmin; j < jmax; j++) {
02368                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
02369                 double yval = table[table_index_y];
02370                 int pd_index = j * image_size + imin;
02371                 for (int i = imin; i < imax; i++, pd_index++) {
02372                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
02373                         double xval = table[table_index_x];
02374                         pd[pd_index] -= (float)(fval * yval * xval);
02375                 }
02376         }
02377 
02378         // add the new point
02379         gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
02380         if (gbox <= 0)
02381                 gbox = 1;
02382         pd = proj->get_data();
02383         s = ind;
02384         xc = vec[0] / apix + image_size / 2;
02385         yc = vec[1] / apix + image_size / 2;
02386         fval = amp;
02387         imin = int (xc) - gbox, imax = int (xc) + gbox;
02388         jmin = int (yc) - gbox, jmax = int (yc) + gbox;
02389 
02390         if (imin < 0) imin = 0;
02391         if (jmin < 0) jmin = 0;
02392         if (imax > image_size) imax = image_size;
02393         if (jmax > image_size) jmax = image_size;
02394 
02395         for (int j = jmin; j < jmax; j++) {
02396                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
02397                 double yval = table[table_index_y];
02398                 int pd_index = j * image_size + imin;
02399                 for (int i = imin; i < imax; i++, pd_index++) {
02400                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
02401                         double xval = table[table_index_x];
02402                         pd[pd_index] -= (float)(fval * yval * xval);
02403                 }
02404         }
02405 
02406         proj->update();
02407         return;
02408 }
02409 
02410 
02411 EMData *PointArray::pdb2mrc_by_nfft(int , float , float )
02412 {
02413 #if defined NFFT
02414         nfft_3D_plan my_plan;           // plan for the nfft
02415 
02417         nfft_3D_init(&my_plan, map_size, get_number_points());
02418 
02420         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
02421                 // FFTW and nfft use row major array layout, EMAN uses column major
02422                 my_plan.v[3 * j + 2] = (fftw_real) (points[i] / (apix * map_size));
02423                 my_plan.v[3 * j + 1] = (fftw_real) (points[i + 1] / (apix * map_size));
02424                 my_plan.v[3 * j] = (fftw_real) (points[i + 2] / (apix * map_size));
02425                 my_plan.f[j].re = (fftw_real) (points[i + 3]);
02426                 my_plan.f[j].im = 0.0;
02427         }
02428 
02430         if (my_plan.nfft_flags & PRE_PSI) {
02431                 nfft_3D_precompute_psi(&my_plan);
02432         }
02433 
02434         // compute the uniform Fourier transform
02435         nfft_3D_transpose(&my_plan);
02436 
02437         // copy the Fourier transform to EMData data array
02438         EMData *fft = new EMData();
02439         fft->set_size(map_size + 2, map_size, map_size);
02440         fft->set_complex(true);
02441         fft->set_ri(true);
02442         fft->to_zero();
02443         float *data = fft->get_data();
02444         double norm = 1.0 / (map_size * map_size * map_size);
02445         for (int k = 0; k < map_size; k++) {
02446                 for (int j = 0; j < map_size; j++) {
02447                         for (int i = 0; i < map_size / 2; i++) {
02448                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
02449                                         (float) (my_plan.
02450                                                          f_hat[k * map_size * map_size + j * map_size + i +
02451                                                                    map_size / 2].re) * norm;
02452                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
02453                                         (float) (my_plan.
02454                                                          f_hat[k * map_size * map_size + j * map_size + i +
02455                                                                    map_size / 2].im) * norm * (-1.0);
02456                         }
02457                 }
02458         }
02460         nfft_3D_finalize(&my_plan);
02461 
02462         // low pass processor
02463         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
02464         int index = 0;
02465         for (int k = 0; k < map_size; k++) {
02466                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
02467                 for (int j = 0; j < map_size; j++) {
02468                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
02469                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
02470                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
02471                                 data[index] *= val;
02472                                 data[index + 1] *= val;
02473                         }
02474                 }
02475         }
02476         fft->update();
02477         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
02478 
02479         fft->process_inplace("xform.phaseorigin.tocorner");     // move phase origin to center of image map_size, instead of at corner
02480         EMData *map = fft->do_ift();
02481         map->set_attr("apix_x", apix);
02482         map->set_attr("apix_y", apix);
02483         map->set_attr("apix_z", apix);
02484         map->set_attr("origin_x", -map_size/2*apix);
02485         map->set_attr("origin_y", -map_size/2*apix);
02486         map->set_attr("origin_z", -map_size/2*apix);
02487         if( fft )
02488         {
02489                 delete fft;
02490                 fft = 0;
02491         }
02492         return map;
02493 #elif defined NFFT2
02494         nfft_plan my_plan;                      // plan for the nfft
02495 
02497         nfft_init_3d(&my_plan, map_size, map_size, map_size, get_number_points());
02498 
02500         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
02501                 // FFTW and nfft use row major array layout, EMAN uses column major
02502                 my_plan.x[3 * j + 2] = (double) (points[i] / (apix * map_size));
02503                 my_plan.x[3 * j + 1] = (double) (points[i + 1] / (apix * map_size));
02504                 my_plan.x[3 * j] = (double) (points[i + 2] / (apix * map_size));
02505                 my_plan.f[j][0] = (double) (points[i + 3]);
02506                 my_plan.f[j][1] = 0.0;
02507         }
02508 
02510         if (my_plan.nfft_flags & PRE_PSI) {
02511                 nfft_precompute_psi(&my_plan);
02512                 if (my_plan.nfft_flags & PRE_FULL_PSI)
02513                         nfft_full_psi(&my_plan, pow(10, -10));
02514         }
02515 
02516         // compute the uniform Fourier transform
02517         nfft_transposed(&my_plan);
02518 
02519         // copy the Fourier transform to EMData data array
02520         EMData *fft = new EMData();
02521         fft->set_size(map_size + 2, map_size, map_size);
02522         fft->set_complex(true);
02523         fft->set_ri(true);
02524         fft->to_zero();
02525         float *data = fft->get_data();
02526         double norm = 1.0 / (map_size * map_size * map_size);
02527         for (int k = 0; k < map_size; k++) {
02528                 for (int j = 0; j < map_size; j++) {
02529                         for (int i = 0; i < map_size / 2; i++) {
02530                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
02531                                         (float) (my_plan.
02532                                                          f_hat[k * map_size * map_size + j * map_size + i +
02533                                                                    map_size / 2][0]) * norm;
02534                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
02535                                         (float) (my_plan.
02536                                                          f_hat[k * map_size * map_size + j * map_size + i +
02537                                                                    map_size / 2][1]) * norm;
02538                         }
02539                 }
02540         }
02542         nfft_finalize(&my_plan);
02543 
02544         // low pass processor
02545         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
02546         int index = 0;
02547         for (int k = 0; k < map_size; k++) {
02548                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
02549                 for (int j = 0; j < map_size; j++) {
02550                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
02551                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
02552                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
02553                                 data[index] *= val;
02554                                 data[index + 1] *= val;
02555                         }
02556                 }
02557         }
02558         fft->update();
02559         //fft->process_inplace(".filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
02560 
02561         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image map_size, instead of at corner
02562         EMData *map = fft->do_ift();
02563         map->set_attr("apix_x", apix);
02564         map->set_attr("apix_y", apix);
02565         map->set_attr("apix_z", apix);
02566         map->set_attr("origin_x", -map_size / 2 * apix);
02567         map->set_attr("origin_y", -map_size / 2 * apix);
02568         map->set_attr("origin_z", -map_size / 2 * apix);
02569         if( fft )
02570         {
02571                 delete fft;
02572                 fft = 0;
02573         }
02574         return map;
02575 #else
02576         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
02577         return 0;
02578 #endif
02579 }
02580 
02581 EMData *PointArray::projection_by_nfft(int , float , float )
02582 {
02583 #if defined NFFT
02584         nfft_2D_plan my_plan;           // plan for the nfft
02585         int N[2], n[2];
02586         N[0] = image_size;
02587         n[0] = next_power_of_2(2 * image_size);
02588         N[1] = image_size;
02589         n[1] = next_power_of_2(2 * image_size);
02590 
02592         nfft_2D_init(&my_plan, image_size, get_number_points());
02593         //nfft_2D_init_specific(&my_plan, N, get_number_points(), n, 3,
02594         //                 PRE_PHI_HUT | PRE_PSI,
02595         //                 FFTW_ESTIMATE | FFTW_OUT_OF_PLACE);
02597         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
02598                 // FFTW and nfft use row major array layout, EMAN uses column major
02599                 my_plan.v[2 * j + 1] = (fftw_real) (points[i] / (apix * image_size));
02600                 my_plan.v[2 * j] = (fftw_real) (points[i + 1] / (apix * image_size));
02601                 my_plan.f[j].re = (fftw_real) points[i + 3];
02602                 my_plan.f[j].im = 0.0;
02603         }
02604 
02606         if (my_plan.nfft_flags & PRE_PSI) {
02607                 nfft_2D_precompute_psi(&my_plan);
02608         }
02609 
02610         // compute the uniform Fourier transform
02611         nfft_2D_transpose(&my_plan);
02612 
02613         // copy the Fourier transform to EMData data array
02614         EMData *fft = new EMData();
02615         fft->set_size(image_size + 2, image_size, 1);
02616         fft->set_complex(true);
02617         fft->set_ri(true);
02618         fft->to_zero();
02619         float *data = fft->get_data();
02620         double norm = 1.0 / (image_size * image_size);
02621         for (int j = 0; j < image_size; j++) {
02622                 for (int i = 0; i < image_size / 2; i++) {
02623                         data[j * (image_size + 2) + 2 * i] =
02624                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].re) * norm;
02625                         data[j * (image_size + 2) + 2 * i + 1] =
02626                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].im) * norm * (-1.0);
02627                 }
02628         }
02630         nfft_2D_finalize(&my_plan);
02631 
02632         if (res > 0) {
02633                 // Gaussian low pass processor
02634                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
02635                 int index = 0;
02636                 for (int j = 0; j < image_size; j++) {
02637                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
02638                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
02639                                 double val = exp(-(i * i + RY2) / sigma2);
02640                                 data[index] *= val;
02641                                 data[index + 1] *= val;
02642                         }
02643                 }
02644         }
02645         fft->update();
02646         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
02647 
02648         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
02649 
02650         return fft;
02651 #elif defined NFFT2
02652         nfft_plan my_plan;                      // plan for the nfft
02653         int N[2], n[2];
02654         N[0] = image_size;
02655         n[0] = next_power_of_2(2 * image_size);
02656         N[1] = image_size;
02657         n[1] = next_power_of_2(2 * image_size);
02658 
02660         //nfft_init_2d(&my_plan,image_size,image_size,get_number_points());
02661         nfft_init_specific(&my_plan, 2, N, get_number_points(), n, 3,
02662                                            PRE_PHI_HUT | PRE_PSI |
02663                                            MALLOC_X | MALLOC_F_HAT | MALLOC_F, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
02665         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
02666                 // FFTW and nfft use row major array layout, EMAN uses column major
02667                 my_plan.x[2 * j + 1] = (double) (points[i] / (apix * image_size));
02668                 my_plan.x[2 * j] = (double) (points[i + 1] / (apix * image_size));
02669                 my_plan.f[j][0] = (double) points[i + 3];
02670                 my_plan.f[j][1] = 0.0;
02671         }
02672 
02674         if (my_plan.nfft_flags & PRE_PSI) {
02675                 nfft_precompute_psi(&my_plan);
02676                 if (my_plan.nfft_flags & PRE_FULL_PSI)
02677                         nfft_full_psi(&my_plan, pow(10, -6));
02678         }
02679 
02680         // compute the uniform Fourier transform
02681         nfft_transposed(&my_plan);
02682 
02683         // copy the Fourier transform to EMData data array
02684         EMData *fft = new EMData();
02685         fft->set_size(image_size + 2, image_size, 1);
02686         fft->set_complex(true);
02687         fft->set_ri(true);
02688         fft->to_zero();
02689         float *data = fft->get_data();
02690         double norm = 1.0 / (image_size * image_size);
02691         for (int j = 0; j < image_size; j++) {
02692                 for (int i = 0; i < image_size / 2; i++) {
02693                         data[j * (image_size + 2) + 2 * i] =
02694                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][0]) * norm;
02695                         data[j * (image_size + 2) + 2 * i + 1] =
02696                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][1]) * norm;
02697                 }
02698         }
02700         nfft_finalize(&my_plan);
02701 
02702         if (res > 0) {
02703                 // Gaussian low pass process
02704                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
02705                 int index = 0;
02706                 for (int j = 0; j < image_size; j++) {
02707                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
02708                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
02709                                 double val = exp(-(i * i + RY2) / sigma2);
02710                                 data[index] *= val;
02711                                 data[index + 1] *= val;
02712                         }
02713                 }
02714         }
02715         fft->update();
02716         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
02717 
02718         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
02719 
02720         return fft;
02721 #else
02722         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
02723         return 0;
02724 #endif
02725 }
02726 
02727 #ifdef OPTPP
02728 #include "NLF.h"
02729 #include "BoundConstraint.h"
02730 #include "OptCG.h"
02731 //#include "OptNewton.h"
02732 #include "newmatap.h"
02733 
02734 vector<EMData*> optdata;
02735 PointArray *optobj;
02736 float optpixres;
02737 
02738 void init_opt_proj(int ndim, ColumnVector& x)
02739 {
02740 int i;
02741 double *data=optobj->get_points_array();
02742 
02743 for (i=0; i<ndim; i++) x(i+1)=data[i];
02744 }
02745 
02746 void calc_opt_proj(int n, const ColumnVector& x, double& fx, int& result)
02747 {
02748         int i;
02749         PointArray pa;
02750         Transform xform;
02751         int size=optdata[0]->get_xsize();
02752         fx=0;
02753 
02754         for (i=0; i<optdata.size(); i++) {
02755                 xform=(optdata[i]->get_transform());
02756                 pa.set_from((double *)x.nric()+1,n/4,std::string("c1"),&xform);
02757                 EMData *p=pa.projection_by_summation(size,1.0,optpixres);
02758                 p->process_inplace("normalize.unitlen");
02759                 fx-=sqrt(p->cmp("dot",EMObject(optdata[i]),Dict()));
02760         }
02761 
02762         result=NLPFunction;
02763 
02764         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));
02765 }
02766 
02767 void calc_opt_projd(int mode,int n, const ColumnVector& x, double& fx, ColumnVector& gx, int& result)
02768 {
02769 
02770 if (mode & NLPFunction) {
02771 
02772 }
02773 
02774 if (mode & NLPGradient) {
02775 
02776 }
02777 
02778 }
02779 #endif
02780 
02781 void PointArray::opt_from_proj(const vector<EMData*> & proj,float pixres) {
02782 #ifdef OPTPP
02783         optdata=proj;
02784         optobj=this;
02785         optpixres=pixres;
02786 
02787         FDNLF1 nlf(get_number_points()*4,calc_opt_proj,init_opt_proj);
02788 //      NLF1 nlf(get_number_points()*4,init_opt_proj,calc_opt_projd);
02789         nlf.initFcn();
02790 
02791         OptCG opt(&nlf);
02792 //      OptQNewton opt(&nlf);
02793         opt.setMaxFeval(2000);
02794         opt.optimize();
02795         opt.printStatus("Done");
02796 #else
02797         (void)proj;             //suppress warning message
02798         (void)pixres;   //suppress warning message
02799         LOGWARN("OPT++ support not enabled.\n");
02800         return;
02801 #endif
02802 }
02803 
02804 Vec3f PointArray::get_vector_at(int i)
02805 {
02806 return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
02807 }
02808 
02809 double PointArray::get_value_at(int i)
02810 {
02811 return points[i*4+3];
02812 }
02813 
02814 void PointArray::set_vector_at(int i,Vec3f vec,double value)
02815 {
02816 points[i*4]=vec[0];
02817 points[i*4+1]=vec[1];
02818 points[i*4+2]=vec[2];
02819 points[i*4+3]=value;
02820 }
02821 
02822 void PointArray::set_vector_at(int i,vector<double> v)
02823 {
02824 points[i*4]  =v[0];
02825 points[i*4+1]=v[1];
02826 points[i*4+2]=v[2];
02827 if (v.size()>=4) points[i*4+3]=v[3];
02828 }