EMAN2
spidutil.cpp
Go to the documentation of this file.
00001 /*
00002  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00003  * Copyright (c) 2000-2006 Baylor College of Medicine
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00029  *
00030  */
00031 
00032 #include "ctf.h"
00033 #include "log.h"
00034 #include "emdata.h"
00035 #include "xydata.h"
00036 #include "Assert.h"
00037 #include "projector.h"
00038 
00039 #include <sys/types.h>
00040 #include <sys/times.h>
00041 #include <time.h>
00042 #include <sys/time.h>
00043 
00044 #include <iostream>
00045 
00046 
00047 using namespace EMAN;
00048 
00049 using std::cout;
00050 using std::cin;
00051 using std::endl;
00052 using std::string;
00053 
00054 #include "spidutil.h"
00055 #include "spidfft.h"
00056 
00057 int  aprings(int  nimg, int nring, float  *imgstk, int nsam, 
00058              int  nrow, int *numr, float *refcstk, int lcirc, 
00059              char mode)
00060 {
00061     int j, status;
00062     float *wr;
00063 
00064     status = 0;
00065     wr = (float*) calloc(nring, sizeof(float));
00066     if (!wr) {
00067         fprintf(stderr, "aprings: failed to allocate wr!\n");
00068         status = -1;
00069         goto EXIT;
00070     }
00071 
00072 
00073     // 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 }
00086 
00087 //-------------------------------------------------------------------
00088 
00089 int  apring1(float *sqimg, int nsam , int  nrow, float *imgcirc, int lcirc, 
00090              int nring   , char mode, int *numr, float *wr)
00091 {
00092    int    status = 0;
00093    int    nsb, nse, nrb, nre, maxrin, ltf, lt, lt2, lt3, ns2, nr2;
00094    int    inr, kcirc, jt, i, j;
00095    float  fnr2, fns2, yq, fi, dpi, x, y;
00096    double dfi;
00097 
00098    dpi    = 2.0*atan(1.0); 
00099 
00100    // 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 }
00184     
00185 // -----------------------------------------------------------------
00186 float quadri(float xx, float yy, int nxdata, int nydata, float *fdata)
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 }
00308 
00309 int alprbs(int *numr, int nring, int *lcirc, char mode)
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 }
00355 
00356 // -----------------------------------------------------------------
00357 void ringwe(float *wr, int *numr, int nring)
00358 {
00359    int i, maxrin;
00360    float dpi;
00361 
00362    dpi = 8.0*atan(1.0);
00363    maxrin = numr(3,nring); 
00364    for  (i=1;i<=nring;i++) {
00365       wr(i)=(float)(numr(1,i))*dpi/(float)(numr(3,i))
00366            *(float)(maxrin)/(float)(numr(3,i));
00367            cout << numr(1,i)<<"  "<< numr(2,i)<<"  "<< numr(3,i)<<"  " << wr(i)<<"  "<<endl;
00368    }
00369 }
00370 
00371 // -----------------------------------------------------------------
00372 int setnring(int mr, int nr, int iskip)
00373 {
00374     int nring = 0, i;
00375 
00376     i = mr;
00377     while (i<=nr) {
00378        nring = nring+1;
00379        i = i + iskip;
00380     } 
00381     return nring;
00382 }
00383 
00384 // -----------------------------------------------------------------
00385 void numrinit(int mr, int nr, int iskip, int *numr)
00386 {
00387     int nring = 0;
00388     int i;
00389 
00390     i = mr;
00391     while (i<=nr) {
00392       nring++;
00393       numr(1,nring) = i;
00394       i=i+iskip;
00395     }
00396 }
00397 
00398 // -----------------------------------------------------------------
00399 void normass(float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2, 
00400              int ir1, int ir2)
00401 {
00402 /*
00403 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 }
00460 
00461 // -----------------------------------------------------------------
00462 
00463 void frngs(float *circ, int *numr, int nring)
00464 {
00465    int i, l; 
00466  
00467    for (i=1; i<=nring;i++) {
00468      l=(int)(log2(numr(3,i)));
00469      fftr_q(&circ(numr(2,i)),l);
00470    }
00471 }
00472 
00473 // -----------------------------------------------------------------
00474 int  applyws(float *circ, int lcirc, int *numr, float *wr,
00475              int nring)
00476 {
00477    int   maxrin, numr3i, numr2i;
00478    float w;
00479    int   i, j, jc;
00480    int   status = 0;
00481 
00482    maxrin = numr(3,nring);
00483  
00484    for (i=1;i<=nring;i++) {
00485       numr3i=numr(3,i);
00486       numr2i=numr(2,i);
00487       w=wr(i);
00488       circ(numr2i)=circ(numr2i)*w;
00489       if (numr3i == maxrin) {
00490          circ(numr2i+1)=circ(numr2i+1)*w;
00491       }
00492       else { 
00493          circ(numr2i+1)=circ(numr2i+1)*0.5*w;
00494       }
00495       for (j=3;j<=numr3i;j++) {
00496          jc=j+numr2i-1;
00497          circ(jc)=circ(jc)*w;
00498       }
00499    } 
00500    if (jc >= lcirc) status = -1;
00501    return status; 
00502 }
00503 
00504 //--------------------------------------------------
00505 void normas(float *sqimg, int nsam, int ns1, int ns2, int nr1, int nr2, 
00506             int ir1, int ir2)
00507 {
00508 /*
00509 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 }
00556 
00557 //-------------------------------------------------------
00558 void alrq(float *xim,  int nsam , int nrow , int *numr,
00559           float *circ, int lcirc, int nring, char mode)
00560 {
00561 /* 
00562 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 }
00635 
00636 //---------------------------------------------------
00637 int crosrng_ms(float *circ1, float *circ2, int  lcirc, int  nring,
00638                int   maxrin, int   *numr , double *qn, float *tot,
00639                double   *qm, double *tmt)
00640 {
00641 /*
00642 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 }
00785 
00786 //---------------------------------------------------
00787 void prb1d(double *b, int npoint, float *pos)
00788 {
00789    double  c2,c3;
00790    int     nhalf;
00791 
00792    nhalf = npoint/2 + 1;
00793    *pos  = 0.0;
00794 
00795    if (npoint == 7) {
00796       c2 = 49.*b(1) + 6.*b(2) - 21.*b(3) - 32.*b(4) - 27.*b(5)
00797          - 6.*b(6) + 31.*b(7);
00798       c3 = 5.*b(1) - 3.*b(3) - 4.*b(4) - 3.*b(5) + 5.*b(7);
00799    } 
00800    else if (npoint == 5) {
00801       c2 = (74.*b(1) - 23.*b(2) - 60.*b(3) - 37.*b(4)
00802          + 46.*b(5) ) / (-70.);
00803       c3 = (2.*b(1) - b(2) - 2.*b(3) - b(4) + 2.*b(5) ) / 14.0;
00804    }
00805    else if (npoint == 3) {
00806       c2 = (5.*b(1) - 8.*b(2) + 3.*b(3) ) / (-2.0);
00807       c3 = (b(1) - 2.*b(2) + b(3) ) / 2.0;
00808    }
00809    else if (npoint == 9) {
00810       c2 = (1708.*b(1) + 581.*b(2) - 246.*b(3) - 773.*b(4)
00811          - 1000.*b(5) - 927.*b(6) - 554.*b(7) + 119.*b(8)
00812          + 1092.*b(9) ) / (-4620.);
00813       c3 = (28.*b(1) + 7.*b(2) - 8.*b(3) - 17.*b(4) - 20.*b(5)
00814          - 17.*b(6) - 8.*b(7) + 7.*b(8) + 28.*b(9) ) / 924.0;
00815    }
00816    if (c3 != 0.0)  *pos = c2/(2.0*c3) - nhalf;
00817 }
00818 // -----------------------------------------------
00819 
00820 void apmaster_1(char mode, float *divas, int nr, int *numth,
00821                 int  lsam, int lrow, int *nsam, int *nrow)
00822 {
00823 /*
00824 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 }
00859 //----------------------------------------------------------------
00860 void win_resize(float *imgfrom, float *imgto, int lsam, int lrow, 
00861                 int nsam, int nrow, int lr1, int lr2, int ls1, int ls2)
00862 {
00863 /*
00864 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 }
00899 
00900 //--------------------------------------------------
00901 int apmd(EMData *refprj, Dict refparams, EMData *expimg, APMDopt options,
00902          float *fangles)
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 }
01097 
01098 //----------------------------------------------------------------
01099 int aprq2d(float   *sqimg, float     *bfc, int     *numr, int      nsam, 
01100            int       nrow, int   ishrange, int     istep, int       nsb, 
01101            int        nse, int        nrb, int       nre, int     lcirc, 
01102            int      nring, int     maxrin, int      nima, char     mode, 
01103            float *refdirs, float  *expdir, float   range, float  *diref, 
01104            float   *ccrot, float *rangnew, float *xshsum, float *yshsum,  
01105            int   *nimalcg, int   ckmirror, int limitrange)
01106 /*
01107 c 
01108 c  parameters:
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 }
01403 
01404 
01405 //-----------------------------------------------
01406 float ang_n(float rkk, char mode, int maxrin)
01407 {
01408     float ang; 
01409 
01410     if (mode == 'H' || mode == 'h') {
01411         ang = fmod(((rkk-1.0) / maxrin+1.0)*180.0, 180.0);
01412     }
01413     else if ( mode == 'F' || mode == 'f') {
01414         ang = fmod(((rkk-1.0) / maxrin+1.0)*360.0, 360.0);
01415     }
01416     return ang;
01417 }
01418 //-----------------------------------------------
01419 void alrq_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
01420              int  *numr, float *circ, int lcirc, int  nring, char  mode)
01421 {
01422    double dpi, dfi;
01423    int    it, jt, inr, l, nsim, kcirc, lt;
01424    float  yq, xold, yold, fi, x, y;
01425 
01426    //     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 }
01490 //-----------------------------------------------
01491 void parabld(double *z33, float *xsh, float *ysh, double *peakv)
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 }
01541 //-----------------------------------------------
01542 void crosrng_e(float *circ1, float *circ2, int lcirc,
01543                int    nring, int   maxrin, int *numr,
01544                double *qn, float *tot, int neg)
01545 {
01546 /*
01547 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 }
01637 //-----------------------------------------------
01638 int apmq(EMData *refprj, Dict refparams, EMData *expimg, APMQopt options,
01639          float  *angles, float *shifts)
01640 {
01641     double quadpi = 3.14159265358979323846;
01642     double dgr_to_rad =  quadpi/180;
01643 
01644     int    mr, nr, iskip;
01645     float  range; 
01646     char   mode;
01647 
01648     int    nsam, nrow, numth, limitrange, maxrin, nimalcg, nring, status;
01649     int    nx, ny, nxr, nyr, nidi, nima, lcirc, nwsam, nwrow, nsb, nse,
01650            nrb, nre, it, ishrange, istep, ckmirror, nrad, nang;
01651     int    iref, ireft, mirrornew, imgref, ngotpar, mirrorold;
01652     float  ccrot, rangnew, xshnew, yshnew, peakv; 
01653     float  rangout, angdif, rangold, xshold, yshold, c, s;
01654  
01655     float  *refstk, *imgstk, *bfc, *imgwindow, *refdirs, *expdirs, 
01656            *angexps, *expdir;
01657 
01658     vector <float> angrefs;
01659     float  *dlist;
01660     int    ldd = 7;
01661     int    *numr;
01662 
01663     float  divas;
01664 
01665     status = 0;
01666 
01667     // 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 }