EMAN2
Classes | Defines | Functions
spidutil.h File Reference
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 Documentation

#define angexps (   i,
 
)    angexps [((j)-1)*3 + (i)-1]

Definition at line 35 of file spidutil.h.

Referenced by apmq().

#define angles (   i,
 
)    angles [((j)-1)*3 + (i)-1]

Definition at line 32 of file spidutil.h.

#define angrefs (   i,
 
)    angrefs [((j)-1)*3 + (i)-1]

Definition at line 34 of file spidutil.h.

Referenced by apmd(), and apmq().

#define b (   i)    b [(i)-1]

Definition at line 56 of file spidutil.h.

#define bfc (   i,
 
)    bfc[((j)-1)*lcirc + (i) - 1]

Definition at line 73 of file spidutil.h.

Referenced by apmq(), and aprq2d().

#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,
 
)    dlist[((j)-1)*ldd + (i)-1]

Definition at line 67 of file spidutil.h.

Referenced by apmd(), and apmq().

#define expdir (   i)    expdir[(i)-1]

Definition at line 68 of file spidutil.h.

Referenced by apmq(), and aprq2d().

#define expdirs (   i,
 
)    expdirs[((j)-1)*3 + (i)-1]

Definition at line 69 of file spidutil.h.

Referenced by apmq().

#define fangles (   i,
 
)    fangles [((j)-1)*3 + (i)-1]

Definition at line 37 of file spidutil.h.

Referenced by apmd().

#define fdata (   i,
 
)    fdata [((j)-1)*nxdata + (i)-1]

Definition at line 43 of file spidutil.h.

#define fit (   i,
 
)    fit[((j)+istep)*(2*istep+1) + (i) + istep]
#define fitp (   i,
 
)    fitp[ ((j)+1)*3 + (i) + 1]

Definition at line 74 of file spidutil.h.

Referenced by aprq2d().

#define imgcirc (   i)    imgcirc[(i)-1]

Definition at line 39 of file spidutil.h.

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

#define imgfrom (   i,
 
)    imgfrom[((j)-1)*lsam + (i)-1]

Definition at line 58 of file spidutil.h.

Referenced by win_resize().

#define imgstk (   i,
  j,
 
)    imgstk[((k)-1)*nsam*nrow + ((j)-1)*nsam + (i)-1]

Definition at line 60 of file spidutil.h.

Referenced by apmd(), apmq(), and aprings().

#define imgto (   i,
 
)    imgto [((j)-1)*nsam + (i)-1]

Definition at line 59 of file spidutil.h.

Referenced by win_resize().

#define imgwindow (   i,
 
)    imgwindow [((j)-1)*nwsam + (i)-1]

Definition at line 62 of file spidutil.h.

Referenced by apmd(), and apmq().

#define lcg (   i)    lcg[(i)-1]

Definition at line 72 of file spidutil.h.

Referenced by aprq2d().

#define max (   a,
 
)    ((a) <= (b) ? (b) : (a))

Definition at line 46 of file spidutil.h.

#define min (   a,
 
)    ((a) >= (b) ? (b) : (a))

Definition at line 45 of file spidutil.h.

#define min0 (   a,
 
)    ((a) >= (b) ? (b) : (a))

Definition at line 44 of file spidutil.h.

Referenced by alprbs().

#define newangles (   i,
 
)    newangles [((j)-1)*3 + (i)-1]

Definition at line 33 of file spidutil.h.

Referenced by main().

#define numr (   i,
 
)    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,
 
)    refcstk[((j)-1)*lcirc + (i) - 1]

Definition at line 61 of file spidutil.h.

Referenced by apmd(), and aprings().

#define refdir (   i)    refdir[(i)-1]

Definition at line 71 of file spidutil.h.

#define refdirs (   i,
 
)    refdirs[((j)-1)*3 + (i)-1]

Definition at line 70 of file spidutil.h.

Referenced by apmq(), and aprq2d().

#define rotmp (   i,
 
)    rotmp[((j)+istep)*(2*istep+1) + (i) + istep]

Definition at line 76 of file spidutil.h.

Referenced by aprq2d().

#define shifts (   i,
 
)    shifts[((j)-1)*2 + (i)-1]
#define sqimg (   i,
 
)    sqimg [((j)-1)*nsam + (i)-1]

Definition at line 41 of file spidutil.h.

Referenced by apring1(), normas(), and normass().

#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]
#define tot (   i)    tot[(i)-1]
#define totmin (   i)    totmin[(i)-1]

Definition at line 63 of file spidutil.h.

Referenced by apmd().

#define totmir (   i)    totmir[(i)-1]

Definition at line 64 of file spidutil.h.

Referenced by apmd().

#define wr (   i)    wr [(i)-1]

Definition at line 40 of file spidutil.h.

Referenced by ali3d_d(), applyws(), apring1(), aprings(), ringwe(), trplot_(), and vrplot_().

#define xcmplx (   i,
 
)    xcmplx [((j)-1)*2 + (i)-1]

Definition at line 48 of file spidutil.h.

#define xim (   i,
 
)    xim [((j)-1)*nsam + (i)-1]

Definition at line 42 of file spidutil.h.

#define z33 (   i,
 
)    z33[((j)-1)*3 + (i)-1]

Definition at line 77 of file spidutil.h.

Referenced by parabld().


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().

{
/*
c  purpose: appears to circular rings, postitioned
c           in a linear array that holds rings concatenated together.
c           output is dependent on number of rings 
c                                                                      *
c  parameters:   numr(1,i) - ring number                      (sent)
c                numr(2,i) - beginning in circ                (ret.)
c                numr(3,i) - length in circ                   (ret.)
c                nring                                        (sent)
c                lcirc - total length of circ.                (ret.)
c
c image_processing_routine
c
*/
    int i, jp, ip;
    double dpi;
    int status = 0; 
    // hardwire for right now
    int MAXFFT = 32768;

    dpi = pi;
    if (mode == 'f' || mode == 'F') dpi=2*pi;

    *lcirc = 0;
    for (i=1;i<=nring;i++) {
       jp = (int)(dpi*numr(1,i));
       // original fortran code ip = 2**log2(jp), log2(jp) rounds up. 
       ip = (int)( pow(2,ceil(log2(jp))) );
       if (i < nring && jp > ip+ip/2)  ip=min0(MAXFFT,2*ip);

       //  last ring should be oversampled to allow higher accuracy
       //  of peak location (?).
       if (i == nring && jp > ip+ip/5) ip=min0(MAXFFT,2*ip);
       numr(3,i) = ip;
       if (i == 1) {
          numr(2,1) = 1;
       }
       else {
          numr(2,i) = numr(2,i-1)+numr(3,i-1);
       }
       *lcirc = *lcirc + ip;
    }
    return status;
}
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().

{
/* 
c                                                                     
c  purpose:                                                          
c                                                                   
c  parameters: convert to polar coordinates
c                                                                  
*/
   //  dimension         xim(nsam,nrow),circ(lcirc)
   //  integer           numr(3,nring)

   double dfi, dpi;
   int    ns2, nr2, i, inr, l, nsim, kcirc, lt, j;
   float  yq, xold, yold, fi, x, y;

   ns2 = nsam/2+1;
   nr2 = nrow/2+1;
   dpi = 2.0*atan(1.0);

//#pragma omp   parallel do private(i,j,inr,yq,l,lt,nsim,dfi,kcirc,
//#pragma omp&  xold,yold,fi,x,y)
   for (i=1;i<=nring;i++) {
     // radius of the ring
     inr = numr(1,i);
     yq  = inr;
     l   = numr(3,i);
     if (mode == 'h' || mode == 'H') {
        lt = l/2;
     }
     else if (mode == 'f' || mode == 'F' ) {
        lt = l/4;
     }

     nsim           = lt-1;
     dfi            = dpi/(nsim+1);
     kcirc          = numr(2,i);
     xold           = 0.0;
     yold           = inr;
     circ(kcirc)    = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
     xold           = inr;
     yold           = 0.0;
     circ(lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);

     if (mode == 'f' || mode == 'F') {
        xold              = 0.0;
        yold              = -inr;
        circ(lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
        xold              = -inr;
        yold              = 0.0;
        circ(lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
     }

     for (j=1;j<=nsim;j++) {
        fi               = dfi*j;
        x                = sin(fi)*yq;
        y                = cos(fi)*yq;
        xold             = x;
        yold             = y;
        circ(j+kcirc)    = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
        xold             =  y;
        yold             = -x;
        circ(j+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);

        if (mode == 'f' || mode == 'F')  {
           xold                = -x;
           yold                = -y;
           circ(j+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
           xold                = -y;
           yold                =  x;
           circ(j+lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
        };
     }
   }
//#pragma omp   end parallel do 
}
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().

{
   double dpi, dfi;
   int    it, jt, inr, l, nsim, kcirc, lt;
   float  yq, xold, yold, fi, x, y;

   //     cns2 and cnr2 are predefined centers
   //     no need to set to zero, all elements are defined

   dpi = 2*atan(1.0);
   for (it=1;it<=nring;it++) {
      // radius of the ring
      inr = numr(1,it);
      yq  = inr;

      l = numr(3,it);
      if ( mode == 'h' || mode == 'H' ) { 
         lt = l / 2;
      }
      else if ( mode == 'f' || mode == 'F' ) {
         lt = l / 4;
      } 

      nsim  = lt - 1;
      dfi   = dpi / (nsim+1);
      kcirc = numr(2,it);
      xold  = 0.0;
      yold  = inr;

      circ(kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

      xold  = inr;
      yold  = 0.0;
      circ(lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

      if ( mode == 'f' || mode == 'F' ) {
         xold = 0.0;
         yold = -inr;
         circ(lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

         xold = -inr;
         yold = 0.0;
         circ(lt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
      }
      
      for (jt=1;jt<=nsim;jt++) {
         fi   = dfi * jt;
         x    = sin(fi) * yq;
         y    = cos(fi) * yq;

         xold = x;
         yold = y;
         circ(jt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

         xold = y;
         yold = -x;
         circ(jt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

         if ( mode == 'f' || mode == 'F' ) {
            xold = -x;
            yold = -y;
            circ(jt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);

            xold = -y;
            yold = x;
            circ(jt+lt+lt+lt+kcirc) = quadri(xold+cns2,yold+cnr2,nsam,nrow,xim);
         }
      } // end for jt
   } //end for it
}
float ang_n ( float  rkk,
char  mode,
int  maxrin 
)

Definition at line 1406 of file spidutil.cpp.

{
    float ang; 

    if (mode == 'H' || mode == 'h') {
        ang = fmod(((rkk-1.0) / maxrin+1.0)*180.0, 180.0);
    }
    else if ( mode == 'F' || mode == 'f') {
        ang = fmod(((rkk-1.0) / maxrin+1.0)*360.0, 360.0);
    }
    return ang;
}
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().

{
/*
c parameters:
c       mode                degree mode                       (input)
c       divas               degrees                           (output)
c       numth               degrees                           (output)
c       lsam                orig size                         (input)
c       lrow                orig size                         (input)
c       nsam                new size                          (output)
c       nrow                new size                          (output)
*/
   int nra;

   if ( mode == 'h') {
      *divas = 180.0;
   }
   else {
      *divas = 360.0;
   }

   *numth = 1;
#ifdef sp_mp
//       find number of omp threads
//        call getthreads(numth)
#endif

   //  calculation of actual dimension of an image to be interpolated
   //  2*(no.of rings)+(0'th element)+2*(margin of 1)

   nra  = ((lsam-1)/2)*2+1;
   if ( ((lrow-1)/2)*2+1 < nra ) nra = ((lrow-1)/2)*2+1;
   if ( 2*nr+3 < nra ) nra = 2*nr+3;

   //  returns circular reduced nsam, nrow
   *nsam = nra;
   *nrow = nra;
}
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().

{
   int   maxrin, numr3i, numr2i;
   float w;
   int   i, j, jc;
   int   status = 0;

   maxrin = numr(3,nring);
 
   for (i=1;i<=nring;i++) {
      numr3i=numr(3,i);
      numr2i=numr(2,i);
      w=wr(i);
      circ(numr2i)=circ(numr2i)*w;
      if (numr3i == maxrin) {
         circ(numr2i+1)=circ(numr2i+1)*w;
      }
      else { 
         circ(numr2i+1)=circ(numr2i+1)*0.5*w;
      }
      for (j=3;j<=numr3i;j++) {
         jc=j+numr2i-1;
         circ(jc)=circ(jc)*w;
      }
   } 
   if (jc >= lcirc) status = -1;
   return status; 
}
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().

{
   int    status = 0;
   int    nsb, nse, nrb, nre, maxrin, ltf, lt, lt2, lt3, ns2, nr2;
   int    inr, kcirc, jt, i, j;
   float  fnr2, fns2, yq, fi, dpi, x, y;
   double dfi;

   dpi    = 2.0*atan(1.0); 

   // calculate dimensions for normass
   nsb = -nsam/2;
   nse = -nsb-1+(nsam%2);
   nrb = -nrow/2;
   nre = -nrb-1+(nrow%2);

   //  normalize under the mask,  tried doing this on the
   //  polar rings but it gives different answers. al

   normass(sqimg,nsam,nsb,nse,nrb,nre,numr(1,1),numr(1,nring));

   maxrin = numr(3,nring);
   ns2    = nsam / 2 + 1;
   nr2    = nrow / 2 + 1;
   fnr2   = nr2;
   fns2   = ns2;

   // convert window from image into polar coordinates
        
   if (mode == 'f' || mode == 'F') {
      ltf = 4;
      for (i=1;i<=nring;i++) {
         //  radius of the ring
         inr             = numr(1,i);
         yq              = inr;
         lt              = numr(3,i) / ltf;
         lt2             = lt  + lt;
         lt3             = lt2 + lt;
         dfi             = dpi / lt;
         kcirc           = numr(2,i);
        
         imgcirc(kcirc)     = sqimg(ns2,     nr2+inr);
         imgcirc(lt+kcirc)  = sqimg(ns2+inr, nr2);
         imgcirc(lt2+kcirc) = sqimg(ns2,     nr2-inr);
         imgcirc(lt3+kcirc) = sqimg(ns2-inr, nr2);

         for (j=1;j<=lt - 1;j++) {
            fi = dfi     * j;
            x  = sin(fi) * yq;
            y  = cos(fi) * yq;
            jt = j + kcirc;

            imgcirc(jt)     = quadri(fns2+x,fnr2+y,nsam,nrow,sqimg);
            imgcirc(jt+lt)  = quadri(fns2+y,fnr2-x,nsam,nrow,sqimg);
            imgcirc(jt+lt2) = quadri(fns2-x,fnr2-y,nsam,nrow,sqimg);
            imgcirc(jt+lt3) = quadri(fns2-y,fnr2+x,nsam,nrow,sqimg);
         }
      }
   }
   else if (mode == 'h' || mode == 'H') {
      ltf = 2;
      for (i=1; i<=nring;i++) {
         // radius of the ring
         inr            = numr(1,i);
         yq             = inr;
         lt             = numr(3,i) / ltf;
         dfi            = dpi / lt;
         kcirc          = numr(2,i);

         imgcirc(kcirc)    = sqimg(ns2,     nr2+inr);
         imgcirc(lt+kcirc) = sqimg(ns2+inr, nr2);
 
         for (j=1;j<=lt - 1;j++) {
            fi          = dfi * j;
            x           = sin(fi) * yq;
            y           = cos(fi) * yq;
            jt          = j + kcirc;

            imgcirc(jt)    = quadri(fns2+x,fnr2+y,nsam,nrow,sqimg);
            imgcirc(jt+lt) = quadri(fns2+y,fnr2-x,nsam,nrow,sqimg);
         }
      }
   }

   //  fourier of circ 
   frngs(imgcirc,numr,nring);

   //    weight circ  using wr
   if (wr(1) > 0.0) {
      status = applyws(imgcirc,lcirc,numr,wr,nring);
   }

   return status;
}
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().

{
    int j, status;
    float *wr;

    status = 0;
    wr = (float*) calloc(nring, sizeof(float));
    if (!wr) {
        fprintf(stderr, "aprings: failed to allocate wr!\n");
        status = -1;
        goto EXIT;
    }


    // get wr weights
    ringwe(wr, numr, nring);
    if ( mode == 'H' || mode == 'h' )
        for (j=1;j<=nring;j++) wr(j) = wr(j)/2.0;
    for (j = 1; j<=nimg; j++) {
       apring1(&imgstk(1,1,j), nsam, nrow, &refcstk(1,j),
               lcirc, nring, mode, numr, wr);
    }

 EXIT:
    if (wr) free(wr);
    return status;
}
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().

{
   float  *imgcirc;
   int    *lcg; 

   double ccrotd,peak,tota,tmta,tmt;
   int    mirrored;

   double quadpi=3.14159265358979323846;
   double dgr_to_rad = quadpi/180.0;

   int    imi, iend, mwant, jtma, itma;
   float  dt, dtabs, rangnewt;
   int    jt, it, irr, ir, ibe, isx, isy, status, idis;
   float  cnr2, cns2, co, so, afit, sx, sy, tot;

   double fitp[9], *fit;
   float  *rotmp;
   int    fitsize;

   status = 0;
 
   imgcirc = (float*)calloc(lcirc,sizeof(float));
   if (imgcirc == NULL) {
       fprintf(stderr,"aprq2d: failed to allocate imgcirc\n");
       status = -1;
       goto EXIT;
   }

   fitsize = (2*istep+1)*(2*istep+1);
   if ( istep >= 1) {
       fit   = (double*)calloc(fitsize,sizeof(double));
       rotmp = (float*) calloc(fitsize,sizeof(float));
   }
   else {
       status = -2;
       goto EXIT;
   }

   peak = 0.0;
   iend  = nima;

   if (limitrange) {
      // restricted range search
      lcg  = (int*) calloc(nima, sizeof(int));
      if (!lcg) {
         mwant = nima;
         status = -1;
         fprintf(stderr,"lcg: %d\n", mwant);
         fprintf(stderr, "range: %g, nima: %d\n", range, nima);
         goto EXIT;
      }

      *nimalcg = 0;
      for (imi=1; imi<=nima; imi++) {
         // dt near 1.0 = not-mirrored, dt near -1.0 = mirrored
         dt    = (expdir(1) * refdirs(1,imi) + 
                  expdir(2) * refdirs(2,imi) +
                  expdir(3) * refdirs(3,imi));
         dtabs = fabs(dt);

         if (dtabs >= range) {
            // mirored or non-mirrored is within range
            *nimalcg++;
            lcg(*nimalcg) = imi;
            if (dt < 0) lcg(*nimalcg) = -imi;
         }
      }

      if (*nimalcg <= 0) {
         // there is no reference within search range
         *xshsum  = 0;
         *yshsum  = 0;
         *diref   = 0;
         *rangnew = 0;
         *ccrot   = -1.0;
         goto EXIT; 
      }
      iend = *nimalcg;
      // end of restricted range search
   }

   ccrotd = -1.0e23;

   // go through centers for shift alignment
   for (jt=-ishrange;jt<=ishrange;jt=jt+istep) {
      cnr2 = nrow/2+1+jt;
      for (it=-ishrange;it<=ishrange;it=it+istep) {
         cns2 = nsam/2+1+it;

         // normalize under the mask
         //'normalize' image values over variance range
         normass(sqimg,nsam,nsb-it,nse-it,nrb-jt,nre-jt,numr(1,1),
                 numr(1,nring));

         // interpolation into polar coordinates
         // creates imgcirc (exp. image circles) for this position

         alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,nring,mode);

         // creates fourier of: a_circ
         frngs(imgcirc,numr,nring);

         // compare exp. image with all reference images
         for (irr=1;irr<=iend;irr++) {
            ir = irr;
            if (limitrange) ir = abs(lcg(irr));
              
            if (ckmirror) {
               if (limitrange) {
                  mirrored = 0;
                  if (lcg(irr) < 0) mirrored = 1;
                  // check either mirrored or non-mirrored position 
                  crosrng_e(&bfc(1,ir),imgcirc,lcirc,nring,
                            maxrin,numr,&tota,&tot,mirrored);
               }
               else {
                  // check both non-mirrored & mirrored positions 
                  status = crosrng_ms(&bfc(1,ir),imgcirc,lcirc,nring,
                                      maxrin,numr,&tota,&tot,&tmta,&tmt);
               }
            }
            else {
               // do not check mirrored position
                mirrored = 0;
                crosrng_e(&bfc(1,ir),imgcirc,lcirc,nring,
                          maxrin,numr,&tota,&tot,mirrored);
            }

            if (tota >= ccrotd) {
               // good match with tota (mirrored or not)  position 
               ccrotd  = tota;
               ibe     = ir;
               isx     = it;
               isy     = jt;
               *rangnew = ang_n(tot,mode,maxrin);
               idis    = ir;
               if (limitrange && lcg(irr) < 0) idis = -ir;
            }

            if (ckmirror && !limitrange) {
               // have to compare with mirrored position 
               if (tmta >= ccrotd) {
                  // good match with mirrored position 
                  ccrotd  = tmta;
                  ibe     = ir;
                  isx     = it;
                  isy     = jt;
                  *rangnew =  ang_n(tmt,mode,maxrin);
                  idis    = -ir;
               }
            }
         } // endfor irr 
      } // endfor it
   } //  endfor jt

   // try to interpolate
   *ccrot = ccrotd;
   sx     = isx;
   sy     = isy;
   *diref = idis;

   // do not interpolate for point on the edge
   if ( abs(isx) != ishrange && abs(isy) != ishrange) {
      // have to find neighbouring values
      fit(0,0)   = ccrotd;
      rotmp(0,0) = *rangnew;

      for (jt=-istep;jt<=istep;jt++) {
         for (it=-istep;it<=istep;it++) {
            if (it !=0 || jt != 0) {
               cnr2 = nrow/2+1+jt+isy;
               cns2 = nsam/2+1+it+isx;

               normass(sqimg, nsam, nsb-(it+isx),nse-(it+isx),
                       nrb-(jt+isy),nre-(jt+isy), numr(1,1), numr(1,nring));

               alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,
                       nring,mode);

               frngs(imgcirc,numr,nring);

               //  if (idis .lt. 0)  check mirrored only
               mirrored = 0;
               if (idis < 0) mirrored = 1;
               crosrng_e(&bfc(1,ibe),imgcirc,lcirc,nring,
                         maxrin,numr,&fit(it,jt),&rotmp(it,jt),
                         mirrored);
               rotmp(it,jt) = ang_n(rotmp(it,jt),mode,maxrin);
            }
         } // endfor it
      } //endfor jt 

      //  find the maximum within +/-istep
      //  maximum cannot be on the edge, i.e., it,jt/=istep
      afit     = fit(0,0);
      jtma     = 0;
      itma     = 0;
      rangnewt = rotmp(0,0);
      if ( istep > 1) {
         for (jt=-istep+1;jt<=istep-1;jt++) {
            for (it=-istep+1;it<=istep-1;it++) {
               if (fit(it,jt) > afit) {
                  afit     = fit(it,jt);
                  rangnewt = rotmp(it,jt);
                  itma     = it;
                  jtma     = jt;
               }
            }
         } 
      }
      //  temp variable overcomes compiler bug on opt 64 pgi 6.0
      *rangnew = rangnewt;

      //  copy values around the peak.
      for (jt=-1;jt<=1;jt++) 
         for (it=-1;it<=1;it++)
            fitp(it,jt) = fit(itma+it,jtma+jt);

      //  update location of the peak
      ccrotd = afit;
      isx    = isx+itma;
      isy    = isy+jtma;
      parabld(fitp,&sx,&sy,&peak);

      //  check whether interpolation is ok.
      if (fabs(sx) < 1.0 && fabs(sy) < 1.0) {
         //  not on edge of 3x3 area
         sx   = sx+isx;
         sy   = sy+isy;
         cnr2 = nrow/2+1+sy;
         cns2 = nsam/2+1+sx;

         normass(sqimg,nsam,nsb-isx,nse-isx,nrb-isy,nre-isy,numr(1,1),
                 numr(1,nring));

         alrq_ms(sqimg,nsam,nrow,cns2,cnr2,numr,imgcirc,lcirc,nring,mode);

         frngs(imgcirc,numr,nring);

         mirrored = 0;
         if (idis < 0) mirrored = 1;

         crosrng_e(&bfc(1,ibe),imgcirc,lcirc,nring,
                   maxrin,numr,&ccrotd,rangnew,mirrored);

         *ccrot   = ccrotd;
         *rangnew = ang_n(*rangnew,mode,maxrin);
      } 
      else {
         //  not on edge of 3x3 area
         sx = isx;
         sy = isy;
      }
   }

   sx = -sx;
   sy = -sy;

   // now have to change order of shift & rotation.
   // in this program image is shifted first, rotated second.
   // in 'rt sq' it is rotation first, shift second.
   // this part corresponds to 'sa p'.
   co      =  cos((*rangnew) * dgr_to_rad);
   so      = -sin((*rangnew) * dgr_to_rad);
   *xshsum = sx*co - sy*so;
   *yshsum = sx*so + sy*co;

   free(fit);
   free(rotmp);
   free(imgcirc);
   if (limitrange) free(lcg);

EXIT:
   return status;
}
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().

{
/*
c checks single position, neg is flag for checking mirrored position
c
c  input - fourier transforms of rings!
c  first set is conjugated (mirrored) if neg
c  circ1 already multiplied by weights!
c       automatic arrays
        dimension         t(maxrin+2)
        double precision  q(maxrin+2)
        double precision  t7(-3:3)
*/
   float *t;
   double t7[7], *q;
   int    i, j, k, ip, jc, numr3i, numr2i, jtot;
   float  pos; 

   ip = maxrin;
   q = (double*)calloc(maxrin+2, sizeof(double));
   t = (float*)calloc(maxrin+2, sizeof(float));
     
   for (i=1;i<=nring;i++) {
      numr3i = numr(3,i);
      numr2i = numr(2,i);

      t(1) = (circ1(numr2i)) * circ2(numr2i);

      if (numr3i != maxrin) {
         // test .ne. first for speed on some compilers
         t(numr3i+1) = circ1(numr2i+1) * circ2(numr2i+1);
         t(2)        = 0.0;

         if (neg) {
            // first set is conjugated (mirrored)
            for (j=3;j<=numr3i;j=j+2) {
               jc = j+numr2i-1;
               t(j) =(circ1(jc))*circ2(jc)-(circ1(jc+1))*circ2(jc+1);
               t(j+1) = -(circ1(jc))*circ2(jc+1)-(circ1(jc+1))*circ2(jc);
            } 
         } 
         else {
            for (j=3;j<=numr3i;j=j+2) {
               jc = j+numr2i-1;
               t(j) = (circ1(jc))*circ2(jc) + (circ1(jc+1))*circ2(jc+1);
               t(j+1) = -(circ1(jc))*circ2(jc+1) + (circ1(jc+1))*circ2(jc);
            }
         } 
         for (j=1;j<=numr3i+1;j++) q(j) = q(j) + t(j);
      }
      else {
         t(2) = circ1(numr2i+1) * circ2(numr2i+1);
         if (neg) {
            // first set is conjugated (mirrored)
            for (j=3;j<=maxrin;j=j+2) {
               jc = j+numr2i-1;
               t(j) = (circ1(jc))*circ2(jc) - (circ1(jc+1))*circ2(jc+1);
               t(j+1) = -(circ1(jc))*circ2(jc+1) - (circ1(jc+1))*circ2(jc);
            }
         }
         else {
            for (j=3;j<=maxrin;j=j+2) {
               jc = j+numr2i-1;
               t(j) = (circ1(jc))*circ2(jc) + (circ1(jc+1))*circ2(jc+1);
               t(j+1) = -(circ1(jc))*circ2(jc+1) + (circ1(jc+1))*circ2(jc);
            } 
         }
         for (j = 1; j <= maxrin+2; j++) q(j) = q(j) + t(j);
      }
   }

   fftr_d(q,ip);

   *qn = -1.0e20;
   for (j=1;j<=maxrin;j++) {
      if (q(j) >= *qn) {
         *qn = q(j);
         jtot = j;
      }
   } 

   for (k=-3;k<=3;k++) {
      j = (jtot+k+maxrin-1)%maxrin + 1;
      t7(k+4) = q(j);
   }

   prb1d(t7,7,&pos);

   *tot = (float)jtot + pos;

   if (q) free(q);
   if (t) free(t);
}
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.

{
   //  dimension  br(1),bi(1)

   int b3,b4,b5,b6,b7,b56;
   int n, k, l, j, i, ix0, ix1; 
   float rni, tr1, ti1, tr2, ti2, cc, c, ss, s, t, x2, x3, x4, x5, sgn;
   float tab1[15]; 
   int status=0;

   tab1(1)=9.58737990959775e-5;
   tab1(2)=1.91747597310703e-4;
   tab1(3)=3.83495187571395e-4;
   tab1(4)=7.66990318742704e-4;
   tab1(5)=1.53398018628476e-3;
   tab1(6)=3.06795676296598e-3;
   tab1(7)=6.13588464915449e-3;
   tab1(8)=1.22715382857199e-2;
   tab1(9)=2.45412285229123e-2;
   tab1(10)=4.90676743274181e-2;
   tab1(11)=9.80171403295604e-2;
   tab1(12)=1.95090322016128e-1;
   tab1(13)=3.82683432365090e-1;
   tab1(14)=7.07106781186546e-1;
   tab1(15)=1.00000000000000;

   n=(int)pow(2,ln);

   k=abs(ks);
   l=16-ln;
   b3=n*k;
   b6=b3;
   b7=k;
   if( ks > 0 ) {
      sgn=1.0;
   } 
   else {
      sgn=-1.0;
      rni=1.0/(float)n;
      j=1;
      for (i=1; i<=n;i++) {
         br(j)=br(j)*rni;
         bi(j)=bi(j)*rni;
         j=j+k;
      }
   }
L12:
   b6=b6/2;
   b5=b6;
   b4=2*b6;
   b56=b5-b6;
L14:
   tr1=br(b5+1);
   ti1=bi(b5+1);
   tr2=br(b56+1);
   ti2=bi(b56+1);

   br(b5+1)=tr2-tr1;
   bi(b5+1)=ti2-ti1;
   br(b56+1)=tr1+tr2;
   bi(b56+1)=ti1+ti2;

   b5=b5+b4;
   b56=b5-b6;
   if (b5 <= b3)  goto  L14;
   if (b6 == b7)  goto  L20;

   b4=b7;
   cc=2.0*pow(tab1(l),2);
   c=1.0-cc;
   l=l+1;
   ss=sgn*tab1(l);
   s=ss;
L16: 
   b5=b6+b4;
   b4=2*b6;
   b56=b5-b6;
L18:
   tr1=br(b5+1);
   ti1=bi(b5+1);
   tr2=br(b56+1);
   ti2=bi(b56+1);
   br(b5+1)=c*(tr2-tr1)-s*(ti2-ti1);
   bi(b5+1)=s*(tr2-tr1)+c*(ti2-ti1);
   br(b56+1)=tr1+tr2;
   bi(b56+1)=ti1+ti2;

   b5=b5+b4;
   b56=b5-b6;
   if(b5 <= b3)  goto L18;
   b4=b5-b6;
   b5=b4-b3;
   c=-c;
   b4=b6-b5;
   if(b5 < b4)  goto  L16;
   b4=b4+b7;
   if(b4 >= b5) goto  L12;

   t=c-cc*c-ss*s;
   s=s+ss*c-cc*s;
   c=t;
   goto  L16;
L20:
   ix0=b3/2;
   b3=b3-b7;
   b4=0;
   b5=0;
   b6=ix0;
   ix1=0;
   if ( b6 == b7) goto EXIT;
L22:
   b4=b3-b4;
   b5=b3-b5;
   x2=br(b4+1);
   x3=br(b5+1);
   x4=bi(b4+1);
   x5=bi(b5+1);
   br(b4+1)=x3;
   br(b5+1)=x2;
   bi(b4+1)=x5;
   bi(b5+1)=x4;
   if (b6 < b4) goto  L22;
L24:
   b4=b4+b7;
   b5=b6+b5;
   x2=br(b4+1);
   x3=br(b5+1);
   x4=bi(b4+1);
   x5=bi(b5+1);
   br(b4+1)=x3;
   br(b5+1)=x2;
   bi(b4+1)=x5;
   bi(b5+1)=x4;
   ix0=b6;
L26:
   ix0=ix0/2;
   ix1=ix1-ix0;
   if(ix1 >= 0)  goto  L26;

   ix0=2*ix0;
   b4=b4+b7;
   ix1=ix1+ix0;
   b5=ix1;
   if (b5 >= b4)  goto  L22;
   if (b4 < b6)   goto  L24;
EXIT:
   status = 0; 
}
void fftr_q ( float *  xcmplx,
int  nv 
)

Definition at line 53 of file spidfft.cpp.

References abs, fftc_q(), t, tab1, and xcmplx.

{
   // dimension xcmplx(2,1); xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary

   float tab1[15];
   int nu, inv, nu1, n, isub, n2, i1, i2, i;
   float ss, cc, c, s, tr, ti, tr1, tr2, ti1, ti2, t;

   tab1(1)=9.58737990959775e-5;
   tab1(2)=1.91747597310703e-4;
   tab1(3)=3.83495187571395e-4;
   tab1(4)=7.66990318742704e-4;
   tab1(5)=1.53398018628476e-3;
   tab1(6)=3.06795676296598e-3;
   tab1(7)=6.13588464915449e-3;
   tab1(8)=1.22715382857199e-2;
   tab1(9)=2.45412285229123e-2;
   tab1(10)=4.90676743274181e-2;
   tab1(11)=9.80171403295604e-2;
   tab1(12)=1.95090322016128e-1;
   tab1(13)=3.82683432365090e-1;
   tab1(14)=7.07106781186546e-1;
   tab1(15)=1.00000000000000;

   nu=abs(nv);
   inv=nv/nu;
   nu1=nu-1;
   n=(int)pow(2,nu1);
   isub=16-nu1;

   ss=-tab1(isub);
   cc=-2.0*pow(tab1(isub-1),2);
   c=1.0;
   s=0.0;
   n2=n/2;
   if ( inv > 0) {
      fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,2);
      tr=xcmplx(1,1);
      ti=xcmplx(2,1);
      xcmplx(1,1)=tr+ti;
      xcmplx(2,1)=tr-ti;
      for (i=1;i<=n2;i++) {
         i1=i+1;
         i2=n-i+1;
         tr1=xcmplx(1,i1);
         tr2=xcmplx(1,i2);
         ti1=xcmplx(2,i1);
         ti2=xcmplx(2,i2);
         t=(cc*c-ss*s)+c;
         s=(cc*s+ss*c)+s;
         c=t;
         xcmplx(1,i1)=0.5*((tr1+tr2)+(ti1+ti2)*c-(tr1-tr2)*s);
         xcmplx(1,i2)=0.5*((tr1+tr2)-(ti1+ti2)*c+(tr1-tr2)*s);
         xcmplx(2,i1)=0.5*((ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c);
         xcmplx(2,i2)=0.5*(-(ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c);
     }
   }
   else {
     tr=xcmplx(1,1);
     ti=xcmplx(2,1);
     xcmplx(1,1)=0.5*(tr+ti);
     xcmplx(2,1)=0.5*(tr-ti);
     for (i=1; i<=n2; i++) {
        i1=i+1;
        i2=n-i+1;
        tr1=xcmplx(1,i1);
        tr2=xcmplx(1,i2);
        ti1=xcmplx(2,i1);
        ti2=xcmplx(2,i2);
        t=(cc*c-ss*s)+c;
        s=(cc*s+ss*c)+s;
        c=t;
        xcmplx(1,i1)=0.5*((tr1+tr2)-(tr1-tr2)*s-(ti1+ti2)*c);
        xcmplx(1,i2)=0.5*((tr1+tr2)+(tr1-tr2)*s+(ti1+ti2)*c);
        xcmplx(2,i1)=0.5*((ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s);
        xcmplx(2,i2)=0.5*(-(ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s);
     }
     fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,-2);
   }
}
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().

{
   int i, l; 
 
   for (i=1; i<=nring;i++) {
     l=(int)(log2(numr(3,i)));
     fftr_q(&circ(numr(2,i)),l);
   }
}
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().

{
/*
c  purpose:    normalizes ring data.  covered area is: ir1....ir2      *
c                                                                      *
c  parameters:                                                         *
c
c  note   :    i think this is for parallel use only, because normass
c              is quicker for non_parallel use!! al sept 01
c                                                                      *
*/
    //  dimension  sqimg(ns1:ns2,nr1:nr2)

    double     av,vr;
    int        i1sq, i2sq, n, i, j, ir, j2, irow, jcol;

    i1sq = ir1 * ir1;
    i2sq = ir2 * ir2;

    av   = 0.0;
    vr   = 0.0;
    n    = 0;

    for (j=nr1;j<=nr2;j++) {
       j2 = j*j;
       for (i=ns1;i<=ns2;i++) {
          ir = j2 + i*i;
          jcol = j-nr1+1;
          if (ir >= i1sq && ir <= i2sq) {
             n++;
             irow = i-ns1+1;
             av = av + sqimg(irow,jcol);
             vr = vr + sqimg(irow,jcol)*sqimg(irow,jcol);
          }
       }
    } 

    av = av / n;

    //   multiplication is faster
    vr = 1.0 / (sqrt((vr-n*av*av) / (n-1)));

    for (j=nr1; j<=nr2; j++) {
       jcol = j-nr1+1;
       for (i=ns1;i<=ns2;i++) {
          irow = i-ns1+1;
          sqimg(irow,jcol) = (sqimg(irow,jcol) - av ) * vr;
       } 
    }
}
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().

{
/*
c serially normalizes x by variance
c
c  note   :    for parallel use normas instead al sept 01
c              difficult to add error flag due to use inside
c              parrallel region aug 05 al
*/
   // this is how sqimg is declared in spider
   // dimension  sqimg(ns1:ns2,nr1:nr2)

   double   av,vr,vrinv,dtemp;
   int      i1sq, i2sq, n, ir, jsq, i, j, irow, jcol;

   i1sq = ir1 * ir1;
   i2sq = ir2 * ir2;
   av   = 0.0;
   vr   = 0.0;
   n    = 0;

   for (j=nr1; j<=nr2; j++) {
      jsq = j * j;
      for (i=ns1;i<=ns2;i++) {
         ir = jsq + i * i;
         if (ir >= i1sq && ir <= i2sq) {
            n  = n  + 1;
            irow = i-ns1+1;
            jcol = j-nr1+1; 
            av = av + sqimg(irow,jcol);
            vr = vr + sqimg(irow,jcol)*sqimg(irow,jcol);
         }
      }
   }

   av   = av / n;
   dtemp = (vr - n * av * av);
   if (dtemp > 0) {
      vr    = sqrt(dtemp / (n-1));
      vrinv = 1.0 / vr;

      //  array operation on x
      for ( j = nr1; j<=nr2;j++) 
         for (i = ns1; i <= ns2; i++) {
            irow = i - ns1 + 1; 
            jcol = j - nr1 + 1; 
            sqimg(irow,jcol) = (sqimg(irow,jcol) - av) * vrinv;
         }
   }
   else {
      // trap for blank image area
      // array operation on x
      for ( j = nr1; j<=nr2;j++) 
         for (i = ns1; i <= ns2; i++) {
            irow = i - ns1 + 1; 
            jcol = j - nr1 + 1; 
            sqimg(irow,jcol) = 0.0;
      }
   }
}
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().

{
    int nring = 0;
    int i;

    i = mr;
    while (i<=nr) {
      nring++;
      numr(1,nring) = i;
      i=i+iskip;
    }
}
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().

{
/*
c parabld  9/25/81 : parabolic fit to 3 by 3 peak neighborhood
c double precision version 
c
c the formula for paraboloid to be fiited into the nine points is:
c
c       f = c1 + c2*y + c3*y**2 + c4*x + c5*xy + c6*x**2
c
c the values of the coefficients c1 - c6 on the basis of the
c nine points around the peak, as evaluated by altran:
*/
   double c1,c2,c3,c4,c5,c6,denom;
   float  xmin, ymin;

   c1 = (26.*z33(1,1) - z33(1,2) + 2*z33(1,3) - z33(2,1) - 19.*z33(2,2)
         -7.*z33(2,3) + 2.*z33(3,1) - 7.*z33(3,2) + 14.*z33(3,3))/9.0;

   c2 = (8.* z33(1,1) - 8.*z33(1,2) + 5.*z33(2,1) - 8.*z33(2,2) + 3.*z33(2,3)
        +2.*z33(3,1) - 8.*z33(3,2) + 6.*z33(3,3))/(-6.);

   c3 = (z33(1,1) - 2.*z33(1,2) + z33(1,3) + z33(2,1) -2.*z33(2,2)
        + z33(2,3) + z33(3,1) - 2.*z33(3,2) + z33(3,3))/6.0;

   c4 = (8.*z33(1,1) + 5.*z33(1,2) + 2.*z33(1,3) -8.*z33(2,1) -8.*z33(2,2)
       - 8.*z33(2,3) + 3.*z33(3,2) + 6.*z33(3,3))/(-6.0);

   c5 = (z33(1,1) - z33(1,3) - z33(3,1) + z33(3,3))/4.0;

   c6 = (z33(1,1) + z33(1,2) + z33(1,3) - 2.*z33(2,1) - 2.*z33(2,2)
        -2.*z33(2,3) + z33(3,1) + z33(3,2) + z33(3,3))/6.0;

   // the peak coordinates of the paraboloid can now be evaluated as:

   *ysh=0.0;
   *xsh=0.0;
   denom=4.*c3*c6 - c5*c5;
   if (denom != 0.0) {
      *ysh=(c4*c5 - 2.*c2*c6) /denom-2.0;
      *xsh=(c2*c5 - 2.*c4*c3) /denom-2.0;
      *peakv= 4.*c1*c3*c6 - c1*c5*c5 -c2*c2*c6 + c2*c4*c5 - c4*c4*c3;
      *peakv= *peakv/denom;
      // limit interplation to +/- 1. range
      xmin = min(*xsh,1.0);
      ymin = min(*ysh,1.0);
      *xsh=max(xmin,-1.0);
      *ysh=max(ymin,-1.0);
   } 
}
void prb1d ( double *  b,
int  npoint,
float *  pos 
)

Definition at line 787 of file spidutil.cpp.

References b.

{
   double  c2,c3;
   int     nhalf;

   nhalf = npoint/2 + 1;
   *pos  = 0.0;

   if (npoint == 7) {
      c2 = 49.*b(1) + 6.*b(2) - 21.*b(3) - 32.*b(4) - 27.*b(5)
         - 6.*b(6) + 31.*b(7);
      c3 = 5.*b(1) - 3.*b(3) - 4.*b(4) - 3.*b(5) + 5.*b(7);
   } 
   else if (npoint == 5) {
      c2 = (74.*b(1) - 23.*b(2) - 60.*b(3) - 37.*b(4)
         + 46.*b(5) ) / (-70.);
      c3 = (2.*b(1) - b(2) - 2.*b(3) - b(4) + 2.*b(5) ) / 14.0;
   }
   else if (npoint == 3) {
      c2 = (5.*b(1) - 8.*b(2) + 3.*b(3) ) / (-2.0);
      c3 = (b(1) - 2.*b(2) + b(3) ) / 2.0;
   }
   else if (npoint == 9) {
      c2 = (1708.*b(1) + 581.*b(2) - 246.*b(3) - 773.*b(4)
         - 1000.*b(5) - 927.*b(6) - 554.*b(7) + 119.*b(8)
         + 1092.*b(9) ) / (-4620.);
      c3 = (28.*b(1) + 7.*b(2) - 8.*b(3) - 17.*b(4) - 20.*b(5)
         - 17.*b(6) - 8.*b(7) + 7.*b(8) + 28.*b(9) ) / 924.0;
   }
   if (c3 != 0.0)  *pos = c2/(2.0*c3) - nhalf;
}
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:

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:

  • f0 = (x0, y0)
  • f1 = (xb, y0)
  • f2 = (xa, y0)
  • f3 = (x0, yb)
  • f4 = (x0, ya)
  • fc = (xc, yc)

Mesh spacings:

  • hxa -- x- mesh spacing to the left of f0
  • hxb -- x- mesh spacing to the right of f0
  • hyb -- y- mesh spacing above f0
  • hya -- y- mesh spacing below f0

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]xx-coord value
[in]yy-coord value
nx
ny
[in]imageImage object (pointer)
Returns:
Interpolated value

Definition at line 186 of file spidutil.cpp.

References fdata, quadri(), x, and y.

{
/*
c  purpose: quadratic interpolation 
c 
c  parameters:       xx,yy treated as circularly closed.
c                    fdata - image 1..nxdata, 1..nydata
c
c                    f3    fc       f0, f1, f2, f3 are the values
c                     +             at the grid points.  x is the
c                     + x           point at which the function
c              f2++++f0++++f1       is to be estimated. (it need
c                     +             not be in the first quadrant).
c                     +             fc - the outer corner point
c                    f4             nearest x.
c
c                                   f0 is the value of the fdata at
c                                   fdata(i,j), it is the interior mesh
c                                   point nearest  x.
c                                   the coordinates of f0 are (x0,y0),
c                                   the coordinates of f1 are (xb,y0),
c                                   the coordinates of f2 are (xa,y0),
c                                   the coordinates of f3 are (x0,yb),
c                                   the coordinates of f4 are (x0,ya),
c                                   the coordinates of fc are (xc,yc),
c
c                   o               hxa, hxb are the mesh spacings
c                   +               in the x-direction to the left
c                  hyb              and right of the center point.
c                   +
c            ++hxa++o++hxb++o       hyb, hya are the mesh spacings
c                   +               in the y-direction.
c                  hya
c                   +               hxc equals either  hxb  or  hxa
c                   o               depending on where the corner
c                                   point is located.
c
c                                   construct the interpolant
c                                   f = f0 + c1*(x-x0) +
c                                       c2*(x-x0)*(x-x1) +
c                                       c3*(y-y0) + c4*(y-y0)*(y-y1)
c                                       + c5*(x-x0)*(y-y0)
c
c
*/
    float x, y, dx0, dy0, f0, c1, c2, c3, c4, c5, dxb, dyb;
    float quadri;
    int   i, j, ip1, im1, jp1, jm1, ic, jc, hxc, hyc;

    x = xx;
    y = yy;

    // circular closure
    if (x < 1.0)               x = x+(1 - floor(x) / nxdata) * nxdata;
    if (x > (float)nxdata+0.5) x = fmod(x-1.0,(float)nxdata) + 1.0;
    if (y < 1.0)               y = y+(1 - floor(y) / nydata) * nydata;
    if (y > (float)nydata+0.5) y = fmod(y-1.0,(float)nydata) + 1.0;


    i   = (int) floor(x);
    j   = (int) floor(y);

    dx0 = x - i;
    dy0 = y - j;

    ip1 = i + 1;
    im1 = i - 1;
    jp1 = j + 1;
    jm1 = j - 1;

    if (ip1 > nxdata) ip1 = ip1 - nxdata;
    if (im1 < 1)      im1 = im1 + nxdata;
    if (jp1 > nydata) jp1 = jp1 - nydata;
    if (jm1 < 1)      jm1 = jm1 + nydata;

    f0  = fdata(i,j);
    c1  = fdata(ip1,j) - f0;
    c2  = (c1 - f0 + fdata(im1,j)) * 0.5;
    c3  = fdata(i,jp1) - f0;
    c4  = (c3 - f0 + fdata(i,jm1)) * 0.5;

    dxb = dx0 - 1;
    dyb = dy0 - 1;

    // hxc & hyc are either 1 or -1
    if (dx0 >= 0) {
       hxc = 1;
    }
    else {
       hxc = -1;
    }
    if (dy0 >= 0) {
       hyc = 1;
    }
    else {
       hyc = -1;
    }
 
    ic  = i + hxc;
    jc  = j + hyc;

    if (ic > nxdata) {
       ic = ic - nxdata;
    }
    else if (ic < 1) {
       ic = ic + nxdata;
    }

    if (jc > nydata) {
       jc = jc - nydata;
    }
    else if (jc < 1) {
       jc = jc + nydata;
    }

    c5  =  ( (fdata(ic,jc) - f0 - hxc * c1 - (hxc * (hxc - 1.0)) * c2 
            - hyc * c3 - (hyc * (hyc - 1.0)) * c4) * (hxc * hyc));

    quadri = f0 + dx0 * (c1 + dxb * c2 + dy0 * c5) + dy0 * (c3 + dyb * c4);

    return quadri; 
}
void ringwe ( float *  wr,
int *  numr,
int  nring 
)

Definition at line 357 of file spidutil.cpp.

References numr, and wr.

{
   int i, maxrin;
   float dpi;

   dpi = 8.0*atan(1.0);
   maxrin = numr(3,nring); 
   for  (i=1;i<=nring;i++) {
      wr(i)=(float)(numr(1,i))*dpi/(float)(numr(3,i))
           *(float)(maxrin)/(float)(numr(3,i));
           cout << numr(1,i)<<"  "<< numr(2,i)<<"  "<< numr(3,i)<<"  " << wr(i)<<"  "<<endl;
   }
}
int setnring ( int  mr,
int  nr,
int  iskip 
)

Definition at line 372 of file spidutil.cpp.

Referenced by apmd(), and apmq().

{
    int nring = 0, i;

    i = mr;
    while (i<=nr) {
       nring = nring+1;
       i = i + iskip;
    } 
    return nring;
}
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().

{
/*
c adpated from SPIDER ap_getdat.f
c purpose:       read read windowed image date into array x for 'ap' ops.
c
c parameters:
c       lsam,lrow           image dimensions                  (input)
c       nsam,nrow           output image dimensions           (input)
c       lr1,lr2,ls1,ls2     output image window               (input)
c       imgfrom             input image                       (input)
c       imgto               output output                     (output)
c
*/
   int window;

   int k3, k2, kt;

   //     real, dimension(lsam)                        :: bufin

    window = 0;
    if (lr1 != 1 || ls1 != 1 || lr2 != lrow || ls2 != lsam) window = 1;

    if (window) {
      // window from the whole image
      for (k2=lr1;k2<=lr2;k2++) {
         kt = k2-lr1+1;
         for (k3=ls1;k3<=ls2;k3++) 
            imgto(k3-ls1+1,kt) = imgfrom(k3,k2);
      }
    }
    else {
      // do a copy
      for (k2=1;k2<=lrow;k2++) 
         for (k3=1;k3<=lsam;k3++) 
            imgto(k3,k2) = imgfrom(k3,k2);
    }
}