This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Classes | |
struct | APMQopt |
struct | APMDopt |
Defines | |
#define | angles(i, j) angles [((j)-1)*3 + (i)-1] |
#define | newangles(i, j) newangles [((j)-1)*3 + (i)-1] |
#define | angrefs(i, j) angrefs [((j)-1)*3 + (i)-1] |
#define | angexps(i, j) angexps [((j)-1)*3 + (i)-1] |
#define | shifts(i, j) shifts[((j)-1)*2 + (i)-1] |
#define | fangles(i, j) fangles [((j)-1)*3 + (i)-1] |
#define | numr(i, j) numr [((j)-1)*3 + (i)-1] |
#define | imgcirc(i) imgcirc[(i)-1] |
#define | wr(i) wr [(i)-1] |
#define | sqimg(i, j) sqimg [((j)-1)*nsam + (i)-1] |
#define | xim(i, j) xim [((j)-1)*nsam + (i)-1] |
#define | fdata(i, j) fdata [((j)-1)*nxdata + (i)-1] |
#define | min0(a, b) ((a) >= (b) ? (b) : (a)) |
#define | min(a, b) ((a) >= (b) ? (b) : (a)) |
#define | max(a, b) ((a) <= (b) ? (b) : (a)) |
#define | tab1(i) tab1 [(i)-1] |
#define | xcmplx(i, j) xcmplx [((j)-1)*2 + (i)-1] |
#define | br(i) br [(i)-1] |
#define | bi(i) bi [(i)-1] |
#define | circ(i) circ [(i)-1] |
#define | circ1(i) circ1 [(i)-1] |
#define | circ2(i) circ2 [(i)-1] |
#define | t(i) t [(i)-1] |
#define | q(i) q [(i)-1] |
#define | b(i) b [(i)-1] |
#define | t7(i) t7 [(i)-1] |
#define | imgfrom(i, j) imgfrom[((j)-1)*lsam + (i)-1] |
#define | imgto(i, j) imgto [((j)-1)*nsam + (i)-1] |
#define | imgstk(i, j, k) imgstk[((k)-1)*nsam*nrow + ((j)-1)*nsam + (i)-1] |
#define | refcstk(i, j) refcstk[((j)-1)*lcirc + (i) - 1] |
#define | imgwindow(i, j) imgwindow [((j)-1)*nwsam + (i)-1] |
#define | totmin(i) totmin[(i)-1] |
#define | totmir(i) totmir[(i)-1] |
#define | tot(i) tot[(i)-1] |
#define | tmt(i) tmt[(i)-1] |
#define | dlist(i, j) dlist[((j)-1)*ldd + (i)-1] |
#define | expdir(i) expdir[(i)-1] |
#define | expdirs(i, j) expdirs[((j)-1)*3 + (i)-1] |
#define | refdirs(i, j) refdirs[((j)-1)*3 + (i)-1] |
#define | refdir(i) refdir[(i)-1] |
#define | lcg(i) lcg[(i)-1] |
#define | bfc(i, j) bfc[((j)-1)*lcirc + (i) - 1] |
#define | fitp(i, j) fitp[ ((j)+1)*3 + (i) + 1] |
#define | fit(i, j) fit[((j)+istep)*(2*istep+1) + (i) + istep] |
#define | rotmp(i, j) rotmp[((j)+istep)*(2*istep+1) + (i) + istep] |
#define | z33(i, j) z33[((j)-1)*3 + (i)-1] |
Functions | |
int | aprings (int nimg, int nring, float *imgstk, int nsam, int nrow, int *numr, float *refcstk, int lcirc, char mode) |
int | apring1 (float *sqimg, int nsam, int nrow, float *imgcirc, int lcirc, int nring, char mode, int *numr, float *wr) |
float | quadri (float xx, float yy, int nxdata, int nydata, float *fdata) |
Quadratic interpolation (2D). | |
void | ringwe (float *wr, int *numr, int nring) |
int | alprbs (int *numr, int nring, int *lcirc, char mode) |
int | setnring (int mr, int nr, int iskip) |
void | numrinit (int mr, int nr, int iskip, int *numr) |
void | normass (float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2, int ir1, int ir2) |
void | normas (float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2, int ir1, int ir2) |
void | frngs (float *circ, int *numr, int nring) |
void | fftr_q (float *xcmplx, int nv) |
void | fftc_q (float *br, float *bi, int ln, int ks) |
int | applyws (float *circ, int lcirc, int *numr, float *wr, int nring) |
void | alrq (float *xim, int nsam, int nrow, int *numr, float *circ, int lcirc, int nring, char mode) |
void | alrq_ms (float *xim, int nsam, int nrow, float cns2, float cnr2, int *numr, float *circ, int lcirc, int nring, char mode) |
int | crosrng_ms (float *circ1, float *circ2, int lcirc, int nring, int maxrin, int *numr, double *qn, double *tot, double *qm, double *tmt) |
void | prb1d (double *b, int npoint, float *pos) |
void | apmaster_1 (char mode, float *divas, int nr, int *numth, int lsam, int lrow, int *nsam, int *nrow) |
void | win_resize (float *imgfrom, float *imgto, int lsam, int lrow, int nsam, int nrow, int lr1, int lr2, int ls1, int ls2) |
int | apmd (EMData *refprj, Dict refparam, EMData *expimg, APMDopt options, float *fangle) |
float | ang_n (float rkk, char mode, int maxrin) |
void | parabld (double *z33, float *xsh, float *ysh, double *peakv) |
void | crosrng_e (float *circ1, float *circ2, int lcirc, int nring, int maxrin, int *numr, double *qn, float *tot, int neg) |
int | apmq (EMData *refprj, Dict refparams, EMData *expimg, APMQopt options, float *angles, float *shifts) |
int | aprq2d (float *sqimg, float *bfc, int *numr, int nsam, int nrow, int ishrange, int istep, int nsb, int nse, int nrb, int nre, int lcirc, int nring, int maxrin, int nima, char mode, float *refdir, float *expdir, float range, float *diref, float *ccrot, float *rangnew, float *xshsum, float *yshsum, int *nimalcg, int ckmirror, int limitrange) |
#define angexps | ( | i, | |||
j | ) | angexps [((j)-1)*3 + (i)-1] |
Definition at line 32 of file spidutil.h.
Referenced by apmq(), EMAN::Util::even_angles(), EMAN::OptimumOrientationGenerator::gen_orientations(), and main().
#define angrefs | ( | i, | |||
j | ) | angrefs [((j)-1)*3 + (i)-1] |
#define b | ( | i | ) | b [(i)-1] |
Definition at line 56 of file spidutil.h.
#define bfc | ( | i, | |||
j | ) | bfc[((j)-1)*lcirc + (i) - 1] |
#define bi | ( | i | ) | bi [(i)-1] |
Definition at line 50 of file spidutil.h.
#define br | ( | i | ) | br [(i)-1] |
Definition at line 49 of file spidutil.h.
#define circ | ( | i | ) | circ [(i)-1] |
Definition at line 51 of file spidutil.h.
#define circ1 | ( | i | ) | circ1 [(i)-1] |
Definition at line 52 of file spidutil.h.
#define circ2 | ( | i | ) | circ2 [(i)-1] |
Definition at line 53 of file spidutil.h.
#define dlist | ( | i, | |||
j | ) | dlist[((j)-1)*ldd + (i)-1] |
#define expdir | ( | i | ) | expdir[(i)-1] |
#define expdirs | ( | i, | |||
j | ) | expdirs[((j)-1)*3 + (i)-1] |
#define fangles | ( | i, | |||
j | ) | fangles [((j)-1)*3 + (i)-1] |
#define fdata | ( | i, | |||
j | ) | fdata [((j)-1)*nxdata + (i)-1] |
Definition at line 43 of file spidutil.h.
#define fit | ( | i, | |||
j | ) | fit[((j)+istep)*(2*istep+1) + (i) + istep] |
Definition at line 75 of file spidutil.h.
Referenced by aprq2d(), EMAN::Symmetry3D::cache_au_planes(), and EMAN::Symmetry3D::delete_au_planes().
#define fitp | ( | i, | |||
j | ) | fitp[ ((j)+1)*3 + (i) + 1] |
#define imgcirc | ( | i | ) | imgcirc[(i)-1] |
#define imgfrom | ( | i, | |||
j | ) | imgfrom[((j)-1)*lsam + (i)-1] |
#define imgstk | ( | i, | |||
j, | |||||
k | ) | imgstk[((k)-1)*nsam*nrow + ((j)-1)*nsam + (i)-1] |
#define imgto | ( | i, | |||
j | ) | imgto [((j)-1)*nsam + (i)-1] |
#define imgwindow | ( | i, | |||
j | ) | imgwindow [((j)-1)*nwsam + (i)-1] |
#define lcg | ( | i | ) | lcg[(i)-1] |
#define max | ( | a, | |||
b | ) | ((a) <= (b) ? (b) : (a)) |
Definition at line 46 of file spidutil.h.
#define min | ( | a, | |||
b | ) | ((a) >= (b) ? (b) : (a)) |
Definition at line 45 of file spidutil.h.
#define min0 | ( | a, | |||
b | ) | ((a) >= (b) ? (b) : (a)) |
#define newangles | ( | i, | |||
j | ) | newangles [((j)-1)*3 + (i)-1] |
#define numr | ( | i, | |||
j | ) | numr [((j)-1)*3 + (i)-1] |
Definition at line 38 of file spidutil.h.
#define q | ( | i | ) | q [(i)-1] |
Definition at line 55 of file spidutil.h.
#define refcstk | ( | i, | |||
j | ) | refcstk[((j)-1)*lcirc + (i) - 1] |
#define refdir | ( | i | ) | refdir[(i)-1] |
Definition at line 71 of file spidutil.h.
#define refdirs | ( | i, | |||
j | ) | refdirs[((j)-1)*3 + (i)-1] |
#define rotmp | ( | i, | |||
j | ) | rotmp[((j)+istep)*(2*istep+1) + (i) + istep] |
#define shifts | ( | i, | |||
j | ) | shifts[((j)-1)*2 + (i)-1] |
#define sqimg | ( | i, | |||
j | ) | sqimg [((j)-1)*nsam + (i)-1] |
#define t | ( | i | ) | t [(i)-1] |
Definition at line 54 of file spidutil.h.
#define t7 | ( | i | ) | t7 [(i)-1] |
Definition at line 57 of file spidutil.h.
#define tab1 | ( | i | ) | tab1 [(i)-1] |
Definition at line 47 of file spidutil.h.
#define tmt | ( | i | ) | tmt[(i)-1] |
Definition at line 66 of file spidutil.h.
Referenced by apmd(), aprq2d(), EMAN::Util::Crosrng_ms(), EMAN::Util::Crosrng_ms_delta(), and EMAN::Util::Crosrng_psi().
#define tot | ( | i | ) | tot[(i)-1] |
Definition at line 65 of file spidutil.h.
Referenced by apmd(), aprq2d(), wustl_mm::SkeletonMaker::Volume::components26(), wustl_mm::SkeletonMaker::Volume::components6(), EMAN::Util::constrained_helix(), EMAN::Util::constrained_helix_test(), EMAN::Util::Crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), EMAN::Util::Crosrng_ms_delta(), EMAN::Util::Crosrng_ns(), EMAN::Util::Crosrng_psi(), EMAN::Util::Crosrng_sm_psi(), wustl_mm::SkeletonMaker::Volume::hasCompleteSheet(), and wustl_mm::SkeletonMaker::Volume::isSheetEnd().
#define totmin | ( | i | ) | totmin[(i)-1] |
#define totmir | ( | i | ) | totmir[(i)-1] |
#define wr | ( | i | ) | wr [(i)-1] |
#define xcmplx | ( | i, | |||
j | ) | xcmplx [((j)-1)*2 + (i)-1] |
Definition at line 48 of file spidutil.h.
#define xim | ( | i, | |||
j | ) | xim [((j)-1)*nsam + (i)-1] |
Definition at line 42 of file spidutil.h.
#define z33 | ( | i, | |||
j | ) | z33[((j)-1)*3 + (i)-1] |
int alprbs | ( | int * | numr, | |
int | nring, | |||
int * | lcirc, | |||
char | mode | |||
) |
Definition at line 309 of file spidutil.cpp.
References EMAN::MAXFFT, min0, numr, pi, and status.
Referenced by apmd(), and apmq().
00310 { 00311 /* 00312 c purpose: appears to circular rings, postitioned 00313 c in a linear array that holds rings concatenated together. 00314 c output is dependent on number of rings 00315 c * 00316 c parameters: numr(1,i) - ring number (sent) 00317 c numr(2,i) - beginning in circ (ret.) 00318 c numr(3,i) - length in circ (ret.) 00319 c nring (sent) 00320 c lcirc - total length of circ. (ret.) 00321 c 00322 c image_processing_routine 00323 c 00324 */ 00325 int i, jp, ip; 00326 double dpi; 00327 int status = 0; 00328 // hardwire for right now 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 // original fortran code ip = 2**log2(jp), log2(jp) rounds up. 00338 ip = (int)( pow(2,ceil(log2(jp))) ); 00339 if (i < nring && jp > ip+ip/2) ip=min0(MAXFFT,2*ip); 00340 00341 // last ring should be oversampled to allow higher accuracy 00342 // of peak location (?). 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 }
void alrq | ( | float * | xim, | |
int | nsam, | |||
int | nrow, | |||
int * | numr, | |||
float * | circ, | |||
int | lcirc, | |||
int | nring, | |||
char | mode | |||
) |
Definition at line 558 of file spidutil.cpp.
References circ, numr, quadri(), x, and y.
Referenced by apmd().
00560 { 00561 /* 00562 c 00563 c purpose: 00564 c 00565 c parameters: convert to polar coordinates 00566 c 00567 */ 00568 // dimension xim(nsam,nrow),circ(lcirc) 00569 // integer numr(3,nring) 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 //#pragma omp parallel do private(i,j,inr,yq,l,lt,nsim,dfi,kcirc, 00580 //#pragma omp& xold,yold,fi,x,y) 00581 for (i=1;i<=nring;i++) { 00582 // radius of the ring 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 //#pragma omp end parallel do 00634 }
void alrq_ms | ( | float * | xim, | |
int | nsam, | |||
int | nrow, | |||
float | cns2, | |||
float | cnr2, | |||
int * | numr, | |||
float * | circ, | |||
int | lcirc, | |||
int | nring, | |||
char | mode | |||
) |
Definition at line 1419 of file spidutil.cpp.
References circ, numr, quadri(), x, and y.
Referenced by aprq2d().
01421 { 01422 double dpi, dfi; 01423 int it, jt, inr, l, nsim, kcirc, lt; 01424 float yq, xold, yold, fi, x, y; 01425 01426 // cns2 and cnr2 are predefined centers 01427 // no need to set to zero, all elements are defined 01428 01429 dpi = 2*atan(1.0); 01430 for (it=1;it<=nring;it++) { 01431 // radius of the ring 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 } // end for jt 01488 } //end for it 01489 }
float ang_n | ( | float | rkk, | |
char | mode, | |||
int | maxrin | |||
) |
Definition at line 1406 of file spidutil.cpp.
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 }
void apmaster_1 | ( | char | mode, | |
float * | divas, | |||
int | nr, | |||
int * | numth, | |||
int | lsam, | |||
int | lrow, | |||
int * | nsam, | |||
int * | nrow | |||
) |
Definition at line 820 of file spidutil.cpp.
Referenced by apmd(), and apmq().
00822 { 00823 /* 00824 c parameters: 00825 c mode degree mode (input) 00826 c divas degrees (output) 00827 c numth degrees (output) 00828 c lsam orig size (input) 00829 c lrow orig size (input) 00830 c nsam new size (output) 00831 c nrow new size (output) 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 // find number of omp threads 00845 // call getthreads(numth) 00846 #endif 00847 00848 // calculation of actual dimension of an image to be interpolated 00849 // 2*(no.of rings)+(0'th element)+2*(margin of 1) 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 // returns circular reduced nsam, nrow 00856 *nsam = nra; 00857 *nrow = nra; 00858 }
int apmd | ( | EMData * | refprj, | |
Dict | refparam, | |||
EMData * | expimg, | |||
APMDopt | options, | |||
float * | fangle | |||
) |
int apmq | ( | EMData * | refprj, | |
Dict | refparams, | |||
EMData * | expimg, | |||
APMQopt | options, | |||
float * | angles, | |||
float * | shifts | |||
) |
int applyws | ( | float * | circ, | |
int | lcirc, | |||
int * | numr, | |||
float * | wr, | |||
int | nring | |||
) |
Definition at line 474 of file spidutil.cpp.
References circ, numr, status, and wr.
Referenced by apring1().
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 }
int apring1 | ( | float * | sqimg, | |
int | nsam, | |||
int | nrow, | |||
float * | imgcirc, | |||
int | lcirc, | |||
int | nring, | |||
char | mode, | |||
int * | numr, | |||
float * | wr | |||
) |
Definition at line 89 of file spidutil.cpp.
References applyws(), frngs(), imgcirc, normass(), numr, quadri(), sqimg, status, wr, x, and y.
Referenced by aprings().
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 // calculate dimensions for normass 00101 nsb = -nsam/2; 00102 nse = -nsb-1+(nsam%2); 00103 nrb = -nrow/2; 00104 nre = -nrb-1+(nrow%2); 00105 00106 // normalize under the mask, tried doing this on the 00107 // polar rings but it gives different answers. al 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 // convert window from image into polar coordinates 00118 00119 if (mode == 'f' || mode == 'F') { 00120 ltf = 4; 00121 for (i=1;i<=nring;i++) { 00122 // radius of the ring 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 // radius of the ring 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 // fourier of circ 00175 frngs(imgcirc,numr,nring); 00176 00177 // weight circ using wr 00178 if (wr(1) > 0.0) { 00179 status = applyws(imgcirc,lcirc,numr,wr,nring); 00180 } 00181 00182 return status; 00183 }
int aprings | ( | int | nimg, | |
int | nring, | |||
float * | imgstk, | |||
int | nsam, | |||
int | nrow, | |||
int * | numr, | |||
float * | refcstk, | |||
int | lcirc, | |||
char | mode | |||
) |
Definition at line 57 of file spidutil.cpp.
References apring1(), imgstk, refcstk, ringwe(), status, and wr.
Referenced by apmd(), and apmq().
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 // get wr weights 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 }
int aprq2d | ( | float * | sqimg, | |
float * | bfc, | |||
int * | numr, | |||
int | nsam, | |||
int | nrow, | |||
int | ishrange, | |||
int | istep, | |||
int | nsb, | |||
int | nse, | |||
int | nrb, | |||
int | nre, | |||
int | lcirc, | |||
int | nring, | |||
int | maxrin, | |||
int | nima, | |||
char | mode, | |||
float * | refdir, | |||
float * | expdir, | |||
float | range, | |||
float * | diref, | |||
float * | ccrot, | |||
float * | rangnew, | |||
float * | xshsum, | |||
float * | yshsum, | |||
int * | nimalcg, | |||
int | ckmirror, | |||
int | limitrange | |||
) |
Definition at line 1099 of file spidutil.cpp.
References abs, alrq_ms(), ang_n(), bfc, crosrng_e(), crosrng_ms(), dgr_to_rad, dt, expdir, fit, fitp, frngs(), imgcirc, lcg, normass(), numr, parabld(), quadpi, refdirs, rotmp, status, tmt, and tot.
Referenced by apmq().
01108 : 01109 c diref number of most similar ref. proj. (output) 01110 c (negative if mirrored) 01111 c ccrot corr coeff. (output) 01112 c rangnew inplane angle (output) 01113 c xshsum shift (output) 01114 c yshsum shift (output) 01115 c nimalcg (output) 01116 c 01117 dimension a(nsam,nrow),bfc(lcirc,nima),numr(3,nring) 01118 double precision fitp(-1:1,-1:1) 01119 double precision, dimension(*) :: tt 01120 real, dimension(3,nima) :: refdirs 01121 real, dimension(3) :: expdir 01122 automatic arrays 01123 double precision fit(-istep:istep,-istep:istep) 01124 dimension rotmp(-istep:istep,-istep:istep) 01125 real, dimension(lcirc) :: a_circ 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 // restricted range search 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 // dt near 1.0 = not-mirrored, dt near -1.0 = mirrored 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 // mirored or non-mirrored is within range 01190 *nimalcg++; 01191 lcg(*nimalcg) = imi; 01192 if (dt < 0) lcg(*nimalcg) = -imi; 01193 } 01194 } 01195 01196 if (*nimalcg <= 0) { 01197 // there is no reference within search range 01198 *xshsum = 0; 01199 *yshsum = 0; 01200 *diref = 0; 01201 *rangnew = 0; 01202 *ccrot = -1.0; 01203 goto EXIT; 01204 } 01205 iend = *nimalcg; 01206 // end of restricted range search 01207 } 01208 01209 ccrotd = -1.0e23; 01210 01211 // go through centers for shift alignment 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 // normalize under the mask 01218 //'normalize' image values over variance range 01219 normass(sqimg,nsam,nsb-it,nse-it,nrb-jt,nre-jt,numr(1,1), 01220 numr(1,nring)); 01221 01222 // interpolation into polar coordinates 01223 // creates imgcirc (exp. image circles) for this position 01224 01225 alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,nring,mode); 01226 01227 // creates fourier of: a_circ 01228 frngs(imgcirc,numr,nring); 01229 01230 // compare exp. image with all reference images 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 // check either mirrored or non-mirrored position 01240 crosrng_e(&bfc(1,ir),imgcirc,lcirc,nring, 01241 maxrin,numr,&tota,&tot,mirrored); 01242 } 01243 else { 01244 // check both non-mirrored & mirrored positions 01245 status = crosrng_ms(&bfc(1,ir),imgcirc,lcirc,nring, 01246 maxrin,numr,&tota,&tot,&tmta,&tmt); 01247 } 01248 } 01249 else { 01250 // do not check mirrored position 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 // good match with tota (mirrored or not) position 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 // have to compare with mirrored position 01269 if (tmta >= ccrotd) { 01270 // good match with mirrored position 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 } // endfor irr 01280 } // endfor it 01281 } // endfor jt 01282 01283 // try to interpolate 01284 *ccrot = ccrotd; 01285 sx = isx; 01286 sy = isy; 01287 *diref = idis; 01288 01289 // do not interpolate for point on the edge 01290 if ( abs(isx) != ishrange && abs(isy) != ishrange) { 01291 // have to find neighbouring values 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 // if (idis .lt. 0) check mirrored only 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 } // endfor it 01318 } //endfor jt 01319 01320 // find the maximum within +/-istep 01321 // maximum cannot be on the edge, i.e., it,jt/=istep 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 // temp variable overcomes compiler bug on opt 64 pgi 6.0 01339 *rangnew = rangnewt; 01340 01341 // copy values around the peak. 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 // update location of the peak 01347 ccrotd = afit; 01348 isx = isx+itma; 01349 isy = isy+jtma; 01350 parabld(fitp,&sx,&sy,&peak); 01351 01352 // check whether interpolation is ok. 01353 if (fabs(sx) < 1.0 && fabs(sy) < 1.0) { 01354 // not on edge of 3x3 area 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 // not on edge of 3x3 area 01378 sx = isx; 01379 sy = isy; 01380 } 01381 } 01382 01383 sx = -sx; 01384 sy = -sy; 01385 01386 // now have to change order of shift & rotation. 01387 // in this program image is shifted first, rotated second. 01388 // in 'rt sq' it is rotation first, shift second. 01389 // this part corresponds to 'sa p'. 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 }
void crosrng_e | ( | float * | circ1, | |
float * | circ2, | |||
int | lcirc, | |||
int | nring, | |||
int | maxrin, | |||
int * | numr, | |||
double * | qn, | |||
float * | tot, | |||
int | neg | |||
) |
Definition at line 1542 of file spidutil.cpp.
References circ1, circ2, fftr_d(), numr, prb1d(), q, t, and t7.
Referenced by aprq2d().
01545 { 01546 /* 01547 c checks single position, neg is flag for checking mirrored position 01548 c 01549 c input - fourier transforms of rings! 01550 c first set is conjugated (mirrored) if neg 01551 c circ1 already multiplied by weights! 01552 c automatic arrays 01553 dimension t(maxrin+2) 01554 double precision q(maxrin+2) 01555 double precision t7(-3:3) 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 // test .ne. first for speed on some compilers 01574 t(numr3i+1) = circ1(numr2i+1) * circ2(numr2i+1); 01575 t(2) = 0.0; 01576 01577 if (neg) { 01578 // first set is conjugated (mirrored) 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 // first set is conjugated (mirrored) 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 }
int crosrng_ms | ( | float * | circ1, | |
float * | circ2, | |||
int | lcirc, | |||
int | nring, | |||
int | maxrin, | |||
int * | numr, | |||
double * | qn, | |||
double * | tot, | |||
double * | qm, | |||
double * | tmt | |||
) |
void fftc_q | ( | float * | br, | |
float * | bi, | |||
int | ln, | |||
int | ks | |||
) |
Definition at line 135 of file spidfft.cpp.
References abs, bi, br, status, t, and tab1.
00136 { 00137 // dimension br(1),bi(1) 00138 00139 int b3,b4,b5,b6,b7,b56; 00140 int n, k, l, j, i, ix0, ix1; 00141 float rni, tr1, ti1, tr2, ti2, cc, c, ss, s, t, x2, x3, x4, x5, sgn; 00142 float tab1[15]; 00143 int status=0; 00144 00145 tab1(1)=9.58737990959775e-5; 00146 tab1(2)=1.91747597310703e-4; 00147 tab1(3)=3.83495187571395e-4; 00148 tab1(4)=7.66990318742704e-4; 00149 tab1(5)=1.53398018628476e-3; 00150 tab1(6)=3.06795676296598e-3; 00151 tab1(7)=6.13588464915449e-3; 00152 tab1(8)=1.22715382857199e-2; 00153 tab1(9)=2.45412285229123e-2; 00154 tab1(10)=4.90676743274181e-2; 00155 tab1(11)=9.80171403295604e-2; 00156 tab1(12)=1.95090322016128e-1; 00157 tab1(13)=3.82683432365090e-1; 00158 tab1(14)=7.07106781186546e-1; 00159 tab1(15)=1.00000000000000; 00160 00161 n=(int)pow(2,ln); 00162 00163 k=abs(ks); 00164 l=16-ln; 00165 b3=n*k; 00166 b6=b3; 00167 b7=k; 00168 if( ks > 0 ) { 00169 sgn=1.0; 00170 } 00171 else { 00172 sgn=-1.0; 00173 rni=1.0/(float)n; 00174 j=1; 00175 for (i=1; i<=n;i++) { 00176 br(j)=br(j)*rni; 00177 bi(j)=bi(j)*rni; 00178 j=j+k; 00179 } 00180 } 00181 L12: 00182 b6=b6/2; 00183 b5=b6; 00184 b4=2*b6; 00185 b56=b5-b6; 00186 L14: 00187 tr1=br(b5+1); 00188 ti1=bi(b5+1); 00189 tr2=br(b56+1); 00190 ti2=bi(b56+1); 00191 00192 br(b5+1)=tr2-tr1; 00193 bi(b5+1)=ti2-ti1; 00194 br(b56+1)=tr1+tr2; 00195 bi(b56+1)=ti1+ti2; 00196 00197 b5=b5+b4; 00198 b56=b5-b6; 00199 if (b5 <= b3) goto L14; 00200 if (b6 == b7) goto L20; 00201 00202 b4=b7; 00203 cc=2.0*pow(tab1(l),2); 00204 c=1.0-cc; 00205 l=l+1; 00206 ss=sgn*tab1(l); 00207 s=ss; 00208 L16: 00209 b5=b6+b4; 00210 b4=2*b6; 00211 b56=b5-b6; 00212 L18: 00213 tr1=br(b5+1); 00214 ti1=bi(b5+1); 00215 tr2=br(b56+1); 00216 ti2=bi(b56+1); 00217 br(b5+1)=c*(tr2-tr1)-s*(ti2-ti1); 00218 bi(b5+1)=s*(tr2-tr1)+c*(ti2-ti1); 00219 br(b56+1)=tr1+tr2; 00220 bi(b56+1)=ti1+ti2; 00221 00222 b5=b5+b4; 00223 b56=b5-b6; 00224 if(b5 <= b3) goto L18; 00225 b4=b5-b6; 00226 b5=b4-b3; 00227 c=-c; 00228 b4=b6-b5; 00229 if(b5 < b4) goto L16; 00230 b4=b4+b7; 00231 if(b4 >= b5) goto L12; 00232 00233 t=c-cc*c-ss*s; 00234 s=s+ss*c-cc*s; 00235 c=t; 00236 goto L16; 00237 L20: 00238 ix0=b3/2; 00239 b3=b3-b7; 00240 b4=0; 00241 b5=0; 00242 b6=ix0; 00243 ix1=0; 00244 if ( b6 == b7) goto EXIT; 00245 L22: 00246 b4=b3-b4; 00247 b5=b3-b5; 00248 x2=br(b4+1); 00249 x3=br(b5+1); 00250 x4=bi(b4+1); 00251 x5=bi(b5+1); 00252 br(b4+1)=x3; 00253 br(b5+1)=x2; 00254 bi(b4+1)=x5; 00255 bi(b5+1)=x4; 00256 if (b6 < b4) goto L22; 00257 L24: 00258 b4=b4+b7; 00259 b5=b6+b5; 00260 x2=br(b4+1); 00261 x3=br(b5+1); 00262 x4=bi(b4+1); 00263 x5=bi(b5+1); 00264 br(b4+1)=x3; 00265 br(b5+1)=x2; 00266 bi(b4+1)=x5; 00267 bi(b5+1)=x4; 00268 ix0=b6; 00269 L26: 00270 ix0=ix0/2; 00271 ix1=ix1-ix0; 00272 if(ix1 >= 0) goto L26; 00273 00274 ix0=2*ix0; 00275 b4=b4+b7; 00276 ix1=ix1+ix0; 00277 b5=ix1; 00278 if (b5 >= b4) goto L22; 00279 if (b4 < b6) goto L24; 00280 EXIT: 00281 status = 0; 00282 }
void fftr_q | ( | float * | xcmplx, | |
int | nv | |||
) |
Definition at line 53 of file spidfft.cpp.
References abs, fftc_q(), t, tab1, and xcmplx.
00054 { 00055 // dimension xcmplx(2,1); xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary 00056 00057 float tab1[15]; 00058 int nu, inv, nu1, n, isub, n2, i1, i2, i; 00059 float ss, cc, c, s, tr, ti, tr1, tr2, ti1, ti2, t; 00060 00061 tab1(1)=9.58737990959775e-5; 00062 tab1(2)=1.91747597310703e-4; 00063 tab1(3)=3.83495187571395e-4; 00064 tab1(4)=7.66990318742704e-4; 00065 tab1(5)=1.53398018628476e-3; 00066 tab1(6)=3.06795676296598e-3; 00067 tab1(7)=6.13588464915449e-3; 00068 tab1(8)=1.22715382857199e-2; 00069 tab1(9)=2.45412285229123e-2; 00070 tab1(10)=4.90676743274181e-2; 00071 tab1(11)=9.80171403295604e-2; 00072 tab1(12)=1.95090322016128e-1; 00073 tab1(13)=3.82683432365090e-1; 00074 tab1(14)=7.07106781186546e-1; 00075 tab1(15)=1.00000000000000; 00076 00077 nu=abs(nv); 00078 inv=nv/nu; 00079 nu1=nu-1; 00080 n=(int)pow(2,nu1); 00081 isub=16-nu1; 00082 00083 ss=-tab1(isub); 00084 cc=-2.0*pow(tab1(isub-1),2); 00085 c=1.0; 00086 s=0.0; 00087 n2=n/2; 00088 if ( inv > 0) { 00089 fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,2); 00090 tr=xcmplx(1,1); 00091 ti=xcmplx(2,1); 00092 xcmplx(1,1)=tr+ti; 00093 xcmplx(2,1)=tr-ti; 00094 for (i=1;i<=n2;i++) { 00095 i1=i+1; 00096 i2=n-i+1; 00097 tr1=xcmplx(1,i1); 00098 tr2=xcmplx(1,i2); 00099 ti1=xcmplx(2,i1); 00100 ti2=xcmplx(2,i2); 00101 t=(cc*c-ss*s)+c; 00102 s=(cc*s+ss*c)+s; 00103 c=t; 00104 xcmplx(1,i1)=0.5*((tr1+tr2)+(ti1+ti2)*c-(tr1-tr2)*s); 00105 xcmplx(1,i2)=0.5*((tr1+tr2)-(ti1+ti2)*c+(tr1-tr2)*s); 00106 xcmplx(2,i1)=0.5*((ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00107 xcmplx(2,i2)=0.5*(-(ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00108 } 00109 } 00110 else { 00111 tr=xcmplx(1,1); 00112 ti=xcmplx(2,1); 00113 xcmplx(1,1)=0.5*(tr+ti); 00114 xcmplx(2,1)=0.5*(tr-ti); 00115 for (i=1; i<=n2; i++) { 00116 i1=i+1; 00117 i2=n-i+1; 00118 tr1=xcmplx(1,i1); 00119 tr2=xcmplx(1,i2); 00120 ti1=xcmplx(2,i1); 00121 ti2=xcmplx(2,i2); 00122 t=(cc*c-ss*s)+c; 00123 s=(cc*s+ss*c)+s; 00124 c=t; 00125 xcmplx(1,i1)=0.5*((tr1+tr2)-(tr1-tr2)*s-(ti1+ti2)*c); 00126 xcmplx(1,i2)=0.5*((tr1+tr2)+(tr1-tr2)*s+(ti1+ti2)*c); 00127 xcmplx(2,i1)=0.5*((ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00128 xcmplx(2,i2)=0.5*(-(ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00129 } 00130 fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,-2); 00131 } 00132 }
void frngs | ( | float * | circ, | |
int * | numr, | |||
int | nring | |||
) |
void normas | ( | float * | sqimg, | |
int | nsam, | |||
int | ns1, | |||
int | ns2, | |||
int | nr1, | |||
int | nr2, | |||
int | ir1, | |||
int | ir2 | |||
) |
Definition at line 505 of file spidutil.cpp.
Referenced by apmd().
00507 { 00508 /* 00509 c purpose: normalizes ring data. covered area is: ir1....ir2 * 00510 c * 00511 c parameters: * 00512 c 00513 c note : i think this is for parallel use only, because normass 00514 c is quicker for non_parallel use!! al sept 01 00515 c * 00516 */ 00517 // dimension sqimg(ns1:ns2,nr1:nr2) 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 // multiplication is faster 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 }
void normass | ( | float * | sqimg, | |
int | nsam, | |||
int | ns1, | |||
int | ns2, | |||
int | nr1, | |||
int | nr2, | |||
int | ir1, | |||
int | ir2 | |||
) |
Definition at line 399 of file spidutil.cpp.
Referenced by apring1(), and aprq2d().
00401 { 00402 /* 00403 c serially normalizes x by variance 00404 c 00405 c note : for parallel use normas instead al sept 01 00406 c difficult to add error flag due to use inside 00407 c parrallel region aug 05 al 00408 */ 00409 // this is how sqimg is declared in spider 00410 // dimension sqimg(ns1:ns2,nr1:nr2) 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 // array operation on x 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 // trap for blank image area 00451 // array operation on x 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 }
void numrinit | ( | int | mr, | |
int | nr, | |||
int | iskip, | |||
int * | numr | |||
) |
void parabld | ( | double * | z33, | |
float * | xsh, | |||
float * | ysh, | |||
double * | peakv | |||
) |
Definition at line 1491 of file spidutil.cpp.
Referenced by aprq2d().
01492 { 01493 /* 01494 c parabld 9/25/81 : parabolic fit to 3 by 3 peak neighborhood 01495 c double precision version 01496 c 01497 c the formula for paraboloid to be fiited into the nine points is: 01498 c 01499 c f = c1 + c2*y + c3*y**2 + c4*x + c5*xy + c6*x**2 01500 c 01501 c the values of the coefficients c1 - c6 on the basis of the 01502 c nine points around the peak, as evaluated by altran: 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 // the peak coordinates of the paraboloid can now be evaluated as: 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 // limit interplation to +/- 1. range 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 }
void prb1d | ( | double * | b, | |
int | npoint, | |||
float * | pos | |||
) |
Definition at line 787 of file spidutil.cpp.
References b.
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 }
float quadri | ( | float | xx, | |
float | yy, | |||
int | nxdata, | |||
int | nydata, | |||
float * | fdata | |||
) |
Quadratic interpolation (2D).
Note: This routine starts counting from 1, not 0!
This routine uses six image points for interpolation:
f3 fc | | x f2-----f0----f1 | | f4 *
f0 - f4 are image values near the interpolated point X. f0 is the interior mesh point nearest x.
Coords:
[in] | x | x-coord value |
[in] | y | y-coord value |
nx | ||
ny | ||
[in] | image | Image object (pointer) |
Definition at line 186 of file spidutil.cpp.
References fdata, quadri(), x, and y.
00187 { 00188 /* 00189 c purpose: quadratic interpolation 00190 c 00191 c parameters: xx,yy treated as circularly closed. 00192 c fdata - image 1..nxdata, 1..nydata 00193 c 00194 c f3 fc f0, f1, f2, f3 are the values 00195 c + at the grid points. x is the 00196 c + x point at which the function 00197 c f2++++f0++++f1 is to be estimated. (it need 00198 c + not be in the first quadrant). 00199 c + fc - the outer corner point 00200 c f4 nearest x. 00201 c 00202 c f0 is the value of the fdata at 00203 c fdata(i,j), it is the interior mesh 00204 c point nearest x. 00205 c the coordinates of f0 are (x0,y0), 00206 c the coordinates of f1 are (xb,y0), 00207 c the coordinates of f2 are (xa,y0), 00208 c the coordinates of f3 are (x0,yb), 00209 c the coordinates of f4 are (x0,ya), 00210 c the coordinates of fc are (xc,yc), 00211 c 00212 c o hxa, hxb are the mesh spacings 00213 c + in the x-direction to the left 00214 c hyb and right of the center point. 00215 c + 00216 c ++hxa++o++hxb++o hyb, hya are the mesh spacings 00217 c + in the y-direction. 00218 c hya 00219 c + hxc equals either hxb or hxa 00220 c o depending on where the corner 00221 c point is located. 00222 c 00223 c construct the interpolant 00224 c f = f0 + c1*(x-x0) + 00225 c c2*(x-x0)*(x-x1) + 00226 c c3*(y-y0) + c4*(y-y0)*(y-y1) 00227 c + c5*(x-x0)*(y-y0) 00228 c 00229 c 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 // circular closure 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 // hxc & hyc are either 1 or -1 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 }
void ringwe | ( | float * | wr, | |
int * | numr, | |||
int | nring | |||
) |
Definition at line 357 of file spidutil.cpp.
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 }
int setnring | ( | int | mr, | |
int | nr, | |||
int | iskip | |||
) |
Definition at line 372 of file spidutil.cpp.
Referenced by apmd(), and apmq().
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 }
void win_resize | ( | float * | imgfrom, | |
float * | imgto, | |||
int | lsam, | |||
int | lrow, | |||
int | nsam, | |||
int | nrow, | |||
int | lr1, | |||
int | lr2, | |||
int | ls1, | |||
int | ls2 | |||
) |
Definition at line 860 of file spidutil.cpp.
References imgfrom, imgto, and window().
Referenced by apmd(), and apmq().
00862 { 00863 /* 00864 c adpated from SPIDER ap_getdat.f 00865 c purpose: read read windowed image date into array x for 'ap' ops. 00866 c 00867 c parameters: 00868 c lsam,lrow image dimensions (input) 00869 c nsam,nrow output image dimensions (input) 00870 c lr1,lr2,ls1,ls2 output image window (input) 00871 c imgfrom input image (input) 00872 c imgto output output (output) 00873 c 00874 */ 00875 int window; 00876 00877 int k3, k2, kt; 00878 00879 // real, dimension(lsam) :: bufin 00880 00881 window = 0; 00882 if (lr1 != 1 || ls1 != 1 || lr2 != lrow || ls2 != lsam) window = 1; 00883 00884 if (window) { 00885 // window from the whole image 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 // do a copy 00894 for (k2=1;k2<=lrow;k2++) 00895 for (k3=1;k3<=lsam;k3++) 00896 imgto(k3,k2) = imgfrom(k3,k2); 00897 } 00898 }