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: imageinterp.c,v $";

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;
}
Revision:	1.4
Committed:	Fri May 3 20:47:15 2013 UTC (10 years, 4 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	globalhs_version_5, Ver_8-7, Ver_8-5, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_8-3, globalhs_version_8, globalhs_version_9, globalhs_version_0, globalhs_version_1, globalhs_version_2, globalhs_version_3, globalhs_version_4, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-2, Ver_8-10, Ver_8-1, Ver_8-6, Ver_9-1, Ver_8-4, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.3:	+12 -42 lines
Log Message:	start using projection.h
#	User	Rev	Content
1	tplarson	1.1
2			/*
3			* imageinterp (modified from obs2helio)
4			*
5			* Purpose:
6			* Interpolate CCD data to estimate value that would
7			* be obtained at specified equal incrememts of heliographic
8			* longitude and sine of latitude.
9			*
10			* Description:
11			*
12			* Employs bi-linear or cubic convolution interpolation to map a
13			* portion of the solar disk to a rectangular array "cols" by "rows" in
14			* width and height respectively.
15			*
16			* Reference:
17			*
18			* Green, Robin M. "Spherical Astronomy," Cambridge University
19			* Press, Cambridge, 1985.
20			*
21			* Responsible: Kay Leibrand KLeibrand@solar.Stanford.EDU
22			*
23			* Restrictions:
24			* The number of output map rows must be even.
25			* The number of output map columns must be less than MAXCOLS.
26			* Assumes orientation is a 4 character string with one of the following
27			* 8 legal values: SESW SENE SWSE SWNW NWNE NWSW NENW NESE
28			*
29			* Bugs:
30			* Possible division by 0
31			*
32			* Planned updates:
33			* Optimize.
34			*
35			* Revision history is at end of file.
36			*/
37
38			#include <math.h>
39	tplarson	1.4	#include "projection.h"
40	tplarson	1.1
41	tplarson	1.4	char *cvsinfo_imageinterp = "cvsinfo: $Header: imageinterp.c,v $";
42	tplarson	1.1
43			const int kMaxCols = 8192;
44			const double kOmegaCarr = (360 / 27.27527); /* degrees/day - synodic Carrington rotation rate */
45			const double kRadsPerDegree = M_PI/180.;
46
47			static void Ccker(double *u, double s);
48			static double Ccintd(double *f, int nx, int ny, double x, double y);
49			static double Linintd(double *f, int nx, int ny, double x, double y);
50			static void Distort(double x,
51			double y,
52			double rsun,
53			int sign,
54			double *xd,
55			double *yd,
56	tplarson	1.4	const LIBPROJECTION_Dist_t *distpars);
57	tplarson	1.1
58			/* Calculate the interpolation kernel. */
59			void Ccker(double *u, double s)
60			{
61			double s2, s3;
62
63			s2= s * s;
64			s3= s2 * s;
65			u[0] = s2 - 0.5 * (s3 + s);
66			u[1] = 1.5s3 - 2.5s2 + 1.0;
67			u[2] = -1.5s3 + 2.0s2 + 0.5*s;
68			u[3] = 0.5 * (s3 - s2);
69			}
70
71			/* Cubic convolution interpolation */
72			/*
73			* Cubic convolution interpolation, based on GONG software, received Apr-88
74			* Interpolation by cubic convolution in 2 dimensions with 32-bit double data
75			*
76			* Reference:
77			* R.G. Keys in IEEE Trans. Acoustics, Speech, and Signal Processing,
78			* Vol. 29, 1981, pp. 1153-1160.
79			*
80			* Usage:
81			* double ccintd (double *f, int nx, int ny, double x, double y)
82			*
83			* Bugs:
84			* Currently interpolation within one pixel of edge is not
85			* implemented. If x or y is in this range or off the picture, the
86			* function value returned is zero.
87			* Only works for double data.
88			*/
89			double Ccintd(double *f, int nx, int ny, double x, double y)
90			{
91			double ux[4], uy[4];
92			/* Sum changed to double to speed up things */
93			double sum;
94
95			int ix, iy, ix1, iy1, i, j;
96
97			if (x < 1. \|\| x >= (double)(nx-2) \|\|
98			y < 1. \|\| y >= (double)(ny-2))
99			return 0.0;
100
101			ix = (int)x;
102			iy = (int)y;
103			Ccker(ux, x - (double)ix);
104			Ccker(uy, y - (double)iy);
105
106			ix1 = ix - 1;
107			iy1 = iy - 1;
108			sum = 0.;
109			for (i = 0; i < 4; i++)
110			for (j = 0; j < 4; j++)
111			sum = sum + f[(iy1+i) * nx + ix1 + j] * uy[i] * ux[j];
112			return sum;
113			}
114
115			/*
116			* Bilinear interpolation, based on pipeLNU remap.
117			*
118			* Reference:
119			* Abramowitz, Milton, and Stegun, Irene, "Handbook of Mathematical
120			* Functions," Dover Publications, New York, 1964.
121			*
122			* Usage:
123			* double linintd (double *f, int nx, int ny, double x, double y)
124			*
125			* Bugs:
126			* Only works for double data.
127			*/
128
129			double Linintd(double *f, int nx, int ny, double x, double y)
130			{
131			double p0, p1, q0, q1; /* weights */
132			int ilow, jlow; /* selected CCD pixels */
133			double fptr; / input array temp */
134			double u;
135
136			ilow = (int)floor(x);
137			jlow = (int)floor(y);
138			if (ilow < 0) ilow = 0;
139			if (ilow+1 >= nx) ilow -= 1;
140			if (jlow < 0) jlow = 0;
141			if (jlow+1 >= ny) jlow -= 1;
142			p1 = x - ilow;
143			p0 = 1.0 - p1;
144			q1 = y - jlow;
145			q0 = 1.0 - q1;
146
147			/* Bilinear interpolation from four data points. */
148			u = 0.0;
149			fptr = f + jlow*nx + ilow;
150			u += p0 * q0 * *fptr++;
151			u += p1 * q0 * *fptr;
152			fptr += nx - 1;
153			u += p0 * q1 * *fptr++;
154			u += p1 * q1 * *fptr;
155			return u;
156			}
157
158
159			int SetDistort(int dist,
160			double cubic,
161			double alpha,
162			double beta,
163			double feff,
164	tplarson	1.4	LIBPROJECTION_Dist_t *dOut)
165	tplarson	1.1	{
166			int error = 0;
167
168			if (dOut != NULL)
169			{
170			int status = 0;
171
172			dOut->disttype = dist;
173
174			/* disttype == 0 is 'no distortion' and is not an error */
175			if (dOut->disttype != 0)
176			{
177			if (dOut->disttype == 1)
178			{
179			/* FD */
180			dOut->xc = 511.5;
181			dOut->yc = 511.5;
182			dOut->scale = 1.0;
183			}
184			else if (dOut->disttype == 2)
185			{
186			/* vw */
187			dOut->xc = 95.1;
188			dOut->yc = 95.1;
189			dOut->scale = 5.0;
190			}
191			else
192			{
193			error = 1;
194			}
195
196			if (!error)
197			{
198			dOut->feff = feff;
199			dOut->alpha = alpha * kRadsPerDegree;
200			beta *= kRadsPerDegree;
201			dOut->cosbeta = cos(beta);
202			dOut->sinbeta = sin(beta);
203			dOut->cdist = cubic;
204
205			if (status != CMDPARAMS_SUCCESS)
206			{
207			error = 1;
208			}
209
210			}
211
212			} /* disttype != 0*/
213			}
214			else
215			{
216			error = 1;
217			}
218
219			return error;
220			}
221
222
223			void Distort(double x,
224			double y,
225			double rsun,
226			int sign,
227			double *xd,
228			double *yd,
229	tplarson	1.4	const LIBPROJECTION_Dist_t *distpars)
230	tplarson	1.1	{
231			/*
232			For sign=+1 transforms ideal coordinates (x,y) to distorted
233			coordinates (xd,yd). rsun needed to correct for scale error.
234			For sign=-1 ideal and distorted coordinates are swapped.
235			Units are pixels from lower left corner.
236			Adapted from /home/schou/calib/dist/distort.pro
237			Note that rsun must be passed in FD pixels!
238			*/
239
240			double xx,yy,xxl,yyl,x0,y0,x00,y00,x0l,y0l,epsx,epsy,rdist;
241
242			if (distpars->disttype == 0)
243			{
244			*xd=x;
245			*yd=y;
246			return;
247			}
248			/* Convert fo FD pixel scale*/
249
250			rsun=rsun*distpars->scale;
251
252			/* Don't do anyway since rsun already corrected. Don't get me started... */
253			/* above comment no longer applies, rsun is now untouched in o2helio (v2helio) */
254
255			/* Shift zero point to nominal image center and convert to FD pixel scale */
256
257			x00=distpars->scale*(x-distpars->xc);
258			y00=distpars->scale*(y-distpars->yc);
259
260			/* Radial distortion */
261
262			rdist=distpars->cdist(x00x00+y00y00-rsunrsun);
263
264			x0=x00+rdist*x00;
265			y0=y00+rdist*y00;
266
267			/* Distortion due to CCD tilt */
268
269			/* First rotate coordinate system */
270			x0l=x0distpars->cosbeta+y0distpars->sinbeta;
271			y0l=-x0distpars->sinbeta+y0distpars->cosbeta;
272
273			/* Apply distortion */
274			xxl=x0l(1.-distpars->alphadistpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
275			yyl=y0l(1.+distpars->alphadistpars->alpha/4.+y0l/distpars->feff*distpars->alpha) -
276			rsunrsun/distpars->feffdistpars->alpha;
277
278			/* Rotate coordinates back */
279			xx=xxldistpars->cosbeta-yyldistpars->sinbeta;
280			yy=xxldistpars->sinbeta+yyldistpars->cosbeta;
281
282			/* Total distortion */
283
284			epsx=xx-x00;
285			epsy=yy-y00;
286
287			/* Return distorted values. Cheating with inverse transform. */
288
289			xd=x+signepsx/distpars->scale;
290			yd=y+signepsy/distpars->scale;
291			}
292
293
294			int imageinterp(
295			float V, / input velocity array */
296			/* assumed declared V[ypixels][xpixels] */
297			float U, / output projected array */
298
299	tplarson	1.2	int xpixels, /* x width of input array */
300			int ypixels, /* y width of input array */
301			double x0, /* x pixel address of disk center */
302			double y0, /* y pixel address of disk center */
303			double P, /* angle between CCD y-axis and solar vertical */
304			/* positive to the east (CCW) */
305	tplarson	1.1	double rsun, /* pixels */
306	tplarson	1.3	double Rmax, /* maximum disk radius to use (e.g. 0.95) */
307			int NaN_beyond_rmax,
308	tplarson	1.1	int interpolation, /* option */
309			int cols, /* width of output array */
310			int rows, /* height of output array */
311			double vsign,
312	tplarson	1.4	const LIBPROJECTION_Dist_t *distpars)
313	tplarson	1.1	{
314
315			float u; /* output temp */
316			double x, y; /* CCD location of desired point */
317			int row, col; /* index into output array */
318			int rowoffset;
319	tplarson	1.3	double xratio, yratio, Rmax2;
320	tplarson	1.1
321			long i;
322			double *vdp;
323			float *vp;
324	tplarson	1.2	double xtmp, ytmp;
325	tplarson	1.1
326			double *VD;
327			VD=(double )(malloc(xpixelsypixels*sizeof(double)));
328			vdp=VD;
329			vp=V;
330			for (i=0;i<xpixelsypixels;i++) vdp++=(double)*vp++;
331	tplarson	1.3
332			Rmax2 = RmaxRmaxrsun*rsun;
333	tplarson	1.1
334			if (cols > kMaxCols)
335			{
336	tplarson	1.4	return kLIBPROJECTION_DimensionMismatch;
337	tplarson	1.1	}
338
339
340	tplarson	1.4	if (interpolation != kLIBPROJECTION_InterCubic && interpolation != kLIBPROJECTION_InterBilinear)
341	tplarson	1.1	{
342	tplarson	1.4	return kLIBPROJECTION_UnsupportedInterp;
343	tplarson	1.1	}
344
345			xratio=(double)xpixels/cols;
346			yratio=(double)ypixels/rows;
347
348			for (row = 0; row < rows; row++) {
349
350			rowoffset = row*cols;
351
352			for (col = 0; col < cols; col++) {
353
354	tplarson	1.3	xtmp = (col + 0.5)*xratio - 0.5 - x0;
355			ytmp = (row + 0.5)*yratio - 0.5 - y0;
356			if ((xtmpxtmp + ytmpytmp) >= Rmax2)
357			{
358			*(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
359			continue;
360			}
361			x= x0 + xtmpcos(P) - ytmpsin(P);
362			y= y0 + xtmpsin(P) + ytmpcos(P);
363	tplarson	1.1
364			Distort(x, y, rsun, +1, &x, &y, distpars);
365
366	tplarson	1.4	if (interpolation == kLIBPROJECTION_InterCubic)
367	tplarson	1.1	{
368			u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
369			}
370	tplarson	1.4	else if (interpolation == kLIBPROJECTION_InterBilinear)
371	tplarson	1.1	{
372			u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
373			}
374
375			*(U + rowoffset + col) = u;
376			}
377			}
378			free(VD);
379	tplarson	1.4	return kLIBPROJECTION_Success;
380	tplarson	1.1	}