00001
00002
00003
00004
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 #include "ctf.h"
00033 #include "log.h"
00034 #include "emdata.h"
00035 #include "xydata.h"
00036 #include "Assert.h"
00037 #include "projector.h"
00038
00039 #include <sys/types.h>
00040 #include <sys/times.h>
00041 #include <time.h>
00042 #include <sys/time.h>
00043
00044 #include <iostream>
00045
00046
00047 using namespace EMAN;
00048
00049 using std::cout;
00050 using std::cin;
00051 using std::endl;
00052 using std::string;
00053
00054 #include "spidutil.h"
00055 #include "spidfft.h"
00056
00057 int aprings(int nimg, int nring, float *imgstk, int nsam,
00058 int nrow, int *numr, float *refcstk, int lcirc,
00059 char mode)
00060 {
00061 int j, status;
00062 float *wr;
00063
00064 status = 0;
00065 wr = (float*) calloc(nring, sizeof(float));
00066 if (!wr) {
00067 fprintf(stderr, "aprings: failed to allocate wr!\n");
00068 status = -1;
00069 goto EXIT;
00070 }
00071
00072
00073
00074 ringwe(wr, numr, nring);
00075 if ( mode == 'H' || mode == 'h' )
00076 for (j=1;j<=nring;j++) wr(j) = wr(j)/2.0;
00077 for (j = 1; j<=nimg; j++) {
00078 apring1(&imgstk(1,1,j), nsam, nrow, &refcstk(1,j),
00079 lcirc, nring, mode, numr, wr);
00080 }
00081
00082 EXIT:
00083 if (wr) free(wr);
00084 return status;
00085 }
00086
00087
00088
00089 int apring1(float *sqimg, int nsam , int nrow, float *imgcirc, int lcirc,
00090 int nring , char mode, int *numr, float *wr)
00091 {
00092 int status = 0;
00093 int nsb, nse, nrb, nre, maxrin, ltf, lt, lt2, lt3, ns2, nr2;
00094 int inr, kcirc, jt, i, j;
00095 float fnr2, fns2, yq, fi, dpi, x, y;
00096 double dfi;
00097
00098 dpi = 2.0*atan(1.0);
00099
00100
00101 nsb = -nsam/2;
00102 nse = -nsb-1+(nsam%2);
00103 nrb = -nrow/2;
00104 nre = -nrb-1+(nrow%2);
00105
00106
00107
00108
00109 normass(sqimg,nsam,nsb,nse,nrb,nre,numr(1,1),numr(1,nring));
00110
00111 maxrin = numr(3,nring);
00112 ns2 = nsam / 2 + 1;
00113 nr2 = nrow / 2 + 1;
00114 fnr2 = nr2;
00115 fns2 = ns2;
00116
00117
00118
00119 if (mode == 'f' || mode == 'F') {
00120 ltf = 4;
00121 for (i=1;i<=nring;i++) {
00122
00123 inr = numr(1,i);
00124 yq = inr;
00125 lt = numr(3,i) / ltf;
00126 lt2 = lt + lt;
00127 lt3 = lt2 + lt;
00128 dfi = dpi / lt;
00129 kcirc = numr(2,i);
00130
00131 imgcirc(kcirc) = sqimg(ns2, nr2+inr);
00132 imgcirc(lt+kcirc) = sqimg(ns2+inr, nr2);
00133 imgcirc(lt2+kcirc) = sqimg(ns2, nr2-inr);
00134 imgcirc(lt3+kcirc) = sqimg(ns2-inr, nr2);
00135
00136 for (j=1;j<=lt - 1;j++) {
00137 fi = dfi * j;
00138 x = sin(fi) * yq;
00139 y = cos(fi) * yq;
00140 jt = j + kcirc;
00141
00142 imgcirc(jt) = quadri(fns2+x,fnr2+y,nsam,nrow,sqimg);
00143 imgcirc(jt+lt) = quadri(fns2+y,fnr2-x,nsam,nrow,sqimg);
00144 imgcirc(jt+lt2) = quadri(fns2-x,fnr2-y,nsam,nrow,sqimg);
00145 imgcirc(jt+lt3) = quadri(fns2-y,fnr2+x,nsam,nrow,sqimg);
00146 }
00147 }
00148 }
00149 else if (mode == 'h' || mode == 'H') {
00150 ltf = 2;
00151 for (i=1; i<=nring;i++) {
00152
00153 inr = numr(1,i);
00154 yq = inr;
00155 lt = numr(3,i) / ltf;
00156 dfi = dpi / lt;
00157 kcirc = numr(2,i);
00158
00159 imgcirc(kcirc) = sqimg(ns2, nr2+inr);
00160 imgcirc(lt+kcirc) = sqimg(ns2+inr, nr2);
00161
00162 for (j=1;j<=lt - 1;j++) {
00163 fi = dfi * j;
00164 x = sin(fi) * yq;
00165 y = cos(fi) * yq;
00166 jt = j + kcirc;
00167
00168 imgcirc(jt) = quadri(fns2+x,fnr2+y,nsam,nrow,sqimg);
00169 imgcirc(jt+lt) = quadri(fns2+y,fnr2-x,nsam,nrow,sqimg);
00170 }
00171 }
00172 }
00173
00174
00175 frngs(imgcirc,numr,nring);
00176
00177
00178 if (wr(1) > 0.0) {
00179 status = applyws(imgcirc,lcirc,numr,wr,nring);
00180 }
00181
00182 return status;
00183 }
00184
00185
00186 float quadri(float xx, float yy, int nxdata, int nydata, float *fdata)
00187 {
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 float x, y, dx0, dy0, f0, c1, c2, c3, c4, c5, dxb, dyb;
00232 float quadri;
00233 int i, j, ip1, im1, jp1, jm1, ic, jc, hxc, hyc;
00234
00235 x = xx;
00236 y = yy;
00237
00238
00239 if (x < 1.0) x = x+(1 - floor(x) / nxdata) * nxdata;
00240 if (x > (float)nxdata+0.5) x = fmod(x-1.0,(float)nxdata) + 1.0;
00241 if (y < 1.0) y = y+(1 - floor(y) / nydata) * nydata;
00242 if (y > (float)nydata+0.5) y = fmod(y-1.0,(float)nydata) + 1.0;
00243
00244
00245 i = (int) floor(x);
00246 j = (int) floor(y);
00247
00248 dx0 = x - i;
00249 dy0 = y - j;
00250
00251 ip1 = i + 1;
00252 im1 = i - 1;
00253 jp1 = j + 1;
00254 jm1 = j - 1;
00255
00256 if (ip1 > nxdata) ip1 = ip1 - nxdata;
00257 if (im1 < 1) im1 = im1 + nxdata;
00258 if (jp1 > nydata) jp1 = jp1 - nydata;
00259 if (jm1 < 1) jm1 = jm1 + nydata;
00260
00261 f0 = fdata(i,j);
00262 c1 = fdata(ip1,j) - f0;
00263 c2 = (c1 - f0 + fdata(im1,j)) * 0.5;
00264 c3 = fdata(i,jp1) - f0;
00265 c4 = (c3 - f0 + fdata(i,jm1)) * 0.5;
00266
00267 dxb = dx0 - 1;
00268 dyb = dy0 - 1;
00269
00270
00271 if (dx0 >= 0) {
00272 hxc = 1;
00273 }
00274 else {
00275 hxc = -1;
00276 }
00277 if (dy0 >= 0) {
00278 hyc = 1;
00279 }
00280 else {
00281 hyc = -1;
00282 }
00283
00284 ic = i + hxc;
00285 jc = j + hyc;
00286
00287 if (ic > nxdata) {
00288 ic = ic - nxdata;
00289 }
00290 else if (ic < 1) {
00291 ic = ic + nxdata;
00292 }
00293
00294 if (jc > nydata) {
00295 jc = jc - nydata;
00296 }
00297 else if (jc < 1) {
00298 jc = jc + nydata;
00299 }
00300
00301 c5 = ( (fdata(ic,jc) - f0 - hxc * c1 - (hxc * (hxc - 1.0)) * c2
00302 - hyc * c3 - (hyc * (hyc - 1.0)) * c4) * (hxc * hyc));
00303
00304 quadri = f0 + dx0 * (c1 + dxb * c2 + dy0 * c5) + dy0 * (c3 + dyb * c4);
00305
00306 return quadri;
00307 }
00308
00309 int alprbs(int *numr, int nring, int *lcirc, char mode)
00310 {
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 int i, jp, ip;
00326 double dpi;
00327 int status = 0;
00328
00329 int MAXFFT = 32768;
00330
00331 dpi = pi;
00332 if (mode == 'f' || mode == 'F') dpi=2*pi;
00333
00334 *lcirc = 0;
00335 for (i=1;i<=nring;i++) {
00336 jp = (int)(dpi*numr(1,i));
00337
00338 ip = (int)( pow(2,ceil(log2(jp))) );
00339 if (i < nring && jp > ip+ip/2) ip=min0(MAXFFT,2*ip);
00340
00341
00342
00343 if (i == nring && jp > ip+ip/5) ip=min0(MAXFFT,2*ip);
00344 numr(3,i) = ip;
00345 if (i == 1) {
00346 numr(2,1) = 1;
00347 }
00348 else {
00349 numr(2,i) = numr(2,i-1)+numr(3,i-1);
00350 }
00351 *lcirc = *lcirc + ip;
00352 }
00353 return status;
00354 }
00355
00356
00357 void ringwe(float *wr, int *numr, int nring)
00358 {
00359 int i, maxrin;
00360 float dpi;
00361
00362 dpi = 8.0*atan(1.0);
00363 maxrin = numr(3,nring);
00364 for (i=1;i<=nring;i++) {
00365 wr(i)=(float)(numr(1,i))*dpi/(float)(numr(3,i))
00366 *(float)(maxrin)/(float)(numr(3,i));
00367 cout << numr(1,i)<<" "<< numr(2,i)<<" "<< numr(3,i)<<" " << wr(i)<<" "<<endl;
00368 }
00369 }
00370
00371
00372 int setnring(int mr, int nr, int iskip)
00373 {
00374 int nring = 0, i;
00375
00376 i = mr;
00377 while (i<=nr) {
00378 nring = nring+1;
00379 i = i + iskip;
00380 }
00381 return nring;
00382 }
00383
00384
00385 void numrinit(int mr, int nr, int iskip, int *numr)
00386 {
00387 int nring = 0;
00388 int i;
00389
00390 i = mr;
00391 while (i<=nr) {
00392 nring++;
00393 numr(1,nring) = i;
00394 i=i+iskip;
00395 }
00396 }
00397
00398
00399 void normass(float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2,
00400 int ir1, int ir2)
00401 {
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 double av,vr,vrinv,dtemp;
00413 int i1sq, i2sq, n, ir, jsq, i, j, irow, jcol;
00414
00415 i1sq = ir1 * ir1;
00416 i2sq = ir2 * ir2;
00417 av = 0.0;
00418 vr = 0.0;
00419 n = 0;
00420
00421 for (j=nr1; j<=nr2; j++) {
00422 jsq = j * j;
00423 for (i=ns1;i<=ns2;i++) {
00424 ir = jsq + i * i;
00425 if (ir >= i1sq && ir <= i2sq) {
00426 n = n + 1;
00427 irow = i-ns1+1;
00428 jcol = j-nr1+1;
00429 av = av + sqimg(irow,jcol);
00430 vr = vr + sqimg(irow,jcol)*sqimg(irow,jcol);
00431 }
00432 }
00433 }
00434
00435 av = av / n;
00436 dtemp = (vr - n * av * av);
00437 if (dtemp > 0) {
00438 vr = sqrt(dtemp / (n-1));
00439 vrinv = 1.0 / vr;
00440
00441
00442 for ( j = nr1; j<=nr2;j++)
00443 for (i = ns1; i <= ns2; i++) {
00444 irow = i - ns1 + 1;
00445 jcol = j - nr1 + 1;
00446 sqimg(irow,jcol) = (sqimg(irow,jcol) - av) * vrinv;
00447 }
00448 }
00449 else {
00450
00451
00452 for ( j = nr1; j<=nr2;j++)
00453 for (i = ns1; i <= ns2; i++) {
00454 irow = i - ns1 + 1;
00455 jcol = j - nr1 + 1;
00456 sqimg(irow,jcol) = 0.0;
00457 }
00458 }
00459 }
00460
00461
00462
00463 void frngs(float *circ, int *numr, int nring)
00464 {
00465 int i, l;
00466
00467 for (i=1; i<=nring;i++) {
00468 l=(int)(log2(numr(3,i)));
00469 fftr_q(&circ(numr(2,i)),l);
00470 }
00471 }
00472
00473
00474 int applyws(float *circ, int lcirc, int *numr, float *wr,
00475 int nring)
00476 {
00477 int maxrin, numr3i, numr2i;
00478 float w;
00479 int i, j, jc;
00480 int status = 0;
00481
00482 maxrin = numr(3,nring);
00483
00484 for (i=1;i<=nring;i++) {
00485 numr3i=numr(3,i);
00486 numr2i=numr(2,i);
00487 w=wr(i);
00488 circ(numr2i)=circ(numr2i)*w;
00489 if (numr3i == maxrin) {
00490 circ(numr2i+1)=circ(numr2i+1)*w;
00491 }
00492 else {
00493 circ(numr2i+1)=circ(numr2i+1)*0.5*w;
00494 }
00495 for (j=3;j<=numr3i;j++) {
00496 jc=j+numr2i-1;
00497 circ(jc)=circ(jc)*w;
00498 }
00499 }
00500 if (jc >= lcirc) status = -1;
00501 return status;
00502 }
00503
00504
00505 void normas(float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2,
00506 int ir1, int ir2)
00507 {
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 double av,vr;
00520 int i1sq, i2sq, n, i, j, ir, j2, irow, jcol;
00521
00522 i1sq = ir1 * ir1;
00523 i2sq = ir2 * ir2;
00524
00525 av = 0.0;
00526 vr = 0.0;
00527 n = 0;
00528
00529 for (j=nr1;j<=nr2;j++) {
00530 j2 = j*j;
00531 for (i=ns1;i<=ns2;i++) {
00532 ir = j2 + i*i;
00533 jcol = j-nr1+1;
00534 if (ir >= i1sq && ir <= i2sq) {
00535 n++;
00536 irow = i-ns1+1;
00537 av = av + sqimg(irow,jcol);
00538 vr = vr + sqimg(irow,jcol)*sqimg(irow,jcol);
00539 }
00540 }
00541 }
00542
00543 av = av / n;
00544
00545
00546 vr = 1.0 / (sqrt((vr-n*av*av) / (n-1)));
00547
00548 for (j=nr1; j<=nr2; j++) {
00549 jcol = j-nr1+1;
00550 for (i=ns1;i<=ns2;i++) {
00551 irow = i-ns1+1;
00552 sqimg(irow,jcol) = (sqimg(irow,jcol) - av ) * vr;
00553 }
00554 }
00555 }
00556
00557
00558 void alrq(float *xim, int nsam , int nrow , int *numr,
00559 float *circ, int lcirc, int nring, char mode)
00560 {
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 double dfi, dpi;
00572 int ns2, nr2, i, inr, l, nsim, kcirc, lt, j;
00573 float yq, xold, yold, fi, x, y;
00574
00575 ns2 = nsam/2+1;
00576 nr2 = nrow/2+1;
00577 dpi = 2.0*atan(1.0);
00578
00579
00580
00581 for (i=1;i<=nring;i++) {
00582
00583 inr = numr(1,i);
00584 yq = inr;
00585 l = numr(3,i);
00586 if (mode == 'h' || mode == 'H') {
00587 lt = l/2;
00588 }
00589 else if (mode == 'f' || mode == 'F' ) {
00590 lt = l/4;
00591 }
00592
00593 nsim = lt-1;
00594 dfi = dpi/(nsim+1);
00595 kcirc = numr(2,i);
00596 xold = 0.0;
00597 yold = inr;
00598 circ(kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00599 xold = inr;
00600 yold = 0.0;
00601 circ(lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00602
00603 if (mode == 'f' || mode == 'F') {
00604 xold = 0.0;
00605 yold = -inr;
00606 circ(lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00607 xold = -inr;
00608 yold = 0.0;
00609 circ(lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00610 }
00611
00612 for (j=1;j<=nsim;j++) {
00613 fi = dfi*j;
00614 x = sin(fi)*yq;
00615 y = cos(fi)*yq;
00616 xold = x;
00617 yold = y;
00618 circ(j+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00619 xold = y;
00620 yold = -x;
00621 circ(j+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00622
00623 if (mode == 'f' || mode == 'F') {
00624 xold = -x;
00625 yold = -y;
00626 circ(j+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00627 xold = -y;
00628 yold = x;
00629 circ(j+lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
00630 };
00631 }
00632 }
00633
00634 }
00635
00636
00637 int crosrng_ms(float *circ1, float *circ2, int lcirc, int nring,
00638 int maxrin, int *numr , double *qn, float *tot,
00639 double *qm, double *tmt)
00640 {
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658 double *t, *q, t7[7];
00659
00660 int ip, jc, numr3i, numr2i, i, j, k, jtot;
00661 float t1, t2, t3, t4, c1, c2, d1, d2, pos;
00662
00663 int status = 0;
00664
00665 *qn = 0.0;
00666 *qm = 0.0;
00667 *tot = 0.0;
00668 *tmt = 0.0;
00669
00670 ip = -(int)(log2(maxrin));
00671
00672
00673
00674
00675 q = (double*)calloc(maxrin+2,sizeof(double));
00676 if (!q) {
00677 status = -1;
00678 return status;
00679 }
00680
00681
00682
00683 t = (double*)calloc(maxrin+2,sizeof(double));
00684 if (!t) {
00685 status = -1;
00686 return status;
00687 }
00688
00689
00690
00691 for (i=1;i<=nring;i++) {
00692
00693 numr3i = numr(3,i);
00694 numr2i = numr(2,i);
00695
00696 t1 = circ1(numr2i) * circ2(numr2i);
00697 q(1) = q(1)+t1;
00698 t(1) = t(1)+t1;
00699
00700 if (numr3i == maxrin) {
00701 t1 = circ1(numr2i+1) * circ2(numr2i+1);
00702 q(2) = q(2)+t1;
00703 t(2) = t(2)+t1;
00704 }
00705 else {
00706 t1 = circ1(numr2i+1) * circ2(numr2i+1);
00707 q(numr3i+1) = q(numr3i+1)+t1;
00708 }
00709
00710 for (j=3;j<=numr3i;j=j+2) {
00711 jc = j+numr2i-1;
00712
00713 c1 = circ1(jc);
00714 c2 = circ1(jc+1);
00715 d1 = circ2(jc);
00716 d2 = circ2(jc+1);
00717
00718 t1 = c1 * d1;
00719 t3 = c1 * d2;
00720 t2 = c2 * d2;
00721 t4 = c2 * d1;
00722
00723 q(j) = q(j) + t1 + t2;
00724 q(j+1) = q(j+1) - t3 + t4;
00725 t(j) = t(j) + t1 - t2;
00726 t(j+1) = t(j+1) - t3 - t4;
00727 }
00728 }
00729
00730 fftr_d(q,ip);
00731
00732 jtot = 0;
00733 *qn = -1.0e20;
00734 for (j=1; j<=maxrin; j++) {
00735 if (q(j) >= *qn) {
00736 *qn = q(j);
00737 jtot = j;
00738 }
00739 }
00740
00741 if (jtot <= 0) {
00742
00743 fprintf(stderr, "no max in crosrng_ms, compiler error\n");
00744 status = -1;
00745 return status;
00746 }
00747
00748 for (k=-3;k<=3;k++) {
00749 j = ((jtot+k+maxrin-1)%maxrin)+1;
00750 t7(k+4) = q(j);
00751 }
00752
00753
00754 prb1d(t7,7,&pos);
00755 *tot = (float)(jtot)+pos;
00756
00757
00758 fftr_d(t,ip);
00759
00760
00761 *qm = -1.0e20;
00762 for (j=1; j<=maxrin;j++) {
00763 if ( t(j) >= *qm ) {
00764 *qm = t(j);
00765 jtot = j;
00766 }
00767 }
00768
00769
00770 for (k=-3;k<=3;k++) {
00771 j = ((jtot+k+maxrin-1)%maxrin) + 1;
00772 t7(k+4) = t(j);
00773 }
00774
00775
00776
00777 prb1d(t7,7,&pos);
00778 *tmt = float(jtot) + pos;
00779
00780 free(t);
00781 free(q);
00782
00783 return status;
00784 }
00785
00786
00787 void prb1d(double *b, int npoint, float *pos)
00788 {
00789 double c2,c3;
00790 int nhalf;
00791
00792 nhalf = npoint/2 + 1;
00793 *pos = 0.0;
00794
00795 if (npoint == 7) {
00796 c2 = 49.*b(1) + 6.*b(2) - 21.*b(3) - 32.*b(4) - 27.*b(5)
00797 - 6.*b(6) + 31.*b(7);
00798 c3 = 5.*b(1) - 3.*b(3) - 4.*b(4) - 3.*b(5) + 5.*b(7);
00799 }
00800 else if (npoint == 5) {
00801 c2 = (74.*b(1) - 23.*b(2) - 60.*b(3) - 37.*b(4)
00802 + 46.*b(5) ) / (-70.);
00803 c3 = (2.*b(1) - b(2) - 2.*b(3) - b(4) + 2.*b(5) ) / 14.0;
00804 }
00805 else if (npoint == 3) {
00806 c2 = (5.*b(1) - 8.*b(2) + 3.*b(3) ) / (-2.0);
00807 c3 = (b(1) - 2.*b(2) + b(3) ) / 2.0;
00808 }
00809 else if (npoint == 9) {
00810 c2 = (1708.*b(1) + 581.*b(2) - 246.*b(3) - 773.*b(4)
00811 - 1000.*b(5) - 927.*b(6) - 554.*b(7) + 119.*b(8)
00812 + 1092.*b(9) ) / (-4620.);
00813 c3 = (28.*b(1) + 7.*b(2) - 8.*b(3) - 17.*b(4) - 20.*b(5)
00814 - 17.*b(6) - 8.*b(7) + 7.*b(8) + 28.*b(9) ) / 924.0;
00815 }
00816 if (c3 != 0.0) *pos = c2/(2.0*c3) - nhalf;
00817 }
00818
00819
00820 void apmaster_1(char mode, float *divas, int nr, int *numth,
00821 int lsam, int lrow, int *nsam, int *nrow)
00822 {
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833 int nra;
00834
00835 if ( mode == 'h') {
00836 *divas = 180.0;
00837 }
00838 else {
00839 *divas = 360.0;
00840 }
00841
00842 *numth = 1;
00843 #ifdef sp_mp
00844
00845
00846 #endif
00847
00848
00849
00850
00851 nra = ((lsam-1)/2)*2+1;
00852 if ( ((lrow-1)/2)*2+1 < nra ) nra = ((lrow-1)/2)*2+1;
00853 if ( 2*nr+3 < nra ) nra = 2*nr+3;
00854
00855
00856 *nsam = nra;
00857 *nrow = nra;
00858 }
00859
00860 void win_resize(float *imgfrom, float *imgto, int lsam, int lrow,
00861 int nsam, int nrow, int lr1, int lr2, int ls1, int ls2)
00862 {
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875 int window;
00876
00877 int k3, k2, kt;
00878
00879
00880
00881 window = 0;
00882 if (lr1 != 1 || ls1 != 1 || lr2 != lrow || ls2 != lsam) window = 1;
00883
00884 if (window) {
00885
00886 for (k2=lr1;k2<=lr2;k2++) {
00887 kt = k2-lr1+1;
00888 for (k3=ls1;k3<=ls2;k3++)
00889 imgto(k3-ls1+1,kt) = imgfrom(k3,k2);
00890 }
00891 }
00892 else {
00893
00894 for (k2=1;k2<=lrow;k2++)
00895 for (k3=1;k3<=lsam;k3++)
00896 imgto(k3,k2) = imgfrom(k3,k2);
00897 }
00898 }
00899
00900
00901 int apmd(EMData *refprj, Dict refparams, EMData *expimg, APMDopt options,
00902 float *fangles)
00903 {
00904
00905
00906
00907
00908
00909
00910
00911
00912 int mr, nr, iskip;
00913 char mode;
00914 int nring, lcirc, nref, nimg, maxrin, status;
00915 int nx, ny, nxr, nyr, nsam, nrow;
00916 int nwsam, nwrow, ns1, ns2, nr1, nr2;
00917 int lq, lr1, lr2, ls1, ls2;
00918 float *refstk, *refcstk, *imgstk, *imgcirc,*imgwindow;
00919 int *numr;
00920 double *totmin, *totmir, *tmt;
00921 float *tot;
00922 float divas;
00923 int j, k, numth, idi, ldd, nang, mm;
00924 float *dlist;
00925 double eav;
00926 float rang;
00927 vector <float> angrefs;
00928
00929 status = 0;
00930
00931
00932 mr = options.mr;
00933 nr = options.nr;
00934 iskip = options.iskip;
00935 mode = options.mode;
00936
00937 nx = expimg->get_xsize();
00938 ny = expimg->get_ysize();
00939 nimg = expimg->get_zsize();
00940
00941 nxr = refprj->get_xsize();
00942 nyr = refprj->get_ysize();
00943 nref = refprj->get_zsize();
00944
00945 if (nx != nxr || ny != nyr) {
00946 status = -2;
00947 goto EXIT;
00948 }
00949
00950 nsam = nx;
00951 nrow = ny;
00952
00953
00954 refstk = refprj->get_data();
00955
00956 nring = setnring(mr, nr, iskip);
00957
00958 numr = (int*) calloc(3*nring,sizeof(int));
00959
00960 numrinit(mr, nr, iskip, numr);
00961
00962 status = alprbs(numr, nring, &lcirc, mode);
00963 if (status != 0) goto EXIT;
00964
00965
00966 refcstk = (float*)calloc(lcirc*nref,sizeof(float));
00967
00968
00969 status = aprings(nref, nring, refstk, nsam, nrow, numr,
00970 refcstk, lcirc, mode);
00971 if (status != 0) goto EXIT;
00972
00973
00974 imgstk = expimg->get_data();
00975
00976
00977
00978 apmaster_1(mode, &divas, nr, &numth, nsam, nrow, &nwsam, &nwrow);
00979
00980
00981 imgwindow = (float*)calloc(nwsam*nwrow,sizeof(float));
00982 imgcirc = (float*)calloc(lcirc, sizeof(float));
00983
00984 totmin = (double*)calloc(nref,sizeof(double));
00985 totmir = (double*)calloc(nref,sizeof(double));
00986 tot = (float*)calloc(nref,sizeof(float));
00987 tmt = (double*)calloc(nref,sizeof(double));
00988 if (totmin == NULL || totmir == NULL || tot == NULL || tmt == NULL) {
00989 fprintf(stderr,"apmd: failed to allocate totmin, totmir, tot, tmt\n");
00990 status = -1;
00991 goto EXIT;
00992 }
00993
00994 maxrin = numr(3,nring);
00995 ldd = 3;
00996 dlist = (float*)calloc(ldd*nimg, sizeof(float));
00997 if ( !dlist ) {
00998 status = -1;
00999 goto EXIT;
01000 }
01001
01002
01003 lq =nsam/2+1;
01004 lr1=(nwrow-1)/2;
01005 lr2=lq+lr1;
01006 lr1=lq-lr1;
01007 lq =nrow/2+1;
01008 ls1=(nwsam-1)/2;
01009 ls2=lq+ls1;
01010 ls1=lq-ls1;
01011
01012 for (j = 1; j<=nimg; j++) {
01013 win_resize(&imgstk(1,1,j), imgwindow, nsam, nrow, nwsam, nwrow,
01014 lr1, lr2, ls1, ls2);
01015
01016
01017 ns1 = -nwsam/2;
01018 ns2 = nwsam/2;
01019 nr1 = -nwrow/2;
01020 nr2 = nwrow/2;
01021 normas(imgwindow, nwsam, ns1, ns2, nr1, nr2, numr(1,1), numr(1,nring));
01022
01023
01024 alrq(imgwindow,nwsam,nwrow,numr,imgcirc,lcirc,nring,mode);
01025
01026
01027 frngs(imgcirc, numr, nring);
01028
01029 for (k = 1; k<=nref; k++) {
01030 status = crosrng_ms(&refcstk(1,k), imgcirc, lcirc , nring,
01031 maxrin , numr , &totmin(k), &tot(k),
01032 &totmir(k) , &tmt(k));
01033 }
01034
01035 eav = 1.0e-20;
01036 for (k=1;k<=nref;k++) {
01037 if ( totmin(k) >= eav && totmin(k) > totmir(k) ) {
01038 eav = totmin(k);
01039 idi = k;
01040 rang = tot(k);
01041 }
01042 else if ( totmir(k) >= eav ) {
01043 eav = totmir(k);
01044 idi = -k;
01045 rang = tmt(k);
01046 }
01047 }
01048
01049 dlist(1,j) = idi;
01050 dlist(2,j) = eav;
01051 rang = (rang-1)/maxrin*divas;
01052 dlist(3,j) = rang;
01053
01054 }
01055
01056
01057 angrefs = refparams["anglelist"];
01058 nang = angrefs.size()/3;
01059 if (nang != nref) {
01060 fprintf(stderr, "apmd: nang = %d, nref = %d\n", nang, nref);
01061 status = -3;
01062 goto EXIT;
01063 }
01064
01065 for (j = 1; j<=nimg; j++) {
01066 mm = (int) dlist(1,j);
01067 if( mm != 0) {
01068 fangles(1,j)=-dlist(3,j)+360.0;
01069 if ( mm < 0) {
01070 mm = -mm;
01071 fangles(1,j)=fangles(1,j)+180.0;
01072 fangles(2,j)=180.0-angrefs(2,mm);
01073 fangles(3,j)=angrefs(1,mm)+180.0;
01074 if ( fangles(1,j) >= 360.0) fangles(1,j)=fangles(1,j)-360.0;
01075 if ( fangles(3,j) >= 360.0) fangles(3,j)=fangles(3,j)-360.0;
01076 }
01077 else {
01078 fangles(2,j)=angrefs(2,mm);
01079 fangles(3,j)=angrefs(1,mm);
01080 }
01081 }
01082 }
01083
01084 EXIT:
01085 if (numr) free(numr);
01086 if (refcstk) free(refcstk);
01087 if (imgwindow) free(imgwindow);
01088 if (imgcirc) free(imgcirc);
01089 if (totmin) free(totmin);
01090 if (totmir) free(totmir);
01091 if (tot) free(tot);
01092 if (tmt) free(tmt);
01093 if (dlist) free(dlist);
01094
01095 return status;
01096 }
01097
01098
01099 int aprq2d(float *sqimg, float *bfc, int *numr, int nsam,
01100 int nrow, int ishrange, int istep, int nsb,
01101 int nse, int nrb, int nre, int lcirc,
01102 int nring, int maxrin, int nima, char mode,
01103 float *refdirs, float *expdir, float range, float *diref,
01104 float *ccrot, float *rangnew, float *xshsum, float *yshsum,
01105 int *nimalcg, int ckmirror, int limitrange)
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128 {
01129 float *imgcirc;
01130 int *lcg;
01131
01132 double ccrotd,peak,tota,tmta,tmt;
01133 int mirrored;
01134
01135 double quadpi=3.14159265358979323846;
01136 double dgr_to_rad = quadpi/180.0;
01137
01138 int imi, iend, mwant, jtma, itma;
01139 float dt, dtabs, rangnewt;
01140 int jt, it, irr, ir, ibe, isx, isy, status, idis;
01141 float cnr2, cns2, co, so, afit, sx, sy, tot;
01142
01143 double fitp[9], *fit;
01144 float *rotmp;
01145 int fitsize;
01146
01147 status = 0;
01148
01149 imgcirc = (float*)calloc(lcirc,sizeof(float));
01150 if (imgcirc == NULL) {
01151 fprintf(stderr,"aprq2d: failed to allocate imgcirc\n");
01152 status = -1;
01153 goto EXIT;
01154 }
01155
01156 fitsize = (2*istep+1)*(2*istep+1);
01157 if ( istep >= 1) {
01158 fit = (double*)calloc(fitsize,sizeof(double));
01159 rotmp = (float*) calloc(fitsize,sizeof(float));
01160 }
01161 else {
01162 status = -2;
01163 goto EXIT;
01164 }
01165
01166 peak = 0.0;
01167 iend = nima;
01168
01169 if (limitrange) {
01170
01171 lcg = (int*) calloc(nima, sizeof(int));
01172 if (!lcg) {
01173 mwant = nima;
01174 status = -1;
01175 fprintf(stderr,"lcg: %d\n", mwant);
01176 fprintf(stderr, "range: %g, nima: %d\n", range, nima);
01177 goto EXIT;
01178 }
01179
01180 *nimalcg = 0;
01181 for (imi=1; imi<=nima; imi++) {
01182
01183 dt = (expdir(1) * refdirs(1,imi) +
01184 expdir(2) * refdirs(2,imi) +
01185 expdir(3) * refdirs(3,imi));
01186 dtabs = fabs(dt);
01187
01188 if (dtabs >= range) {
01189
01190 *nimalcg++;
01191 lcg(*nimalcg) = imi;
01192 if (dt < 0) lcg(*nimalcg) = -imi;
01193 }
01194 }
01195
01196 if (*nimalcg <= 0) {
01197
01198 *xshsum = 0;
01199 *yshsum = 0;
01200 *diref = 0;
01201 *rangnew = 0;
01202 *ccrot = -1.0;
01203 goto EXIT;
01204 }
01205 iend = *nimalcg;
01206
01207 }
01208
01209 ccrotd = -1.0e23;
01210
01211
01212 for (jt=-ishrange;jt<=ishrange;jt=jt+istep) {
01213 cnr2 = nrow/2+1+jt;
01214 for (it=-ishrange;it<=ishrange;it=it+istep) {
01215 cns2 = nsam/2+1+it;
01216
01217
01218
01219 normass(sqimg,nsam,nsb-it,nse-it,nrb-jt,nre-jt,numr(1,1),
01220 numr(1,nring));
01221
01222
01223
01224
01225 alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,nring,mode);
01226
01227
01228 frngs(imgcirc,numr,nring);
01229
01230
01231 for (irr=1;irr<=iend;irr++) {
01232 ir = irr;
01233 if (limitrange) ir = abs(lcg(irr));
01234
01235 if (ckmirror) {
01236 if (limitrange) {
01237 mirrored = 0;
01238 if (lcg(irr) < 0) mirrored = 1;
01239
01240 crosrng_e(&bfc(1,ir),imgcirc,lcirc,nring,
01241 maxrin,numr,&tota,&tot,mirrored);
01242 }
01243 else {
01244
01245 status = crosrng_ms(&bfc(1,ir),imgcirc,lcirc,nring,
01246 maxrin,numr,&tota,&tot,&tmta,&tmt);
01247 }
01248 }
01249 else {
01250
01251 mirrored = 0;
01252 crosrng_e(&bfc(1,ir),imgcirc,lcirc,nring,
01253 maxrin,numr,&tota,&tot,mirrored);
01254 }
01255
01256 if (tota >= ccrotd) {
01257
01258 ccrotd = tota;
01259 ibe = ir;
01260 isx = it;
01261 isy = jt;
01262 *rangnew = ang_n(tot,mode,maxrin);
01263 idis = ir;
01264 if (limitrange && lcg(irr) < 0) idis = -ir;
01265 }
01266
01267 if (ckmirror && !limitrange) {
01268
01269 if (tmta >= ccrotd) {
01270
01271 ccrotd = tmta;
01272 ibe = ir;
01273 isx = it;
01274 isy = jt;
01275 *rangnew = ang_n(tmt,mode,maxrin);
01276 idis = -ir;
01277 }
01278 }
01279 }
01280 }
01281 }
01282
01283
01284 *ccrot = ccrotd;
01285 sx = isx;
01286 sy = isy;
01287 *diref = idis;
01288
01289
01290 if ( abs(isx) != ishrange && abs(isy) != ishrange) {
01291
01292 fit(0,0) = ccrotd;
01293 rotmp(0,0) = *rangnew;
01294
01295 for (jt=-istep;jt<=istep;jt++) {
01296 for (it=-istep;it<=istep;it++) {
01297 if (it !=0 || jt != 0) {
01298 cnr2 = nrow/2+1+jt+isy;
01299 cns2 = nsam/2+1+it+isx;
01300
01301 normass(sqimg, nsam, nsb-(it+isx),nse-(it+isx),
01302 nrb-(jt+isy),nre-(jt+isy), numr(1,1), numr(1,nring));
01303
01304 alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,
01305 nring,mode);
01306
01307 frngs(imgcirc,numr,nring);
01308
01309
01310 mirrored = 0;
01311 if (idis < 0) mirrored = 1;
01312 crosrng_e(&bfc(1,ibe),imgcirc,lcirc,nring,
01313 maxrin,numr,&fit(it,jt),&rotmp(it,jt),
01314 mirrored);
01315 rotmp(it,jt) = ang_n(rotmp(it,jt),mode,maxrin);
01316 }
01317 }
01318 }
01319
01320
01321
01322 afit = fit(0,0);
01323 jtma = 0;
01324 itma = 0;
01325 rangnewt = rotmp(0,0);
01326 if ( istep > 1) {
01327 for (jt=-istep+1;jt<=istep-1;jt++) {
01328 for (it=-istep+1;it<=istep-1;it++) {
01329 if (fit(it,jt) > afit) {
01330 afit = fit(it,jt);
01331 rangnewt = rotmp(it,jt);
01332 itma = it;
01333 jtma = jt;
01334 }
01335 }
01336 }
01337 }
01338
01339 *rangnew = rangnewt;
01340
01341
01342 for (jt=-1;jt<=1;jt++)
01343 for (it=-1;it<=1;it++)
01344 fitp(it,jt) = fit(itma+it,jtma+jt);
01345
01346
01347 ccrotd = afit;
01348 isx = isx+itma;
01349 isy = isy+jtma;
01350 parabld(fitp,&sx,&sy,&peak);
01351
01352
01353 if (fabs(sx) < 1.0 && fabs(sy) < 1.0) {
01354
01355 sx = sx+isx;
01356 sy = sy+isy;
01357 cnr2 = nrow/2+1+sy;
01358 cns2 = nsam/2+1+sx;
01359
01360 normass(sqimg,nsam,nsb-isx,nse-isx,nrb-isy,nre-isy,numr(1,1),
01361 numr(1,nring));
01362
01363 alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,nring,mode);
01364
01365 frngs(imgcirc,numr,nring);
01366
01367 mirrored = 0;
01368 if (idis < 0) mirrored = 1;
01369
01370 crosrng_e(&bfc(1,ibe),imgcirc,lcirc,nring,
01371 maxrin,numr,&ccrotd,rangnew,mirrored);
01372
01373 *ccrot = ccrotd;
01374 *rangnew = ang_n(*rangnew,mode,maxrin);
01375 }
01376 else {
01377
01378 sx = isx;
01379 sy = isy;
01380 }
01381 }
01382
01383 sx = -sx;
01384 sy = -sy;
01385
01386
01387
01388
01389
01390 co = cos((*rangnew) * dgr_to_rad);
01391 so = -sin((*rangnew) * dgr_to_rad);
01392 *xshsum = sx*co - sy*so;
01393 *yshsum = sx*so + sy*co;
01394
01395 free(fit);
01396 free(rotmp);
01397 free(imgcirc);
01398 if (limitrange) free(lcg);
01399
01400 EXIT:
01401 return status;
01402 }
01403
01404
01405
01406 float ang_n(float rkk, char mode, int maxrin)
01407 {
01408 float ang;
01409
01410 if (mode == 'H' || mode == 'h') {
01411 ang = fmod(((rkk-1.0) / maxrin+1.0)*180.0, 180.0);
01412 }
01413 else if ( mode == 'F' || mode == 'f') {
01414 ang = fmod(((rkk-1.0) / maxrin+1.0)*360.0, 360.0);
01415 }
01416 return ang;
01417 }
01418
01419 void alrq_ms(float *xim, int nsam, int nrow, float cns2, float cnr2,
01420 int *numr, float *circ, int lcirc, int nring, char mode)
01421 {
01422 double dpi, dfi;
01423 int it, jt, inr, l, nsim, kcirc, lt;
01424 float yq, xold, yold, fi, x, y;
01425
01426
01427
01428
01429 dpi = 2*atan(1.0);
01430 for (it=1;it<=nring;it++) {
01431
01432 inr = numr(1,it);
01433 yq = inr;
01434
01435 l = numr(3,it);
01436 if ( mode == 'h' || mode == 'H' ) {
01437 lt = l / 2;
01438 }
01439 else if ( mode == 'f' || mode == 'F' ) {
01440 lt = l / 4;
01441 }
01442
01443 nsim = lt - 1;
01444 dfi = dpi / (nsim+1);
01445 kcirc = numr(2,it);
01446 xold = 0.0;
01447 yold = inr;
01448
01449 circ(kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01450
01451 xold = inr;
01452 yold = 0.0;
01453 circ(lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01454
01455 if ( mode == 'f' || mode == 'F' ) {
01456 xold = 0.0;
01457 yold = -inr;
01458 circ(lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01459
01460 xold = -inr;
01461 yold = 0.0;
01462 circ(lt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01463 }
01464
01465 for (jt=1;jt<=nsim;jt++) {
01466 fi = dfi * jt;
01467 x = sin(fi) * yq;
01468 y = cos(fi) * yq;
01469
01470 xold = x;
01471 yold = y;
01472 circ(jt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01473
01474 xold = y;
01475 yold = -x;
01476 circ(jt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01477
01478 if ( mode == 'f' || mode == 'F' ) {
01479 xold = -x;
01480 yold = -y;
01481 circ(jt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01482
01483 xold = -y;
01484 yold = x;
01485 circ(jt+lt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
01486 }
01487 }
01488 }
01489 }
01490
01491 void parabld(double *z33, float *xsh, float *ysh, double *peakv)
01492 {
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504 double c1,c2,c3,c4,c5,c6,denom;
01505 float xmin, ymin;
01506
01507 c1 = (26.*z33(1,1) - z33(1,2) + 2*z33(1,3) - z33(2,1) - 19.*z33(2,2)
01508 -7.*z33(2,3) + 2.*z33(3,1) - 7.*z33(3,2) + 14.*z33(3,3))/9.0;
01509
01510 c2 = (8.* z33(1,1) - 8.*z33(1,2) + 5.*z33(2,1) - 8.*z33(2,2) + 3.*z33(2,3)
01511 +2.*z33(3,1) - 8.*z33(3,2) + 6.*z33(3,3))/(-6.);
01512
01513 c3 = (z33(1,1) - 2.*z33(1,2) + z33(1,3) + z33(2,1) -2.*z33(2,2)
01514 + z33(2,3) + z33(3,1) - 2.*z33(3,2) + z33(3,3))/6.0;
01515
01516 c4 = (8.*z33(1,1) + 5.*z33(1,2) + 2.*z33(1,3) -8.*z33(2,1) -8.*z33(2,2)
01517 - 8.*z33(2,3) + 3.*z33(3,2) + 6.*z33(3,3))/(-6.0);
01518
01519 c5 = (z33(1,1) - z33(1,3) - z33(3,1) + z33(3,3))/4.0;
01520
01521 c6 = (z33(1,1) + z33(1,2) + z33(1,3) - 2.*z33(2,1) - 2.*z33(2,2)
01522 -2.*z33(2,3) + z33(3,1) + z33(3,2) + z33(3,3))/6.0;
01523
01524
01525
01526 *ysh=0.0;
01527 *xsh=0.0;
01528 denom=4.*c3*c6 - c5*c5;
01529 if (denom != 0.0) {
01530 *ysh=(c4*c5 - 2.*c2*c6) /denom-2.0;
01531 *xsh=(c2*c5 - 2.*c4*c3) /denom-2.0;
01532 *peakv= 4.*c1*c3*c6 - c1*c5*c5 -c2*c2*c6 + c2*c4*c5 - c4*c4*c3;
01533 *peakv= *peakv/denom;
01534
01535 xmin = min(*xsh,1.0);
01536 ymin = min(*ysh,1.0);
01537 *xsh=max(xmin,-1.0);
01538 *ysh=max(ymin,-1.0);
01539 }
01540 }
01541
01542 void crosrng_e(float *circ1, float *circ2, int lcirc,
01543 int nring, int maxrin, int *numr,
01544 double *qn, float *tot, int neg)
01545 {
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557 float *t;
01558 double t7[7], *q;
01559 int i, j, k, ip, jc, numr3i, numr2i, jtot;
01560 float pos;
01561
01562 ip = maxrin;
01563 q = (double*)calloc(maxrin+2, sizeof(double));
01564 t = (float*)calloc(maxrin+2, sizeof(float));
01565
01566 for (i=1;i<=nring;i++) {
01567 numr3i = numr(3,i);
01568 numr2i = numr(2,i);
01569
01570 t(1) = (circ1(numr2i)) * circ2(numr2i);
01571
01572 if (numr3i != maxrin) {
01573
01574 t(numr3i+1) = circ1(numr2i+1) * circ2(numr2i+1);
01575 t(2) = 0.0;
01576
01577 if (neg) {
01578
01579 for (j=3;j<=numr3i;j=j+2) {
01580 jc = j+numr2i-1;
01581 t(j) =(circ1(jc))*circ2(jc)-(circ1(jc+1))*circ2(jc+1);
01582 t(j+1) = -(circ1(jc))*circ2(jc+1)-(circ1(jc+1))*circ2(jc);
01583 }
01584 }
01585 else {
01586 for (j=3;j<=numr3i;j=j+2) {
01587 jc = j+numr2i-1;
01588 t(j) = (circ1(jc))*circ2(jc) + (circ1(jc+1))*circ2(jc+1);
01589 t(j+1) = -(circ1(jc))*circ2(jc+1) + (circ1(jc+1))*circ2(jc);
01590 }
01591 }
01592 for (j=1;j<=numr3i+1;j++) q(j) = q(j) + t(j);
01593 }
01594 else {
01595 t(2) = circ1(numr2i+1) * circ2(numr2i+1);
01596 if (neg) {
01597
01598 for (j=3;j<=maxrin;j=j+2) {
01599 jc = j+numr2i-1;
01600 t(j) = (circ1(jc))*circ2(jc) - (circ1(jc+1))*circ2(jc+1);
01601 t(j+1) = -(circ1(jc))*circ2(jc+1) - (circ1(jc+1))*circ2(jc);
01602 }
01603 }
01604 else {
01605 for (j=3;j<=maxrin;j=j+2) {
01606 jc = j+numr2i-1;
01607 t(j) = (circ1(jc))*circ2(jc) + (circ1(jc+1))*circ2(jc+1);
01608 t(j+1) = -(circ1(jc))*circ2(jc+1) + (circ1(jc+1))*circ2(jc);
01609 }
01610 }
01611 for (j = 1; j <= maxrin+2; j++) q(j) = q(j) + t(j);
01612 }
01613 }
01614
01615 fftr_d(q,ip);
01616
01617 *qn = -1.0e20;
01618 for (j=1;j<=maxrin;j++) {
01619 if (q(j) >= *qn) {
01620 *qn = q(j);
01621 jtot = j;
01622 }
01623 }
01624
01625 for (k=-3;k<=3;k++) {
01626 j = (jtot+k+maxrin-1)%maxrin + 1;
01627 t7(k+4) = q(j);
01628 }
01629
01630 prb1d(t7,7,&pos);
01631
01632 *tot = (float)jtot + pos;
01633
01634 if (q) free(q);
01635 if (t) free(t);
01636 }
01637
01638 int apmq(EMData *refprj, Dict refparams, EMData *expimg, APMQopt options,
01639 float *angles, float *shifts)
01640 {
01641 double quadpi = 3.14159265358979323846;
01642 double dgr_to_rad = quadpi/180;
01643
01644 int mr, nr, iskip;
01645 float range;
01646 char mode;
01647
01648 int nsam, nrow, numth, limitrange, maxrin, nimalcg, nring, status;
01649 int nx, ny, nxr, nyr, nidi, nima, lcirc, nwsam, nwrow, nsb, nse,
01650 nrb, nre, it, ishrange, istep, ckmirror, nrad, nang;
01651 int iref, ireft, mirrornew, imgref, ngotpar, mirrorold;
01652 float ccrot, rangnew, xshnew, yshnew, peakv;
01653 float rangout, angdif, rangold, xshold, yshold, c, s;
01654
01655 float *refstk, *imgstk, *bfc, *imgwindow, *refdirs, *expdirs,
01656 *angexps, *expdir;
01657
01658 vector <float> angrefs;
01659 float *dlist;
01660 int ldd = 7;
01661 int *numr;
01662
01663 float divas;
01664
01665 status = 0;
01666
01667
01668 mr = options.mr;
01669 nr = options.nr;
01670 ishrange = options.shrange;
01671 istep = options.istep;
01672 mode = options.mode;
01673 range = options.range;
01674
01675 angexps = options.angexps;
01676
01677 nx = expimg->get_xsize();
01678 ny = expimg->get_ysize();
01679 nidi = expimg->get_zsize();
01680
01681 nxr = refprj->get_xsize();
01682 nyr = refprj->get_ysize();
01683 nima = refprj->get_zsize();
01684
01685 if (nx != nxr || ny != nyr) {
01686 status = -2;
01687 goto EXIT;
01688 }
01689
01690 nsam = nx;
01691 nrow = ny;
01692
01693 nrad = min(nsam/2-1, nrow/2-1);
01694 if ( mr <=0 ) {
01695 fprintf(stderr,"first ring must be > 0\n");
01696 status = 10;
01697 goto EXIT;
01698 }
01699 if ( nr >= nrad) {
01700 fprintf(stderr,"last ring must be < %d\n", nrad);
01701 status = 10;
01702 goto EXIT;
01703 }
01704 if ( (ishrange+nr) > (nrad-1)) {
01705 fprintf(stderr,"last ring + translation must be < %d\n", nrad);
01706 status = 10;
01707 goto EXIT;
01708 }
01709
01710 ngotpar = 0;
01711
01712
01713 angrefs = refparams["anglelist"];
01714 nang = angrefs.size()/3;
01715 if (nang != nima) {
01716 fprintf(stderr, "apmd: nang = %d, nima = %d\n", nang, nima);
01717 status = -3;
01718 goto EXIT;
01719 }
01720
01721
01722 refstk = refprj->get_data();
01723
01724 iskip = 1;
01725 nring = setnring(mr, nr, iskip);
01726
01727 numr = (int*)calloc(3*nring,sizeof(int));
01728 if (!numr) {
01729 fprintf(stderr,"apmq: failed to allocate numr\n");
01730 status = -1;
01731 goto EXIT;
01732 }
01733
01734 numrinit(mr, nr, iskip, numr);
01735
01736
01737 status = alprbs(numr, nring, &lcirc, mode);
01738 if (status != 0) goto EXIT;
01739
01740 maxrin = numr(3,nring);
01741
01742
01743 bfc = (float*)calloc(lcirc*nima,sizeof(float));
01744
01745
01746 status = aprings(nima, nring, refstk, nsam, nrow, numr, bfc, lcirc, mode);
01747 if (status !=0 ) goto EXIT;
01748
01749
01750 imgstk = expimg->get_data();
01751
01752
01753 apmaster_1(mode,&divas,nr,&numth,nsam,nrow,&nwsam,&nwrow);
01754
01755
01756
01757
01758
01759 nwsam = nsam;
01760 nwrow = nrow;
01761 imgwindow = (float*)calloc(nwsam*nwrow,sizeof(float));
01762
01763
01764 nsb = -nsam/2;
01765 nse = -nsb-1+(nsam%2);
01766 nrb = -nrow/2;
01767 nre = -nrb-1+(nrow%2);
01768
01769 limitrange = 0;
01770 ckmirror = 1;
01771 if (range > 0.0 && range < 360.0) limitrange = 1;
01772 range = cos(range*dgr_to_rad);
01773 nimalcg = 1;
01774
01775 if (limitrange) {
01776
01777 refdirs = (float*)calloc(3*nima,sizeof(float));
01778 if (!refdirs) {
01779 fprintf(stderr, "apmq: failed to allocate refdirs!\n");
01780 status = -1;
01781 goto EXIT;
01782 }
01783
01784
01785
01786
01787
01788 angexps = options.angexps;
01789 if (!angexps) {
01790 fprintf(stderr, "apmq: no angexps available !\n");
01791 status = -4;
01792 goto EXIT;
01793 }
01794
01795 expdirs = (float*)calloc(3*nidi, sizeof(float));
01796 if (!expdirs) {
01797 fprintf(stderr, "apmq: failed to allocate expdirs!\n");
01798 status = -1;
01799 goto EXIT;
01800 }
01801
01802
01803
01804 }
01805
01806 dlist = (float*)calloc(ldd*nidi,sizeof(float));
01807 if (dlist == NULL || expdirs == NULL) {
01808 fprintf(stderr, "apmq: failed to allocated dlist or expdirs...\n");
01809 status = -1;
01810 goto EXIT;
01811 }
01812
01813
01814 for (it=1;it<=nidi;it++) {
01815
01816 win_resize(&imgstk(1,1,it), imgwindow, nsam, nrow, nwsam, nwrow,
01817 1, nrow, 1, nsam);
01818 printf("it = %d\n", it);
01819
01820 if (limitrange) expdir = &expdirs(1,it);
01821
01822 status = aprq2d(imgwindow , bfc , numr , nsam ,
01823 nrow , ishrange , istep, nsb ,
01824 nse , nrb , nre , lcirc ,
01825 nring , maxrin , nima , mode ,
01826 refdirs , expdir , range, &dlist(2,it),
01827 &dlist(3,it), &dlist(4,it), &dlist(5,it),
01828 &dlist(6,it), &nimalcg , ckmirror ,
01829 limitrange);
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847 iref = (int)dlist(2,it);
01848 if (iref < 0) {
01849
01850 imgref = -iref;
01851
01852
01853 ireft = -iref;
01854 mirrornew = 1;
01855 }
01856 else if (iref == 0) {
01857
01858 imgref = 0;
01859
01860 ireft = 1;
01861 mirrornew = 0;
01862 }
01863 else {
01864 imgref = iref;
01865
01866 ireft = iref;
01867 mirrornew = 0;
01868 }
01869
01870 ccrot = dlist(3,it);
01871 rangnew = dlist(4,it);
01872 xshnew = dlist(5,it);
01873 yshnew = dlist(6,it);
01874 peakv = 1.0;
01875
01876
01877 if (imgref <= 0) {
01878
01879 ccrot = -1.0;
01880 peakv = 0.0;
01881 }
01882
01883
01884 angles(1,it) = 0.0;
01885 angles(2,it) = 0.0;
01886 angles(2,it) = 0.0;
01887
01888 if (imgref > 0) {
01889
01890 angles(1,it) = angrefs(1,iref);
01891 angles(2,it) = angrefs(2,iref);
01892 angles(3,it) = angrefs(3,iref);
01893
01894 if (mirrornew) {
01895
01896 angles(1,it) = -angles(1,it);
01897 angles(2,it) = 180+angles(2,it);
01898 }
01899 }
01900 else if (ngotpar >= 3) {
01901
01902 angles(1,it) = angexps(1,it);
01903 angles(2,it) = angexps(2,it);
01904 angles(3,it) = angexps(3,it);
01905 }
01906
01907 rangold = 0.0;
01908 xshold = 0.0;
01909 yshold = 0.0;
01910
01911 if (ngotpar >= 7 && ishrange > 0) {
01912
01913 rangold = angexps(4,it);
01914 xshold = angexps(5,it);
01915 yshold = angexps(6,it);
01916
01917 mirrorold = 0;
01918 if (angexps(7,it) > 0) mirrorold = 1;
01919 if (mirrorold) {
01920 status = 5;
01921 fprintf(stderr, "mirrorold = %d\n", mirrorold);
01922 goto EXIT;
01923 }
01924 }
01925
01926
01927
01928 c = cos(rangnew * dgr_to_rad);
01929 s = -sin(rangnew * dgr_to_rad);
01930
01931 shifts(1,it) = xshnew + xshold*c - yshold*s;
01932 shifts(2,it) = yshnew + xshold*s + yshold*c;
01933 rangout = rangold + rangnew;
01934
01935
01936 while ( rangout < 0.0) {
01937 rangout = rangout + 360.0;
01938 }
01939 while(rangout >= 360.0) {
01940 rangout = rangout - 360.0;
01941 }
01942
01943
01944 angdif = -1.0;
01945
01946 if (imgref <= 0) {
01947
01948 angdif = 0.0;
01949 }
01950 else if (ngotpar >= 3) {
01951
01952 angdif = fabs(expdirs(1,it) * refdirs(1,iref) +
01953 expdirs(2,it) * refdirs(2,iref) +
01954 expdirs(3,it) * refdirs(3,iref));
01955 angdif = min(1.0,angdif);
01956 angdif = acos(angdif) / dgr_to_rad;
01957 }
01958 }
01959
01960 EXIT:
01961 free(numr);
01962 if (limitrange) {
01963 free(refdirs);
01964 free(expdirs);
01965 }
01966 free(dlist);
01967 free(bfc);
01968 free(imgwindow);
01969 return status;
01970 }