util_sparx.cpp

Go to the documentation of this file.

/*
 * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
 * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

#ifdef _WIN32
        #pragma warning(disable:4819)
#endif  //_WIN32

#include <cstring>
#include <ctime>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <boost/format.hpp>
#include "emdata.h"
#include "util.h"
#include "fundamentals.h"
#include "lapackblas.h"
#include "lbfgsb.h"
using namespace EMAN;
#include "steepest.h"
#include "emassert.h"
#include "randnum.h"

#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_bessel.h>
#include <cmath>
using namespace std;
using std::complex;

vector<float> Util::infomask(EMData* Vol, EMData* mask, bool flip = false)
//  flip true:  find statistics under the mask (mask >0.5)
//  flip false: find statistics ourside the mask (mask <0.5)
{
        ENTERFUNC;
        vector<float> stats;
        float *Volptr, *maskptr,MAX,MIN;
        long double Sum1,Sum2;
        long count;

        MAX = -FLT_MAX;
        MIN =  FLT_MAX;
        count = 0L;
        Sum1 = 0.0L;
        Sum2 = 0.0L;

        if (mask == NULL) {
           //Vol->update_stat();
           stats.push_back(Vol->get_attr("mean"));
           stats.push_back(Vol->get_attr("sigma"));
           stats.push_back(Vol->get_attr("minimum"));
           stats.push_back(Vol->get_attr("maximum"));
           return stats;
        }

        /* Check if the sizes of the mask and image are same */

        size_t nx = Vol->get_xsize();
        size_t ny = Vol->get_ysize();
        size_t nz = Vol->get_zsize();

        size_t mask_nx = mask->get_xsize();
        size_t mask_ny = mask->get_ysize();
        size_t mask_nz = mask->get_zsize();

        if  (nx != mask_nx || ny != mask_ny || nz != mask_nz )
                throw ImageDimensionException("The dimension of the image does not match the dimension of the mask!");

 /*       if (nx != mask_nx ||
            ny != mask_ny ||
            nz != mask_nz  ) {
           // should throw an exception here!!! (will clean it up later CY)
           fprintf(stderr, "The dimension of the image does not match the dimension of the mask!\n");
           fprintf(stderr, " nx = %d, mask_nx = %d\n", nx, mask_nx);
           fprintf(stderr, " ny = %d, mask_ny = %d\n", ny, mask_ny);
           fprintf(stderr, " nz = %d, mask_nz = %d\n", nz, mask_nz);
           exit(1);
        }
 */
        Volptr = Vol->get_data();
        maskptr = mask->get_data();

        for (size_t i = 0; i < nx*ny*nz; i++) {
              if (maskptr[i]>0.5f == flip) {
                Sum1 += Volptr[i];
                Sum2 += Volptr[i]*Volptr[i];
                MAX = (MAX < Volptr[i])?Volptr[i]:MAX;
                MIN = (MIN > Volptr[i])?Volptr[i]:MIN;
                count++;
              }
        }

        if (count == 0) {
                LOGERR("Invalid mask");
                throw ImageFormatException( "Invalid mask");
        }

       float avg = static_cast<float>(Sum1/count);
       float sig2 = static_cast<float>(Sum2 - count*avg*avg)/(count-1);
       float sig = sqrt(sig2);

       stats.push_back(avg);
       stats.push_back(sig);
       stats.push_back(MIN);
       stats.push_back(MAX);

       return stats;
}


//----------------------------------------------------------------------------------------------------------

Dict Util::im_diff(EMData* V1, EMData* V2, EMData* mask)
{
        ENTERFUNC;

        if (!EMUtil::is_same_size(V1, V2)) {
                LOGERR("images not same size");
                throw ImageFormatException( "images not same size");
        }

        size_t nx = V1->get_xsize();
        size_t ny = V1->get_ysize();
        size_t nz = V1->get_zsize();
        size_t size = nx*ny*nz;

        EMData *BD = new EMData();
        BD->set_size(nx, ny, nz);

        float *params = new float[2];

        float *V1ptr, *V2ptr, *MASKptr, *BDptr, A, B;
        long double S1=0.L,S2=0.L,S3=0.L,S4=0.L;
        int nvox = 0L;

        V1ptr = V1->get_data();
        V2ptr = V2->get_data();
        BDptr = BD->get_data();


        if(!mask){
                EMData * Mask = new EMData();
                Mask->set_size(nx,ny,nz);
                Mask->to_one();
                MASKptr = Mask->get_data();
        } else {
                if (!EMUtil::is_same_size(V1, mask)) {
                        LOGERR("mask not same size");
                        throw ImageFormatException( "mask not same size");
                }

                MASKptr = mask->get_data();
        }



//       calculation of S1,S2,S3,S3,nvox

        for (size_t i = 0L;i < size; i++) {
              if (MASKptr[i]>0.5f) {
               S1 += V1ptr[i]*V2ptr[i];
               S2 += V1ptr[i]*V1ptr[i];
               S3 += V2ptr[i];
               S4 += V1ptr[i];
               nvox ++;
              }
        }

        if ((nvox*S1 - S3*S4) == 0. || (nvox*S2 - S4*S4) == 0) {
                A =1.0f ;
        } else {
                A = static_cast<float>( (nvox*S1 - S3*S4)/(nvox*S2 - S4*S4) );
        }
        B = static_cast<float> (A*S4  -  S3)/nvox;

        // calculation of the difference image

        for (size_t i = 0L;i < size; i++) {
             if (MASKptr[i]>0.5f) {
               BDptr[i] = A*V1ptr[i] -  B  - V2ptr[i];
             }  else  {
               BDptr[i] = 0.f;
             }
        }

        BD->update();

        params[0] = A;
        params[1] = B;

        Dict BDnParams;
        BDnParams["imdiff"] = BD;
        BDnParams["A"] = params[0];
        BDnParams["B"] = params[1];

        EXITFUNC;
        return BDnParams;
 }

//----------------------------------------------------------------------------------------------------------



EMData *Util::TwoDTestFunc(int Size, float p, float q,  float a, float b, int flag, float alphaDeg) //PRB
{
        ENTERFUNC;
        int Mid= (Size+1)/2;

        if (flag==0) { // This is the real function
                EMData* ImBW = new EMData();
                ImBW->set_size(Size,Size,1);
                ImBW->to_zero();

                float tempIm;
                float x,y;

                for (int ix=(1-Mid);  ix<Mid; ix++){
                        for (int iy=(1-Mid);  iy<Mid; iy++){
                                x = (float)ix;
                                y = (float)iy;
                        tempIm= static_cast<float>( (1/(2*M_PI)) * cos(p*x)* cos(q*y) * exp(-.5*x*x/(a*a))* exp(-.5*y*y/(b*b)) );
                                (*ImBW)(ix+Mid-1,iy+Mid-1) = tempIm * exp(.5f*p*p*a*a)* exp(.5f*q*q*b*b);
                        }
                }
                ImBW->update();
                ImBW->set_complex(false);
                ImBW->set_ri(true);

                return ImBW;
        }
        else if (flag==1) {  // This is the Fourier Transform
                EMData* ImBWFFT = new EMData();
                ImBWFFT ->set_size(2*Size,Size,1);
                ImBWFFT ->to_zero();

                float r,s;

                for (int ir=(1-Mid);  ir<Mid; ir++){
                        for (int is=(1-Mid);  is<Mid; is++){
                                r = (float)ir;
                                s = (float)is;
                        (*ImBWFFT)(2*(ir+Mid-1),is+Mid-1)= cosh(p*r*a*a) * cosh(q*s*b*b) *
                                exp(-.5f*r*r*a*a)* exp(-.5f*s*s*b*b);
                        }
                }
                ImBWFFT->update();
                ImBWFFT->set_complex(true);
                ImBWFFT->set_ri(true);
                ImBWFFT->set_shuffled(true);
                ImBWFFT->set_fftodd(true);

                return ImBWFFT;
        }
        else if (flag==2 || flag==3) { //   This is the projection in Real Space
                float alpha = static_cast<float>( alphaDeg*M_PI/180.0 );
                float C=cos(alpha);
                float S=sin(alpha);
                float D= sqrt(S*S*b*b + C*C*a*a);
                //float D2 = D*D;   PAP - to get rid of warning

                float P = p * C *a*a/D ;
                float Q = q * S *b*b/D ;

                if (flag==2) {
                        EMData* pofalpha = new EMData();
                        pofalpha ->set_size(Size,1,1);
                        pofalpha ->to_zero();

                        float Norm0 =  D*(float)sqrt(2*pi);
                        float Norm1 =  exp( .5f*(P+Q)*(P+Q)) / Norm0 ;
                        float Norm2 =  exp( .5f*(P-Q)*(P-Q)) / Norm0 ;
                        float sD;

                        for (int is=(1-Mid);  is<Mid; is++){
                                sD = is/D ;
                                (*pofalpha)(is+Mid-1) =  Norm1 * exp(-.5f*sD*sD)*cos(sD*(P+Q))
                         + Norm2 * exp(-.5f*sD*sD)*cos(sD*(P-Q));
                        }
                        pofalpha-> update();
                        pofalpha-> set_complex(false);
                        pofalpha-> set_ri(true);

                        return pofalpha;
                }
                if (flag==3) { // This is the projection in Fourier Space
                        float vD;

                        EMData* pofalphak = new EMData();
                        pofalphak ->set_size(2*Size,1,1);
                        pofalphak ->to_zero();

                        for (int iv=(1-Mid);  iv<Mid; iv++){
                                vD = iv*D ;
                                (*pofalphak)(2*(iv+Mid-1)) =  exp(-.5f*vD*vD)*(cosh(vD*(P+Q)) + cosh(vD*(P-Q)) );
                        }
                        pofalphak-> update();
                        pofalphak-> set_complex(false);
                        pofalphak-> set_ri(true);

                        return pofalphak;
                }
        }
        else if (flag==4) {
                cout <<" FH under construction";
                EMData* OutFT= TwoDTestFunc(Size, p, q, a, b, 1);
                EMData* TryFH= OutFT -> real2FH(4.0);
                return TryFH;
        } else {
                cout <<" flag must be 0,1,2,3, or 4";
        }

        EXITFUNC;
        return 0;
}


void Util::spline_mat(float *x, float *y, int n,  float *xq, float *yq, int m) //PRB
{

        float x0= x[0];
        float x1= x[1];
        float x2= x[2];
        float y0= y[0];
        float y1= y[1];
        float y2= y[2];
        float yp1 =  (y1-y0)/(x1-x0) +  (y2-y0)/(x2-x0) - (y2-y1)/(x2-x1)  ;
        float xn  = x[n];
        float xnm1= x[n-1];
        float xnm2= x[n-2];
        float yn  = y[n];
        float ynm1= y[n-1];
        float ynm2= y[n-2];
        float ypn=  (yn-ynm1)/(xn-xnm1) +  (yn-ynm2)/(xn-xnm2) - (ynm1-ynm2)/(xnm1-xnm2) ;
        float *y2d = new float[n];
        Util::spline(x,y,n,yp1,ypn,y2d);
        Util::splint(x,y,y2d,n,xq,yq,m); //PRB
        delete [] y2d;
        return;
}


void Util::spline(float *x, float *y, int n, float yp1, float ypn, float *y2) //PRB
{
        int i,k;
        float p, qn, sig, un, *u;
        u = new float[n-1];

        if (yp1 > .99e30){
                y2[0]=u[0]=0.0;
        } else {
                y2[0]=-.5f;
                u[0] =(3.0f/ (x[1] -x[0]))*( (y[1]-y[0])/(x[1]-x[0]) -yp1);
        }

        for (i=1; i < n-1; i++) {
                sig= (x[i] - x[i-1])/(x[i+1] - x[i-1]);
                p = sig*y2[i-1] + 2.0f;
                y2[i]  = (sig-1.0f)/p;
                u[i] = (y[i+1] - y[i] )/(x[i+1]-x[i] ) -  (y[i] - y[i-1] )/(x[i] -x[i-1]);
                u[i] = (6.0f*u[i]/ (x[i+1]-x[i-1]) - sig*u[i-1])/p;
        }

        if (ypn>.99e30){
                qn=0; un=0;
        } else {
                qn= .5f;
                un= (3.0f/(x[n-1] -x[n-2])) * (ypn -  (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
        }
        y2[n-1]= (un - qn*u[n-2])/(qn*y2[n-2]+1.0f);
        for (k=n-2; k>=0; k--){
                y2[k]=y2[k]*y2[k+1]+u[k];
        }
        delete [] u;
}


void Util::splint( float *xa, float *ya, float *y2a, int n,  float *xq, float *yq, int m) //PRB
{
        int klo, khi, k;
        float h, b, a;

//      klo=0; // can try to put here
        for (int j=0; j<m;j++){
                klo=0;
                khi=n-1;
                while (khi-klo >1) {
                        k=(khi+klo) >>1;
                        if  (xa[k]>xq[j]){ khi=k;}
                        else { klo=k;}
                }
                h=xa[khi]- xa[klo];
                if (h==0.0) printf("Bad XA input to routine SPLINT \n");
                a =(xa[khi]-xq[j])/h;
                b=(xq[j]-xa[klo])/h;
                yq[j]=a*ya[klo] + b*ya[khi]
                        + ((a*a*a-a)*y2a[klo]
                             +(b*b*b-b)*y2a[khi]) *(h*h)/6.0f;
        }
//      printf("h=%f, a = %f, b=%f, ya[klo]=%f, ya[khi]=%f , yq=%f\n",h, a, b, ya[klo], ya[khi],yq[0]);
}


void Util::Radialize(int *PermMatTr, float *kValsSorted,   // PRB
               float *weightofkValsSorted, int Size, int *SizeReturned)
{
        int iMax = (int) floor( (Size-1.0)/2 +.01);
        int CountMax = (iMax+2)*(iMax+1)/2;
        int Count=-1;
        float *kVals     = new float[CountMax];
        float *weightMat = new float[CountMax];
        int *PermMat     = new   int[CountMax];
        SizeReturned[0] = CountMax;

//      printf("Aa \n");        fflush(stdout);
        for (int jkx=0; jkx< iMax+1; jkx++) {
                for (int jky=0; jky< jkx+1; jky++) {
                        Count++;
                        kVals[Count] = sqrtf((float) (jkx*jkx +jky*jky));
                        weightMat[Count]=  1.0;
                        if (jkx!=0)  { weightMat[Count] *=2;}
                        if (jky!=0)  { weightMat[Count] *=2;}
                        if (jkx!=jky){ weightMat[Count] *=2;}
                        PermMat[Count]=Count+1;
                }
        }

        int lkVals = Count+1;
//      printf("Cc \n");fflush(stdout);

        sort_mat(&kVals[0],&kVals[Count],
             &PermMat[0],  &PermMat[Count]);  //PermMat is
                                //also returned as well as kValsSorted
        fflush(stdout);

        int newInd;

        for (int iP=0; iP < lkVals ; iP++ ) {
                newInd =  PermMat[iP];
                PermMatTr[newInd-1] = iP+1;
        }

//      printf("Ee \n"); fflush(stdout);

        int CountA=-1;
        int CountB=-1;

        while (CountB< (CountMax-1)) {
                CountA++;
                CountB++;
//              printf("CountA=%d ; CountB=%d \n", CountA,CountB);fflush(stdout);
                kValsSorted[CountA] = kVals[CountB] ;
                if (CountB<(CountMax-1) ) {
                        while (fabs(kVals[CountB] -kVals[CountB+1])<.0000001  ) {
                                SizeReturned[0]--;
                                for (int iP=0; iP < lkVals; iP++){
//                                      printf("iP=%d \n", iP);fflush(stdout);
                                        if  (PermMatTr[iP]>CountA+1) {
                                                PermMatTr[iP]--;
                                        }
                                }
                                CountB++;
                        }
                }
        }


        for (int CountD=0; CountD < CountMax; CountD++) {
            newInd = PermMatTr[CountD];
            weightofkValsSorted[newInd-1] += weightMat[CountD];
        }

}


vector<float>
Util::even_angles(float delta, float t1, float t2, float p1, float p2)
{
        vector<float> angles;
        float psi = 0.0;
        if ((0.0 == t1)&&(0.0 == t2)||(t1 >= t2)) {
                t1 = 0.0f;
                t2 = 90.0f;
        }
        if ((0.0 == p1)&&(0.0 == p2)||(p1 >= p2)) {
                p1 = 0.0f;
                p2 = 359.9f;
        }
        bool skip = ((t1 < 90.0)&&(90.0 == t2)&&(0.0 == p1)&&(p2 > 180.0));
        for (float theta = t1; theta <= t2; theta += delta) {
                float detphi;
                int lt;
                if ((0.0 == theta)||(180.0 == theta)) {
                        detphi = 360.0f;
                        lt = 1;
                } else {
                        detphi = delta/sin(theta*static_cast<float>(dgr_to_rad));
                        lt = int((p2 - p1)/detphi)-1;
                        if (lt < 1) lt = 1;
                        detphi = (p2 - p1)/lt;
                }
                for (int i = 0; i < lt; i++) {
                        float phi = p1 + i*detphi;
                        if (skip&&(90.0 == theta)&&(phi > 180.0)) continue;
                        angles.push_back(phi);
                        angles.push_back(theta);
                        angles.push_back(psi);
                }
        }
        return angles;
}


#define  fdata(i,j)      fdata[ i-1 + (j-1)*nxdata ]
/*float Util::quadri(float xx, float yy, int nxdata, int nydata, float* fdata)
{

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

    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
        while ( x < 1.0 ) x += nxdata;
        while ( x >= (float)(nxdata+1) )  x -= nxdata;
        while ( y < 1.0 ) y += nydata;
        while ( y >= (float)(nydata+1) )  y -= nydata;


    i   = (int) x;
    j   = (int) 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;
}*/
float Util::quadri(float xx, float yy, int nxdata, int nydata, float* fdata)
{
//  purpose: quadratic interpolation
//  Optimized for speed, circular closer removed, checking of ranges removed
        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;

        //     any xx and yy
        while ( x < 1.0 )                 x += nxdata;
        while ( x >= (float)(nxdata+1) )  x -= nxdata;
        while ( y < 1.0 )                 y += nydata;
        while ( y >= (float)(nydata+1) )  y -= nydata;

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

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

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

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

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

        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 -= nxdata;  else if (ic < 1) ic += nxdata;
        if (jc > nydata) jc -= nydata;  else if (jc < 1) jc += nydata;

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


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

        return quadri;
}

#undef fdata

#define  fdata(i,j)      fdata[ i-1 + (j-1)*nxdata ]
float Util::quadri_background(float xx, float yy, int nxdata, int nydata, float* fdata, int xnew, int ynew)
{
//  purpose: quadratic interpolation
//  Optimized for speed, circular closer removed, checking of ranges removed
        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;

        // wrap around is not done circulantly; if (x,y) is not in the image, then x = xnew and y = ynew
        if ( (x < 1.0) || ( x >= (float)(nxdata+1) ) || ( y < 1.0 ) || ( y >= (float)(nydata+1) )){
              x = xnew;
                  y = ynew;
     }
            

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

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

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

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

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

        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 -= nxdata;  else if (ic < 1) ic += nxdata;
        if (jc > nydata) jc -= nydata;  else if (jc < 1) jc += nydata;

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


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

        return quadri;
}

#undef fdata


float  Util::get_pixel_conv_new(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb) {
// Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]

/* Commented by Zhengfan Yang on 04/20/07
This function is written to replace get_pixel_conv(), which is too slow to use in practice.
I made the following changes to get_pixel_conv():
1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
3. Unfold the 'for' loop
4. Reduce the usage of multiplications through some bracketing (from 98 times to 57 times in 2D case, from 1029 times to 400 times in 3D case)

The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
size, say N=5, you can easily modify it by referring my code.
*/
        int K = kb.get_window_size();
        int kbmin = -K/2;
        int kbmax = -kbmin;
        int kbc = kbmax+1;

        float pixel =0.0f;
        float w=0.0f;

        delx = restrict1(delx, nx);
        int inxold = int(round(delx));
        if ( ny < 2 ) {  //1D
                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                int x1, x2, x3, x4, x5, x6, x7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;
                }

                pixel = data[x1]*tablex1 + data[x2]*tablex2 + data[x3]*tablex3 +
                        data[x4]*tablex4 + data[x5]*tablex5 + data[x6]*tablex6 +
                        data[x7]*tablex7 ;

                w = tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7;
        } else if ( nz < 2 ) {  // 2D
                dely = restrict1(dely, ny);
                int inyold = int(round(dely));
                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                float tabley1 = kb.i0win_tab(dely-inyold+3);
                float tabley2 = kb.i0win_tab(dely-inyold+2);
                float tabley3 = kb.i0win_tab(dely-inyold+1);
                float tabley4 = kb.i0win_tab(dely-inyold);
                float tabley5 = kb.i0win_tab(dely-inyold-1);
                float tabley6 = kb.i0win_tab(dely-inyold-2);
                float tabley7 = kb.i0win_tab(dely-inyold-3);

                int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;

                        y1 = ((inyold-3+ny)%ny)*nx;
                        y2 = ((inyold-2+ny)%ny)*nx;
                        y3 = ((inyold-1+ny)%ny)*nx;
                        y4 = ((inyold  +ny)%ny)*nx;
                        y5 = ((inyold+1+ny)%ny)*nx;
                        y6 = ((inyold+2+ny)%ny)*nx;
                        y7 = ((inyold+3+ny)%ny)*nx;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;

                        y1 = (inyold-3)*nx;
                        y2 = (inyold-2)*nx;
                        y3 = (inyold-1)*nx;
                        y4 = inyold*nx;
                        y5 = (inyold+1)*nx;
                        y6 = (inyold+2)*nx;
                        y7 = (inyold+3)*nx;
                }

                pixel    = ( data[x1+y1]*tablex1 + data[x2+y1]*tablex2 + data[x3+y1]*tablex3 +
                             data[x4+y1]*tablex4 + data[x5+y1]*tablex5 + data[x6+y1]*tablex6 +
                             data[x7+y1]*tablex7 ) * tabley1 +
                           ( data[x1+y2]*tablex1 + data[x2+y2]*tablex2 + data[x3+y2]*tablex3 +
                             data[x4+y2]*tablex4 + data[x5+y2]*tablex5 + data[x6+y2]*tablex6 +
                             data[x7+y2]*tablex7 ) * tabley2 +
                           ( data[x1+y3]*tablex1 + data[x2+y3]*tablex2 + data[x3+y3]*tablex3 +
                             data[x4+y3]*tablex4 + data[x5+y3]*tablex5 + data[x6+y3]*tablex6 +
                             data[x7+y3]*tablex7 ) * tabley3 +
                           ( data[x1+y4]*tablex1 + data[x2+y4]*tablex2 + data[x3+y4]*tablex3 +
                             data[x4+y4]*tablex4 + data[x5+y4]*tablex5 + data[x6+y4]*tablex6 +
                             data[x7+y4]*tablex7 ) * tabley4 +
                           ( data[x1+y5]*tablex1 + data[x2+y5]*tablex2 + data[x3+y5]*tablex3 +
                             data[x4+y5]*tablex4 + data[x5+y5]*tablex5 + data[x6+y5]*tablex6 +
                             data[x7+y5]*tablex7 ) * tabley5 +
                           ( data[x1+y6]*tablex1 + data[x2+y6]*tablex2 + data[x3+y6]*tablex3 +
                             data[x4+y6]*tablex4 + data[x5+y6]*tablex5 + data[x6+y6]*tablex6 +
                             data[x7+y6]*tablex7 ) * tabley6 +
                           ( data[x1+y7]*tablex1 + data[x2+y7]*tablex2 + data[x3+y7]*tablex3 +
                             data[x4+y7]*tablex4 + data[x5+y7]*tablex5 + data[x6+y7]*tablex6 +
                             data[x7+y7]*tablex7 ) * tabley7;

                w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
                    (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
        } else {  //  3D
                dely = restrict1(dely, ny);
                int inyold = int(Util::round(dely));
                delz = restrict1(delz, nz);
                int inzold = int(Util::round(delz));

                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                float tabley1 = kb.i0win_tab(dely-inyold+3);
                float tabley2 = kb.i0win_tab(dely-inyold+2);
                float tabley3 = kb.i0win_tab(dely-inyold+1);
                float tabley4 = kb.i0win_tab(dely-inyold);
                float tabley5 = kb.i0win_tab(dely-inyold-1);
                float tabley6 = kb.i0win_tab(dely-inyold-2);
                float tabley7 = kb.i0win_tab(dely-inyold-3);

                float tablez1 = kb.i0win_tab(delz-inzold+3);
                float tablez2 = kb.i0win_tab(delz-inzold+2);
                float tablez3 = kb.i0win_tab(delz-inzold+1);
                float tablez4 = kb.i0win_tab(delz-inzold);
                float tablez5 = kb.i0win_tab(delz-inzold-1);
                float tablez6 = kb.i0win_tab(delz-inzold-2);
                float tablez7 = kb.i0win_tab(delz-inzold-3);

                int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7, z1, z2, z3, z4, z5, z6, z7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >= nz-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;

                        y1 = ((inyold-3+ny)%ny)*nx;
                        y2 = ((inyold-2+ny)%ny)*nx;
                        y3 = ((inyold-1+ny)%ny)*nx;
                        y4 = ((inyold  +ny)%ny)*nx;
                        y5 = ((inyold+1+ny)%ny)*nx;
                        y6 = ((inyold+2+ny)%ny)*nx;
                        y7 = ((inyold+3+ny)%ny)*nx;

                        z1 = ((inzold-3+nz)%nz)*nx*ny;
                        z2 = ((inzold-2+nz)%nz)*nx*ny;
                        z3 = ((inzold-1+nz)%nz)*nx*ny;
                        z4 = ((inzold  +nz)%nz)*nx*ny;
                        z5 = ((inzold+1+nz)%nz)*nx*ny;
                        z6 = ((inzold+2+nz)%nz)*nx*ny;
                        z7 = ((inzold+3+nz)%nz)*nx*ny;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;

                        y1 = (inyold-3)*nx;
                        y2 = (inyold-2)*nx;
                        y3 = (inyold-1)*nx;
                        y4 = inyold*nx;
                        y5 = (inyold+1)*nx;
                        y6 = (inyold+2)*nx;
                        y7 = (inyold+3)*nx;

                        z1 = (inzold-3)*nx*ny;
                        z2 = (inzold-2)*nx*ny;
                        z3 = (inzold-1)*nx*ny;
                        z4 = inzold*nx*ny;
                        z5 = (inzold+1)*nx*ny;
                        z6 = (inzold+2)*nx*ny;
                        z7 = (inzold+3)*nx*ny;
                }

                pixel  = ( ( data[x1+y1+z1]*tablex1 + data[x2+y1+z1]*tablex2 + data[x3+y1+z1]*tablex3 +
                             data[x4+y1+z1]*tablex4 + data[x5+y1+z1]*tablex5 + data[x6+y1+z1]*tablex6 +
                             data[x7+y1+z1]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z1]*tablex1 + data[x2+y2+z1]*tablex2 + data[x3+y2+z1]*tablex3 +
                             data[x4+y2+z1]*tablex4 + data[x5+y2+z1]*tablex5 + data[x6+y2+z1]*tablex6 +
                             data[x7+y2+z1]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z1]*tablex1 + data[x2+y3+z1]*tablex2 + data[x3+y3+z1]*tablex3 +
                             data[x4+y3+z1]*tablex4 + data[x5+y3+z1]*tablex5 + data[x6+y3+z1]*tablex6 +
                             data[x7+y3+z1]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z1]*tablex1 + data[x2+y4+z1]*tablex2 + data[x3+y4+z1]*tablex3 +
                             data[x4+y4+z1]*tablex4 + data[x5+y4+z1]*tablex5 + data[x6+y4+z1]*tablex6 +
                             data[x7+y4+z1]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z1]*tablex1 + data[x2+y5+z1]*tablex2 + data[x3+y5+z1]*tablex3 +
                             data[x4+y5+z1]*tablex4 + data[x5+y5+z1]*tablex5 + data[x6+y5+z1]*tablex6 +
                             data[x7+y5+z1]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z1]*tablex1 + data[x2+y6+z1]*tablex2 + data[x3+y6+z1]*tablex3 +
                             data[x4+y6+z1]*tablex4 + data[x5+y6+z1]*tablex5 + data[x6+y6+z1]*tablex6 +
                             data[x7+y6+z1]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z1]*tablex1 + data[x2+y7+z1]*tablex2 + data[x3+y7+z1]*tablex3 +
                             data[x4+y7+z1]*tablex4 + data[x5+y7+z1]*tablex5 + data[x6+y7+z1]*tablex6 +
                             data[x7+y7+z1]*tablex7 ) * tabley7 ) *tablez1 +
                         ( ( data[x1+y1+z2]*tablex1 + data[x2+y1+z2]*tablex2 + data[x3+y1+z2]*tablex3 +
                             data[x4+y1+z2]*tablex4 + data[x5+y1+z2]*tablex5 + data[x6+y1+z2]*tablex6 +
                             data[x7+y1+z2]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z2]*tablex1 + data[x2+y2+z2]*tablex2 + data[x3+y2+z2]*tablex3 +
                             data[x4+y2+z2]*tablex4 + data[x5+y2+z2]*tablex5 + data[x6+y2+z2]*tablex6 +
                             data[x7+y2+z2]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z2]*tablex1 + data[x2+y3+z2]*tablex2 + data[x3+y3+z2]*tablex3 +
                             data[x4+y3+z2]*tablex4 + data[x5+y3+z2]*tablex5 + data[x6+y3+z2]*tablex6 +
                             data[x7+y3+z2]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z2]*tablex1 + data[x2+y4+z2]*tablex2 + data[x3+y4+z2]*tablex3 +
                             data[x4+y4+z2]*tablex4 + data[x5+y4+z2]*tablex5 + data[x6+y4+z2]*tablex6 +
                             data[x7+y4+z2]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z2]*tablex1 + data[x2+y5+z2]*tablex2 + data[x3+y5+z2]*tablex3 +
                             data[x4+y5+z2]*tablex4 + data[x5+y5+z2]*tablex5 + data[x6+y5+z2]*tablex6 +
                             data[x7+y5+z2]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z2]*tablex1 + data[x2+y6+z2]*tablex2 + data[x3+y6+z2]*tablex3 +
                             data[x4+y6+z2]*tablex4 + data[x5+y6+z2]*tablex5 + data[x6+y6+z2]*tablex6 +
                             data[x7+y6+z2]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z2]*tablex1 + data[x2+y7+z2]*tablex2 + data[x3+y7+z2]*tablex3 +
                             data[x4+y7+z2]*tablex4 + data[x5+y7+z2]*tablex5 + data[x6+y7+z2]*tablex6 +
                             data[x7+y7+z2]*tablex7 ) * tabley7 ) *tablez2 +
                         ( ( data[x1+y1+z3]*tablex1 + data[x2+y1+z3]*tablex2 + data[x3+y1+z3]*tablex3 +
                             data[x4+y1+z3]*tablex4 + data[x5+y1+z3]*tablex5 + data[x6+y1+z3]*tablex6 +
                             data[x7+y1+z3]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z3]*tablex1 + data[x2+y2+z3]*tablex2 + data[x3+y2+z3]*tablex3 +
                             data[x4+y2+z3]*tablex4 + data[x5+y2+z3]*tablex5 + data[x6+y2+z3]*tablex6 +
                             data[x7+y2+z3]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z3]*tablex1 + data[x2+y3+z3]*tablex2 + data[x3+y3+z3]*tablex3 +
                             data[x4+y3+z3]*tablex4 + data[x5+y3+z3]*tablex5 + data[x6+y3+z3]*tablex6 +
                             data[x7+y3+z3]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z3]*tablex1 + data[x2+y4+z3]*tablex2 + data[x3+y4+z3]*tablex3 +
                             data[x4+y4+z3]*tablex4 + data[x5+y4+z3]*tablex5 + data[x6+y4+z3]*tablex6 +
                             data[x7+y4+z3]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z3]*tablex1 + data[x2+y5+z3]*tablex2 + data[x3+y5+z3]*tablex3 +
                             data[x4+y5+z3]*tablex4 + data[x5+y5+z3]*tablex5 + data[x6+y5+z3]*tablex6 +
                             data[x7+y5+z3]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z3]*tablex1 + data[x2+y6+z3]*tablex2 + data[x3+y6+z3]*tablex3 +
                             data[x4+y6+z3]*tablex4 + data[x5+y6+z3]*tablex5 + data[x6+y6+z3]*tablex6 +
                             data[x7+y6+z3]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z3]*tablex1 + data[x2+y7+z3]*tablex2 + data[x3+y7+z3]*tablex3 +
                             data[x4+y7+z3]*tablex4 + data[x5+y7+z3]*tablex5 + data[x6+y7+z3]*tablex6 +
                             data[x7+y7+z3]*tablex7 ) * tabley7 ) *tablez3 +
                         ( ( data[x1+y1+z4]*tablex1 + data[x2+y1+z4]*tablex2 + data[x3+y1+z4]*tablex3 +
                             data[x4+y1+z4]*tablex4 + data[x5+y1+z4]*tablex5 + data[x6+y1+z4]*tablex6 +
                             data[x7+y1+z4]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z4]*tablex1 + data[x2+y2+z4]*tablex2 + data[x3+y2+z4]*tablex3 +
                             data[x4+y2+z4]*tablex4 + data[x5+y2+z4]*tablex5 + data[x6+y2+z4]*tablex6 +
                             data[x7+y2+z4]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z4]*tablex1 + data[x2+y3+z4]*tablex2 + data[x3+y3+z4]*tablex3 +
                             data[x4+y3+z4]*tablex4 + data[x5+y3+z4]*tablex5 + data[x6+y3+z4]*tablex6 +
                             data[x7+y3+z4]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z4]*tablex1 + data[x2+y4+z4]*tablex2 + data[x3+y4+z4]*tablex3 +
                             data[x4+y4+z4]*tablex4 + data[x5+y4+z4]*tablex5 + data[x6+y4+z4]*tablex6 +
                             data[x7+y4+z4]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z4]*tablex1 + data[x2+y5+z4]*tablex2 + data[x3+y5+z4]*tablex3 +
                             data[x4+y5+z4]*tablex4 + data[x5+y5+z4]*tablex5 + data[x6+y5+z4]*tablex6 +
                             data[x7+y5+z4]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z4]*tablex1 + data[x2+y6+z4]*tablex2 + data[x3+y6+z4]*tablex3 +
                             data[x4+y6+z4]*tablex4 + data[x5+y6+z4]*tablex5 + data[x6+y6+z4]*tablex6 +
                             data[x7+y6+z4]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z4]*tablex1 + data[x2+y7+z4]*tablex2 + data[x3+y7+z4]*tablex3 +
                             data[x4+y7+z4]*tablex4 + data[x5+y7+z4]*tablex5 + data[x6+y7+z4]*tablex6 +
                             data[x7+y7+z4]*tablex7 ) * tabley7 ) *tablez4 +
                         ( ( data[x1+y1+z5]*tablex1 + data[x2+y1+z5]*tablex2 + data[x3+y1+z5]*tablex3 +
                             data[x4+y1+z5]*tablex4 + data[x5+y1+z5]*tablex5 + data[x6+y1+z5]*tablex6 +
                             data[x7+y1+z5]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z5]*tablex1 + data[x2+y2+z5]*tablex2 + data[x3+y2+z5]*tablex3 +
                             data[x4+y2+z5]*tablex4 + data[x5+y2+z5]*tablex5 + data[x6+y2+z5]*tablex6 +
                             data[x7+y2+z5]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z5]*tablex1 + data[x2+y3+z5]*tablex2 + data[x3+y3+z5]*tablex3 +
                             data[x4+y3+z5]*tablex4 + data[x5+y3+z5]*tablex5 + data[x6+y3+z5]*tablex6 +
                             data[x7+y3+z5]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z5]*tablex1 + data[x2+y4+z5]*tablex2 + data[x3+y4+z5]*tablex3 +
                             data[x4+y4+z5]*tablex4 + data[x5+y4+z5]*tablex5 + data[x6+y4+z5]*tablex6 +
                             data[x7+y4+z5]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z5]*tablex1 + data[x2+y5+z5]*tablex2 + data[x3+y5+z5]*tablex3 +
                             data[x4+y5+z5]*tablex4 + data[x5+y5+z5]*tablex5 + data[x6+y5+z5]*tablex6 +
                             data[x7+y5+z5]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z5]*tablex1 + data[x2+y6+z5]*tablex2 + data[x3+y6+z5]*tablex3 +
                             data[x4+y6+z5]*tablex4 + data[x5+y6+z5]*tablex5 + data[x6+y6+z5]*tablex6 +
                             data[x7+y6+z5]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z5]*tablex1 + data[x2+y7+z5]*tablex2 + data[x3+y7+z5]*tablex3 +
                             data[x4+y7+z5]*tablex4 + data[x5+y7+z5]*tablex5 + data[x6+y7+z5]*tablex6 +
                             data[x7+y7+z5]*tablex7 ) * tabley7 ) *tablez5 +
                         ( ( data[x1+y1+z6]*tablex1 + data[x2+y1+z6]*tablex2 + data[x3+y1+z6]*tablex3 +
                             data[x4+y1+z6]*tablex4 + data[x5+y1+z6]*tablex5 + data[x6+y1+z6]*tablex6 +
                             data[x7+y1+z6]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z6]*tablex1 + data[x2+y2+z6]*tablex2 + data[x3+y2+z6]*tablex3 +
                             data[x4+y2+z6]*tablex4 + data[x5+y2+z6]*tablex5 + data[x6+y2+z6]*tablex6 +
                             data[x7+y2+z6]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z6]*tablex1 + data[x2+y3+z6]*tablex2 + data[x3+y3+z6]*tablex3 +
                             data[x4+y3+z6]*tablex4 + data[x5+y3+z6]*tablex5 + data[x6+y3+z6]*tablex6 +
                             data[x7+y3+z6]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z6]*tablex1 + data[x2+y4+z6]*tablex2 + data[x3+y4+z6]*tablex3 +
                             data[x4+y4+z6]*tablex4 + data[x5+y4+z6]*tablex5 + data[x6+y4+z6]*tablex6 +
                             data[x7+y4+z6]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z6]*tablex1 + data[x2+y5+z6]*tablex2 + data[x3+y5+z6]*tablex3 +
                             data[x4+y5+z6]*tablex4 + data[x5+y5+z6]*tablex5 + data[x6+y5+z6]*tablex6 +
                             data[x7+y5+z6]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z6]*tablex1 + data[x2+y6+z6]*tablex2 + data[x3+y6+z6]*tablex3 +
                             data[x4+y6+z6]*tablex4 + data[x5+y6+z6]*tablex5 + data[x6+y6+z6]*tablex6 +
                             data[x7+y6+z6]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z6]*tablex1 + data[x2+y7+z6]*tablex2 + data[x3+y7+z6]*tablex3 +
                             data[x4+y7+z6]*tablex4 + data[x5+y7+z6]*tablex5 + data[x6+y7+z6]*tablex6 +
                             data[x7+y7+z6]*tablex7 ) * tabley7 ) *tablez6 +
                         ( ( data[x1+y1+z7]*tablex1 + data[x2+y1+z7]*tablex2 + data[x3+y1+z7]*tablex3 +
                             data[x4+y1+z7]*tablex4 + data[x5+y1+z7]*tablex5 + data[x6+y1+z7]*tablex6 +
                             data[x7+y1+z7]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z7]*tablex1 + data[x2+y2+z7]*tablex2 + data[x3+y2+z7]*tablex3 +
                             data[x4+y2+z7]*tablex4 + data[x5+y2+z7]*tablex5 + data[x6+y2+z7]*tablex6 +
                             data[x7+y2+z7]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z7]*tablex1 + data[x2+y3+z7]*tablex2 + data[x3+y3+z7]*tablex3 +
                             data[x4+y3+z7]*tablex4 + data[x5+y3+z7]*tablex5 + data[x6+y3+z7]*tablex6 +
                             data[x7+y3+z7]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z7]*tablex1 + data[x2+y4+z7]*tablex2 + data[x3+y4+z7]*tablex3 +
                             data[x4+y4+z7]*tablex4 + data[x5+y4+z7]*tablex5 + data[x6+y4+z7]*tablex6 +
                             data[x7+y4+z7]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z7]*tablex1 + data[x2+y5+z7]*tablex2 + data[x3+y5+z7]*tablex3 +
                             data[x4+y5+z7]*tablex4 + data[x5+y5+z7]*tablex5 + data[x6+y5+z7]*tablex6 +
                             data[x7+y5+z7]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z7]*tablex1 + data[x2+y6+z7]*tablex2 + data[x3+y6+z7]*tablex3 +
                             data[x4+y6+z7]*tablex4 + data[x5+y6+z7]*tablex5 + data[x6+y6+z7]*tablex6 +
                             data[x7+y6+z7]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z7]*tablex1 + data[x2+y7+z7]*tablex2 + data[x3+y7+z7]*tablex3 +
                             data[x4+y7+z7]*tablex4 + data[x5+y7+z7]*tablex5 + data[x6+y7+z7]*tablex6 +
                             data[x7+y7+z7]*tablex7 ) * tabley7 ) *tablez7;

                w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
                    (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7) *
                    (tablez1+tablez2+tablez3+tablez4+tablez5+tablez6+tablez7);
        }
        return pixel/w;
}

float  Util::get_pixel_conv_new_background(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb, int xnew, int ynew) {
// Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]

/* Commented by Zhengfan Yang on 04/20/07
This function is written to replace get_pixel_conv(), which is too slow to use in practice.
I made the following changes to get_pixel_conv():
1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
3. Unfold the 'for' loop
4. Reduce the usage of multiplications through some bracketing (from 98 times to 57 times in 2D case, from 1029 times to 400 times in 3D case)

The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
size, say N=5, you can easily modify it by referring my code.
*/
        int K = kb.get_window_size();
        int kbmin = -K/2;
        int kbmax = -kbmin;
        int kbc = kbmax+1;

        float pixel =0.0f;
        float w=0.0f;

    float argdelx = delx; // adding this for 2D case where the wrap around is not done circulantly using restrict1.
        delx = restrict1(delx, nx);
        int inxold = int(round(delx));
        if ( ny < 2 ) {  //1D
                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                int x1, x2, x3, x4, x5, x6, x7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;
                }

                pixel = data[x1]*tablex1 + data[x2]*tablex2 + data[x3]*tablex3 +
                        data[x4]*tablex4 + data[x5]*tablex5 + data[x6]*tablex6 +
                        data[x7]*tablex7 ;

                w = tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7;
        } else if ( nz < 2 ) {  // 2D
                
                delx = argdelx;
                // the wrap around is not done circulantly for 2D case; if (argdelx, argdely) is not in the image, then make them (xnew, ynew) which is definitely in the image
                if ((delx < 0.0f) || (delx >= (float) (nx)) || (dely < 0.0f) || (dely >= (float) (ny)) ){
                delx = (float)xnew*2.0;
                dely = (float)ynew*2.0;
                }
                
                int inxold = int(round(delx));
                int inyold = int(round(dely));
                
                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                float tabley1 = kb.i0win_tab(dely-inyold+3);
                float tabley2 = kb.i0win_tab(dely-inyold+2);
                float tabley3 = kb.i0win_tab(dely-inyold+1);
                float tabley4 = kb.i0win_tab(dely-inyold);
                float tabley5 = kb.i0win_tab(dely-inyold-1);
                float tabley6 = kb.i0win_tab(dely-inyold-2);
                float tabley7 = kb.i0win_tab(dely-inyold-3);

                int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;

                        y1 = ((inyold-3+ny)%ny)*nx;
                        y2 = ((inyold-2+ny)%ny)*nx;
                        y3 = ((inyold-1+ny)%ny)*nx;
                        y4 = ((inyold  +ny)%ny)*nx;
                        y5 = ((inyold+1+ny)%ny)*nx;
                        y6 = ((inyold+2+ny)%ny)*nx;
                        y7 = ((inyold+3+ny)%ny)*nx;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;

                        y1 = (inyold-3)*nx;
                        y2 = (inyold-2)*nx;
                        y3 = (inyold-1)*nx;
                        y4 = inyold*nx;
                        y5 = (inyold+1)*nx;
                        y6 = (inyold+2)*nx;
                        y7 = (inyold+3)*nx;
                }

                pixel    = ( data[x1+y1]*tablex1 + data[x2+y1]*tablex2 + data[x3+y1]*tablex3 +
                             data[x4+y1]*tablex4 + data[x5+y1]*tablex5 + data[x6+y1]*tablex6 +
                             data[x7+y1]*tablex7 ) * tabley1 +
                           ( data[x1+y2]*tablex1 + data[x2+y2]*tablex2 + data[x3+y2]*tablex3 +
                             data[x4+y2]*tablex4 + data[x5+y2]*tablex5 + data[x6+y2]*tablex6 +
                             data[x7+y2]*tablex7 ) * tabley2 +
                           ( data[x1+y3]*tablex1 + data[x2+y3]*tablex2 + data[x3+y3]*tablex3 +
                             data[x4+y3]*tablex4 + data[x5+y3]*tablex5 + data[x6+y3]*tablex6 +
                             data[x7+y3]*tablex7 ) * tabley3 +
                           ( data[x1+y4]*tablex1 + data[x2+y4]*tablex2 + data[x3+y4]*tablex3 +
                             data[x4+y4]*tablex4 + data[x5+y4]*tablex5 + data[x6+y4]*tablex6 +
                             data[x7+y4]*tablex7 ) * tabley4 +
                           ( data[x1+y5]*tablex1 + data[x2+y5]*tablex2 + data[x3+y5]*tablex3 +
                             data[x4+y5]*tablex4 + data[x5+y5]*tablex5 + data[x6+y5]*tablex6 +
                             data[x7+y5]*tablex7 ) * tabley5 +
                           ( data[x1+y6]*tablex1 + data[x2+y6]*tablex2 + data[x3+y6]*tablex3 +
                             data[x4+y6]*tablex4 + data[x5+y6]*tablex5 + data[x6+y6]*tablex6 +
                             data[x7+y6]*tablex7 ) * tabley6 +
                           ( data[x1+y7]*tablex1 + data[x2+y7]*tablex2 + data[x3+y7]*tablex3 +
                             data[x4+y7]*tablex4 + data[x5+y7]*tablex5 + data[x6+y7]*tablex6 +
                             data[x7+y7]*tablex7 ) * tabley7;

                w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
                    (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
        } else {  //  3D
                dely = restrict1(dely, ny);
                int inyold = int(Util::round(dely));
                delz = restrict1(delz, nz);
                int inzold = int(Util::round(delz));

                float tablex1 = kb.i0win_tab(delx-inxold+3);
                float tablex2 = kb.i0win_tab(delx-inxold+2);
                float tablex3 = kb.i0win_tab(delx-inxold+1);
                float tablex4 = kb.i0win_tab(delx-inxold);
                float tablex5 = kb.i0win_tab(delx-inxold-1);
                float tablex6 = kb.i0win_tab(delx-inxold-2);
                float tablex7 = kb.i0win_tab(delx-inxold-3);

                float tabley1 = kb.i0win_tab(dely-inyold+3);
                float tabley2 = kb.i0win_tab(dely-inyold+2);
                float tabley3 = kb.i0win_tab(dely-inyold+1);
                float tabley4 = kb.i0win_tab(dely-inyold);
                float tabley5 = kb.i0win_tab(dely-inyold-1);
                float tabley6 = kb.i0win_tab(dely-inyold-2);
                float tabley7 = kb.i0win_tab(dely-inyold-3);

                float tablez1 = kb.i0win_tab(delz-inzold+3);
                float tablez2 = kb.i0win_tab(delz-inzold+2);
                float tablez3 = kb.i0win_tab(delz-inzold+1);
                float tablez4 = kb.i0win_tab(delz-inzold);
                float tablez5 = kb.i0win_tab(delz-inzold-1);
                float tablez6 = kb.i0win_tab(delz-inzold-2);
                float tablez7 = kb.i0win_tab(delz-inzold-3);

                int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7, z1, z2, z3, z4, z5, z6, z7;

                if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >= nz-kbc-2 )  {
                        x1 = (inxold-3+nx)%nx;
                        x2 = (inxold-2+nx)%nx;
                        x3 = (inxold-1+nx)%nx;
                        x4 = (inxold  +nx)%nx;
                        x5 = (inxold+1+nx)%nx;
                        x6 = (inxold+2+nx)%nx;
                        x7 = (inxold+3+nx)%nx;

                        y1 = ((inyold-3+ny)%ny)*nx;
                        y2 = ((inyold-2+ny)%ny)*nx;
                        y3 = ((inyold-1+ny)%ny)*nx;
                        y4 = ((inyold  +ny)%ny)*nx;
                        y5 = ((inyold+1+ny)%ny)*nx;
                        y6 = ((inyold+2+ny)%ny)*nx;
                        y7 = ((inyold+3+ny)%ny)*nx;

                        z1 = ((inzold-3+nz)%nz)*nx*ny;
                        z2 = ((inzold-2+nz)%nz)*nx*ny;
                        z3 = ((inzold-1+nz)%nz)*nx*ny;
                        z4 = ((inzold  +nz)%nz)*nx*ny;
                        z5 = ((inzold+1+nz)%nz)*nx*ny;
                        z6 = ((inzold+2+nz)%nz)*nx*ny;
                        z7 = ((inzold+3+nz)%nz)*nx*ny;
                } else {
                        x1 = inxold-3;
                        x2 = inxold-2;
                        x3 = inxold-1;
                        x4 = inxold;
                        x5 = inxold+1;
                        x6 = inxold+2;
                        x7 = inxold+3;

                        y1 = (inyold-3)*nx;
                        y2 = (inyold-2)*nx;
                        y3 = (inyold-1)*nx;
                        y4 = inyold*nx;
                        y5 = (inyold+1)*nx;
                        y6 = (inyold+2)*nx;
                        y7 = (inyold+3)*nx;

                        z1 = (inzold-3)*nx*ny;
                        z2 = (inzold-2)*nx*ny;
                        z3 = (inzold-1)*nx*ny;
                        z4 = inzold*nx*ny;
                        z5 = (inzold+1)*nx*ny;
                        z6 = (inzold+2)*nx*ny;
                        z7 = (inzold+3)*nx*ny;
                }

                pixel  = ( ( data[x1+y1+z1]*tablex1 + data[x2+y1+z1]*tablex2 + data[x3+y1+z1]*tablex3 +
                             data[x4+y1+z1]*tablex4 + data[x5+y1+z1]*tablex5 + data[x6+y1+z1]*tablex6 +
                             data[x7+y1+z1]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z1]*tablex1 + data[x2+y2+z1]*tablex2 + data[x3+y2+z1]*tablex3 +
                             data[x4+y2+z1]*tablex4 + data[x5+y2+z1]*tablex5 + data[x6+y2+z1]*tablex6 +
                             data[x7+y2+z1]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z1]*tablex1 + data[x2+y3+z1]*tablex2 + data[x3+y3+z1]*tablex3 +
                             data[x4+y3+z1]*tablex4 + data[x5+y3+z1]*tablex5 + data[x6+y3+z1]*tablex6 +
                             data[x7+y3+z1]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z1]*tablex1 + data[x2+y4+z1]*tablex2 + data[x3+y4+z1]*tablex3 +
                             data[x4+y4+z1]*tablex4 + data[x5+y4+z1]*tablex5 + data[x6+y4+z1]*tablex6 +
                             data[x7+y4+z1]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z1]*tablex1 + data[x2+y5+z1]*tablex2 + data[x3+y5+z1]*tablex3 +
                             data[x4+y5+z1]*tablex4 + data[x5+y5+z1]*tablex5 + data[x6+y5+z1]*tablex6 +
                             data[x7+y5+z1]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z1]*tablex1 + data[x2+y6+z1]*tablex2 + data[x3+y6+z1]*tablex3 +
                             data[x4+y6+z1]*tablex4 + data[x5+y6+z1]*tablex5 + data[x6+y6+z1]*tablex6 +
                             data[x7+y6+z1]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z1]*tablex1 + data[x2+y7+z1]*tablex2 + data[x3+y7+z1]*tablex3 +
                             data[x4+y7+z1]*tablex4 + data[x5+y7+z1]*tablex5 + data[x6+y7+z1]*tablex6 +
                             data[x7+y7+z1]*tablex7 ) * tabley7 ) *tablez1 +
                         ( ( data[x1+y1+z2]*tablex1 + data[x2+y1+z2]*tablex2 + data[x3+y1+z2]*tablex3 +
                             data[x4+y1+z2]*tablex4 + data[x5+y1+z2]*tablex5 + data[x6+y1+z2]*tablex6 +
                             data[x7+y1+z2]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z2]*tablex1 + data[x2+y2+z2]*tablex2 + data[x3+y2+z2]*tablex3 +
                             data[x4+y2+z2]*tablex4 + data[x5+y2+z2]*tablex5 + data[x6+y2+z2]*tablex6 +
                             data[x7+y2+z2]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z2]*tablex1 + data[x2+y3+z2]*tablex2 + data[x3+y3+z2]*tablex3 +
                             data[x4+y3+z2]*tablex4 + data[x5+y3+z2]*tablex5 + data[x6+y3+z2]*tablex6 +
                             data[x7+y3+z2]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z2]*tablex1 + data[x2+y4+z2]*tablex2 + data[x3+y4+z2]*tablex3 +
                             data[x4+y4+z2]*tablex4 + data[x5+y4+z2]*tablex5 + data[x6+y4+z2]*tablex6 +
                             data[x7+y4+z2]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z2]*tablex1 + data[x2+y5+z2]*tablex2 + data[x3+y5+z2]*tablex3 +
                             data[x4+y5+z2]*tablex4 + data[x5+y5+z2]*tablex5 + data[x6+y5+z2]*tablex6 +
                             data[x7+y5+z2]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z2]*tablex1 + data[x2+y6+z2]*tablex2 + data[x3+y6+z2]*tablex3 +
                             data[x4+y6+z2]*tablex4 + data[x5+y6+z2]*tablex5 + data[x6+y6+z2]*tablex6 +
                             data[x7+y6+z2]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z2]*tablex1 + data[x2+y7+z2]*tablex2 + data[x3+y7+z2]*tablex3 +
                             data[x4+y7+z2]*tablex4 + data[x5+y7+z2]*tablex5 + data[x6+y7+z2]*tablex6 +
                             data[x7+y7+z2]*tablex7 ) * tabley7 ) *tablez2 +
                         ( ( data[x1+y1+z3]*tablex1 + data[x2+y1+z3]*tablex2 + data[x3+y1+z3]*tablex3 +
                             data[x4+y1+z3]*tablex4 + data[x5+y1+z3]*tablex5 + data[x6+y1+z3]*tablex6 +
                             data[x7+y1+z3]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z3]*tablex1 + data[x2+y2+z3]*tablex2 + data[x3+y2+z3]*tablex3 +
                             data[x4+y2+z3]*tablex4 + data[x5+y2+z3]*tablex5 + data[x6+y2+z3]*tablex6 +
                             data[x7+y2+z3]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z3]*tablex1 + data[x2+y3+z3]*tablex2 + data[x3+y3+z3]*tablex3 +
                             data[x4+y3+z3]*tablex4 + data[x5+y3+z3]*tablex5 + data[x6+y3+z3]*tablex6 +
                             data[x7+y3+z3]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z3]*tablex1 + data[x2+y4+z3]*tablex2 + data[x3+y4+z3]*tablex3 +
                             data[x4+y4+z3]*tablex4 + data[x5+y4+z3]*tablex5 + data[x6+y4+z3]*tablex6 +
                             data[x7+y4+z3]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z3]*tablex1 + data[x2+y5+z3]*tablex2 + data[x3+y5+z3]*tablex3 +
                             data[x4+y5+z3]*tablex4 + data[x5+y5+z3]*tablex5 + data[x6+y5+z3]*tablex6 +
                             data[x7+y5+z3]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z3]*tablex1 + data[x2+y6+z3]*tablex2 + data[x3+y6+z3]*tablex3 +
                             data[x4+y6+z3]*tablex4 + data[x5+y6+z3]*tablex5 + data[x6+y6+z3]*tablex6 +
                             data[x7+y6+z3]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z3]*tablex1 + data[x2+y7+z3]*tablex2 + data[x3+y7+z3]*tablex3 +
                             data[x4+y7+z3]*tablex4 + data[x5+y7+z3]*tablex5 + data[x6+y7+z3]*tablex6 +
                             data[x7+y7+z3]*tablex7 ) * tabley7 ) *tablez3 +
                         ( ( data[x1+y1+z4]*tablex1 + data[x2+y1+z4]*tablex2 + data[x3+y1+z4]*tablex3 +
                             data[x4+y1+z4]*tablex4 + data[x5+y1+z4]*tablex5 + data[x6+y1+z4]*tablex6 +
                             data[x7+y1+z4]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z4]*tablex1 + data[x2+y2+z4]*tablex2 + data[x3+y2+z4]*tablex3 +
                             data[x4+y2+z4]*tablex4 + data[x5+y2+z4]*tablex5 + data[x6+y2+z4]*tablex6 +
                             data[x7+y2+z4]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z4]*tablex1 + data[x2+y3+z4]*tablex2 + data[x3+y3+z4]*tablex3 +
                             data[x4+y3+z4]*tablex4 + data[x5+y3+z4]*tablex5 + data[x6+y3+z4]*tablex6 +
                             data[x7+y3+z4]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z4]*tablex1 + data[x2+y4+z4]*tablex2 + data[x3+y4+z4]*tablex3 +
                             data[x4+y4+z4]*tablex4 + data[x5+y4+z4]*tablex5 + data[x6+y4+z4]*tablex6 +
                             data[x7+y4+z4]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z4]*tablex1 + data[x2+y5+z4]*tablex2 + data[x3+y5+z4]*tablex3 +
                             data[x4+y5+z4]*tablex4 + data[x5+y5+z4]*tablex5 + data[x6+y5+z4]*tablex6 +
                             data[x7+y5+z4]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z4]*tablex1 + data[x2+y6+z4]*tablex2 + data[x3+y6+z4]*tablex3 +
                             data[x4+y6+z4]*tablex4 + data[x5+y6+z4]*tablex5 + data[x6+y6+z4]*tablex6 +
                             data[x7+y6+z4]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z4]*tablex1 + data[x2+y7+z4]*tablex2 + data[x3+y7+z4]*tablex3 +
                             data[x4+y7+z4]*tablex4 + data[x5+y7+z4]*tablex5 + data[x6+y7+z4]*tablex6 +
                             data[x7+y7+z4]*tablex7 ) * tabley7 ) *tablez4 +
                         ( ( data[x1+y1+z5]*tablex1 + data[x2+y1+z5]*tablex2 + data[x3+y1+z5]*tablex3 +
                             data[x4+y1+z5]*tablex4 + data[x5+y1+z5]*tablex5 + data[x6+y1+z5]*tablex6 +
                             data[x7+y1+z5]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z5]*tablex1 + data[x2+y2+z5]*tablex2 + data[x3+y2+z5]*tablex3 +
                             data[x4+y2+z5]*tablex4 + data[x5+y2+z5]*tablex5 + data[x6+y2+z5]*tablex6 +
                             data[x7+y2+z5]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z5]*tablex1 + data[x2+y3+z5]*tablex2 + data[x3+y3+z5]*tablex3 +
                             data[x4+y3+z5]*tablex4 + data[x5+y3+z5]*tablex5 + data[x6+y3+z5]*tablex6 +
                             data[x7+y3+z5]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z5]*tablex1 + data[x2+y4+z5]*tablex2 + data[x3+y4+z5]*tablex3 +
                             data[x4+y4+z5]*tablex4 + data[x5+y4+z5]*tablex5 + data[x6+y4+z5]*tablex6 +
                             data[x7+y4+z5]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z5]*tablex1 + data[x2+y5+z5]*tablex2 + data[x3+y5+z5]*tablex3 +
                             data[x4+y5+z5]*tablex4 + data[x5+y5+z5]*tablex5 + data[x6+y5+z5]*tablex6 +
                             data[x7+y5+z5]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z5]*tablex1 + data[x2+y6+z5]*tablex2 + data[x3+y6+z5]*tablex3 +
                             data[x4+y6+z5]*tablex4 + data[x5+y6+z5]*tablex5 + data[x6+y6+z5]*tablex6 +
                             data[x7+y6+z5]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z5]*tablex1 + data[x2+y7+z5]*tablex2 + data[x3+y7+z5]*tablex3 +
                             data[x4+y7+z5]*tablex4 + data[x5+y7+z5]*tablex5 + data[x6+y7+z5]*tablex6 +
                             data[x7+y7+z5]*tablex7 ) * tabley7 ) *tablez5 +
                         ( ( data[x1+y1+z6]*tablex1 + data[x2+y1+z6]*tablex2 + data[x3+y1+z6]*tablex3 +
                             data[x4+y1+z6]*tablex4 + data[x5+y1+z6]*tablex5 + data[x6+y1+z6]*tablex6 +
                             data[x7+y1+z6]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z6]*tablex1 + data[x2+y2+z6]*tablex2 + data[x3+y2+z6]*tablex3 +
                             data[x4+y2+z6]*tablex4 + data[x5+y2+z6]*tablex5 + data[x6+y2+z6]*tablex6 +
                             data[x7+y2+z6]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z6]*tablex1 + data[x2+y3+z6]*tablex2 + data[x3+y3+z6]*tablex3 +
                             data[x4+y3+z6]*tablex4 + data[x5+y3+z6]*tablex5 + data[x6+y3+z6]*tablex6 +
                             data[x7+y3+z6]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z6]*tablex1 + data[x2+y4+z6]*tablex2 + data[x3+y4+z6]*tablex3 +
                             data[x4+y4+z6]*tablex4 + data[x5+y4+z6]*tablex5 + data[x6+y4+z6]*tablex6 +
                             data[x7+y4+z6]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z6]*tablex1 + data[x2+y5+z6]*tablex2 + data[x3+y5+z6]*tablex3 +
                             data[x4+y5+z6]*tablex4 + data[x5+y5+z6]*tablex5 + data[x6+y5+z6]*tablex6 +
                             data[x7+y5+z6]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z6]*tablex1 + data[x2+y6+z6]*tablex2 + data[x3+y6+z6]*tablex3 +
                             data[x4+y6+z6]*tablex4 + data[x5+y6+z6]*tablex5 + data[x6+y6+z6]*tablex6 +
                             data[x7+y6+z6]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z6]*tablex1 + data[x2+y7+z6]*tablex2 + data[x3+y7+z6]*tablex3 +
                             data[x4+y7+z6]*tablex4 + data[x5+y7+z6]*tablex5 + data[x6+y7+z6]*tablex6 +
                             data[x7+y7+z6]*tablex7 ) * tabley7 ) *tablez6 +
                         ( ( data[x1+y1+z7]*tablex1 + data[x2+y1+z7]*tablex2 + data[x3+y1+z7]*tablex3 +
                             data[x4+y1+z7]*tablex4 + data[x5+y1+z7]*tablex5 + data[x6+y1+z7]*tablex6 +
                             data[x7+y1+z7]*tablex7 ) * tabley1 +
                           ( data[x1+y2+z7]*tablex1 + data[x2+y2+z7]*tablex2 + data[x3+y2+z7]*tablex3 +
                             data[x4+y2+z7]*tablex4 + data[x5+y2+z7]*tablex5 + data[x6+y2+z7]*tablex6 +
                             data[x7+y2+z7]*tablex7 ) * tabley2 +
                           ( data[x1+y3+z7]*tablex1 + data[x2+y3+z7]*tablex2 + data[x3+y3+z7]*tablex3 +
                             data[x4+y3+z7]*tablex4 + data[x5+y3+z7]*tablex5 + data[x6+y3+z7]*tablex6 +
                             data[x7+y3+z7]*tablex7 ) * tabley3 +
                           ( data[x1+y4+z7]*tablex1 + data[x2+y4+z7]*tablex2 + data[x3+y4+z7]*tablex3 +
                             data[x4+y4+z7]*tablex4 + data[x5+y4+z7]*tablex5 + data[x6+y4+z7]*tablex6 +
                             data[x7+y4+z7]*tablex7 ) * tabley4 +
                           ( data[x1+y5+z7]*tablex1 + data[x2+y5+z7]*tablex2 + data[x3+y5+z7]*tablex3 +
                             data[x4+y5+z7]*tablex4 + data[x5+y5+z7]*tablex5 + data[x6+y5+z7]*tablex6 +
                             data[x7+y5+z7]*tablex7 ) * tabley5 +
                           ( data[x1+y6+z7]*tablex1 + data[x2+y6+z7]*tablex2 + data[x3+y6+z7]*tablex3 +
                             data[x4+y6+z7]*tablex4 + data[x5+y6+z7]*tablex5 + data[x6+y6+z7]*tablex6 +
                             data[x7+y6+z7]*tablex7 ) * tabley6 +
                           ( data[x1+y7+z7]*tablex1 + data[x2+y7+z7]*tablex2 + data[x3+y7+z7]*tablex3 +
                             data[x4+y7+z7]*tablex4 + data[x5+y7+z7]*tablex5 + data[x6+y7+z7]*tablex6 +
                             data[x7+y7+z7]*tablex7 ) * tabley7 ) *tablez7;

                w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
                    (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7) *
                    (tablez1+tablez2+tablez3+tablez4+tablez5+tablez6+tablez7);
        }
        return pixel/w;
}

/*
complex<float> Util::extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb) {

        int nxreal = nx - 2;
        if (nxreal != ny)
                throw ImageDimensionException("extractpoint requires ny == nx");
        int nhalf = nxreal/2;
        int kbsize = kb.get_window_size();
        int kbmin = -kbsize/2;
        int kbmax = -kbmin;
        bool flip = (nuxnew < 0.f);
        if (flip) {
                nuxnew *= -1;
                nuynew *= -1;
        }
        // put (xnew,ynew) on a grid.  The indices will be wrong for
        // the Fourier elements in the image, but the grid sizing will
        // be correct.
        int ixn = int(Util::round(nuxnew));
        int iyn = int(Util::round(nuynew));
        // set up some temporary weighting arrays
        float* wy0 = new float[kbmax - kbmin + 1];
        float* wy = wy0 - kbmin; // wy[kbmin:kbmax]
        float* wx0 = new float[kbmax - kbmin + 1];
        float* wx = wx0 - kbmin;
        for (int i = kbmin; i <= kbmax; i++) {
                        int iyp = iyn + i;
                        wy[i] = kb.i0win_tab(nuynew - iyp);
                        int ixp = ixn + i;
                        wx[i] = kb.i0win_tab(nuxnew - ixp);
        }
        // restrict loops to non-zero elements
        int iymin = 0;
        for (int iy = kbmin; iy <= -1; iy++) {
                if (wy[iy] != 0.f) {
                        iymin = iy;
                        break;
                }
        }
        int iymax = 0;
        for (int iy = kbmax; iy >= 1; iy--) {
                if (wy[iy] != 0.f) {
                        iymax = iy;
                        break;
                }
        }
        int ixmin = 0;
        for (int ix = kbmin; ix <= -1; ix++) {
                if (wx[ix] != 0.f) {
                        ixmin = ix;
                        break;
                }
        }
        int ixmax = 0;
        for (int ix = kbmax; ix >= 1; ix--) {
                if (wx[ix] != 0.f) {
                        ixmax = ix;
                        break;
                }
        }
        float wsum = 0.0f;
        for (int iy = iymin; iy <= iymax; iy++)
                for (int ix = ixmin; ix <= ixmax; ix++)
                        wsum += wx[ix]*wy[iy];

        complex<float> result(0.f,0.f);
        if ((ixn >= -kbmin) && (ixn <= nhalf-1-kbmax) && (iyn >= -nhalf-kbmin) && (iyn <= nhalf-1-kbmax)) {
                // (xin,yin) not within window border from the edge
                for (int iy = iymin; iy <= iymax; iy++) {
                        int iyp = iyn + iy;
                        for (int ix = ixmin; ix <= ixmax; ix++) {
                                int ixp = ixn + ix;
                                float w = wx[ix]*wy[iy];
                                complex<float> val = fimage->cmplx(ixp,iyp);
                                result += val*w;
                        }
                }
        } else {
                // points that "stick out"
                for (int iy = iymin; iy <= iymax; iy++) {
                        int iyp = iyn + iy;
                        for (int ix = ixmin; ix <= ixmax; ix++) {
                                int ixp = ixn + ix;
                                bool mirror = false;
                                int ixt= ixp, iyt= iyp;
                                if (ixt < 0) {
                                        ixt = -ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (ixt > nhalf) {
                                        ixt = nxreal - ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (iyt > nhalf-1)  iyt -= nxreal;
                                if (iyt < -nhalf)   iyt += nxreal;
                                float w = wx[ix]*wy[iy];
                                complex<float> val = fimage->cmplx(ixt,iyt);
                                if (mirror)  result += conj(val)*w;
                                else         result += val*w;
                        }
                }
        }
        if (flip)  result = conj(result)/wsum;
        else result /= wsum;
        delete [] wx0;
        delete [] wy0;
        return result;
}*/

/*
complex<float> Util::extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb) {

        int nxreal = nx - 2;
        if (nxreal != ny)
                throw ImageDimensionException("extractpoint requires ny == nx");
        int nhalf = nxreal/2;
        bool flip = false;
        if (nuxnew < 0.f) {
                nuxnew *= -1;
                nuynew *= -1;
                flip = true;
        }
        if (nuynew >= nhalf-0.5)  {
                nuynew -= nxreal;
        } else if (nuynew < -nhalf-0.5) {
                nuynew += nxreal;
        }

        // put (xnew,ynew) on a grid.  The indices will be wrong for
        // the Fourier elements in the image, but the grid sizing will
        // be correct.
        int ixn = int(Util::round(nuxnew));
        int iyn = int(Util::round(nuynew));

        // set up some temporary weighting arrays
        static float wy[7];
        static float wx[7];

        float iynn = nuynew - iyn;
        wy[0] = kb.i0win_tab(iynn+3);
        wy[1] = kb.i0win_tab(iynn+2);
        wy[2] = kb.i0win_tab(iynn+1);
        wy[3] = kb.i0win_tab(iynn);
        wy[4] = kb.i0win_tab(iynn-1);
        wy[5] = kb.i0win_tab(iynn-2);
        wy[6] = kb.i0win_tab(iynn-3);

        float ixnn = nuxnew - ixn;
        wx[0] = kb.i0win_tab(ixnn+3);
        wx[1] = kb.i0win_tab(ixnn+2);
        wx[2] = kb.i0win_tab(ixnn+1);
        wx[3] = kb.i0win_tab(ixnn);
        wx[4] = kb.i0win_tab(ixnn-1);
        wx[5] = kb.i0win_tab(ixnn-2);
        wx[6] = kb.i0win_tab(ixnn-3);

        float wsum = (wx[0]+wx[1]+wx[2]+wx[3]+wx[4]+wx[5]+wx[6])*(wy[0]+wy[1]+wy[2]+wy[3]+wy[4]+wy[5]+wy[6]);

        complex<float> result(0.f,0.f);
        for (int iy = 0; iy < 7; iy++) {
                int iyp = iyn + iy - 3 ;
                for (int ix = 0; ix < 7; ix++) {
                        int ixp = ixn + ix - 3;
                        float w = wx[ix]*wy[iy];
                        complex<float> val = fimage->cmplx(ixp,iyp);
                        result += val*w;
                }
        }

        if (flip)  result = conj(result)/wsum;
        else result /= wsum;

        return result;
}*/


complex<float> Util::extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb) {

        int nxreal = nx - 2;
        if (nxreal != ny)
                throw ImageDimensionException("extractpoint requires ny == nx");
        int nhalf = nxreal/2;
        bool flip = (nuxnew < 0.f);
        if (flip) {
                nuxnew *= -1;
                nuynew *= -1;
        }
        if (nuynew >= nhalf-0.5)  {
                nuynew -= nxreal;
        } else if (nuynew < -nhalf-0.5) {
                nuynew += nxreal;
        }

        // put (xnew,ynew) on a grid.  The indices will be wrong for
        // the Fourier elements in the image, but the grid sizing will
        // be correct.
        int ixn = int(Util::round(nuxnew));
        int iyn = int(Util::round(nuynew));

        // set up some temporary weighting arrays
        static float wy[7];
        static float wx[7];

        float iynn = nuynew - iyn;
        wy[0] = kb.i0win_tab(iynn+3);
        wy[1] = kb.i0win_tab(iynn+2);
        wy[2] = kb.i0win_tab(iynn+1);
        wy[3] = kb.i0win_tab(iynn);
        wy[4] = kb.i0win_tab(iynn-1);
        wy[5] = kb.i0win_tab(iynn-2);
        wy[6] = kb.i0win_tab(iynn-3);

        float ixnn = nuxnew - ixn;
        wx[0] = kb.i0win_tab(ixnn+3);
        wx[1] = kb.i0win_tab(ixnn+2);
        wx[2] = kb.i0win_tab(ixnn+1);
        wx[3] = kb.i0win_tab(ixnn);
        wx[4] = kb.i0win_tab(ixnn-1);
        wx[5] = kb.i0win_tab(ixnn-2);
        wx[6] = kb.i0win_tab(ixnn-3);

        float wsum = (wx[0]+wx[1]+wx[2]+wx[3]+wx[4]+wx[5]+wx[6])*(wy[0]+wy[1]+wy[2]+wy[3]+wy[4]+wy[5]+wy[6]);

        complex<float> result(0.f,0.f);
        if ((ixn >= 3) && (ixn <= nhalf-3) && (iyn >= -nhalf+3) && (iyn <= nhalf-4)) {
                // (xin,yin) not within window border from the edge
                for (int iy = 0; iy < 7; iy++) {
                        int iyp = iyn + iy - 3 ;
                        for (int ix = 0; ix < 7; ix++) {
                                int ixp = ixn + ix - 3;
                                float w = wx[ix]*wy[iy];
                                complex<float> val = fimage->cmplx(ixp,iyp);
                                result += val*w;
                        }
                }
        } else {
                // points that "stick out"
                for (int iy = 0; iy < 7; iy++) {
                        int iyp = iyn + iy - 3;
                        for (int ix = 0; ix < 7; ix++) {
                                int ixp = ixn + ix - 3;
                                bool mirror = false;
                                int ixt = ixp, iyt = iyp;
                                if (ixt < 0) {
                                        ixt = -ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (ixt > nhalf) {
                                        ixt = nxreal - ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (iyt > nhalf-1)  iyt -= nxreal;
                                if (iyt < -nhalf)   iyt += nxreal;
                                float w = wx[ix]*wy[iy];
                                complex<float> val = fimage->cmplx(ixt,iyt);
                                if (mirror)  result += conj(val)*w;
                                else         result += val*w;
                        }
                }
        }
        if (flip)  result = conj(result)/wsum;
        else result /= wsum;
        return result;
}

/*
complex<float> Util::extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb) {

        int nxreal = nx - 2;
        if (nxreal != ny)
                throw ImageDimensionException("extractpoint requires ny == nx");
        int nhalf = nxreal/2;
        bool flip = (nuxnew < 0.f);
        if (flip) {
                nuxnew *= -1;
                nuynew *= -1;
        }
        // put (xnew,ynew) on a grid.  The indices will be wrong for
        // the Fourier elements in the image, but the grid sizing will
        // be correct.
        int ixn = int(Util::round(nuxnew));
        int iyn = int(Util::round(nuynew));
        // set up some temporary weighting arrays
        static float wy[7];
        static float wx[7];

        float iynn = nuynew - iyn;
        wy[0] = kb.i0win_tab(iynn+3);
        wy[1] = kb.i0win_tab(iynn+2);
        wy[2] = kb.i0win_tab(iynn+1);
        wy[3] = kb.i0win_tab(iynn);
        wy[4] = kb.i0win_tab(iynn-1);
        wy[5] = kb.i0win_tab(iynn-2);
        wy[6] = kb.i0win_tab(iynn-3);

        float ixnn = nuxnew - ixn;
        wx[0] = kb.i0win_tab(ixnn+3);
        wx[1] = kb.i0win_tab(ixnn+2);
        wx[2] = kb.i0win_tab(ixnn+1);
        wx[3] = kb.i0win_tab(ixnn);
        wx[4] = kb.i0win_tab(ixnn-1);
        wx[5] = kb.i0win_tab(ixnn-2);
        wx[6] = kb.i0win_tab(ixnn-3);

        float wsum = (wx[0]+wx[1]+wx[2]+wx[3]+wx[4]+wx[5]+wx[6])*(wy[0]+wy[1]+wy[2]+wy[3]+wy[4]+wy[5]+wy[6]);

        complex<float> result(0.f,0.f);

        if ((ixn >= 3) && (ixn <= nhalf-3) && (iyn >= -nhalf+3) && (iyn <= nhalf-4)) {
                // (xin,yin) not within window border from the edge
                result = ( fimage->cmplx(ixn-3,iyn-3)*wx[0] +
                           fimage->cmplx(ixn-2,iyn-3)*wx[1] +
                           fimage->cmplx(ixn-1,iyn-3)*wx[2] +
                           fimage->cmplx(ixn+0,iyn-3)*wx[3] +
                           fimage->cmplx(ixn+1,iyn-3)*wx[4] +
                           fimage->cmplx(ixn+2,iyn-3)*wx[5] +
                           fimage->cmplx(ixn+3,iyn-3)*wx[6] )*wy[0] +
                           ( fimage->cmplx(ixn-3,iyn-2)*wx[0] +
                           fimage->cmplx(ixn-2,iyn-2)*wx[1] +
                           fimage->cmplx(ixn-1,iyn-2)*wx[2] +
                           fimage->cmplx(ixn+0,iyn-2)*wx[3] +
                           fimage->cmplx(ixn+1,iyn-2)*wx[4] +
                           fimage->cmplx(ixn+2,iyn-2)*wx[5] +
                           fimage->cmplx(ixn+3,iyn-2)*wx[6] )*wy[1] +
                           ( fimage->cmplx(ixn-3,iyn-1)*wx[0] +
                           fimage->cmplx(ixn-2,iyn-1)*wx[1] +
                           fimage->cmplx(ixn-1,iyn-1)*wx[2] +
                           fimage->cmplx(ixn+0,iyn-1)*wx[3] +
                           fimage->cmplx(ixn+1,iyn-1)*wx[4] +
                           fimage->cmplx(ixn+2,iyn-1)*wx[5] +
                           fimage->cmplx(ixn+3,iyn-1)*wx[6] )*wy[2] +
                           ( fimage->cmplx(ixn-3,iyn+0)*wx[0] +
                           fimage->cmplx(ixn-2,iyn+0)*wx[1] +
                           fimage->cmplx(ixn-1,iyn+0)*wx[2] +
                           fimage->cmplx(ixn+0,iyn+0)*wx[3] +
                           fimage->cmplx(ixn+1,iyn+0)*wx[4] +
                           fimage->cmplx(ixn+2,iyn+0)*wx[5] +
                           fimage->cmplx(ixn+3,iyn+0)*wx[6] )*wy[3] +
                           ( fimage->cmplx(ixn-3,iyn+1)*wx[0] +
                           fimage->cmplx(ixn-2,iyn+1)*wx[1] +
                           fimage->cmplx(ixn-1,iyn+1)*wx[2] +
                           fimage->cmplx(ixn+0,iyn+1)*wx[3] +
                           fimage->cmplx(ixn+1,iyn+1)*wx[4] +
                           fimage->cmplx(ixn+2,iyn+1)*wx[5] +
                           fimage->cmplx(ixn+3,iyn+1)*wx[6] )*wy[4] +
                           ( fimage->cmplx(ixn-3,iyn+2)*wx[0] +
                           fimage->cmplx(ixn-2,iyn+2)*wx[1] +
                           fimage->cmplx(ixn-1,iyn+2)*wx[2] +
                           fimage->cmplx(ixn+0,iyn+2)*wx[3] +
                           fimage->cmplx(ixn+1,iyn+2)*wx[4] +
                           fimage->cmplx(ixn+2,iyn+2)*wx[5] +
                           fimage->cmplx(ixn+3,iyn+2)*wx[6] )*wy[5] +
                           ( fimage->cmplx(ixn-3,iyn+3)*wx[0] +
                           fimage->cmplx(ixn-2,iyn+3)*wx[1] +
                           fimage->cmplx(ixn-1,iyn+3)*wx[2] +
                           fimage->cmplx(ixn+0,iyn+3)*wx[3] +
                           fimage->cmplx(ixn+1,iyn+3)*wx[4] +
                           fimage->cmplx(ixn+2,iyn+3)*wx[5] +
                           fimage->cmplx(ixn+3,iyn+3)*wx[6] )*wy[6];

        } else {
                // points that "stick out"
                for (int iy = 0; iy < 7; iy++) {
                        int iyp = iyn + iy - 3;
                        for (int ix = 0; ix < 7; ix++) {
                                int ixp = ixn + ix - 3;
                                bool mirror = false;
                                int ixt= ixp, iyt= iyp;
                                if (ixt < 0) {
                                        ixt = -ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (ixt > nhalf) {
                                        ixt = nxreal - ixt;
                                        iyt = -iyt;
                                        mirror = !mirror;
                                }
                                if (iyt > nhalf-1)  iyt -= nxreal;
                                if (iyt < -nhalf)   iyt += nxreal;
                                float w = wx[ix]*wy[iy];
                                complex<float> val = fimage->cmplx(ixt,iyt);
                                if (mirror)  result += conj(val)*w;
                                else         result += val*w;
                        }
                }
        }
        if (flip)  result = conj(result)/wsum;
        else result /= wsum;
        return result;
}*/


float Util::triquad(float R, float S, float T, float* fdata)
{

    const float C2 = 0.5f;    //1.0 / 2.0;
    const float C4 = 0.25f;   //1.0 / 4.0;
    const float C8 = 0.125f;  //1.0 / 8.0;

    float  RS   = R * S;
    float  ST   = S * T;
    float  RT   = R * T;
    float  RST  = R * ST;

    float  RSQ  = 1-R*R;
    float  SSQ  = 1-S*S;
    float  TSQ  = 1-T*T;

    float  RM1  = (1-R);
    float  SM1  = (1-S);
    float  TM1  = (1-T);

    float  RP1  = (1+R);
    float  SP1  = (1+S);
    float  TP1  = (1+T);

    float triquad =
    (-C8) * RST * RM1  * SM1  * TM1 * fdata[0] +
        ( C4) * ST  * RSQ  * SM1  * TM1 * fdata[1] +
        ( C8) * RST * RP1  * SM1  * TM1 * fdata[2] +
        ( C4) * RT  * RM1  * SSQ  * TM1 * fdata[3] +
        (-C2) * T   * RSQ  * SSQ  * TM1 * fdata[4] +
        (-C4) * RT  * RP1  * SSQ  * TM1 * fdata[5] +
        ( C8) * RST * RM1  * SP1  * TM1 * fdata[6] +
        (-C4) * ST  * RSQ  * SP1  * TM1 * fdata[7] +
        (-C8) * RST * RP1  * SP1  * TM1 * fdata[8] +
//
        ( C4) * RS  * RM1  * SM1  * TSQ * fdata[9]  +
        (-C2) * S   * RSQ  * SM1  * TSQ * fdata[10] +
        (-C4) * RS  * RP1  * SM1  * TSQ * fdata[11] +
        (-C2) * R   * RM1  * SSQ  * TSQ * fdata[12] +
                      RSQ  * SSQ  * TSQ * fdata[13] +
        ( C2) * R   * RP1  * SSQ  * TSQ * fdata[14] +
        (-C4) * RS  * RM1  * SP1  * TSQ * fdata[15] +
        ( C2) * S   * RSQ  * SP1  * TSQ * fdata[16] +
        ( C4) * RS  * RP1  * SP1  * TSQ * fdata[17] +
 //
        ( C8) * RST * RM1  * SM1  * TP1 * fdata[18] +
        (-C4) * ST  * RSQ  * SM1  * TP1 * fdata[19] +
        (-C8) * RST * RP1  * SM1  * TP1 * fdata[20] +
        (-C4) * RT  * RM1  * SSQ  * TP1 * fdata[21] +
        ( C2) * T   * RSQ  * SSQ  * TP1 * fdata[22] +
        ( C4) * RT  * RP1  * SSQ  * TP1 * fdata[23] +
        (-C8) * RST * RM1  * SP1  * TP1 * fdata[24] +
        ( C4) * ST  * RSQ  * SP1  * TP1 * fdata[25] +
        ( C8) * RST * RP1  * SP1  * TP1 * fdata[26]   ;
     return triquad;
}

Util::sincBlackman::sincBlackman(int M_, float fc_, int ntable_)
                : M(M_), fc(fc_), ntable(ntable_) {
        // Sinc-Blackman kernel
        build_sBtable();
}

void Util::sincBlackman::build_sBtable() {
        sBtable.resize(ntable+1);
        int ltab = int(round(float(ntable)/1.25f));
        int M2 = M/2;
        fltb = float(ltab)/M2;
        for (int i=ltab+1; i <= ntable; i++) sBtable[i] = 0.0f;
        float x = 1.0e-7f;
        sBtable[0] = (float)(sin(twopi*fc*x)/x*(0.52-0.5*cos(twopi*(x-M2)/M)+0.08*cos(2*twopi*(x-M2)/M)));
        for (int i=1; i <= ltab; i++) {
                x = float(i)/fltb;
                sBtable[i] = (float)(sin(twopi*fc*x)/x*(0.52-0.5*cos(twopi*(x-M2)/M)+0.08*cos(2*twopi*(x-M2)/M)));
                //cout << "  "<<x<<"  "<<sBtable[i] <<endl;
        }
}

Util::KaiserBessel::KaiserBessel(float alpha_, int K_, float r_, float v_,
                                         int N_, float vtable_, int ntable_)
                : alpha(alpha_), v(v_), r(r_), N(N_), K(K_), vtable(vtable_),
                  ntable(ntable_) {
        // Default values are alpha=1.25, K=6, r=0.5, v = K/2
        if (0.f == v) v = float(K)/2;
        if (0.f == vtable) vtable = v;
        alphar = alpha*r;
        fac = static_cast<float>(twopi)*alphar*v;
        vadjust = 1.0f*v;
        facadj = static_cast<float>(twopi)*alphar*vadjust;
        build_I0table();
}

float Util::KaiserBessel::i0win(float x) const {
        float val0 = float(gsl_sf_bessel_I0(facadj));
        float absx = fabs(x);
        if (absx > vadjust) return 0.f;
        float rt = sqrt(1.f - pow(absx/vadjust, 2));
        float res = static_cast<float>(gsl_sf_bessel_I0(facadj*rt))/val0;
        return res;
}

void Util::KaiserBessel::build_I0table() {
        i0table.resize(ntable+1); // i0table[0:ntable]
        int ltab = int(round(float(ntable)/1.25f));
        fltb = float(ltab)/(K/2);
        float val0 = static_cast<float>(gsl_sf_bessel_I0(facadj));
        for (int i=ltab+1; i <= ntable; i++) i0table[i] = 0.f;
        for (int i=0; i <= ltab; i++) {
                float s = float(i)/fltb/N;
                if (s < vadjust) {
                        float rt = sqrt(1.f - pow(s/vadjust, 2));
                        i0table[i] = static_cast<float>(gsl_sf_bessel_I0(facadj*rt))/val0;
                } else {
                        i0table[i] = 0.f;
                }
//              cout << "  "<<s*N<<"  "<<i0table[i] <<endl;
        }
}

float Util::KaiserBessel::I0table_maxerror() {
        float maxdiff = 0.f;
        for (int i = 1; i <= ntable; i++) {
                float diff = fabs(i0table[i] - i0table[i-1]);
                if (diff > maxdiff) maxdiff = diff;
        }
        return maxdiff;
}

float Util::KaiserBessel::sinhwin(float x) const {
        float val0 = sinh(fac)/fac;
        float absx = fabs(x);
        if (0.0 == x) {
                float res = 1.0f;
                return res;
        } else if (absx == alphar) {
                return 1.0f/val0;
        } else if (absx < alphar) {
                float rt = sqrt(1.0f - pow((x/alphar), 2));
                float facrt = fac*rt;
                float res = (sinh(facrt)/facrt)/val0;
                return res;
        } else {
                float rt = sqrt(pow((x/alphar),2) - 1.f);
                float facrt = fac*rt;
                float res = (sin(facrt)/facrt)/val0;
                return res;
        }
}

float Util::FakeKaiserBessel::i0win(float x) const {
        float val0 = sqrt(facadj)*float(gsl_sf_bessel_I1(facadj));
        float absx = fabs(x);
        if (absx > vadjust) return 0.f;
        float rt = sqrt(1.f - pow(absx/vadjust, 2));
        float res = sqrt(facadj*rt)*float(gsl_sf_bessel_I1(facadj*rt))/val0;
        return res;
}

void Util::FakeKaiserBessel::build_I0table() {
        i0table.resize(ntable+1); // i0table[0:ntable]
        int ltab = int(round(float(ntable)/1.1f));
        fltb = float(ltab)/(K/2);
        float val0 = sqrt(facadj)*static_cast<float>(gsl_sf_bessel_I1(facadj));
        for (int i=ltab+1; i <= ntable; i++) i0table[i] = 0.f;
        for (int i=0; i <= ltab; i++) {
                float s = float(i)/fltb/N;
                if (s < vadjust) {
                        float rt = sqrt(1.f - pow(s/vadjust, 2));
                        i0table[i] = sqrt(facadj*rt)*static_cast<float>(gsl_sf_bessel_I1(facadj*rt))/val0;
                } else {
                        i0table[i] = 0.f;
                }
        }
}

float Util::FakeKaiserBessel::sinhwin(float x) const {
        float val0 = sinh(fac)/fac;
        float absx = fabs(x);
        if (0.0 == x) {
                float res = 1.0f;
                return res;
        } else if (absx == alphar) {
                return 1.0f/val0;
        } else if (absx < alphar) {
                float rt = sqrt(1.0f - pow((x/alphar), 2));
                float facrt = fac*rt;
                float res = (sinh(facrt)/facrt)/val0;
                return res;
        } else {
                float rt = sqrt(pow((x/alphar),2) - 1.f);
                float facrt = fac*rt;
                float res = (sin(facrt)/facrt)/val0;
                return res;
        }
}

#if 0 // 1-st order KB window
float Util::FakeKaiserBessel::sinhwin(float x) const {
        //float val0 = sinh(fac)/fac;
        float prefix = 2*facadj*vadjust/float(gsl_sf_bessel_I1(facadj));
        float val0 = prefix*(cosh(facadj) - sinh(facadj)/facadj);
        float absx = fabs(x);
        if (0.0 == x) {
                //float res = 1.0f;
                float res = val0;
                return res;
        } else if (absx == alphar) {
                //return 1.0f/val0;
                return prefix;
        } else if (absx < alphar) {
                float rt = sqrt(1.0f - pow((x/alphar), 2));
                //float facrt = fac*rt;
                float facrt = facadj*rt;
                //float res = (sinh(facrt)/facrt)/val0;
                float res = prefix*(cosh(facrt) - sinh(facrt)/facrt);
                return res;
        } else {
                float rt = sqrt(pow((x/alphar),2) - 1.f);
                //float facrt = fac*rt;
                float facrt = facadj*rt;
                //float res = (sin(facrt)/facrt)/val0;
                float res = prefix*(sin(facrt)/facrt - cos(facrt));
                return res;
        }
}
#endif // 0



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

EMData* Util::Polar2D(EMData* image, vector<int> numr, string cmode){
        int nsam = image->get_xsize();
        int nrow = image->get_ysize();
        int nring = numr.size()/3;
        int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        EMData* out = new EMData();
        out->set_size(lcirc,1,1);
        char mode = (cmode == "F" || cmode == "f") ? 'f' : 'h';
        float *xim  = image->get_data();
        float *circ = out->get_data();
/*   alrq(image->get_data(), nsam, nrow, &numr[0], out->get_data(), lcirc, nring, cmode);
   return out;
}
void Util::alrq(float *xim,  int nsam , int nrow , int *numr,
          float *circ, int lcirc, int nring, char mode)
{*/
/*
c
c  purpose:
c
c  resmaple 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);

        for (i=1;i<=nring;i++) {
                // radius of the ring
                inr = numr(1,i);
                yq  = static_cast<float>(inr);
                l   = numr(3,i);
                if (mode == 'h' || mode == 'H')  lt = l/2;
                else                             lt = l/4;

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

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

                for (j=1;j<=nsim;j++) {
                        fi               = static_cast<float>(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);
                        }
                }
        }
        return  out;
}

EMData* Util::Polar2Dm(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode){
        int nsam = image->get_xsize();
        int nrow = image->get_ysize();
        int nring = numr.size()/3;
        int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        EMData* out = new EMData();
        out->set_size(lcirc,1,1);
        char mode = (cmode == "F" || cmode == "f") ? 'f' : 'h';
        float *xim  = image->get_data();
        float *circ = out->get_data();
        double dpi, dfi;
        int    it, jt, inr, l, nsim, kcirc, lt;
        float  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);

                // "F" means a full circle interpolation
                // "H" means a half circle interpolation

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

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

                Assert( kcirc <= lcirc );
                circ(kcirc) = quadri(xold,yold,nsam,nrow,xim);    // Sampling on 90 degree

                xold  = inr+cns2;
                yold  = 0.0f+cnr2;
                Assert( lt+kcirc <= lcirc );
                circ(lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on 0 degree

                if ( mode == 'f' || mode == 'F' ) {
                        xold = 0.0f+cns2;
                        yold = -inr+cnr2;
                        Assert( lt+lt+kcirc <= lcirc );
                        circ(lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on 270 degree

                        xold = -inr+cns2;
                        yold = 0.0f+cnr2;
                        Assert(lt+lt+lt+kcirc <= lcirc );
                        circ(lt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim); // Sampling on 180 degree
                }

                for (jt=1; jt<=nsim; jt++) {
                        fi   = static_cast<float>(dfi * jt);
                        x    = sin(fi) * inr;
                        y    = cos(fi) * inr;

                        xold = x+cns2;
                        yold = y+cnr2;

                        Assert( jt+kcirc <= lcirc );
                        circ(jt+kcirc) = quadri(xold,yold,nsam,nrow,xim);      // Sampling on the first quadrant

                        xold = y+cns2;
                        yold = -x+cnr2;

                        Assert( jt+lt+kcirc <= lcirc );
                        circ(jt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);    // Sampling on the fourth quadrant

                        if ( mode == 'f' || mode == 'F' ) {
                                xold = -x+cns2;
                                yold = -y+cnr2;

                                Assert( jt+lt+lt+kcirc <= lcirc );
                                circ(jt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim); // Sampling on the third quadrant

                                xold = -y+cns2;
                                yold = x+cnr2;

                                Assert( jt+lt+lt+lt+kcirc <= lcirc );
                                circ(jt+lt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on the second quadrant
                        }
                } // end for jt
        } //end for it
        return out;
}

float Util::bilinear(float xold, float yold, int nsam, int nrow, float* xim)
{
/*
c  purpose: linear interpolation
  Optimized for speed, circular closer removed, checking of ranges removed
*/
    float bilinear;
    int   ixold, iyold;

/*
        float xdif, ydif, xrem, yrem;
        ixold   = (int) floor(xold);
        iyold   = (int) floor(yold);
        ydif = yold - iyold;
        yrem = 1.0f - ydif;

        //  May want to insert if?
//              IF ((IYOLD .GE. 1 .AND. IYOLD .LE. NROW-1) .AND.
//     &            (IXOLD .GE. 1 .AND. IXOLD .LE. NSAM-1)) THEN
//c                INSIDE BOUNDARIES OF OUTPUT IMAGE
        xdif = xold - ixold;
        xrem = 1.0f- xdif;
//                 RBUF(K) = YDIF*(BUF(NADDR+NSAM)*XREM
//     &                    +BUF(NADDR+NSAM+1)*XDIF)
//     &                    +YREM*(BUF(NADDR)*XREM + BUF(NADDR+1)*XDIF)
        bilinear = ydif*(xim(ixold,iyold+1)*xrem + xim(ixold+1,iyold+1)*xdif) +
                                        yrem*(xim(ixold,iyold)*xrem+xim(ixold+1,iyold)*xdif);

    return bilinear;
}
*/
        float xdif, ydif;

        ixold   = (int) xold;
        iyold   = (int) yold;
        ydif = yold - iyold;

        //  May want to insert it?
//              IF ((IYOLD .GE. 1 .AND. IYOLD .LE. NROW-1) .AND.
//     &            (IXOLD .GE. 1 .AND. IXOLD .LE. NSAM-1)) THEN
//c                INSIDE BOUNDARIES OF OUTPUT IMAGE
        xdif = xold - ixold;
        bilinear = xim(ixold, iyold) + ydif* (xim(ixold, iyold+1) - xim(ixold, iyold)) +
                   xdif* (xim(ixold+1, iyold) - xim(ixold, iyold) +
                           ydif* (xim(ixold+1, iyold+1) - xim(ixold+1, iyold) - xim(ixold, iyold+1) + xim(ixold, iyold)) );

        return bilinear;
}

void Util::alrl_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
             int  *numr, float *circ, int , int  nring, char  mode) {
        double dpi, dfi;
        int    it, jt, inr, l, nsim, kcirc, lt;
        float   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);

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

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


                xold  = 0.0f+cns2;
                yold  = inr+cnr2;

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

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

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

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

                for (jt=1; jt<=nsim; jt++) {
                        fi   = static_cast<float>(dfi * jt);
                        x    = sin(fi) * inr;
                        y    = cos(fi) * inr;

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

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

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

                                xold = -y+cns2;
                                yold = x+cnr2;
                                circ(jt+lt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);
                        }
                } // end for jt
        } //end for it
}
/*
void Util::alrl_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
             int  *numr, float *circ, int lcirc, int  nring, char  mode)
{
   double dpi, dfi;
   int    it, jt, inr, l, nsim, kcirc, lt, xold, yold;
   float  yq, 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 = (int) (0.0+cns2);
        yold = (int) (inr+cnr2);

        circ(kcirc) = xim(xold, yold);

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

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

         xold  = (int) (-inr+cns2);
         yold = (int) (0.0+cnr2);
         circ(lt+lt+lt+kcirc) = xim(xold, yold);
      }

      for (jt=1; jt<=nsim; jt++) {
         fi   = dfi * jt;
         x    = sin(fi) * yq;
         y    = cos(fi) * yq;

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

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

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

            xold  = (int) (-y+cns2);
            yold = (int) (x+cnr2);
            circ(jt+lt+lt+lt+kcirc) = xim(xold, yold);
         }
      } // end for jt
   } //end for it
}
*/
//xim((int) floor(xold), (int) floor(yold))
#undef  xim

EMData* Util::Polar2Dmi(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode, Util::KaiserBessel& kb){
// input image is twice the size of the original image
        int nring = numr.size()/3;
        int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        EMData* out = new EMData();
        out->set_size(lcirc,1,1);
        char mode = (cmode == "F" || cmode == "f") ? 'f' : 'h';
        float *circ = out->get_data();
        float *fimage = image->get_data();
        int nx = image->get_xsize();
        int ny = image->get_ysize();
        int nz = image->get_zsize();
        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  = static_cast<float>(inr);

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

                nsim  = lt - 1;
                dfi   = dpi / (nsim+1);
                kcirc = numr(2,it);
                xold  = 0.0f;
                yold  = static_cast<float>(inr);
                circ(kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//      circ(kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

                xold  = static_cast<float>(inr);
                yold  = 0.0f;
                circ(lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//      circ(lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

        if ( mode == 'f' || mode == 'F' ) {
                xold = 0.0f;
                yold = static_cast<float>(-inr);
                circ(lt+lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//         circ(lt+lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

                xold = static_cast<float>(-inr);
                yold = 0.0f;
                circ(lt+lt+lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//         circ(lt+lt+lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);
        }

        for (jt=1;jt<=nsim;jt++) {
                fi   = static_cast<float>(dfi * jt);
                x    = sin(fi) * yq;
                y    = cos(fi) * yq;

                xold = x;
                yold = y;
                circ(jt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//         circ(jt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

                xold = y;
                yold = -x;
                circ(jt+lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//         circ(jt+lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

        if ( mode == 'f' || mode == 'F' ) {
                xold = -x;
                yold = -y;
                circ(jt+lt+lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//            circ(jt+lt+lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);

                xold = -y;
                yold = x;
                circ(jt+lt+lt+lt+kcirc) = get_pixel_conv_new(nx,ny,nz,2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,fimage,kb);
//            circ(jt+lt+lt+lt+kcirc) = image->get_pixel_conv(2*(xold+cns2-1.0f),2*(yold+cnr2-1.0f),0,kb);
        }
        } // end for jt
        } //end for it
        return  out;
}

/*

        A set of 1-D power-of-two FFTs
        Pawel & Chao 01/20/06

fftr_q(xcmplx,nv)
  single precision

 dimension xcmplx(2,iabs(nv)/2);
 xcmplx(1,1) --- R(0), xcmplx(2,1) --- R(NV/2)
 xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary


fftr_d(xcmplx,nv)
  double precision

 dimension xcmplx(2,iabs(nv)/2);
 xcmplx(1,1) --- R(0), xcmplx(2,1) --- R(NV/2)
 xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary



*/
#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]
//-----------------------------------------
void Util::fftc_d(double *br, double *bi, int ln, int ks)
{
        double rni,sgn,tr1,tr2,ti1,ti2;
        double cc,c,ss,s,t,x2,x3,x4,x5;
        int    b3,b4,b5,b6,b7,b56;
        int    n, k, l, j, i, ix0, ix1, status=0;

        const double tab1[] = {
                9.58737990959775e-5,
                1.91747597310703e-4,
                3.83495187571395e-4,
                7.66990318742704e-4,
                1.53398018628476e-3,
                3.06795676296598e-3,
                6.13588464915449e-3,
                1.22715382857199e-2,
                2.45412285229123e-2,
                4.90676743274181e-2,
                9.80171403295604e-2,
                1.95090322016128e-1,
                3.82683432365090e-1,
                7.07106781186546e-1,
                1.00000000000000,
        };

        n=(int)pow(2.0f,ln);

        k=abs(ks);
        l=16-ln;
        b3=n*k;
        b6=b3;
        b7=k;
        if (ks > 0) {
                sgn=1.0f;
        } else {
                sgn=-1.0f;
                rni=1.0f/(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.0f*pow(tab1(l),2);
   c=1.0f-cc;
   l++;
   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 Util::fftc_q(float *br, float *bi, int ln, int ks)
{
        //  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;
        int status=0;

        const float tab1[] = {
                9.58737990959775e-5f,
                1.91747597310703e-4f,
                3.83495187571395e-4f,
                7.66990318742704e-4f,
                1.53398018628476e-3f,
                3.06795676296598e-3f,
                6.13588464915449e-3f,
                1.22715382857199e-2f,
                2.45412285229123e-2f,
                4.90676743274181e-2f,
                9.80171403295604e-2f,
                1.95090322016128e-1f,
                3.82683432365090e-1f,
                7.07106781186546e-1f,
                1.00000000000000f,
        };

        n=(int)pow(2.0f,ln);

        k=abs(ks);
        l=16-ln;
        b3=n*k;
        b6=b3;
        b7=k;
        if( ks > 0 ) {
                sgn=1.0f;
        } else {
                sgn=-1.0f;
                rni=1.0f/(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.0f*pow(tab1(l),2);
   c=1.0f-cc;
   l++;
   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  Util::fftr_q(float *xcmplx, int nv)
{
   // dimension xcmplx(2,1); xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary

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

        const float tab1[] = {
                9.58737990959775e-5f,
                1.91747597310703e-4f,
                3.83495187571395e-4f,
                7.66990318742704e-4f,
                1.53398018628476e-3f,
                3.06795676296598e-3f,
                6.13588464915449e-3f,
                1.22715382857199e-2f,
                2.45412285229123e-2f,
                4.90676743274181e-2f,
                9.80171403295604e-2f,
                1.95090322016128e-1f,
                3.82683432365090e-1f,
                7.07106781186546e-1f,
                1.00000000000000f,
        };

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

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

// -------------------------------------------
void  Util::fftr_d(double *xcmplx, int nv)
{
        // double precision  x(2,1)
        int    i1, i2,  nu, inv, nu1, n, isub, n2, i;
        double tr1,tr2,ti1,ti2,tr,ti;
        double cc,c,ss,s,t;
        const double tab1[] = {
                9.58737990959775e-5,
                1.91747597310703e-4,
                3.83495187571395e-4,
                7.66990318742704e-4,
                1.53398018628476e-3,
                3.06795676296598e-3,
                6.13588464915449e-3,
                1.22715382857199e-2,
                2.45412285229123e-2,
                4.90676743274181e-2,
                9.80171403295604e-2,
                1.95090322016128e-1,
                3.82683432365090e-1,
                7.07106781186546e-1,
                1.00000000000000,
        };

        nu=abs(nv);
        inv=nv/nu;
        nu1=nu-1;
        n=(int)pow(2.0f,nu1);
        isub=16-nu1;
        ss=-tab1(isub);
        cc=-2.0*pow(tab1(isub-1),2);
        c=1.0f;
        s=0.0f;
        n2=n/2;

        if ( inv > 0 ) {
                fftc_d(&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_d(&xcmplx(1,1),&xcmplx(2,1),nu1,-2);
        }
}
#undef  tab1
#undef  xcmplx
#undef  br
#undef  bi

// This function conducts the Single Precision Fourier Transform for a set of rings
void Util::Frngs(EMData* circp, vector<int> numr){
        int nring = numr.size()/3;
        float *circ = circp->get_data();
        int i, l;
        for (i=1; i<=nring;i++) {

#ifdef _WIN32
                l = (int)( log((float)numr(3,i))/log(2.0f) );
#else
                l=(int)(log2(numr(3,i)));
#endif  //_WIN32

                fftr_q(&circ(numr(2,i)),l);
        }
}
// This function conducts the Single Precision Inverse Fourier Transform for a set of rings
void Util::Frngs_inv(EMData* circp, vector<int> numr){
        int nring = numr.size()/3;
        float *circ = circp->get_data();
        int i, l;
        for (i=1; i<=nring;i++) {

#ifdef _WIN32
                l = (int)( log((float)numr(3,i))/log(2.0f) );
#else
                l=(int)(log2(numr(3,i)));
#endif  //_WIN32

                fftr_q(&circ(numr(2,i)),-l);
        }
}
#undef  circ

#define  b(i)            b[i-1]
void Util::prb1d(double *b, int npoint, float *pos) {
        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) {
        else  { // at least one has to be true!!
                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 = static_cast<float>(c2/(2.0*c3) - nhalf);
}
#undef  b

#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]
Dict Util::Crosrng_e(EMData*  circ1p, EMData* circ2p, vector<int> numr, int neg) {
        //  neg = 0 straight,  neg = 1 mirrored
        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];
        double qn;   float  tot;
        float *circ1 = circ1p->get_data();
        float *circ2 = circ2p->get_data();
/*
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)  removed +2 as it is only needed for other ffts
        double precision  q(maxrin)
        double precision  t7(-3:3)
*/
        float *t;
        double t7[7], *q;
        int    i, j, k, ip, jc, numr3i, numr2i, jtot = 0;
        float  pos;

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int) (log2(maxrin));
#endif  //_WIN32

        q = (double*)calloc(maxrin, sizeof(double));
        t = (float*)calloc(maxrin, sizeof(float));

//   cout << *qn <<"  " <<*tot<<"  "<<ip<<endl;
        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; 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);

        Dict retvals;
        retvals["qn"] = qn;
        retvals["tot"] = tot;
        return  retvals;
}

Dict Util::Crosrng_ew(EMData*  circ1p, EMData* circ2p, vector<int> numr, vector<float> w, int neg) {
   //  neg = 0 straight,  neg = 1 mirrored
        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];
        double qn;   float  tot;
        float *circ1 = circ1p->get_data();
        float *circ2 = circ2p->get_data();
/*
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  multiplication by weights!
c       automatic arrays
        dimension         t(maxrin)  removed +2 as it is only needed for other ffts
        double precision  q(maxrin)
        double precision  t7(-3:3)
*/
        float *t;
        double t7[7], *q;
        int    i, j, k, ip, jc, numr3i, numr2i, jtot = 0;
        float  pos;

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int) (log2(maxrin));
#endif  //_WIN32

        q = (double*)calloc(maxrin, sizeof(double));
        t = (float*)calloc(maxrin, sizeof(float));

//   cout << *qn <<"  " <<*tot<<"  "<<ip<<endl;
        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) += t(j)*w[i-1];
                } 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; j++) q(j) += t(j)*w[i-1];
                }
        }

        fftr_d(q,ip);

        qn = -1.0e20;
        for (j=1;j<=maxrin;j++) {
                //cout << j << "  " << q(j) << endl;
                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);

        Dict retvals;
        //tot = 1;
        //qn = q(1);
        retvals["qn"] = qn;
        retvals["tot"] = tot;

        if (q) free(q);

        return  retvals;
}

Dict Util::Crosrng_ms(EMData* circ1p, EMData* circ2p, vector<int> numr) {
        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];
        double qn; float tot; double qm; float tmt;
        float *circ1 = circ1p->get_data();
        float *circ2 = circ2p->get_data();
/*
c
c  checks both straight & mirrored positions
c
c  input - fourier transforms of rings!!
c  circ1 already multiplied by weights!
c
*/

        // dimension             circ1(lcirc),circ2(lcirc)

        // t(maxrin), q(maxrin), t7(-3:3)  //maxrin+2 removed
        double *t, *q, t7[7];

        int   ip, jc, numr3i, numr2i, i, j, k, jtot = 0;
        float t1, t2, t3, t4, c1, c2, d1, d2, pos;

        qn  = 0.0f;
        qm  = 0.0f;
        tot = 0.0f;
        tmt = 0.0f;
#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32
  //for (j=1; j<=maxrin;j++) cout <<"  "<<j<<"   "<<circ1(j)<<"   "<<circ2(j) <<endl;

        //  c - straight  = circ1 * conjg(circ2)
        //  zero q array

        q = (double*)calloc(maxrin,sizeof(double));

        //   t - mirrored  = conjg(circ1) * conjg(circ2)
        //   zero t array
        t = (double*)calloc(maxrin,sizeof(double));

   //   premultiply  arrays ie( circ12 = circ1 * circ2) much slower
        for (i=1; i<=nring; i++) {

                numr3i = numr(3,i);   // Number of samples of this ring
                numr2i = numr(2,i);   // The beginning point of this ring

                t1   = circ1(numr2i) * circ2(numr2i);
                q(1) += t1;
                t(1) += t1;

                t1   = circ1(numr2i+1) * circ2(numr2i+1);
                if (numr3i == maxrin)  {
                        q(2) += t1;
                        t(2) += t1;
                } else {
                        q(numr3i+1) += t1;
                        t(numr3i+1) += t1;
                }

                for (j=3; j<=numr3i; j += 2) {
                        jc     = j+numr2i-1;

// Here, (c1+c2i)*conj(d1+d2i) = (c1*d1+c2*d2)+(-c1*d2+c2*d1)i
//                                ----- -----    ----- -----
//                                 t1     t2      t3    t4
// Here, conj(c1+c2i)*conj(d1+d2i) = (c1*d1-c2*d2)+(-c1*d2-c2*d1)i
//                                    ----- -----    ----- -----
//                                     t1    t2       t3    t4

                        c1     = circ1(jc);
                        c2     = circ1(jc+1);
                        d1     = circ2(jc);
                        d2     = circ2(jc+1);

                        t1     = c1 * d1;
                        t2     = c2 * d2;
                        t3     = c1 * d2;
                        t4     = c2 * d1;

                        q(j)   += t1 + t2;
                        q(j+1) += -t3 + t4;
                        t(j)   += t1 - t2;
                        t(j+1) += -t3 - t4;
                }
        }
        //for (j=1; j<=maxrin; j++) cout <<"  "<<j<<"   "<<q(j) <<"   "<<t(j) <<endl;
        fftr_d(q,ip);

        qn  = -1.0e20;
        for (j=1; j<=maxrin; j++) {//cout <<"  "<<j<<"   "<<q(j) <<endl;
                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);
        }

        // interpolate
        prb1d(t7,7,&pos);
        tot = (float)(jtot)+pos;
        // Do not interpolate
        //tot = (float)(jtot);

        // mirrored
        fftr_d(t,ip);

        // find angle
        qm = -1.0e20;
        for (j=1; j<=maxrin;j++) {//cout <<"  "<<j<<"   "<<t(j) <<endl;
                if ( t(j) >= qm ) {
                        qm   = t(j);
                        jtot = j;
                }
        }

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

        // interpolate

        prb1d(t7,7,&pos);
        tmt = float(jtot) + pos;
        // Do not interpolate
        //tmt = float(jtot);

        free(t);
        free(q);

        Dict retvals;
        retvals["qn"] = qn;
        retvals["tot"] = tot;
        retvals["qm"] = qm;
        retvals["tmt"] = tmt;
        return retvals;
}


Dict Util::Crosrng_sm_psi(EMData* circ1p, EMData* circ2p, vector<int> numr, float psi, int flag) {
// flag 0 - straignt, 1 - mirror

        int nring = numr.size()/3;
        int maxrin = numr[numr.size()-1];
        double qn; float tot; double qm; float tmt;
        float *circ1 = circ1p->get_data();
        float *circ2 = circ2p->get_data();
/*
c
c  checks both straight & mirrored positions
c
c  input - fourier transforms of rings!!
c  circ1 already multiplied by weights!
c
*/

        double *q, t7[7];

        int   ip, jc, numr3i, numr2i, i, j, k, jtot = 0;
        float t1, t2, t3, t4, c1, c2, d1, d2, pos;

        qn  = 0.0f;
        qm  = 0.0f;
        tot = 0.0f;
        tmt = 0.0f;
#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32

        //  c - straight  = circ1 * conjg(circ2)
        //  zero q array

        q = (double*)calloc(maxrin,sizeof(double));

   //   premultiply  arrays ie( circ12 = circ1 * circ2) much slower
        if (flag==0) {
                for (i=1; i<=nring; i++) {

                        numr3i = numr(3,i);   // Number of samples of this ring
                        numr2i = numr(2,i);   // The beginning point of this ring

                        t1   = circ1(numr2i) * circ2(numr2i);
                        q(1) += t1;

                        t1   = circ1(numr2i+1) * circ2(numr2i+1);
                        if (numr3i == maxrin)  {
                                q(2) += t1;
                        } else {
                                q(numr3i+1) += t1;
                        }

                        for (j=3; j<=numr3i; j += 2) {
                                jc     = j+numr2i-1;

        // Here, (c1+c2i)*conj(d1+d2i) = (c1*d1+c2*d2)+(-c1*d2+c2*d1)i
        //                                ----- -----    ----- -----
        //                                 t1     t2      t3    t4

                                c1     = circ1(jc);
                                c2     = circ1(jc+1);
                                d1     = circ2(jc);
                                d2     = circ2(jc+1);

                                t1     = c1 * d1;
                                t3     = c1 * d2;
                                t2     = c2 * d2;
                                t4     = c2 * d1;

                                q(j)   += t1 + t2;
                                q(j+1) += -t3 + t4;
                        }
                }
        } else {
                for (i=1; i<=nring; i++) {

                        numr3i = numr(3,i);   // Number of samples of this ring
                        numr2i = numr(2,i);   // The beginning point of this ring

                        t1   = circ1(numr2i) * circ2(numr2i);
                        q(1) += t1;

                        t1   = circ1(numr2i+1) * circ2(numr2i+1);
                        if (numr3i == maxrin)  {
                                q(2) += t1;
                        } else {
                                q(numr3i+1) += t1;
                        }

                        for (j=3; j<=numr3i; j += 2) {
                                jc     = j+numr2i-1;

        // Here, (c1+c2i)*conj(d1+d2i) = (c1*d1+c2*d2)+(-c1*d2+c2*d1)i
        //                                ----- -----    ----- -----
        //                                 t1     t2      t3    t4

                                c1     = circ1(jc);
                                c2     = circ1(jc+1);
                                d1     = circ2(jc);
                                d2     = circ2(jc+1);

                                t1     = c1 * d1;
                                t3     = c1 * d2;
                                t2     = c2 * d2;
                                t4     = c2 * d1;

                                q(j)   += t1 - t2;
                                q(j+1) += -t3 - t4;
                        }
                }
        }
        fftr_d(q,ip);

        qn  = -1.0e20;
        int psi_pos = int(psi/360.0*maxrin+0.5);

        for (k=-5; k<=5; k++) {
                j = (psi_pos+maxrin-1)%maxrin+1;
                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);
        }

        // interpolate
        prb1d(t7,7,&pos);
        tot = (float)(jtot)+pos;
        free(q);

        Dict retvals;
        retvals["qn"] = qn;
        retvals["tot"] = tot;
        return retvals;
}


Dict Util::Crosrng_ns(EMData* circ1p, EMData* circ2p, vector<int> numr) {
        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];
        double qn; float tot;
        float *circ1 = circ1p->get_data();
        float *circ2 = circ2p->get_data();
/*
c
c  checks only straight position
c
c  input - fourier transforms of rings!!
c  circ1 already multiplied by weights!
c
*/

        // dimension             circ1(lcirc),circ2(lcirc)

        // q(maxrin), t7(-3:3)  //maxrin+2 removed
        double *q, t7[7];

        int   ip, jc, numr3i, numr2i, i, j, k, jtot = 0;
        float c1, c2, d1, d2, pos;

        qn  = 0.0;
        tot = 0.0;
#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
   ip = -(int)(log2(maxrin));
#endif  //_WIN32
        //for (j=1; j<=maxrin;j++) cout <<"  "<<j<<"   "<<circ1(j)<<"   "<<circ2(j) <<endl;

        //  c - straight  = circ1 * conjg(circ2)
        //  zero q array

        q = (double*)calloc(maxrin,sizeof(double));

                        //   premultiply  arrays ie( circ12 = circ1 * circ2) much slower
        for (i=1; i<=nring; i++) {

                numr3i = numr(3,i);   // Number of samples of this ring
                numr2i = numr(2,i);   // The beginning point of this ring

                q(1) += circ1(numr2i) * circ2(numr2i);

                if (numr3i == maxrin)   q(2) += circ1(numr2i+1) * circ2(numr2i+1);
                else  q(numr3i+1) += circ1(numr2i+1) * circ2(numr2i+1);

                for (j=3; j<=numr3i; j += 2) {
                        jc     = j+numr2i-1;

// Here, (c1+c2i)*conj(d1+d2i) = (c1*d1+c2*d2)+(-c1*d2+c2*d1)i
//                                ----- -----    ----- -----
//                                 t1     t2      t3    t4

                        c1     = circ1(jc);
                        c2     = circ1(jc+1);
                        d1     = circ2(jc);
                        d2     = circ2(jc+1);

                        q(j)   += c1 * d1 + c2 * d2;
                        q(j+1) += -c1 * d2 + c2 * d1;
                }
        }
//for (j=1; j<=maxrin; j++) cout <<"  "<<j<<"   "<<q(j) <<endl;
        fftr_d(q,ip);

        qn  = -1.0e20;
        for (j=1; j<=maxrin; j++) {//cout <<"  "<<j<<"   "<<q(j) <<endl;
                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);
        }

        // interpolate
        prb1d(t7,7,&pos);
        tot = (float)(jtot)+pos;
        // Do not interpolate
        //*tot = (float)(jtot);

        free(q);

        Dict retvals;
        retvals["qn"] = qn;
        retvals["tot"] = tot;
        return retvals;
}

#define  dout(i,j)        dout[i+maxrin*j]
#define  circ1b(i)        circ1b[i-1]
#define  circ2b(i)        circ2b[i-1]

EMData* Util::Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr) {
/*
c
c  checks both straight & mirrored positions
c
c  input - fourier transforms of rings!!
c  circ1 already multiplied by weights!
c
  returns EM object with 1D ccf

*/

   // dimension         circ1(lcirc),circ2(lcirc)

        int   ip, jc, numr3i, numr2i, i, j;
        float t1, t2, t3, t4, c1, c2, d1, d2;

        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];

        float* circ1b = circ1->get_data();
        float* circ2b = circ2->get_data();

        // t(maxrin), q(maxrin)  // removed +2
        double *t, *q;

        q = (double*)calloc(maxrin,sizeof(double));
        t = (double*)calloc(maxrin,sizeof(double));

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32

        //  q - straight  = circ1 * conjg(circ2)

        //   t - mirrored  = conjg(circ1) * conjg(circ2)

        //   premultiply  arrays ie( circ12 = circ1 * circ2) much slower

        for (i=1; i<=nring; i++) {

                numr3i = numr(3,i);
                numr2i = numr(2,i);

                t1   = circ1b(numr2i) * circ2b(numr2i);
                q(1) = q(1)+t1;
                t(1) = t(1)+t1;

                t1   = circ1b(numr2i+1) * circ2b(numr2i+1);
                if (numr3i == maxrin)  {
                        q(2) += t1;
                        t(2) += t1;
                } else {
                        q(numr3i+1) += t1;
                        t(numr3i+1) += t1;
                }

                for (j=3; j<=numr3i; j=j+2) {
                        jc     = j+numr2i-1;

                        c1     = circ1b(jc);
                        c2     = circ1b(jc+1);
                        d1     = circ2b(jc);
                        d2     = circ2b(jc+1);

                        t1     = c1 * d1;
                        t3     = c1 * d2;
                        t2     = c2 * d2;
                        t4     = c2 * d1;

                        q(j)   += t1 + t2;
                        q(j+1) += - t3 + t4;
                        t(j)   += t1 - t2;
                        t(j+1) += - t3 - t4;
                }
        }

        // straight
        fftr_d(q,ip);

        // mirrored
        fftr_d(t,ip);

        EMData* out = new EMData();
        out->set_size(maxrin,2,1);
        float *dout = out->get_data();
        for (int i=0; i<maxrin; i++) {dout(i,0)=static_cast<float>(q[i]); dout(i,1)=static_cast<float>(t[i]);}
        //out->set_size(maxrin,1,1);
        //float *dout = out->get_data();
        //for (int i=0; i<maxrin; i++) {dout(i,0)=q[i];}
        free(t);
        free(q);
        return out;
}


vector<float> Util::Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr ) {

        int maxrin = numr[numr.size()-1];

        vector<float> r(2*maxrin);

        Crosrng_msg_vec( circ1, circ2, numr, &r[0], &r[maxrin] );

        return r;
}

#define  dout(i,j)        dout[i+maxrin*j]
#define  circ1b(i)        circ1b[i-1]
#define  circ2b(i)        circ2b[i-1]

void Util::Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t) {
/*
c
c  checks both straight & mirrored positions
c
c  input - fourier transforms of rings!!
c  circ1 already multiplied by weights!
c
  returns EM object with 1D ccf

*/

   // dimension         circ1(lcirc),circ2(lcirc)

        int   ip, jc, numr3i, numr2i, i, j;
        float t1, t2, t3, t4, c1, c2, d1, d2;

        int nring = numr.size()/3;
        //int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
        int maxrin = numr[numr.size()-1];

        float* circ1b = circ1->get_data();
        float* circ2b = circ2->get_data();

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32
        for (int i=1; i<=maxrin; i++)  {q(i) = 0.0f; t(i) = 0.0f;}

        //  q - straight  = circ1 * conjg(circ2)

        //   t - mirrored  = conjg(circ1) * conjg(circ2)

        for (i=1; i<=nring; i++) {

                numr3i = numr(3,i);
                numr2i = numr(2,i);

                t1   = circ1b(numr2i) * circ2b(numr2i);
                q(1) += t1;
                t(1) += t1;

                t1   = circ1b(numr2i+1) * circ2b(numr2i+1);
                if (numr3i == maxrin)  {
                        q(2) += t1;
                        t(2) += t1;
                } else {
                        q(numr3i+1) += t1;
                        t(numr3i+1) += t1;
                }

                for (j=3; j<=numr3i; j=j+2) {
                        jc     = j+numr2i-1;

                        c1     = circ1b(jc);
                        c2     = circ1b(jc+1);
                        d1     = circ2b(jc);
                        d2     = circ2b(jc+1);

                        t1     = c1 * d1;
                        t3     = c1 * d2;
                        t2     = c2 * d2;
                        t4     = c2 * d1;

                        q(j)   += t1 + t2;
                        q(j+1) += -t3 + t4;
                        t(j)   += t1 - t2;
                        t(j+1) += -t3 - t4;
                }
        }
        // straight
        fftr_q(q,ip);
        //for (int i=0; i<maxrin; i++) cout<<i<<"  B    "<<q[i]<<"       "<<t[i]<<endl;

        // mirrored
        fftr_q(t,ip);
}



EMData* Util::Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr)
{
/*
  This program is half of the Crosrng_msg. It only checks straight position.

  input - fourier transforms of rings!!
  circ1 already multiplied by weights!
  returns EM object with 1D ccf

*/

        int   ip, jc, numr3i, numr2i, i, j;
        float t1, t2, t3, t4, c1, c2, d1, d2;

        int nring = numr.size()/3;
        int maxrin = numr[numr.size()-1];

        float* circ1b = circ1->get_data();
        float* circ2b = circ2->get_data();

        double *q;

        q = (double*)calloc(maxrin,sizeof(double));

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32

         //  q - straight  = circ1 * conjg(circ2)

        for (i=1;i<=nring;i++) {

                numr3i = numr(3,i);
                numr2i = numr(2,i);

                t1   = circ1b(numr2i) * circ2b(numr2i);
                q(1) = q(1)+t1;

                if (numr3i == maxrin)  {
                        t1   = circ1b(numr2i+1) * circ2b(numr2i+1);
                        q(2) = q(2)+t1;
                } else {
                        t1              = circ1b(numr2i+1) * circ2b(numr2i+1);
                        q(numr3i+1) = q(numr3i+1)+t1;
                }

                for (j=3;j<=numr3i;j=j+2) {
                        jc     = j+numr2i-1;

                        c1     = circ1b(jc);
                        c2     = circ1b(jc+1);
                        d1     = circ2b(jc);
                        d2     = circ2b(jc+1);

                        t1     = c1 * d1;
                        t3     = c1 * d2;
                        t2     = c2 * d2;
                        t4     = c2 * d1;

                        q(j)   = q(j)   + t1 + t2;
                        q(j+1) = q(j+1) - t3 + t4;
                }
        }

        // straight
        fftr_d(q,ip);

        EMData* out = new EMData();
        out->set_size(maxrin,1,1);
        float *dout = out->get_data();
        for (int i=0; i<maxrin; i++) dout[i]=static_cast<float>(q[i]);
        free(q);
        return out;

}


EMData* Util::Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr)
{
/*

  This program is half of the Crosrng_msg. It only checks mirrored position.

  input - fourier transforms of rings!!
  circ1 already multiplied by weights!
  returns EM object with 1D ccf

*/

        int   ip, jc, numr3i, numr2i, i, j;
        float t1, t2, t3, t4, c1, c2, d1, d2;

        int nring = numr.size()/3;
        int maxrin = numr[numr.size()-1];

        float* circ1b = circ1->get_data();
        float* circ2b = circ2->get_data();

        double *t;

        t = (double*)calloc(maxrin,sizeof(double));

#ifdef _WIN32
        ip = -(int)(log((float)maxrin)/log(2.0f));
#else
        ip = -(int)(log2(maxrin));
#endif  //_WIN32

         //   t - mirrored  = conjg(circ1) * conjg(circ2)

        for (i=1;i<=nring;i++) {

                numr3i = numr(3,i);
                numr2i = numr(2,i);

                t1   = circ1b(numr2i) * circ2b(numr2i);
                t(1) = t(1)+t1;

                if (numr3i == maxrin)  {
                        t1   = circ1b(numr2i+1) * circ2b(numr2i+1);
                        t(2) = t(2)+t1;
                }

                for (j=3;j<=numr3i;j=j+2) {
                        jc     = j+numr2i-1;

                        c1     = circ1b(jc);
                        c2     = circ1b(jc+1);
                        d1     = circ2b(jc);
                        d2     = circ2b(jc+1);

                        t1     = c1 * d1;
                        t3     = c1 * d2;
                        t2     = c2 * d2;
                        t4     = c2 * d1;

                        t(j)   = t(j)   + t1 - t2;
                        t(j+1) = t(j+1) - t3 - t4;
                }
        }

        // mirrored
        fftr_d(t,ip);

        EMData* out = new EMData();
        out->set_size(maxrin,1,1);
        float *dout = out->get_data();
        for (int i=0; i<maxrin; i++) dout[i]=static_cast<float>(t[i]);
        free(t);
        return out;

}

#undef circ1b
#undef circ2b
#undef dout

#undef  circ1
#undef  circ2
#undef  t
#undef  q
#undef  b
#undef  t7


#define    QUADPI                   3.141592653589793238462643383279502884197
#define    PI2                      2*QUADPI

// helper functions for ali2d_ra

float Util::ener(EMData* ave, vector<int> numr) {
        ENTERFUNC;
        long double ener,en;

        int nring = numr.size()/3;
        float *aveptr = ave->get_data();

        ener = 0.0;
        for (int i=1; i<=nring; i++) {
                int numr3i = numr(3,i);
                int np     = numr(2,i)-1;
                float tq = static_cast<float>(PI2*numr(1,i)/numr3i);
                en = tq*(aveptr[np]*aveptr[np]+aveptr[np+1]*aveptr[np+1])*0.5;
                for (int j=np+2; j<np+numr3i-1; j++) en += tq*aveptr[j]*aveptr[j];
                ener += en/numr3i;
        }
        EXITFUNC;
        return static_cast<float>(ener);
}

float Util::ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot) {
        ENTERFUNC;
        long double ener, en;
        float arg, cs, si;

        int nima = data.size();
        int nring = numr.size()/3;
        int maxrin = numr(3,nring);

        ener = 0.0;
        for (int i=1; i<=nring; i++) {
                int numr3i = numr(3,i);
                int np     = numr(2,i)-1;
                float tq = static_cast<float>(PI2*numr(1,i)/numr3i);
                float temp1 = 0.0, temp2 = 0.0;
                for (int kk=0; kk<nima; kk++) {
                        float *ptr = data[kk]->get_data();
                        temp1 += ptr[np];
                        temp2 += static_cast<float>(ptr[np+1]*cos(PI2*(tot[kk]-1.0f)/2.0f*numr3i/maxrin));
                }
                en = tq*(temp1*temp1+temp2*temp2)*0.5;
                for (int j=2; j<numr3i; j+=2) {
                        float tempr = 0.0, tempi = 0.0;
                        for (int kk=0; kk<nima; kk++) {
                                float *ptr = data[kk]->get_data();
                                arg = static_cast<float>( PI2*(tot[kk]-1.0)*(j/2)/maxrin );
                                cs = cos(arg);
                                si = sin(arg);
                                tempr += ptr[np + j]*cs - ptr[np + j +1]*si;
                                tempi += ptr[np + j]*si + ptr[np + j +1]*cs;
                        }
                        en += tq*(tempr*tempr+tempi*tempi);
                }
                ener += en/numr3i;
        }
        EXITFUNC;
        return static_cast<float>(ener);
}

void Util::update_fav (EMData* avep, EMData* datp, float tot, int mirror, vector<int> numr) {
        int nring = numr.size()/3;
        float *ave = avep->get_data();
        float *dat = datp->get_data();
        int i, j, numr3i, np;
        float  arg, cs, si;
        int maxrin = numr(3,nring);
        if(mirror == 1) { //for mirrored data has to be conjugated
                for (i=1; i<=nring; i++) {
                        numr3i = numr(3,i);
                        np     = numr(2,i)-1;
                        ave[np]   += dat[np];
                        ave[np+1] += static_cast<float>( dat[np+1]*cos(PI2*(tot-1.0f)/2.0f*numr3i/maxrin) );
                        for (j=2; j<numr3i; j=j+2) {
                                arg = static_cast<float>( PI2*(tot-1.)*(j/2)/maxrin );
                                cs = cos(arg);
                                si = sin(arg);
                                //complex(data[np + j],data[np + j +1])*complex(cos(arg),sin(arg))
                                ave[np + j]    += dat[np + j]*cs - dat[np + j +1]*si;
                                ave[np + j +1] -= dat[np + j]*si + dat[np + j +1]*cs;
                        }
                }
        } else {
                for (i=1; i<=nring; i++) {
                        numr3i = numr(3,i);
                        np     = numr(2,i)-1;
                        ave[np]   += dat[np];
                        ave[np+1] += static_cast<float>( dat[np+1]*cos(PI2*(tot-1.0f)/2.0f*numr3i/maxrin) );
                        for (j=2; j<numr3i; j=j+2) {
                                arg = static_cast<float>( PI2*(tot-1.)*(j/2)/maxrin );
                                cs = cos(arg);
                                si = sin(arg);
                                //complex(data[np + j],data[np + j +1])*complex(cos(arg),sin(arg))
                                ave[np + j]    += dat[np + j]*cs - dat[np + j +1]*si;
                                ave[np + j +1] += dat[np + j]*si + dat[np + j +1]*cs;
                        }
                }
        }
        avep->update();
        EXITFUNC;
}

void Util::sub_fav(EMData* avep, EMData* datp, float tot, int mirror, vector<int> numr) {
        int nring = numr.size()/3;
        float *ave = avep->get_data();
        float *dat = datp->get_data();
        int i, j, numr3i, np;
        float  arg, cs, si;
        int maxrin = numr(3,nring);
        if(mirror == 1) { //for mirrored data has to be conjugated
                for (i=1; i<=nring; i++) {
                        numr3i = numr(3,i);
                        np     = numr(2,i)-1;
                        ave[np]   -= dat[np];
                        ave[np+1] -= static_cast<float>( dat[np+1]*cos(PI2*(tot-1.0f)/2.0f*numr3i/maxrin) );
                        for (j=2; j<numr3i; j=j+2) {
                                arg = static_cast<float>( PI2*(tot-1.)*(j/2)/maxrin );
                                cs = cos(arg);
                                si = sin(arg);
                                //complex(data[np + j],data[np + j +1])*complex(cos(arg),sin(arg))
                                ave[np + j]    -= dat[np + j]*cs - dat[np + j +1]*si;
                                ave[np + j +1] += dat[np + j]*si + dat[np + j +1]*cs;
                        }
                }
        } else {
                for (i=1; i<=nring; i++) {
                        numr3i = numr(3,i);
                        np     = numr(2,i)-1;
                        ave[np]   -= dat[np];
                        ave[np+1] -= static_cast<float>( dat[np+1]*cos(PI2*(tot-1.0f)/2.0f*numr3i/maxrin) );
                        for (j=2; j<numr3i; j=j+2) {
                                arg = static_cast<float>( PI2*(tot-1.)*(j/2)/maxrin );
                                cs = cos(arg);
                                si = sin(arg);
                                //complex(data[np + j],data[np + j +1])*complex(cos(arg),sin(arg))
                                ave[np + j]    -= dat[np + j]*cs - dat[np + j +1]*si;
                                ave[np + j +1] -= dat[np + j]*si + dat[np + j +1]*cs;
                        }
                }
        }
        avep->update();
        EXITFUNC;
}


#undef    QUADPI
#undef    PI2

#undef  numr
#undef  circ


#define QUADPI   3.141592653589793238462643383279502884197
#define PI2      QUADPI*2
#define deg_rad  QUADPI/180.0
#define rad_deg  180.0/QUADPI

struct ori_t
{
    int iphi;
    int itht;
    int id;
};


struct cmpang
{
    bool operator()( const ori_t& a, const ori_t& b )
    {
        if( a.itht != b.itht )
        {
            return a.itht < b.itht;
        }

        return a.iphi < b.iphi;
    }
};

//helper function for the weights calculation by Voronoi to Cml
vector<double> Util::cml_weights(const vector<float>& cml){
        static const int NBIN = 100;
        int nline=cml.size()/2;
        vector<double> weights(nline);

        vector<ori_t> angs(nline);
        for( int i=0; i < nline; ++i ) {
                angs[i].iphi = int( NBIN*cml[2*i] );
                angs[i].itht = int( NBIN*cml[2*i+1] );
                if( angs[i].itht == 180*NBIN ) angs[i].itht = 0;
                angs[i].id = i;
        }

        //std::cout << "# of angs: " << angs.size() << std::endl;

        std::sort( angs.begin(), angs.end(), cmpang() );

        vector<float> newphi;
        vector<float> newtht;
        vector< vector<int> > indices;

        int curt_iphi = -1;
        int curt_itht = -1;
        for(unsigned int i=0 ;i < angs.size(); ++i ) {
                if( angs[i].iphi==curt_iphi && angs[i].itht==curt_itht ) {
                        Assert( indices.size() > 0 );
                        indices.back().push_back(angs[i].id);
                } else {
                        curt_iphi = angs[i].iphi;
                        curt_itht = angs[i].itht;

                        newphi.push_back( float(curt_iphi)/NBIN );
                        newtht.push_back( float(curt_itht)/NBIN );
                        indices.push_back( vector<int>(1,angs[i].id) );
                }
        }

        //std::cout << "# of indpendent ang: " << newphi.size() << std::endl;


        int num_agl = newphi.size();

        if(num_agl>2) {
                vector<double> w=Util::vrdg(newphi, newtht);

                Assert( w.size()==newphi.size() );
                Assert( indices.size()==newphi.size() );

                for(unsigned int i=0; i < newphi.size(); ++i ) {
                    /*
                    std::cout << "phi,tht,w,n: ";
                    std::cout << boost::format( "%10.3f" ) % newphi[i] << " ";
                    std::cout << boost::format( "%10.3f" ) % newtht[i] << " ";
                    std::cout << boost::format( "%8.6f"  ) % w[i] << " ";
                    std::cout << indices[i].size() << "(";
                    */

                    for(unsigned int j=0; j < indices[i].size(); ++j ) {
                            int id = indices[i][j];
                            weights[id] = w[i]/indices[i].size();
                            //std::cout << id << " ";
                    }

                    //std::cout << ")" << std::endl;

                }
        } else {
                cout<<"warning in Util.cml_weights"<<endl;
                double val = PI2/float(nline);
                for(int i=0; i<nline; i++)  weights[i]=val;
        }

        return weights;

}

/****************************************************
 * New code for common-lines
 ****************************************************/

/*  This function drop a line (line) to an 2D image (img).
 *  The position of the line to the image is defined by (postline).
 *  The part of the line paste is defined by (offset), the begin position
 *  and (length) the size.
 */
void Util::set_line(