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

pointarray.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Wen Jiang, 08/11/2004 (jiang12@purdue.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "pointarray.h"
00037 #include <vector>
00038 #include <cstring>
00039 
00040 using namespace EMAN;
00041 
00042 int cmp_axis_x(const void *a, const void *b)
00043 {
00044         double diff = ((double *) a)[0] - ((double *) b)[0];
00045         if (diff < 0.0)
00046                 return -1;
00047         else if (diff > 0.0)
00048                 return 1;
00049         else
00050                 return 0;
00051 }
00052 int cmp_axis_y(const void *a, const void *b)
00053 {
00054         double diff = ((double *) a)[1] - ((double *) b)[1];
00055         if (diff < 0.0)
00056                 return -1;
00057         else if (diff > 0.0)
00058                 return 1;
00059         else
00060                 return 0;
00061 }
00062 int cmp_axis_z(const void *a, const void *b)
00063 {
00064         double diff = ((double *) a)[2] - ((double *) b)[2];
00065         if (diff < 0.0)
00066                 return -1;
00067         else if (diff > 0.0)
00068                 return 1;
00069         else
00070                 return 0;
00071 }
00072 int cmp_val(const void *a, const void *b)
00073 {
00074         double diff = ((double *) a)[3] - ((double *) b)[3];
00075         if (diff < 0.0)
00076                 return -1;
00077         else if (diff > 0.0)
00078                 return 1;
00079         else
00080                 return 0;
00081 }
00082 // This will do a sort in descending order
00083 int cmp_float(const void *a, const void *b)
00084 {
00085         double diff = *((float *) a) - *((float *) b);
00086         if (diff < 0.0)
00087                 return 1;
00088         else if (diff > 0.0)
00089                 return -1;
00090         else
00091                 return 0;
00092 }
00093 
00094 
00095 PointArray::PointArray()
00096 {
00097         points = 0;
00098         bfactor = 0;
00099         n = 0;
00100 }
00101 
00102 PointArray::PointArray( int nn)
00103 {
00104         n = nn;
00105         points = (double *) calloc(4 * n, sizeof(double));
00106 }
00107 
00108 PointArray::~PointArray()
00109 {
00110         if( points )
00111         {
00112                 free(points);
00113                 points = 0;
00114         }
00115 }
00116 
00117 void PointArray::zero()
00118 {
00119         memset((void *) points, 0, 4 * n * sizeof(double));
00120 }
00121 
00122 PointArray *PointArray::copy()
00123 {
00124         PointArray *pa2 = new PointArray();
00125         pa2->set_number_points(get_number_points());
00126         double *pa2data = pa2->get_points_array();
00127         memcpy(pa2data, get_points_array(), sizeof(double) * 4 * get_number_points());
00128 
00129         return pa2;
00130 }
00131 
00132 PointArray & PointArray::operator=(PointArray & pa)
00133 {
00134         if (this != &pa) {
00135                 set_number_points(pa.get_number_points());
00136                 memcpy(get_points_array(), pa.get_points_array(), sizeof(double) * 4 * get_number_points());
00137         }
00138         return *this;
00139 }
00140 
00141  size_t PointArray::get_number_points() const
00142 {
00143         return n;
00144 }
00145 
00146 void PointArray::set_number_points(size_t nn)
00147 {       
00148         if (n != nn) {
00149                 n = nn;         
00150                 points = (double *) realloc(points, 4 * n * sizeof(double));
00151                 bfactor = (double *) realloc(bfactor, n * sizeof(double));
00152         }
00153 }
00154 
00155 double *PointArray::get_points_array()
00156 {
00157         return points;
00158 }
00159 
00160 void PointArray::set_points_array(double *p)
00161 {
00162         points = p;
00163 }
00164 
00165 EMData *PointArray::distmx(int sortrows) {
00166 if (n==0) return NULL;
00167 
00168 unsigned int i,j;
00169 
00170 EMData *ret= new EMData(n,n,1);
00171 ret->to_zero();
00172 
00173 for (i=0; i<n; i++) {
00174         for (j=i+1; j<n; j++) {
00175                 float r=(get_vector_at(i)-get_vector_at(j)).length();
00176                 ret->set_value_at(i,j,0,r);
00177                 ret->set_value_at(j,i,0,r);
00178         }
00179 }
00180 
00181 if (sortrows) {
00182         float *data=ret->get_data();
00183         for (i=0; i<n; i++) qsort(&data[i*n],n,sizeof(float),cmp_float);
00184         ret->update();
00185 }
00186 
00187 return ret;
00188 }
00189 
00190 vector<int> PointArray::match_points(PointArray *to,float max_miss) {
00191 EMData *mx0=distmx(1);
00192 EMData *mx1=to->distmx(1);
00193 unsigned int n2=mx1->get_xsize();       // same as get_number_points on to
00194 
00195 if (max_miss<0) max_miss=(float)mx0->get_attr("sigma")/10.0f;
00196 //printf("max error %f\n",max_miss);
00197 
00198 
00199 
00200 vector<int> ret(n,-1);
00201 vector<float> rete(n,0.0);
00202 unsigned int i,j,k,l;
00203 
00204 if (!mx0 || !mx1) {
00205         if (mx0) delete mx0;
00206         if (mx1) delete mx1;
00207         return ret;
00208 }
00209 
00210 // i iterates over elements of 'this', j looks for a match in 'to'
00211 // k and l iterate over the individual distances
00212 for (i=0; i<n; i++) {
00213         int bestn=-1;                   // number of best match in mx1
00214         double bestd=1.0e38;            // residual error distance in best match
00215         for (j=0; j<n2; j++) {
00216                 double d=0;
00217                 int nd=0;
00218                 for (k=l=0; k<n-1 && l<n2-1; k++,l++) {
00219                         float d1=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l,j));
00220                         float d2=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l,j));
00221                         float d3=fabs(mx0->get_value_at(k,i)-mx1->get_value_at(l+1,j));
00222                         float d4=fabs(mx0->get_value_at(k+1,i)-mx1->get_value_at(l+1,j));
00223                         if (d2<d1 && d4>d2) { l--; continue; }
00224                         if (d3<d1 && d4>d3) { k--; continue; }
00225                         d+=d1;
00226                         nd++;
00227                 }
00228                 d/=(float)nd;
00229 //              printf("%d -> %d\t%f\t%d\n",i,j,d,nd);
00230                 if (d<bestd) { bestd=d; bestn=j; }
00231         }
00232         ret[i]=bestn;
00233         rete[i]=static_cast<float>(bestd);
00234 }
00235 
00236 // This will remove duplicate assignments, keeping the best one
00237 // also remove any assignments with large errors
00238 for (i=0; i<n; i++) {
00239         for (j=0; j<n; j++) {
00240                 if (rete[i]>max_miss) { ret[i]=-1; break; }
00241                 if (i==j || ret[i]!=ret[j] || ret[i]==-1) continue;
00242                 if (rete[i]>rete[j]) { ret[i]=-1; break; }
00243         }
00244 }
00245 
00246 delete mx0;
00247 delete mx1;
00248 
00249 return ret;
00250 }
00251 
00252 // uses bilinear least-squares to generate a transformation
00253 // matrix for pairs of points
00254 Transform *PointArray::align_2d(PointArray *to,float max_dist) {
00255 vector<int> match=match_points(to,max_dist);
00256 Transform *ret=new Transform();
00257 
00258 // we use bilinear least squares to get 3/6 matrix components
00259 unsigned int i,j;
00260 
00261 vector<float> pts;
00262 for (i=0; i<match.size(); i++) {
00263         if (match[i]==-1) continue;
00264 
00265 //      printf("%d -> %d\n",i,match[i]);
00266         pts.push_back(get_vector_at(i)[0]);
00267         pts.push_back(get_vector_at(i)[1]);
00268         pts.push_back(to->get_vector_at(match[i])[0]);
00269 }
00270 
00271 Vec3f vx=Util::calc_bilinear_least_square(pts);
00272 
00273 // then we get the other 3/6
00274 for (i=j=0; i<match.size(); i++) {
00275         if (match[i]==-1) continue;
00276         pts[j*3]  =get_vector_at(i)[0];
00277         pts[j*3+1]=get_vector_at(i)[1];
00278         pts[j*3+2]=to->get_vector_at(match[i])[1];
00279         j++;
00280 }
00281 
00282 Vec3f vy=Util::calc_bilinear_least_square(pts);
00283 
00284 //ret->set_rotation(vx[1],vx[2],0.0,vy[1],vy[2],0.0,0.0,0.0,1.0);
00285 //ret->set_rotation(vx[1],vy[1],0.0,vx[2],vy[2],0.0,0.0,0.0,1.0);
00286 //ret->set_pretrans(Vec3f(-vx[0],-vy[0],0));
00287 
00288 ret->set(0, 0, vx[1]);
00289 ret->set(0, 1, vy[1]);
00290 ret->set(0, 2, 0.0f);
00291 ret->set(1, 0, vx[2]);
00292 ret->set(1, 1, vy[2]);
00293 ret->set(1, 2, 0.0f);
00294 ret->set(2, 0, 0.0f);
00295 ret->set(2, 1, 0.0f);
00296 ret->set(2, 2, 1.0f);
00297 ret->set_pre_trans(Vec3f(-vx[0],-vy[0],0));
00298 
00299 return ret;
00300 }
00301 
00302 vector<float> PointArray::align_trans_2d(PointArray *to, int flags, float dxhint,float dyhint) {
00303 printf("Warning, this is old code. Use align_2d.\n");
00304 
00305 // returns (dx,dy,residual error,n points used)
00306 // dxhint,dyhint should translate this->to
00307 // flags : 1 - use hint values, 2 - center by strongest point (highest 'value')
00308 int na=   get_number_points();
00309 int nb=to->get_number_points();
00310 if (na<=0 || nb<=0) return vector<float>(4,0);
00311 
00312 int *a2b = (int *)malloc(na*sizeof(int));
00313 int *b2a = (int *)malloc(nb*sizeof(int));
00314 
00315 
00316 // find unweighted centers
00317 float cax,cay,cbx,cby;
00318 int i,j;
00319 
00320 if (flags&1) {
00321         cbx=dxhint;
00322         cby=dyhint;
00323         cax=cay=0;
00324 }
00325 else if (flags&2) {
00326         // find the 'a' list peak
00327         float hia=0.0f;
00328         int hina=0;
00329         for (i=0; i<na; i++) {
00330                 if (get_value_at(i)>hia) { hia=static_cast<float>(get_value_at(i)); hina=i; }
00331         }
00332         cax=get_vector_at(hina)[0];
00333         cay=get_vector_at(hina)[1];
00334 
00335         // find the 'b' list peak
00336         float hib=0;
00337         int hinb=0;
00338         for (i=0; i<na; i++) {
00339                 if (to->get_value_at(i)>hib) { hib=static_cast<float>(to->get_value_at(i)); hinb=i; }
00340         }
00341         cbx=to->get_vector_at(hinb)[0];
00342         cby=to->get_vector_at(hinb)[1];
00343 }
00344 else {
00345         cax=cay=cbx=cby=0;
00346 
00347         for (i=0; i<na; i++) { cax+=get_vector_at(i)[0]; cay+=get_vector_at(i)[1]; }
00348         cax/=(float)na;
00349         cay/=(float)na;
00350 
00351         for (i=0; i<nb; i++) { cbx+=to->get_vector_at(i)[0]; cby+=to->get_vector_at(i)[1]; }
00352         cbx/=(float)nb;
00353         cby/=(float)nb;
00354 }
00355 
00356 Vec3f offset(cbx-cax,cby-cay,0);
00357 
00358 // find the nearest point for each x point, taking the estimated centers into account
00359 for (i=0; i<na; i++) {
00360         float rmin=1.0e30f;
00361         for (j=0; j<nb; j++) {
00362                 float r=(get_vector_at(i)+offset-to->get_vector_at(j)).length();
00363                 if (r<rmin) { a2b[i]=j; rmin=r; }
00364         }
00365 }
00366 
00367 // find the nearest point for each y point
00368 for (i=0; i<nb; i++) {
00369         float rmin=1.0e30f;
00370         for (j=0; j<na; j++) {
00371                 float r=(get_vector_at(j)+offset-to->get_vector_at(i)).length();
00372                 if (r<rmin) { b2a[i]=j; rmin=r; }
00373         }
00374 }
00375 
00376 // now keep only points where x->y matches y->x
00377 for (i=0; i<na; i++) {
00378         if (a2b[i]<0) continue;
00379         if (b2a[a2b[i]]!=i) { printf(" #%d!=%d# ",b2a[a2b[i]],i);  b2a[a2b[i]]=-1; a2b[i]=-1; }
00380         printf("%d->%d  ",i,a2b[i]);
00381 }
00382 printf("\n");
00383 
00384 for (i=0; i<nb; i++) {
00385         if (b2a[i]<0) continue;
00386         if (a2b[b2a[i]]!=i) { a2b[b2a[i]]=-1; b2a[i]=-1; }
00387         printf("%d->%d  ",i,b2a[i]);
00388 }
00389 printf("\n");
00390 
00391 // Compute the average translation required to align the points
00392 float dx=0,dy=0,dr=0,nr=0;
00393 for (i=0; i<na; i++) {
00394         if (a2b[i]==-1) continue;
00395         dx+=to->get_vector_at(a2b[i])[0]-get_vector_at(i)[0];
00396         dy+=to->get_vector_at(a2b[i])[1]-get_vector_at(i)[1];
00397         nr+=1.0;
00398 }
00399 //printf("%f %f %f\n",dx,dy,nr);
00400 if (nr<2) return vector<float>(4,0);
00401 dx/=nr;
00402 dy/=nr;
00403 
00404 // Compute the residual error
00405 for (i=0; i<na; i++) {
00406         if (i==-1  || a2b[i]==-1) continue;
00407         dr+=(to->get_vector_at(a2b[i])-get_vector_at(i)-Vec3f(dx,dy,0)).length();
00408 }
00409 dr/=nr;
00410 
00411 free(a2b);
00412 free(b2a);
00413 vector<float> ret(4);
00414 ret[0]=dx;
00415 ret[1]=dy;
00416 ret[2]=dr;
00417 ret[3]=(float)nr;
00418 return ret;
00419 }
00420 
00421 
00422 bool PointArray::read_from_pdb(const char *file)
00423 {
00424         struct stat filestat;
00425         stat(file, &filestat);
00426         set_number_points(( int)(filestat.st_size / 80 + 1));
00427 
00428         #ifdef DEBUG
00429         printf("PointArray::read_from_pdb(): try %4lu atoms first\n", get_number_points());
00430         #endif
00431 
00432         FILE *fp = fopen(file, "r");
00433         if(!fp) {
00434                 fprintf(stderr,"ERROR in PointArray::read_from_pdb(): cannot open file %s\n",file);
00435                 throw;
00436         }
00437         char s[200];
00438         size_t count = 0;
00439         
00440         while ((fgets(s, 200, fp) != NULL)) {
00441                 if (strncmp(s, "ENDMDL", 6) == 0)
00442                         break;
00443                 if (strncmp(s, "ATOM", 4) != 0)
00444                         continue;
00445 
00446                 if (s[13] == ' ')
00447                         s[13] = s[14];
00448                 if (s[13] == ' ')
00449                         s[13] = s[15];
00450 
00451                 float e = 0;
00452                 char ctt, ctt2 = ' ';
00453                 if (s[13] == ' ')
00454                         ctt = s[14];
00455                 else if (s[12] == ' ') {
00456                         ctt = s[13];
00457                         ctt2 = s[14];
00458                 }
00459                 else {
00460                         ctt = s[12];
00461                         ctt2 = s[13];
00462                 }
00463 
00464                 switch (ctt) {
00465                 case 'H':
00466                         e = 1.0;
00467                         break;
00468                 case 'C':
00469                         e = 6.0;
00470                         break;
00471                 case 'A':
00472                         if (ctt2 == 'U') {
00473                                 e = 79.0;
00474                                 break;
00475                         }
00476                         // treat 'A'mbiguous atoms as N, not perfect, but good enough
00477                 case 'N':
00478                         e = 7.0;
00479                         break;
00480                 case 'O':
00481                         e = 8.0;
00482                         break;
00483                 case 'P':
00484                         e = 15.0;
00485                         break;
00486                 case 'S':
00487                         e = 16.0;
00488                         break;
00489                 case 'W':
00490                         e = 18.0;
00491                         break;                          // ficticious water 'atom'
00492                 default:
00493                         fprintf(stderr, "Unknown atom %c%c\n", ctt, ctt2);
00494                         e = 0;
00495                 }
00496                 if (e == 0)
00497                         continue;
00498 
00499                 float x, y, z, q;
00500                 sscanf(&s[28], " %f %f %f", &x, &y, &z);
00501                 sscanf(&s[60], " %f", &q);
00502 
00503                 if (count + 1 > get_number_points())
00504                         set_number_points(2 * (count + 1));    //makes sure point array is big enough
00505                 
00506 #ifdef DEBUG
00507                 printf("Atom %4lu: x,y,z = %8g,%8g,%8g\te = %g\n", count, x, y, z, e);
00508 #endif
00509                 points[4 * count] = x;
00510                 points[4 * count + 1] = y;
00511                 points[4 * count + 2] = z;
00512                 points[4 * count + 3] = e;
00513                 bfactor[count] = q;
00514                 count++;
00515         }
00516         fclose(fp);
00517         set_number_points(count);
00518         return true;
00519 }
00520 
00521 
00522 void PointArray::save_to_pdb(const char *file)
00523 {
00524         FILE *fp = fopen(file, "w");
00525         for ( size_t i = 0; i < get_number_points(); i++) {
00526                 fprintf(fp, "ATOM  %5lu  CA  ALA A%4lu    %8.3f%8.3f%8.3f%6.2f%6.2f%8s\n", i, i,
00527                                 points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4 * i + 3], 0.0, " ");
00528         }
00529         fclose(fp);
00530 }
00531 
00532 
00533 FloatPoint PointArray::get_center()
00534 {
00535         double xc, yc, zc;
00536         xc = yc = zc = 0.0;
00537         double norm = 0.0;
00538         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00539                 xc += points[i] * points[i + 3];
00540                 yc += points[i + 1] * points[i + 3];
00541                 zc += points[i + 2] * points[i + 3];
00542                 norm += points[i + 3];
00543         }
00544         if (norm <= 0) {
00545                 fprintf(stderr, "Abnormal total value (%g) for PointArray, it should be >0\n", norm);
00546                 return FloatPoint(0, 0, 0);
00547         }
00548         else
00549                 return FloatPoint(xc / norm, yc / norm, zc / norm);
00550 }
00551 
00552 void PointArray::center_to_zero()
00553 {
00554         FloatPoint center = get_center();
00555         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00556                 points[i] -= center[0];
00557                 points[i + 1] -= center[1];
00558                 points[i + 2] -= center[2];
00559         }
00560 }
00561 
00562 Region PointArray::get_bounding_box()
00563 {
00564         double xmin, ymin, zmin;
00565         double xmax, ymax, zmax;
00566         xmin = xmax = points[0];
00567         ymin = ymax = points[1];
00568         zmin = zmax = points[2];
00569         for ( size_t i = 0; i < 4 * get_number_points(); i += 4) {
00570                 if (points[i] > xmax)
00571                         xmax = points[i];
00572                 if (points[i] < xmin)
00573                         xmin = points[i];
00574                 if (points[i + 1] > ymax)
00575                         ymax = points[i + 1];
00576                 if (points[i + 1] < ymin)
00577                         ymin = points[i + 1];
00578                 if (points[i + 2] > zmax)
00579                         zmax = points[i + 2];
00580                 if (points[i + 2] < zmin)
00581                         zmin = points[i + 2];
00582         }
00583         return Region(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin);
00584 }
00585 
00586 
00587 void PointArray::mask(double rmax, double rmin)
00588 {
00589         double rmax2 = rmax * rmax, rmin2 = rmin * rmin;
00590         PointArray *tmp = this->copy();
00591         double *tmp_points = tmp->get_points_array();
00592          int count = 0;
00593         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00594                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v =
00595                         tmp_points[i + 3];
00596                 double r2 = x * x + y * y + z * z;
00597                 if (r2 >= rmin2 && r2 <= rmax2) {
00598                         points[count * 4] = x;
00599                         points[count * 4 + 1] = y;
00600                         points[count * 4 + 2] = z;
00601                         points[count * 4 + 3] = v;
00602                         ++count;
00603                 }
00604         }
00605         set_number_points(count);
00606         if( tmp )
00607         {
00608                 delete tmp;
00609                 tmp = 0;
00610         }
00611 }
00612 
00613 
00614 void PointArray::mask_asymmetric_unit(const string & sym)
00615 {
00616         if (sym == "c1" || sym == "C1")
00617                 return;                                 // do nothing for C1 symmetry
00618         double alt0 = 0, alt1 = M_PI, alt2 = M_PI;
00619         double az0 = 0, az1 = M_PI;
00620         if (sym[0] == 'c' || sym[0] == 'C') {
00621                 int nsym = atoi(sym.c_str() + 1);
00622                 az1 = 2.0 * M_PI / nsym / 2.0;
00623         }
00624         else if (sym[0] == 'd' || sym[0] == 'D') {
00625                 int nsym = atoi(sym.c_str() + 1);
00626                 alt1 = M_PI / 2.0;
00627                 alt2 = alt1;
00628                 az1 = 2.0 * M_PI / nsym / 2.0;
00629         }
00630         else if (sym == "icos" || sym == "ICOS") {
00631                 alt1 = 0.652358139784368185995; // 5fold to 3fold
00632                 alt2 = 0.55357435889704525151;  // half of edge ie. 5fold to 2fold along the edge
00633                 az1 = 2.0 * M_PI / 5 / 2.0;
00634         }
00635         else {
00636                 LOGERR("PointArray::set_to_asymmetric_unit(): sym = %s is not implemented yet",
00637                            sym.c_str());
00638                 return;
00639         }
00640 #ifdef DEBUG
00641         printf("Sym %s: alt0 = %8g\talt1 = %8g\talt2 = %8g\taz0 = %8g\taz1 = %8g\n", sym.c_str(), alt0*180.0/M_PI, alt1*180.0/M_PI, alt2*180.0/M_PI, az0*180.0/M_PI, az1*180.0/M_PI);
00642 #endif
00643 
00644         PointArray *tmp = this->copy();
00645         double *tmp_points = tmp->get_points_array();
00646          int count = 0;
00647         for ( size_t i = 0; i < 4 * tmp->get_number_points(); i += 4) {
00648                 double x = tmp_points[i], y = tmp_points[i + 1], z = tmp_points[i + 2], v = tmp_points[i + 3];
00649                 double az = atan2(y, x);
00650                 double az_abs = fabs(az - az0);
00651                 if (az_abs < (az1 - az0)) {
00652                         double alt_max = alt1 + (alt2 - alt1) * az_abs / (az1 - az0);
00653                         double alt = acos(z / sqrt(x * x + y * y + z * z));
00654                         if (alt < alt_max && alt >= alt0) {
00655 #ifdef DEBUG
00656                                 printf("Point %3lu: x,y,z = %8g,%8g,%8g\taz = %8g\talt = %8g\n",i/4,x,y,z,az*180.0/M_PI, alt*180.0/M_PI);
00657 #endif
00658                                 points[count * 4] = x;
00659                                 points[count * 4 + 1] = y;
00660                                 points[count * 4 + 2] = z;
00661                                 points[count * 4 + 3] = v;
00662                                 count++;
00663                         }
00664                 }
00665         }
00666         set_number_points(count);
00667         if( tmp )
00668         {
00669                 delete tmp;
00670                 tmp = 0;
00671         }
00672 }
00673 
00674 vector<float> PointArray::get_points() {
00675 vector<float> ret;
00676 for (unsigned int i=0; i<n; i++) {
00677         ret.push_back((float)points[i*4]);
00678         ret.push_back((float)points[i*4+1]);
00679         ret.push_back((float)points[i*4+2]);
00680 }
00681 
00682 return ret;
00683 }
00684 
00685 void PointArray::transform(const Transform& xf) {
00686 
00687         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00688                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00689                 v=v*xf;
00690                 points[i]  =v[0];
00691                 points[i+1]=v[1];
00692                 points[i+2]=v[2];
00693         }
00694 
00695 }
00696 
00697 void PointArray::right_transform(const Transform& transform) {
00698         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00699                 Transform s = transform.transpose();
00700                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00701                 v= s*v;
00702                 points[i]  =v[0];
00703                 points[i+1]=v[1];
00704                 points[i+2]=v[2];
00705         }
00706 
00707 }
00708 void PointArray::set_from(PointArray * source, const string & sym, Transform *transform)
00709 {
00710         set_from(source->get_points_array(), source->get_number_points(), sym, transform);
00711 
00712 }
00713 
00714 void PointArray::set_from(double *src,  int num, const string & sym, Transform *xform)
00715 {
00716          int nsym = xform->get_nsym(sym);
00717 
00718         if (get_number_points() != (size_t)nsym * num)
00719                 set_number_points((size_t)nsym * num);
00720 
00721         double *target = get_points_array();
00722 
00723         for ( int s = 0; s < nsym; s++) {
00724                 int index = s * 4 * num;
00725                 for ( int i = 0; i < 4 * num; i += 4, index += 4) {
00726                         Vec3f v((float)src[i],(float)src[i+1],(float)src[i+2]);
00727                         v=v*xform->get_sym(sym,s);
00728                         target[index]  =v[0];
00729                         target[index+1]=v[1];
00730                         target[index+2]=v[2];
00731                         target[index+3]=src[i+3];
00732                 }
00733         }
00734 }
00735 
00736 void PointArray::set_from(vector<float> pts) {
00737         set_number_points(pts.size()/4);
00738         for (unsigned int i=0; i<pts.size(); i++) points[i]=pts[i];
00739 
00740 }
00741 
00742 void PointArray::set_from_density_map(EMData * map, int num, float thresh, float apix,
00743                                                                           Density2PointsArrayAlgorithm mode)
00744 {
00745         if (mode == PEAKS_SUB || mode == PEAKS_DIV) {
00746                 // find out how many voxels are useful voxels
00747                 int num_voxels = 0;
00748                 int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00749                 size_t size = (size_t)nx * ny * nz;
00750                 EMData *tmp_map = map->copy();
00751                 float *pd = tmp_map->get_data();
00752                 for (size_t i = 0; i < size; ++i) {
00753                         if (pd[i] > thresh)
00754                                 ++num_voxels;
00755                 }
00756 
00757                 double pointvol = double (num_voxels) / double (num);
00758                 double gauss_real_width = pow(pointvol, 1. / 3.);       // in pixels
00759 #ifdef DEBUG
00760                 printf("Average point range radius = %g pixels for %d points from %d used voxels\n",
00761                            gauss_real_width, num, num_voxels);
00762 #endif
00763 
00764                 double min_table_val = 1e-4;
00765                 double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
00766 
00767                 double table_step_size = 1.;    // number of steps for each pixel
00768                 double inv_table_step_size = 1.0 / table_step_size;
00769                 int table_size = int (max_table_x * gauss_real_width / (table_step_size) * 1.25) + 1;
00770                 double *table = (double *) malloc(sizeof(double) * table_size);
00771                 for (int i = 0; i < table_size; i++) {
00772                         double x = i * table_step_size / gauss_real_width;
00773                         table[i] = exp(-x * x);
00774                 }
00775 
00776                 int gbox = int (max_table_x * gauss_real_width);        // local box half size in pixels to consider for each point
00777                 if (gbox <= 0)
00778                         gbox = 1;
00779 
00780                 set_number_points(num);
00781                 for (int count = 0; count < num; count++) {
00782                         float cmax = pd[0];
00783                         int cmaxpos = 0;
00784                         for (size_t i = 0; i < size; ++i) {
00785                                 if (pd[i] > cmax) {
00786                                         cmax = pd[i];
00787                                         cmaxpos = i;
00788                                 }
00789                         }
00790                         int iz = cmaxpos / (nx * ny), iy = (cmaxpos - iz * nx * ny) / nx, ix =
00791                                 cmaxpos - iz * nx * ny - iy * nx;
00792 
00793                         // update coordinates in pixels
00794                         points[4*count  ] = ix;
00795                         points[4*count+1] = iy;
00796                         points[4*count+2] = iz;
00797                         points[4*count+3] = cmax;
00798 #ifdef DEBUG
00799                         printf("Point %d: val = %g\tat  %d, %d, %d\n", count, cmax, ix, iy, iz);
00800 #endif
00801 
00802                         int imin = ix - gbox, imax = ix + gbox;
00803                         int jmin = iy - gbox, jmax = iy + gbox;
00804                         int kmin = iz - gbox, kmax = iz + gbox;
00805                         if (imin < 0)
00806                                 imin = 0;
00807                         if (jmin < 0)
00808                                 jmin = 0;
00809                         if (kmin < 0)
00810                                 kmin = 0;
00811                         if (imax > nx)
00812                                 imax = nx;
00813                         if (jmax > ny)
00814                                 jmax = ny;
00815                         if (kmax > nz)
00816                                 kmax = nz;
00817 
00818                         for (int k = kmin; k < kmax; ++k) {
00819                                 int table_index_z = int (fabs(double (k - iz)) * inv_table_step_size);
00820                                 double zval = table[table_index_z];
00821                                 size_t pd_index_z = k * nx * ny;
00822                                 //printf("k = %8d\tx = %8g\tval = %8g\n", k, float(k-iz), zval);
00823                                 for (int j = jmin; j < jmax; ++j) {
00824                                         int table_index_y = int (fabs(double (j - iy)) * inv_table_step_size);
00825                                         double yval = table[table_index_y];
00826                                         size_t pd_index = pd_index_z + j * nx + imin;
00827                                         for (int i = imin; i < imax; ++i, ++pd_index) {
00828                                                 int table_index_x = int (fabs(double (i - ix)) * inv_table_step_size);
00829                                                 double xval = table[table_index_x];
00830                                                 if (mode == PEAKS_SUB)
00831                                                         pd[pd_index] -= (float)(cmax * zval * yval * xval);
00832                                                 else
00833                                                         pd[pd_index] *= (float)(1.0 - zval * yval * xval);      // mode == PEAKS_DIV
00834                                         }
00835                                 }
00836                         }
00837                 }
00838                 set_number_points(num);
00839                 tmp_map->update();
00840                 if( tmp_map )
00841                 {
00842                         delete tmp_map;
00843                         tmp_map = 0;
00844                 }
00845         }
00846         else if (mode == KMEANS) {
00847                 set_number_points(num);
00848                 zero();
00849 
00850                 PointArray tmp_pa;
00851                 tmp_pa.set_number_points(num);
00852                 tmp_pa.zero();
00853 
00854                  int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00855                 float *pd = map->get_data();
00856 
00857                 // initialize segments with random centers at pixels with values > thresh
00858 #ifdef DEBUG
00859                 printf("Start initial random seeding\n");
00860 #endif
00861                 for ( size_t i = 0; i < get_number_points(); i++) {
00862                          int x, y, z;
00863                         double v;
00864                         do {
00865                                 x = ( int) Util::get_frand(0, nx - 1);
00866                                 y = ( int) Util::get_frand(0, ny - 1);
00867                                 z = ( int) Util::get_frand(0, nz - 1);
00868                                 v = pd[z * nx * ny + y * nx + x];
00869 #ifdef DEBUG
00870                                 printf("Trying Point %lu: val = %g\tat  %d, %d, %d\tfrom map (%d,%d,%d)\n", i, v, x,
00871                                            y, z, nx, ny, nz);
00872 #endif
00873                         } while (v <= thresh);
00874                         points[4 * i] = (double) x;
00875                         points[4 * i + 1] = (double) y;
00876                         points[4 * i + 2] = (double) z;
00877                         points[4 * i + 3] = (double) v;
00878 #ifdef DEBUG
00879                         printf("Point %lu: val = %g\tat  %g, %g, %g\n", i, points[4 * i + 3], points[4 * i],
00880                                    points[4 * i + 1], points[4 * i + 2]);
00881 #endif
00882                 }
00883 
00884                 double min_dcen = 1e0;  // minimal mean segment center shift as convergence criterion
00885                 double dcen = 0.0;
00886                 int iter = 0, max_iter = 100;
00887                 do {
00888 #ifdef DEBUG
00889                         printf("Iteration %3d, start\n", iter);
00890 #endif
00891                         double *tmp_points = tmp_pa.get_points_array();
00892 
00893                         // reassign each pixel to the best segment
00894                         for ( int k = 0; k < nz; k++) {
00895                                 for ( int j = 0; j < ny; j++) {
00896                                         for ( int i = 0; i < nx; i++) {
00897                                                 size_t idx = k * nx * ny + j * nx + i;
00898                                                 if (pd[idx] > thresh) {
00899                                                         double min_dist = 1e60; // just a large distance
00900                                                          int min_s = 0;
00901                                                         for ( size_t s = 0; s < get_number_points(); ++s) {
00902                                                                 double x = points[4 * s];
00903                                                                 double y = points[4 * s + 1];
00904                                                                 double z = points[4 * s + 2];
00905                                                                 double dist =
00906                                                                         (k - z) * (k - z) + (j - y) * (j - y) + (i - x) * (i - x);
00907                                                                 if (dist < min_dist) {
00908                                                                         min_dist = dist;
00909                                                                         min_s = s;
00910                                                                 }
00911                                                         }
00912                                                         tmp_points[4 * min_s] += i;
00913                                                         tmp_points[4 * min_s + 1] += j;
00914                                                         tmp_points[4 * min_s + 2] += k;
00915                                                         tmp_points[4 * min_s + 3] += 1.0;
00916                                                 }
00917                                         }
00918                                 }
00919                         }
00920 #ifdef DEBUG
00921                         printf("Iteration %3d, finished reassigning segments\n", iter);
00922 #endif
00923                         // update each segment's center
00924                         dcen = 0.0;
00925                         for ( size_t s = 0; s < get_number_points(); ++s) {
00926                                 if (tmp_points[4 * s + 3]) {
00927                                         tmp_points[4 * s] /= tmp_points[4 * s + 3];
00928                                         tmp_points[4 * s + 1] /= tmp_points[4 * s + 3];
00929                                         tmp_points[4 * s + 2] /= tmp_points[4 * s + 3];
00930 #ifdef DEBUG
00931                                         printf("Iteration %3d, Point %3lu at %8g, %8g, %8g -> %8g, %8g, %8g\n", iter, s,
00932                                                    points[4 * s], points[4 * s + 1], points[4 * s + 2], tmp_points[4 * s],
00933                                                    tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00934 #endif
00935                                 }
00936                                 else {                  // empty segments are reseeded
00937                                          int x, y, z;
00938                                         double v;
00939                                         do {
00940                                                 x = ( int) Util::get_frand(0, nx - 1);
00941                                                 y = ( int) Util::get_frand(0, ny - 1);
00942                                                 z = ( int) Util::get_frand(0, nz - 1);
00943                                                 v = pd[z * nx * ny + y * nx + x];
00944                                         } while (v <= thresh);
00945                                         tmp_points[4 * s] = (double) x;
00946                                         tmp_points[4 * s + 1] = (double) y;
00947                                         tmp_points[4 * s + 2] = (double) z;
00948                                         tmp_points[4 * s + 3] = (double) v;
00949 #ifdef DEBUG
00950                                         printf
00951                                                 ("Iteration %3d, Point %3lu reseeded from %8g, %8g, %8g -> %8g, %8g, %8g\n",
00952                                                  iter, s, points[4 * s], points[4 * s + 1], points[4 * s + 2],
00953                                                  tmp_points[4 * s], tmp_points[4 * s + 1], tmp_points[4 * s + 2]);
00954 #endif
00955                                 }
00956                                 double dx = tmp_points[4 * s] - points[4 * s];
00957                                 double dy = tmp_points[4 * s + 1] - points[4 * s + 1];
00958                                 double dz = tmp_points[4 * s + 2] - points[4 * s + 2];
00959                                 dcen += dx * dx + dy * dy + dz * dz;
00960                         }
00961                         dcen = sqrt(dcen / get_number_points());
00962                         //swap pointter, faster but risky
00963 #ifdef DEBUG
00964                         printf("before swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
00965                                    (long)tmp_pa.get_points_array());
00966 #endif
00967                         double *tp = get_points_array();
00968                         set_points_array(tmp_points);
00969                         tmp_pa.set_points_array(tp);
00970                         tmp_pa.zero();
00971 #ifdef DEBUG
00972                         printf("after  swap: points = %ld\ttmp_points = %ld\n", (long)get_points_array(),
00973                                    (long)tmp_pa.get_points_array());
00974                         printf("Iteration %3d, finished updating segment centers with dcen = %g pixels\n", iter,
00975                                    dcen);
00976 #endif
00977 
00978                         iter++;
00979                 } while (dcen > min_dcen && iter <= max_iter);
00980                 map->update();
00981 
00982                 sort_by_axis(2);        // x,y,z axes = 0, 1, 2
00983         }
00984         else {
00985                 LOGERR("PointArray::set_from_density_map(): mode = %d is not implemented yet", mode);
00986         }
00987         //update to use apix and origin
00988         int nx = map->get_xsize(), ny = map->get_ysize(), nz = map->get_zsize();
00989         float origx, origy, origz;
00990         try {
00991                 origx = map->get_attr("origin_x");
00992                 origy = map->get_attr("origin_y");
00993                 origz = map->get_attr("origin_z");
00994         }
00995         catch(...) {
00996                 origx = -nx / 2 * apix;
00997                 origy = -ny / 2 * apix;
00998                 origz = -nz / 2 * apix;
00999         }
01000 
01001 #ifdef DEBUG
01002         printf("Apix = %g\torigin x,y,z = %8g,%8g,%8g\n",apix, origx, origy, origz);
01003 #endif
01004 
01005         float *pd = map->get_data();
01006         for ( size_t i = 0; i < get_number_points(); ++i) {
01007 #ifdef DEBUG
01008                 printf("Point %4lu: x,y,z,v = %8g,%8g,%8g,%8g",i, points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
01009 #endif
01010                 points[4 * i + 3] =
01011                         pd[(int) points[4 * i + 2] * nx * ny + (int) points[4 * i + 1] * nx +
01012                            (int) points[4 * i]];
01013                 points[4 * i] = points[4 * i] * apix + origx;
01014                 points[4 * i + 1] = points[4 * i + 1] * apix + origy;
01015                 points[4 * i + 2] = points[4 * i + 2] * apix + origz;
01016 #ifdef DEBUG
01017                 printf("\t->\t%8g,%8g,%8g,%8g\n",points[4 * i],points[4 * i + 1],points[4 * i + 2],points[4 * i + 3]);
01018 #endif
01019         }
01020         map->update();
01021 }
01022 
01023 
01024 
01025 
01026 void PointArray::sort_by_axis(int axis)
01027 {
01028         if (axis == 0)
01029                 qsort(points, n, sizeof(double) * 4, cmp_axis_x);
01030         else if (axis == 1)
01031                 qsort(points, n, sizeof(double) * 4, cmp_axis_y);
01032         else if (axis == 2)
01033                 qsort(points, n, sizeof(double) * 4, cmp_axis_z);
01034         else
01035                 qsort(points, n, sizeof(double) * 4, cmp_val);
01036 }
01037 
01038 
01039 EMData *PointArray::pdb2mrc_by_summation(int map_size, float apix, float res, int addpdbbfactor)
01040 {
01041 #ifdef DEBUG
01042         printf("PointArray::pdb2mrc_by_summation(): %lu points\tmapsize = %4d\tapix = %g\tres = %g\n",get_number_points(),map_size, apix, res);
01043 #endif
01044         //if ( gauss_real_width < apix) LOGERR("PointArray::projection_by_summation(): apix(%g) is too large for resolution (%g Angstrom in Fourier space) with %g pixels of 1/e half width", apix, res, gauss_real_width);
01045 
01046         double min_table_val = 1e-7;
01047         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01048 
01049         double table_step_size = 0.001; // number of steps for each pixel
01050         double inv_table_step_size = 1.0 / table_step_size;
01051         
01052 //      sort_by_axis(2);                        // sort by Z-axis
01053 
01054         EMData *map = new EMData();
01055         map->set_size(map_size, map_size, map_size);
01056         map->to_zero();
01057         float *pd = map->get_data();
01058         
01059         vector<double> table;
01060         double gauss_real_width;
01061         int table_size;
01062         int gbox;
01063         
01064         
01065         for ( size_t s = 0; s < get_number_points(); ++s) {
01066                 double xc = points[4 * s] / apix + map_size / 2;
01067                 double yc = points[4 * s + 1] / apix + map_size / 2;
01068                 double zc = points[4 * s + 2] / apix + map_size / 2;
01069                 double fval = points[4 * s + 3];
01070                 
01071                 //printf("\n bfactor=%f",bfactor[s]);
01072                 
01073                 
01074                 
01075                 
01076                 
01077                 if(addpdbbfactor==-1){
01078                         gauss_real_width = res/M_PI;    // in Angstrom, res is in Angstrom
01079                 }else{
01080                         gauss_real_width = (bfactor[s])/(4*sqrt(2.0)*M_PI);     // in Angstrom, res is in Angstrom
01081                 }
01082                 
01083                 
01084                 table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01085                 table.resize(table_size);
01086                 for (int i = 0; i < table_size; i++){
01087                         table[i] = 0;
01088                 }
01089                 
01090                 for (int i = 0; i < table_size; i++) {                                          
01091                         double x = -i * table_step_size * apix / gauss_real_width;
01092                         if(addpdbbfactor==-1){
01093                                 table[i] =  exp(-x * x);        
01094                         }
01095                         else{
01096                         table[i] =  exp(-x * x)/sqrt(gauss_real_width * gauss_real_width * 2 * M_PI);
01097                         }
01098                 }
01099 
01100                 gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
01101                 if (gbox <= 0)
01102                         gbox = 1;
01103                 
01104                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01105                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01106                 int kmin = int (zc) - gbox, kmax = int (zc) + gbox;
01107                 if (imin < 0)
01108                         imin = 0;
01109                 if (jmin < 0)
01110                         jmin = 0;
01111                 if (kmin < 0)
01112                         kmin = 0;
01113                 if (imax > map_size)
01114                         imax = map_size;
01115                 if (jmax > map_size)
01116                         jmax = map_size;
01117                 if (kmax > map_size)
01118                         kmax = map_size;                
01119                 
01120                 for (int k = kmin; k < kmax; k++) {
01121                         size_t table_index_z = size_t (fabs(k - zc) * inv_table_step_size);
01122                         if ( table_index_z >= table.size() ) continue;
01123                         double zval = table[table_index_z];
01124                         size_t pd_index_z = k * map_size * map_size;
01125                         
01126                         for (int j = jmin; j < jmax; j++) {
01127                                 size_t table_index_y = size_t (fabs(j - yc) * inv_table_step_size);
01128                                 if ( table_index_y >= table.size() ) continue;
01129                                 double yval = table[table_index_y];
01130                                 size_t pd_index = pd_index_z + j * map_size + imin;
01131                                 for (int i = imin; i < imax; i++, pd_index++) {
01132                                         size_t table_index_x = size_t (fabs(i - xc) * inv_table_step_size);
01133                                         if ( table_index_x >= table.size() ) continue;
01134                                         double xval = table[table_index_x];
01135                                         pd[pd_index] += (float) (fval * zval * yval * xval);
01136                                 }
01137                         }
01138                 }
01139         }
01140 
01141         map->update();
01142         map->set_attr("apix_x", apix);
01143         map->set_attr("apix_y", apix);
01144         map->set_attr("apix_z", apix);
01145         map->set_attr("origin_x", -map_size/2*apix);
01146         map->set_attr("origin_y", -map_size/2*apix);
01147         map->set_attr("origin_z", -map_size/2*apix);
01148 
01149         return map;
01150 }
01151 
01152 
01153 EMData *PointArray::projection_by_summation(int image_size, float apix, float res)
01154 {
01155         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01156         //if ( gauss_real_width < apix) LOGERR("PointArray::projection_by_summation(): apix(%g) is too large for resolution (%g Angstrom in Fourier space) with %g pixels of 1/e half width", apix, res, gauss_real_width);
01157 
01158         double min_table_val = 1e-7;
01159         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01160 
01161         //double table_step_size = 0.001;    // number of steps for x=[0,1] in exp(-x*x)
01162         //int table_size = int(max_table_x / table_step_size *1.25);
01163         //double* table = (double*)malloc(sizeof(double) * table_size);
01164         //for(int i=0; i<table_size; i++) table[i]=exp(-i*i*table_step_size*table_step_size);
01165 
01166         double table_step_size = 0.001; // number of steps for each pixel
01167         double inv_table_step_size = 1.0 / table_step_size;
01168         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01169         double *table = (double *) malloc(sizeof(double) * table_size);
01170         for (int i = 0; i < table_size; i++) {
01171                 double x = -i * table_step_size * apix / gauss_real_width;
01172                 table[i] = exp(-x * x);
01173         }
01174 
01175         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01176         if (gbox <= 0)
01177                 gbox = 1;
01178         EMData *proj = new EMData();
01179         proj->set_size(image_size, image_size, 1);
01180         proj->to_zero();
01181         float *pd = proj->get_data();
01182         for ( size_t s = 0; s < get_number_points(); ++s) {
01183                 double xc = points[4 * s] / apix + image_size / 2;
01184                 double yc = points[4 * s + 1] / apix + image_size / 2;
01185                 double fval = points[4 * s + 3];
01186                 int imin = int (xc) - gbox, imax = int (xc) + gbox;
01187                 int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01188                 if (imin < 0)
01189                         imin = 0;
01190                 if (jmin < 0)
01191                         jmin = 0;
01192                 if (imax > image_size)
01193                         imax = image_size;
01194                 if (jmax > image_size)
01195                         jmax = image_size;
01196 
01197                 for (int j = jmin; j < jmax; j++) {
01198                         //int table_index_y = int(fabs(j-yc)*apix/gauss_real_width/table_step_size);
01199                         int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01200                         double yval = table[table_index_y];
01201 #ifdef DEBUG
01202                         //double yval2 = exp( - (j-yc)*(j-yc)*apix*apix/(gauss_real_width*gauss_real_width));
01203                         //if(fabs(yval2-yval)/yval2>1e-2) printf("s=%d\txc,yc=%g,%g\tyval,yval2=%g,%g\tdiff=%g\n",s,xc,yc,yval,yval2,fabs(yval2-yval)/yval2);
01204 #endif
01205                         int pd_index = j * image_size + imin;
01206                         for (int i = imin; i < imax; i++, pd_index++) {
01207                                 //int table_index_x = int(fabs(i-xc)*apix/gauss_real_width/table_step_size);
01208                                 int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01209                                 double xval = table[table_index_x];
01210 #ifdef DEBUG
01211                                 //double xval2 = exp( - (i-xc)*(i-xc)*apix*apix/(gauss_real_width*gauss_real_width));
01212                                 //if(fabs(xval2-xval)/xval2>1e-2) printf("\ts=%d\txc,yc=%g,%g\txval,xval2=%g,%g\tdiff=%g\n",s,xc,yc,xval,xval2,fabs(xval2-xval)/xval2);
01213 #endif
01214                                 pd[pd_index] += (float)(fval * yval * xval);
01215                         }
01216                 }
01217         }
01218         for (int i = 0; i < image_size * image_size; i++)
01219                 pd[i] /= sqrt(M_PI);
01220         proj->update();
01221         return proj;
01222 }
01223 
01224 void PointArray::replace_by_summation(EMData *proj, int ind, Vec3f vec, float amp, float apix, float res)
01225 {
01226         double gauss_real_width = res / (M_PI); // in Angstrom, res is in Angstrom
01227 
01228         double min_table_val = 1e-7;
01229         double max_table_x = sqrt(-log(min_table_val)); // for exp(-x*x)
01230 
01231         double table_step_size = 0.001; // number of steps for each pixel
01232         double inv_table_step_size = 1.0 / table_step_size;
01233         int table_size = int (max_table_x * gauss_real_width / (apix * table_step_size) * 1.25);
01234         double *table = (double *) malloc(sizeof(double) * table_size);
01235         for (int i = 0; i < table_size; i++) {
01236                 double x = -i * table_step_size * apix / gauss_real_width;
01237                 table[i] = exp(-x * x)/pow((float)M_PI,.25f);
01238         }
01239         int image_size=proj->get_xsize();
01240 
01241         // subtract the old point
01242         int gbox = int (max_table_x * gauss_real_width / apix); // local box half size in pixels to consider for each point
01243         if (gbox <= 0)
01244                 gbox = 1;
01245         float *pd = proj->get_data();
01246          int s = ind;
01247         double xc = points[4 * s] / apix + image_size / 2;
01248         double yc = points[4 * s + 1] / apix + image_size / 2;
01249         double fval = points[4 * s + 3];
01250         int imin = int (xc) - gbox, imax = int (xc) + gbox;
01251         int jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01252 
01253         if (imin < 0) imin = 0;
01254         if (jmin < 0) jmin = 0;
01255         if (imax > image_size) imax = image_size;
01256         if (jmax > image_size) jmax = image_size;
01257 
01258         for (int j = jmin; j < jmax; j++) {
01259                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01260                 double yval = table[table_index_y];
01261                 int pd_index = j * image_size + imin;
01262                 for (int i = imin; i < imax; i++, pd_index++) {
01263                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01264                         double xval = table[table_index_x];
01265                         pd[pd_index] -= (float)(fval * yval * xval);
01266                 }
01267         }
01268 
01269         // add the new point
01270         gbox = int (max_table_x * gauss_real_width / apix);     // local box half size in pixels to consider for each point
01271         if (gbox <= 0)
01272                 gbox = 1;
01273         pd = proj->get_data();
01274         s = ind;
01275         xc = vec[0] / apix + image_size / 2;
01276         yc = vec[1] / apix + image_size / 2;
01277         fval = amp;
01278         imin = int (xc) - gbox, imax = int (xc) + gbox;
01279         jmin = int (yc) - gbox, jmax = int (yc) + gbox;
01280 
01281         if (imin < 0) imin = 0;
01282         if (jmin < 0) jmin = 0;
01283         if (imax > image_size) imax = image_size;
01284         if (jmax > image_size) jmax = image_size;
01285 
01286         for (int j = jmin; j < jmax; j++) {
01287                 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01288                 double yval = table[table_index_y];
01289                 int pd_index = j * image_size + imin;
01290                 for (int i = imin; i < imax; i++, pd_index++) {
01291                         int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01292                         double xval = table[table_index_x];
01293                         pd[pd_index] -= (float)(fval * yval * xval);
01294                 }
01295         }
01296 
01297         proj->update();
01298         return;
01299 }
01300 
01301 
01302 EMData *PointArray::pdb2mrc_by_nfft(int , float , float )
01303 {
01304 #if defined NFFT
01305         nfft_3D_plan my_plan;           // plan for the nfft
01306 
01308         nfft_3D_init(&my_plan, map_size, get_number_points());
01309 
01311         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01312                 // FFTW and nfft use row major array layout, EMAN uses column major
01313                 my_plan.v[3 * j + 2] = (fftw_real) (points[i] / (apix * map_size));
01314                 my_plan.v[3 * j + 1] = (fftw_real) (points[i + 1] / (apix * map_size));
01315                 my_plan.v[3 * j] = (fftw_real) (points[i + 2] / (apix * map_size));
01316                 my_plan.f[j].re = (fftw_real) (points[i + 3]);
01317                 my_plan.f[j].im = 0.0;
01318         }
01319 
01321         if (my_plan.nfft_flags & PRE_PSI) {
01322                 nfft_3D_precompute_psi(&my_plan);
01323         }
01324 
01325         // compute the uniform Fourier transform
01326         nfft_3D_transpose(&my_plan);
01327 
01328         // copy the Fourier transform to EMData data array
01329         EMData *fft = new EMData();
01330         fft->set_size(map_size + 2, map_size, map_size);
01331         fft->set_complex(true);
01332         fft->set_ri(true);
01333         fft->to_zero();
01334         float *data = fft->get_data();
01335         double norm = 1.0 / (map_size * map_size * map_size);
01336         for (int k = 0; k < map_size; k++) {
01337                 for (int j = 0; j < map_size; j++) {
01338                         for (int i = 0; i < map_size / 2; i++) {
01339                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01340                                         (float) (my_plan.
01341                                                          f_hat[k * map_size * map_size + j * map_size + i +
01342                                                                    map_size / 2].re) * norm;
01343                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01344                                         (float) (my_plan.
01345                                                          f_hat[k * map_size * map_size + j * map_size + i +
01346                                                                    map_size / 2].im) * norm * (-1.0);
01347                         }
01348                 }
01349         }
01351         nfft_3D_finalize(&my_plan);
01352 
01353         // low pass processor
01354         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01355         int index = 0;
01356         for (int k = 0; k < map_size; k++) {
01357                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01358                 for (int j = 0; j < map_size; j++) {
01359                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01360                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01361                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01362                                 data[index] *= val;
01363                                 data[index + 1] *= val;
01364                         }
01365                 }
01366         }
01367         fft->update();
01368         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
01369 
01370         fft->process_inplace("xform.phaseorigin.tocorner");     // move phase origin to center of image map_size, instead of at corner
01371         EMData *map = fft->do_ift();
01372         map->set_attr("apix_x", apix);
01373         map->set_attr("apix_y", apix);
01374         map->set_attr("apix_z", apix);
01375         map->set_attr("origin_x", -map_size/2*apix);
01376         map->set_attr("origin_y", -map_size/2*apix);
01377         map->set_attr("origin_z", -map_size/2*apix);
01378         if( fft )
01379         {
01380                 delete fft;
01381                 fft = 0;
01382         }
01383         return map;
01384 #elif defined NFFT2
01385         nfft_plan my_plan;                      // plan for the nfft
01386 
01388         nfft_init_3d(&my_plan, map_size, map_size, map_size, get_number_points());
01389 
01391         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01392                 // FFTW and nfft use row major array layout, EMAN uses column major
01393                 my_plan.x[3 * j + 2] = (double) (points[i] / (apix * map_size));
01394                 my_plan.x[3 * j + 1] = (double) (points[i + 1] / (apix * map_size));
01395                 my_plan.x[3 * j] = (double) (points[i + 2] / (apix * map_size));
01396                 my_plan.f[j][0] = (double) (points[i + 3]);
01397                 my_plan.f[j][1] = 0.0;
01398         }
01399 
01401         if (my_plan.nfft_flags & PRE_PSI) {
01402                 nfft_precompute_psi(&my_plan);
01403                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01404                         nfft_full_psi(&my_plan, pow(10, -10));
01405         }
01406 
01407         // compute the uniform Fourier transform
01408         nfft_transposed(&my_plan);
01409 
01410         // copy the Fourier transform to EMData data array
01411         EMData *fft = new EMData();
01412         fft->set_size(map_size + 2, map_size, map_size);
01413         fft->set_complex(true);
01414         fft->set_ri(true);
01415         fft->to_zero();
01416         float *data = fft->get_data();
01417         double norm = 1.0 / (map_size * map_size * map_size);
01418         for (int k = 0; k < map_size; k++) {
01419                 for (int j = 0; j < map_size; j++) {
01420                         for (int i = 0; i < map_size / 2; i++) {
01421                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i] =
01422                                         (float) (my_plan.
01423                                                          f_hat[k * map_size * map_size + j * map_size + i +
01424                                                                    map_size / 2][0]) * norm;
01425                                 data[k * map_size * (map_size + 2) + j * (map_size + 2) + 2 * i + 1] =
01426                                         (float) (my_plan.
01427                                                          f_hat[k * map_size * map_size + j * map_size + i +
01428                                                                    map_size / 2][1]) * norm;
01429                         }
01430                 }
01431         }
01433         nfft_finalize(&my_plan);
01434 
01435         // low pass processor
01436         double sigma2 = (map_size * apix / res) * (map_size * apix / res);
01437         int index = 0;
01438         for (int k = 0; k < map_size; k++) {
01439                 double RZ2 = (k - map_size / 2) * (k - map_size / 2);
01440                 for (int j = 0; j < map_size; j++) {
01441                         double RY2 = (j - map_size / 2) * (j - map_size / 2);
01442                         for (int i = 0; i < map_size / 2 + 1; i++, index += 2) {
01443                                 float val = exp(-(i * i + RY2 + RZ2) / sigma2);
01444                                 data[index] *= val;
01445                                 data[index + 1] *= val;
01446                         }
01447                 }
01448         }
01449         fft->update();
01450         //fft->process_inplace(".filter.lowpass.gauss",Dict("cutoff_abs", map_size*apix/res));
01451 
01452         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image map_size, instead of at corner
01453         EMData *map = fft->do_ift();
01454         map->set_attr("apix_x", apix);
01455         map->set_attr("apix_y", apix);
01456         map->set_attr("apix_z", apix);
01457         map->set_attr("origin_x", -map_size / 2 * apix);
01458         map->set_attr("origin_y", -map_size / 2 * apix);
01459         map->set_attr("origin_z", -map_size / 2 * apix);
01460         if( fft )
01461         {
01462                 delete fft;
01463                 fft = 0;
01464         }
01465         return map;
01466 #else
01467         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01468         return 0;
01469 #endif
01470 }
01471 
01472 EMData *PointArray::projection_by_nfft(int , float , float )
01473 {
01474 #if defined NFFT
01475         nfft_2D_plan my_plan;           // plan for the nfft
01476         int N[2], n[2];
01477         N[0] = image_size;
01478         n[0] = next_power_of_2(2 * image_size);
01479         N[1] = image_size;
01480         n[1] = next_power_of_2(2 * image_size);
01481 
01483         nfft_2D_init(&my_plan, image_size, get_number_points());
01484         //nfft_2D_init_specific(&my_plan, N, get_number_points(), n, 3,
01485         //                 PRE_PHI_HUT | PRE_PSI,
01486         //                 FFTW_ESTIMATE | FFTW_OUT_OF_PLACE);
01488         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01489                 // FFTW and nfft use row major array layout, EMAN uses column major
01490                 my_plan.v[2 * j + 1] = (fftw_real) (points[i] / (apix * image_size));
01491                 my_plan.v[2 * j] = (fftw_real) (points[i + 1] / (apix * image_size));
01492                 my_plan.f[j].re = (fftw_real) points[i + 3];
01493                 my_plan.f[j].im = 0.0;
01494         }
01495 
01497         if (my_plan.nfft_flags & PRE_PSI) {
01498                 nfft_2D_precompute_psi(&my_plan);
01499         }
01500 
01501         // compute the uniform Fourier transform
01502         nfft_2D_transpose(&my_plan);
01503 
01504         // copy the Fourier transform to EMData data array
01505         EMData *fft = new EMData();
01506         fft->set_size(image_size + 2, image_size, 1);
01507         fft->set_complex(true);
01508         fft->set_ri(true);
01509         fft->to_zero();
01510         float *data = fft->get_data();
01511         double norm = 1.0 / (image_size * image_size);
01512         for (int j = 0; j < image_size; j++) {
01513                 for (int i = 0; i < image_size / 2; i++) {
01514                         data[j * (image_size + 2) + 2 * i] =
01515                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].re) * norm;
01516                         data[j * (image_size + 2) + 2 * i + 1] =
01517                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2].im) * norm * (-1.0);
01518                 }
01519         }
01521         nfft_2D_finalize(&my_plan);
01522 
01523         if (res > 0) {
01524                 // Gaussian low pass processor
01525                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01526                 int index = 0;
01527                 for (int j = 0; j < image_size; j++) {
01528                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01529                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01530                                 double val = exp(-(i * i + RY2) / sigma2);
01531                                 data[index] *= val;
01532                                 data[index + 1] *= val;
01533                         }
01534                 }
01535         }
01536         fft->update();
01537         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
01538 
01539         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01540 
01541         return fft;
01542 #elif defined NFFT2
01543         nfft_plan my_plan;                      // plan for the nfft
01544         int N[2], n[2];
01545         N[0] = image_size;
01546         n[0] = next_power_of_2(2 * image_size);
01547         N[1] = image_size;
01548         n[1] = next_power_of_2(2 * image_size);
01549 
01551         //nfft_init_2d(&my_plan,image_size,image_size,get_number_points());
01552         nfft_init_specific(&my_plan, 2, N, get_number_points(), n, 3,
01553                                            PRE_PHI_HUT | PRE_PSI |
01554                                            MALLOC_X | MALLOC_F_HAT | MALLOC_F, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
01556         for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01557                 // FFTW and nfft use row major array layout, EMAN uses column major
01558                 my_plan.x[2 * j + 1] = (double) (points[i] / (apix * image_size));
01559                 my_plan.x[2 * j] = (double) (points[i + 1] / (apix * image_size));
01560                 my_plan.f[j][0] = (double) points[i + 3];
01561                 my_plan.f[j][1] = 0.0;
01562         }
01563 
01565         if (my_plan.nfft_flags & PRE_PSI) {
01566                 nfft_precompute_psi(&my_plan);
01567                 if (my_plan.nfft_flags & PRE_FULL_PSI)
01568                         nfft_full_psi(&my_plan, pow(10, -6));
01569         }
01570 
01571         // compute the uniform Fourier transform
01572         nfft_transposed(&my_plan);
01573 
01574         // copy the Fourier transform to EMData data array
01575         EMData *fft = new EMData();
01576         fft->set_size(image_size + 2, image_size, 1);
01577         fft->set_complex(true);
01578         fft->set_ri(true);
01579         fft->to_zero();
01580         float *data = fft->get_data();
01581         double norm = 1.0 / (image_size * image_size);
01582         for (int j = 0; j < image_size; j++) {
01583                 for (int i = 0; i < image_size / 2; i++) {
01584                         data[j * (image_size + 2) + 2 * i] =
01585                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][0]) * norm;
01586                         data[j * (image_size + 2) + 2 * i + 1] =
01587                                 (float) (my_plan.f_hat[j * image_size + i + image_size / 2][1]) * norm;
01588                 }
01589         }
01591         nfft_finalize(&my_plan);
01592 
01593         if (res > 0) {
01594                 // Gaussian low pass process
01595                 double sigma2 = (image_size * apix / res) * (image_size * apix / res);
01596                 int index = 0;
01597                 for (int j = 0; j < image_size; j++) {
01598                         double RY2 = (j - image_size / 2) * (j - image_size / 2);
01599                         for (int i = 0; i < image_size / 2 + 1; i++, index += 2) {
01600                                 double val = exp(-(i * i + RY2) / sigma2);
01601                                 data[index] *= val;
01602                                 data[index + 1] *= val;
01603                         }
01604                 }
01605         }
01606         fft->update();
01607         //fft->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs", box*apix/res));
01608 
01609         fft->process_inplace("xform.phaseorigin.tocenter");     // move phase origin to center of image box, instead of at corner
01610 
01611         return fft;
01612 #else
01613         LOGWARN("nfft support is not enabled. please recompile with nfft support enabled\n");
01614         return 0;
01615 #endif
01616 }
01617 
01618 #ifdef OPTPP
01619 #include "NLF.h"
01620 #include "BoundConstraint.h"
01621 #include "OptCG.h"
01622 //#include "OptNewton.h"
01623 #include "newmatap.h"
01624 
01625 vector<EMData*> optdata;
01626 PointArray *optobj;
01627 float optpixres;
01628 
01629 void init_opt_proj(int ndim, ColumnVector& x)
01630 {
01631 int i;
01632 double *data=optobj->get_points_array();
01633 
01634 for (i=0; i<ndim; i++) x(i+1)=data[i];
01635 }
01636 
01637 void calc_opt_proj(int n, const ColumnVector& x, double& fx, int& result)
01638 {
01639         int i;
01640         PointArray pa;
01641         Transform xform;
01642         int size=optdata[0]->get_xsize();
01643         fx=0;
01644 
01645         for (i=0; i<optdata.size(); i++) {
01646                 xform=(optdata[i]->get_transform());
01647                 pa.set_from((double *)x.nric()+1,n/4,std::string("c1"),&xform);
01648                 EMData *p=pa.projection_by_summation(size,1.0,optpixres);
01649                 p->process_inplace("normalize.unitlen");
01650                 fx-=sqrt(p->cmp("dot",EMObject(optdata[i]),Dict()));
01651         }
01652 
01653         result=NLPFunction;
01654 
01655         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));
01656 }
01657 
01658 void calc_opt_projd(int mode,int n, const ColumnVector& x, double& fx, ColumnVector& gx, int& result)
01659 {
01660 
01661 if (mode & NLPFunction) {
01662 
01663 }
01664 
01665 if (mode & NLPGradient) {
01666 
01667 }
01668 
01669 }
01670 #endif
01671 
01672 void PointArray::opt_from_proj(const vector<EMData*> & proj,float pixres) {
01673 #ifdef OPTPP
01674         optdata=proj;
01675         optobj=this;
01676         optpixres=pixres;
01677 
01678         FDNLF1 nlf(get_number_points()*4,calc_opt_proj,init_opt_proj);
01679 //      NLF1 nlf(get_number_points()*4,init_opt_proj,calc_opt_projd);
01680         nlf.initFcn();
01681 
01682         OptCG opt(&nlf);
01683 //      OptQNewton opt(&nlf);
01684         opt.setMaxFeval(2000);
01685         opt.optimize();
01686         opt.printStatus("Done");
01687 #else
01688         (void)proj;             //suppress warning message
01689         (void)pixres;   //suppress warning message
01690         LOGWARN("OPT++ support not enabled.\n");
01691         return;
01692 #endif
01693 }
01694 
01695 Vec3f PointArray::get_vector_at(int i)
01696 {
01697 return Vec3f((float)points[i*4],(float)points[i*4+1],(float)points[i*4+2]);
01698 }
01699 
01700 double PointArray::get_value_at(int i)
01701 {
01702 return points[i*4+3];
01703 }
01704 
01705 void PointArray::set_vector_at(int i,Vec3f vec,double value)
01706 {
01707 points[i*4]=vec[0];
01708 points[i*4+1]=vec[1];
01709 points[i*4+2]=vec[2];
01710 points[i*4+3]=value;
01711 }
01712 
01713 void PointArray::set_vector_at(int i,vector<double> v)
01714 {
01715 points[i*4]  =v[0];
01716 points[i*4+1]=v[1];
01717 points[i*4+2]=v[2];
01718 if (v.size()>=4) points[i*4+3]=v[3];
01719 }

Generated on Tue Jun 11 13:40:40 2013 for EMAN2 by  doxygen 1.3.9.1