proj/globalhs/apps/imageinterp.c

/*
 *  imageinterp (modified from obs2helio)                       
 *
 *  Purpose:
 *     Interpolate CCD data to estimate value that would
 *     be obtained at specified equal incrememts of heliographic
 *     longitude and sine of latitude.
 *
 *  Description:
 *
 *     Employs bi-linear or cubic convolution interpolation to map a
 *     portion of the solar disk to a rectangular array "cols" by "rows" in
 *     width and height respectively.
 *
 *  Reference:
 *
 *     Green, Robin M. "Spherical Astronomy," Cambridge University
 *     Press, Cambridge, 1985.
 *
 *  Responsible:  Kay Leibrand                   KLeibrand@solar.Stanford.EDU
 *
 *  Restrictions:
 *     The number of output map rows must be even.
 *     The number of output map columns must be less than MAXCOLS.
 *     Assumes orientation is a 4 character string with one of the following
 *     8 legal values: SESW SENE SWSE SWNW NWNE NWSW NENW NESE 
 *
 *  Bugs:
 *     Possible division by 0 
 *
 *  Planned updates:
 *     Optimize.
 *
 *  Revision history is at end of file.
 */

#include <math.h>
#include "projection.h"

char *cvsinfo_imageinterp = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/imageinterp.c,v 1.4 2013/05/03 20:47:15 tplarson Exp $";

const int kMaxCols = 8192;
const double kOmegaCarr = (360 /  27.27527); /* degrees/day - synodic Carrington rotation rate */
const double kRadsPerDegree = M_PI/180.;

static void Ccker(double *u, double s);
static double Ccintd(double *f, int nx, int ny, double x, double y);
static double Linintd(double *f, int nx, int ny, double x, double y);
static void Distort(double x, 
            double y, 
            double rsun, 
            int sign, 
            double *xd, 
            double *yd, 
            const LIBPROJECTION_Dist_t *distpars);

/* Calculate the interpolation kernel. */
void Ccker(double *u, double s)
{
   double s2, s3;
     
   s2= s * s;
   s3= s2 * s;
   u[0] = s2 - 0.5 * (s3 + s);
   u[1] = 1.5*s3 - 2.5*s2 + 1.0;
   u[2] = -1.5*s3 + 2.0*s2 + 0.5*s;
   u[3] = 0.5 * (s3 - s2);
}

/* Cubic convolution interpolation */
/*
 *  Cubic convolution interpolation, based on GONG software, received Apr-88
 *  Interpolation by cubic convolution in 2 dimensions with 32-bit double data
 *
 *  Reference:
 *    R.G. Keys in IEEE Trans. Acoustics, Speech, and Signal Processing, 
 *    Vol. 29, 1981, pp. 1153-1160.
 *
 *  Usage:
 *    double ccintd (double *f, int nx, int ny, double x, double y)
 *
 *  Bugs:
 *    Currently interpolation within one pixel of edge is not
 *      implemented.  If x or y is in this range or off the picture, the
 *      function value returned is zero.
 *    Only works for double data.
 */
double Ccintd(double *f, int nx, int ny, double x, double y)
{
   double  ux[4], uy[4];
   /* Sum changed to double to speed up things */
   double sum;

   int ix, iy, ix1, iy1, i, j;
     
   if (x < 1. || x >= (double)(nx-2) ||
       y < 1. || y >= (double)(ny-2))
     return 0.0;
     
   ix = (int)x;
   iy = (int)y;
   Ccker(ux,  x - (double)ix);
   Ccker(uy,  y - (double)iy);
     
   ix1 = ix - 1;
   iy1 = iy - 1;
   sum = 0.;
   for (i = 0; i < 4; i++)
     for (j = 0; j < 4; j++)
       sum = sum + f[(iy1+i) * nx + ix1 + j] * uy[i] * ux[j];
   return sum;
}

/*
 *  Bilinear interpolation, based on pipeLNU remap.
 *
 *  Reference: 
 *     Abramowitz, Milton, and Stegun, Irene, "Handbook of Mathematical
 *     Functions," Dover Publications, New York, 1964.
 *
 *  Usage:
 *    double linintd (double *f, int nx, int ny, double x, double y)
 *
 *  Bugs:
 *    Only works for double data.
 */

double Linintd(double *f, int nx, int ny, double x, double y) 
{
   double p0, p1, q0, q1;          /* weights                              */
   int ilow, jlow;   /* selected CCD pixels                  */
   double *fptr;                    /* input array temp                     */
   double u;                        
     
   ilow = (int)floor(x);
   jlow = (int)floor(y);
   if (ilow < 0) ilow = 0;
   if (ilow+1 >= nx) ilow -= 1;
   if (jlow < 0) jlow = 0;
   if (jlow+1 >= ny) jlow -= 1;
   p1 = x - ilow;
   p0 = 1.0 - p1;
   q1 = y - jlow;
   q0 = 1.0 - q1;
     
   /* Bilinear interpolation from four data points. */
   u = 0.0;
   fptr = f + jlow*nx + ilow;
   u += p0 * q0 * *fptr++;
   u += p1 * q0 * *fptr;
   fptr += nx - 1;
   u += p0 * q1 * *fptr++;
   u += p1 * q1 * *fptr;
   return u;
}


int SetDistort(int dist, 
           double cubic, 
           double alpha, 
           double beta, 
           double feff,
           LIBPROJECTION_Dist_t *dOut)
{
   int error = 0;

   if (dOut != NULL)
   {
      int status = 0;

      dOut->disttype = dist;

      /* disttype == 0 is 'no distortion' and is not an error */
      if (dOut->disttype != 0)
      {
     if (dOut->disttype == 1) 
     { 
        /* FD */
        dOut->xc = 511.5;
        dOut->yc = 511.5;
        dOut->scale = 1.0;
     }
     else if (dOut->disttype == 2)
     { 
        /* vw */
        dOut->xc = 95.1;
        dOut->yc = 95.1;
        dOut->scale = 5.0;
     }
     else 
     {
        error = 1;
     }

     if (!error)
     {
        dOut->feff = feff;
        dOut->alpha = alpha * kRadsPerDegree;
        beta *= kRadsPerDegree;
        dOut->cosbeta = cos(beta);
        dOut->sinbeta = sin(beta);
        dOut->cdist = cubic;

        if (status != CMDPARAMS_SUCCESS)
        {
           error = 1;
        }

     }

      } /* disttype != 0*/
   }
   else
   {
      error = 1;
   }

   return error;
}


void Distort(double x, 
         double y, 
         double rsun, 
         int sign, 
         double *xd, 
         double *yd, 
         const LIBPROJECTION_Dist_t *distpars)
{
   /*
     For sign=+1 transforms ideal coordinates (x,y) to distorted
     coordinates (xd,yd). rsun needed to correct for scale error.
     For sign=-1 ideal and distorted coordinates are swapped.
     Units are pixels from lower left corner.
     Adapted from /home/schou/calib/dist/distort.pro
     Note that rsun must be passed in FD pixels!
   */
     
   double xx,yy,xxl,yyl,x0,y0,x00,y00,x0l,y0l,epsx,epsy,rdist;
     
   if (distpars->disttype == 0) 
   {
      *xd=x;
      *yd=y;
      return;
   }
   /* Convert fo FD pixel scale*/

     rsun=rsun*distpars->scale;

   /* Don't do anyway since rsun already corrected. Don't get me started... */
   /* above comment no longer applies, rsun is now untouched in o2helio (v2helio) */  
     
   /* Shift zero point to nominal image center and convert to FD pixel scale */
     
   x00=distpars->scale*(x-distpars->xc);
   y00=distpars->scale*(y-distpars->yc);
     
   /* Radial distortion */
     
   rdist=distpars->cdist*(x00*x00+y00*y00-rsun*rsun);
     
   x0=x00+rdist*x00;
   y0=y00+rdist*y00;
     
   /* Distortion due to CCD tilt */
     
   /* First rotate coordinate system */
   x0l=x0*distpars->cosbeta+y0*distpars->sinbeta;
   y0l=-x0*distpars->sinbeta+y0*distpars->cosbeta;
     
   /* Apply distortion */
   xxl=x0l*(1.-distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
   yyl=y0l*(1.+distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha) - 
     rsun*rsun/distpars->feff*distpars->alpha;
     
   /* Rotate coordinates back */
   xx=xxl*distpars->cosbeta-yyl*distpars->sinbeta;
   yy=xxl*distpars->sinbeta+yyl*distpars->cosbeta;
     
   /* Total distortion */
     
   epsx=xx-x00;
   epsy=yy-y00;
     
   /* Return distorted values. Cheating with inverse transform. */
     
   *xd=x+sign*epsx/distpars->scale;
   *yd=y+sign*epsy/distpars->scale;
}


int imageinterp(
        float *V,  /* input velocity array    */
        /* assumed declared V[ypixels][xpixels]  */
        float *U,  /* output projected array   */
        
        int xpixels, /* x width of input array */
        int ypixels, /* y width of input array */
        double   x0, /* x pixel address of disk center      */
        double   y0, /* y pixel address of disk center      */
        double   P,  /* angle between CCD y-axis and solar vertical */
                     /* positive to the east (CCW) */
        double   rsun,  /* pixels */
        double   Rmax, /* maximum disk radius to use (e.g. 0.95)    */
        int      NaN_beyond_rmax,
        int interpolation, /* option */
        int cols,  /* width of output array   */
        int rows,  /* height of output array   */
        double   vsign,
        const LIBPROJECTION_Dist_t *distpars)
{
     
   float u;   /* output temp    */
   double x, y;   /* CCD location of desired point */
   int row, col;   /* index into output array  */
   int rowoffset;
   double xratio, yratio, Rmax2;

   long i;
   double *vdp;
   float *vp;
   double xtmp, ytmp;
     
   double *VD;
   VD=(double *)(malloc(xpixels*ypixels*sizeof(double)));
   vdp=VD;
   vp=V;
   for (i=0;i<xpixels*ypixels;i++) *vdp++=(double)*vp++;

   Rmax2 = Rmax*Rmax*rsun*rsun;
     
   if (cols > kMaxCols) 
   {
      return kLIBPROJECTION_DimensionMismatch;
   }

     
   if (interpolation != kLIBPROJECTION_InterCubic && interpolation != kLIBPROJECTION_InterBilinear) 
   {     
      return kLIBPROJECTION_UnsupportedInterp;
   }

   xratio=(double)xpixels/cols;
   yratio=(double)ypixels/rows;
     
   for (row = 0; row < rows; row++) {
 
      rowoffset = row*cols;
  
      for (col = 0; col < cols; col++) {

         xtmp = (col + 0.5)*xratio - 0.5 - x0;
         ytmp = (row + 0.5)*yratio - 0.5 - y0;
         if ((xtmp*xtmp + ytmp*ytmp) >= Rmax2)
         {
           *(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
           continue;
         }
         x= x0 + xtmp*cos(P) - ytmp*sin(P);
         y= y0 + xtmp*sin(P) + ytmp*cos(P);

         Distort(x, y, rsun, +1, &x, &y, distpars);

         if (interpolation == kLIBPROJECTION_InterCubic) 
         {
            u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
         }
         else if (interpolation == kLIBPROJECTION_InterBilinear) 
         {
            u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
         }

         *(U + rowoffset + col) = u;
      }
   }
   free(VD);
   return kLIBPROJECTION_Success;
}

---