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