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 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();
00194
00195 if (max_miss<0) max_miss=(float)mx0->get_attr("sigma")/10.0f;
00196
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
00211
00212 for (i=0; i<n; i++) {
00213 int bestn=-1;
00214 double bestd=1.0e38;
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
00230 if (d<bestd) { bestd=d; bestn=j; }
00231 }
00232 ret[i]=bestn;
00233 rete[i]=static_cast<float>(bestd);
00234 }
00235
00236
00237
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
00253
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
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
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
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
00285
00286
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
00306
00307
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
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
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
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
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
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
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
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
00400 if (nr<2) return vector<float>(4,0);
00401 dx/=nr;
00402 dy/=nr;
00403
00404
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
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;
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));
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;
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;
00632 alt2 = 0.55357435889704525151;
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
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.);
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));
00766
00767 double table_step_size = 1.;
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);
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
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
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);
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
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;
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
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;
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
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 {
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
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);
00983 }
00984 else {
00985 LOGERR("PointArray::set_from_density_map(): mode = %d is not implemented yet", mode);
00986 }
00987
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
01045
01046 double min_table_val = 1e-7;
01047 double max_table_x = sqrt(-log(min_table_val));
01048
01049 double table_step_size = 0.001;
01050 double inv_table_step_size = 1.0 / table_step_size;
01051
01052
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
01072
01073
01074
01075
01076
01077 if(addpdbbfactor==-1){
01078 gauss_real_width = res/M_PI;
01079 }else{
01080 gauss_real_width = (bfactor[s])/(4*sqrt(2.0)*M_PI);
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);
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);
01156
01157
01158 double min_table_val = 1e-7;
01159 double max_table_x = sqrt(-log(min_table_val));
01160
01161
01162
01163
01164
01165
01166 double table_step_size = 0.001;
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);
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
01199 int table_index_y = int (fabs(j - yc) * inv_table_step_size);
01200 double yval = table[table_index_y];
01201 #ifdef DEBUG
01202
01203
01204 #endif
01205 int pd_index = j * image_size + imin;
01206 for (int i = imin; i < imax; i++, pd_index++) {
01207
01208 int table_index_x = int (fabs(i - xc) * inv_table_step_size);
01209 double xval = table[table_index_x];
01210 #ifdef DEBUG
01211
01212
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);
01227
01228 double min_table_val = 1e-7;
01229 double max_table_x = sqrt(-log(min_table_val));
01230
01231 double table_step_size = 0.001;
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
01242 int gbox = int (max_table_x * gauss_real_width / apix);
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
01270 gbox = int (max_table_x * gauss_real_width / apix);
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;
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
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
01326 nfft_3D_transpose(&my_plan);
01327
01328
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
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
01369
01370 fft->process_inplace("xform.phaseorigin.tocorner");
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;
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
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
01408 nfft_transposed(&my_plan);
01409
01410
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
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
01451
01452 fft->process_inplace("xform.phaseorigin.tocenter");
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;
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
01485
01486
01488 for (int j = 0, i = 0; j < my_plan.M; j++, i += 4) {
01489
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
01502 nfft_2D_transpose(&my_plan);
01503
01504
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
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
01538
01539 fft->process_inplace("xform.phaseorigin.tocenter");
01540
01541 return fft;
01542 #elif defined NFFT2
01543 nfft_plan my_plan;
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
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
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
01572 nfft_transposed(&my_plan);
01573
01574
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
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
01608
01609 fft->process_inplace("xform.phaseorigin.tocenter");
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
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
01680 nlf.initFcn();
01681
01682 OptCG opt(&nlf);
01683
01684 opt.setMaxFeval(2000);
01685 opt.optimize();
01686 opt.printStatus("Done");
01687 #else
01688 (void)proj;
01689 (void)pixres;
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 }