spidutil.cpp File Reference

#include "ctf.h"
#include "log.h"
#include "emdata.h"
#include "xydata.h"
#include "Assert.h"
#include "projector.h"
#include <sys/types.h>
#include <sys/times.h>
#include <time.h>
#include <sys/time.h>
#include <iostream>
#include "spidutil.h"
#include "spidfft.h"

Include dependency graph for spidutil.cpp:

Go to the source code of this file.

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).
int alprbs (int *numr, int nring, int *lcirc, char mode)
void ringwe (float *wr, int *numr, int nring)
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 frngs (float *circ, int *numr, int nring)
int applyws (float *circ, int lcirc, int *numr, float *wr, int nring)
void normas (float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2, int ir1, int ir2)
void alrq (float *xim, int nsam, int nrow, 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, float *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 refparams, EMData *expimg, APMDopt options, float *fangles)
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 *refdirs, float *expdir, float range, float *diref, float *ccrot, float *rangnew, float *xshsum, float *yshsum, int *nimalcg, int ckmirror, int limitrange)
float ang_n (float rkk, char mode, int maxrin)
void alrq_ms (float *xim, int nsam, int nrow, float cns2, float cnr2, int *numr, float *circ, int lcirc, int nring, char mode)
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)


Function Documentation

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  refparams,
EMData expimg,
APMDopt  options,
float *  fangles 
)

Definition at line 901 of file spidutil.cpp.

References alprbs(), alrq(), angrefs, apmaster_1(), aprings(), crosrng_ms(), dlist, fangles, frngs(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), imgcirc, imgstk, imgwindow, APMDopt::iskip, APMDopt::mode, APMDopt::mr, normas(), APMDopt::nr, numr, numrinit(), nx, ny, refcstk, setnring(), EMAN::Dict::size(), status, tmt, tot, totmin, totmir, and win_resize().

Referenced by main().

00903 {
00904     // align      expimg using refprj as the reference
00905     // mr:        first ring 
00906     // nr:        last ring 
00907     // iskip:     number of rings skipped between the first and the last
00908     // mode:      ?
00909     // refparams: contains the angles used to generate refprj
00910     // fangles:   output angles assigned to expimg.
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     // retrieve APMD options
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     // extract image data from the reference image EM object
00954     refstk = refprj->get_data();
00955     //  find number of reference-rings
00956     nring = setnring(mr, nr, iskip);
00957 
00958     numr = (int*) calloc(3*nring,sizeof(int));
00959 
00960     numrinit(mr, nr, iskip, numr);
00961     // calculates pointers for rings and stores them in numr
00962     status = alprbs(numr, nring, &lcirc, mode);
00963     if (status != 0) goto EXIT;
00964 
00965     // convert reference images to rings
00966     refcstk = (float*)calloc(lcirc*nref,sizeof(float));
00967     
00968     // apply ring weights to rings
00969     status = aprings(nref, nring, refstk, nsam, nrow, numr,
00970                      refcstk, lcirc, mode);
00971     if (status != 0) goto EXIT;
00972 
00973     // extract the image data from experimental EM image object
00974     imgstk = expimg->get_data();
00975 
00976     // set the window size of the exp image: 
00977     // nwsam,nwrow  is the new size of the image
00978     apmaster_1(mode, &divas, nr, &numth, nsam, nrow, &nwsam, &nwrow);
00979 
00980     // allocate work space
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     // calculate window paramters
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        // normalize the image 
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        // turn the image into rings
01024        alrq(imgwindow,nwsam,nwrow,numr,imgcirc,lcirc,nring,mode);
01025 
01026        // should we use frng instead??  frng is ||, but || level was moved up PAP
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        // printf("j = %d, %g %g %g\n", j, dlist(1,j), dlist(2,j), dlist(3,j));
01054     }
01055 
01056     // now turn dlist into output angles (spider VOMD)
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 }

int apmq ( EMData refprj,
Dict  refparams,
EMData expimg,
APMQopt  options,
float *  angles,
float *  shifts 
)

Definition at line 1638 of file spidutil.cpp.

References alprbs(), APMQopt::angexps, angexps, angles, angrefs, apmaster_1(), aprings(), aprq2d(), bfc, dgr_to_rad, dlist, expdir, expdirs, EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), imgstk, imgwindow, APMQopt::istep, min, APMQopt::mode, APMQopt::mr, APMQopt::nr, numr, numrinit(), nx, ny, quadpi, APMQopt::range, refdirs, setnring(), shifts, APMQopt::shrange, EMAN::Dict::size(), status, and win_resize().

Referenced by main().

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     // retrieve APMQ options
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     // reference angles are passed in from refparams
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     // extract image data from the reference image EM object
01722     refstk = refprj->get_data();
01723     //  find number of reference-rings
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     //  Calculate pointers for rings and store them in numr
01737     status = alprbs(numr, nring, &lcirc, mode);
01738     if (status != 0) goto EXIT;
01739 
01740     maxrin = numr(3,nring);
01741 
01742     // convert reference images to rings
01743     bfc = (float*)calloc(lcirc*nima,sizeof(float));
01744 
01745     // read reference images into reference rings (bfc) array 
01746     status = aprings(nima, nring, refstk, nsam, nrow, numr, bfc, lcirc, mode);
01747     if (status !=0 ) goto EXIT;
01748 
01749     // extract the image data from experimental EM image object
01750     imgstk = expimg->get_data();
01751 
01752     // find divas, numth, nsam, & nrow
01753     apmaster_1(mode,&divas,nr,&numth,nsam,nrow,&nwsam,&nwrow);
01754 
01755     // *****STRANGELY, APMQ DOES NOT USE THE WINDOWED IMAGE*****
01756     //  Because it is not necessary - in the old APMD code that was done
01757     //     to save space  PAP  
01758     // allocate work space
01759     nwsam = nsam;
01760     nwrow = nrow; 
01761     imgwindow = (float*)calloc(nwsam*nwrow,sizeof(float));
01762 
01763     // calculate dimensions for normas
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        // retrieve refangles for restricted angular search
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        // convert ref. angles to unitary directional vectors (refdirs).
01784        // call ap_getsata(angrefs,refdirs,3,nima,irtflg)
01785 
01786 
01787        // read previously determined exp. angles into angexps
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        // convert exp. angles to unitary directional vectors(expdirs).
01803        // call ap_getsata(angexp,expdirs,7,nidi,irtflg)
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     // loop over all sets of experimental (sample) images
01814     for (it=1;it<=nidi;it++) {
01815        // APMQ DOES NOT WINDOW, IT JUST COPIES imgstk to imgwindow
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); // otherwise expdir is a dummy
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 //debug
01831 //       printf("dlist(2,it) = %11.3e\n", dlist(2,it));
01832 //       printf("dlist(3,it) = %11.3e\n", dlist(3,it));
01833 //       printf("dlist(4,it) = %11.3e\n", dlist(4,it));
01834 //       printf("dlist(5,it) = %11.3e\n", dlist(5,it));
01835 //       printf("dlist(6,it) = %11.3e\n", dlist(6,it));
01836 /*
01837 c      output (in dlist position is increased by 1, no.1 is the key).
01838 c      2 - number of the most similar reference projection.
01839 c      3 - not-normalized correlation coefficient.
01840 c      4 - psi angle. (in=plane rotation)
01841 c      5 - sx shift
01842 c      6 - sy shift
01843 c      7 - input image number.
01844 */
01845        //  dlist(2,it) is list number of most similar ref. image 
01846        //  (<0 if mirrored, 0 if none )
01847        iref = (int)dlist(2,it);
01848        if (iref < 0) {
01849           //  mirrored reference image
01850           imgref = -iref;
01851 
01852           //      ireft is for refdirs index
01853           ireft     = -iref;
01854           mirrornew = 1;
01855        } 
01856        else if (iref == 0) {
01857           // no nearby reference image
01858           imgref = 0;
01859           // ireft is for refdirs index
01860           ireft  = 1;
01861           mirrornew = 0;
01862        }
01863        else {
01864           imgref = iref;
01865           //  ireft is for refdirs index
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        // imgref is number of most similar ref. image 
01877        if (imgref <= 0) {
01878           // no reference image selected
01879           ccrot = -1.0;
01880           peakv = 0.0;
01881        }
01882 //  THE WHOLE NEXT SECTION IS SUSPICIOUS!  I WOULD REMOVE IT! PAP
01883        // set new projection angles
01884        angles(1,it) = 0.0; // default value
01885        angles(2,it) = 0.0; // 
01886        angles(2,it) = 0.0; //
01887 
01888        if (imgref > 0) {
01889           // use ref. angles as new projection angles
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              // ref. projection must be mirrored
01896              angles(1,it) = -angles(1,it);
01897              angles(2,it) = 180+angles(2,it);
01898           }
01899        }
01900        else if (ngotpar >= 3) {
01901           // keep old exp. proj. angles 
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           // use old inplane rot. & shift  
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        // combine rot. & shift with previous transformation
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        // list angles in range 0...360
01936        while ( rangout < 0.0) {
01937           rangout = rangout + 360.0;
01938        }
01939        while(rangout >= 360.0) {
01940           rangout = rangout - 360.0;
01941        } 
01942 
01943        // set flag for no angdif determined
01944        angdif = -1.0;
01945 
01946        if (imgref <= 0) {
01947           //  no relevant ref. image found
01948           angdif = 0.0;
01949        }
01950        else if (ngotpar >= 3) {
01951           //  can find angdif
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     } // endfor it
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 }

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 *  refdirs,
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,
float *  tot,
double *  qm,
double *  tmt 
)

Definition at line 637 of file spidutil.cpp.

References circ1, circ2, fftr_d(), numr, prb1d(), q, status, t, and t7.

Referenced by apmd(), and aprq2d().

00640 {
00641 /*
00642 c
00643 c  checks both straight & mirrored positions
00644 c
00645 c  input - fourier transforms of rings!!
00646 c  circ1 already multiplied by weights!
00647 c
00648 c  notes: aug 04 attempted speedup using 
00649 c       premultiply  arrays ie( circ12 = circ1 * circ2) much slower
00650 c       various  other attempts  failed to yield improvement
00651 c       this is a very important compute demand in alignmen & refine.
00652 c       optional limit on angular search should be added.
00653 */
00654 
00655    // dimension         circ1(lcirc),circ2(lcirc)
00656 
00657    // t(maxrin+2), q(maxrin+2), t7(-3:3)
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    //  c - straight  = circ1 * conjg(circ2)
00673    //  zero q array
00674   
00675    q = (double*)calloc(maxrin+2,sizeof(double));  
00676    if (!q) {
00677       status = -1;
00678       return status;
00679    }
00680 
00681    //   t - mirrored  = conjg(circ1) * conjg(circ2)
00682    //   zero t array
00683    t = (double*)calloc(maxrin+2,sizeof(double));
00684    if (!t) {
00685       status = -1;
00686       return status;
00687    } 
00688 
00689    //   premultiply  arrays ie( circ12 = circ1 * circ2) much slower
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      // some sort of error (probably compiler on mp on altix)
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   // this appears to interpolate? al
00754   prb1d(t7,7,&pos);
00755   *tot = (float)(jtot)+pos;
00756 
00757   // mirrored
00758   fftr_d(t,ip);
00759 
00760   // find angle
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   // find angle
00770   for (k=-3;k<=3;k++) {
00771     j       = ((jtot+k+maxrin-1)%maxrin) + 1;
00772     t7(k+4) = t(j);
00773   }
00774 
00775   // this appears to interpolate? al
00776 
00777   prb1d(t7,7,&pos);
00778   *tmt = float(jtot) + pos;
00779 
00780   free(t);
00781   free(q);
00782 
00783   return status;
00784 }

void frngs ( float *  circ,
int *  numr,
int  nring 
)

Definition at line 463 of file spidutil.cpp.

References circ, fftr_q(), and numr.

Referenced by apmd(), apring1(), and aprq2d().

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 }

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.

References sqimg, and sqrt().

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.

References sqimg, and sqrt().

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 
)

Definition at line 385 of file spidutil.cpp.

References numr.

Referenced by apmd(), and apmq().

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 }

void parabld ( double *  z33,
float *  xsh,
float *  ysh,
double *  peakv 
)

Definition at line 1491 of file spidutil.cpp.

References max, min, and z33.

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  x,
float  y,
int  nx,
int  ny,
float *  image 
)

Quadratic interpolation (2D).

Note: This routine starts counting from 1, not 0!

This routine uses six image points for interpolation:

See also:
M. Abramowitz & I.E. Stegun, Handbook of Mathematical Functions (Dover, New York, 1964), Sec. 25.2.67. http://www.math.sfu.ca/~cbm/aands/page_882.htm

http://www.cl.cam.ac.uk/users/nad/pubs/quad.pdf

                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:

Mesh spacings: Interpolant: f = f0 + c1*(x-x0) + c2*(x-x0)*(x-x1) + c3*(y-y0) + c4*(y-y0)*(y-y1) + c5*(x-x0)*(y-y0)

Parameters:
[in] x x-coord value
[in] y y-coord value
nx 
ny 
[in] image Image object (pointer)
Returns:
Interpolated value

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.

References numr, and wr.

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 }


Generated on Tue Jun 11 12:41:54 2013 for EMAN2 by  doxygen 1.4.7