(file) Return to imageinterp.c CVS log (file) (dir) Up to [Development] / JSOC / proj / globalhs / apps

  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 tplarson 1.1  *
 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              //#include "astro.h"
 40              
 41              typedef struct LIBASTRO_Dist_Struct 
 42              {
 43 tplarson 1.1   int disttype;
 44                double scale; /* Image scale relative to FD. 5 for vw. */
 45                double xc, yc; /* Nominal image center in pixels. */
 46                double cdist; /* Cubic distortion constant for FD images */
 47                double feff, alpha, cosbeta, sinbeta; /* Tilt constants for FD. */
 48              } LIBASTRO_Dist_t;
 49              
 50              typedef enum 
 51              {
 52                 kLIBASTRO_InterBilinear = 0,
 53                 kLIBASTRO_InterCubic = 1,
 54              } LIBASTRO_Interpolation_t;
 55              
 56              typedef enum
 57              {
 58                 kLIBASTRO_Success = 0,
 59                 kLIBASTRO_BadDimensionality,
 60                 kLIBASTRO_BadData,
 61                 kLIBASTRO_CouldntCreateData,
 62                 kLIBASTRO_DimensionMismatch,
 63                 kLIBASTRO_CantDoOddNumLats,
 64 tplarson 1.1    kLIBASTRO_UnsupportedInterp,
 65                 kLIBASTRO_UnsupportedMCOR,
 66                 kLIBASTRO_UnsupportedVCOR,
 67                 kLIBASTRO_InconsistentConstraints,
 68                 kLIBASTRO_InvalidArgs,
 69                 kLIBASTRO_Interpolation,
 70                 kLIBASTRO_InsufficientData
 71              } LIBASTRO_Error_t;
 72              
 73              const int kMaxCols = 8192;
 74              const double kOmegaCarr = (360 /  27.27527); /* degrees/day - synodic Carrington rotation rate */
 75              const double kRadsPerDegree = M_PI/180.;
 76              
 77              static void Ccker(double *u, double s);
 78              static double Ccintd(double *f, int nx, int ny, double x, double y);
 79              static double Linintd(double *f, int nx, int ny, double x, double y);
 80              static void Distort(double x, 
 81              		    double y, 
 82              		    double rsun, 
 83              		    int sign, 
 84              		    double *xd, 
 85 tplarson 1.1 		    double *yd, 
 86              		    const LIBASTRO_Dist_t *distpars);
 87              
 88              /* Calculate the interpolation kernel. */
 89              void Ccker(double *u, double s)
 90              {
 91                 double s2, s3;
 92                   
 93                 s2= s * s;
 94                 s3= s2 * s;
 95                 u[0] = s2 - 0.5 * (s3 + s);
 96                 u[1] = 1.5*s3 - 2.5*s2 + 1.0;
 97                 u[2] = -1.5*s3 + 2.0*s2 + 0.5*s;
 98                 u[3] = 0.5 * (s3 - s2);
 99              }
100              
101              /* Cubic convolution interpolation */
102              /*
103               *  Cubic convolution interpolation, based on GONG software, received Apr-88
104               *  Interpolation by cubic convolution in 2 dimensions with 32-bit double data
105               *
106 tplarson 1.1  *  Reference:
107               *    R.G. Keys in IEEE Trans. Acoustics, Speech, and Signal Processing, 
108               *    Vol. 29, 1981, pp. 1153-1160.
109               *
110               *  Usage:
111               *    double ccintd (double *f, int nx, int ny, double x, double y)
112               *
113               *  Bugs:
114               *    Currently interpolation within one pixel of edge is not
115               *      implemented.  If x or y is in this range or off the picture, the
116               *      function value returned is zero.
117               *    Only works for double data.
118               */
119              double Ccintd(double *f, int nx, int ny, double x, double y)
120              {
121                 double  ux[4], uy[4];
122                 /* Sum changed to double to speed up things */
123                 double sum;
124              
125                 int ix, iy, ix1, iy1, i, j;
126                   
127 tplarson 1.1    if (x < 1. || x >= (double)(nx-2) ||
128                     y < 1. || y >= (double)(ny-2))
129                   return 0.0;
130                   
131                 ix = (int)x;
132                 iy = (int)y;
133                 Ccker(ux,  x - (double)ix);
134                 Ccker(uy,  y - (double)iy);
135                   
136                 ix1 = ix - 1;
137                 iy1 = iy - 1;
138                 sum = 0.;
139                 for (i = 0; i < 4; i++)
140                   for (j = 0; j < 4; j++)
141                     sum = sum + f[(iy1+i) * nx + ix1 + j] * uy[i] * ux[j];
142                 return sum;
143              }
144              
145              /*
146               *  Bilinear interpolation, based on pipeLNU remap.
147               *
148 tplarson 1.1  *  Reference: 
149               *     Abramowitz, Milton, and Stegun, Irene, "Handbook of Mathematical
150               *     Functions," Dover Publications, New York, 1964.
151               *
152               *  Usage:
153               *    double linintd (double *f, int nx, int ny, double x, double y)
154               *
155               *  Bugs:
156               *    Only works for double data.
157               */
158              
159              double Linintd(double *f, int nx, int ny, double x, double y) 
160              {
161                 double p0, p1, q0, q1;          /* weights                              */
162                 int ilow, jlow;   /* selected CCD pixels                  */
163                 double *fptr;                    /* input array temp                     */
164                 double u;                        
165                   
166                 ilow = (int)floor(x);
167                 jlow = (int)floor(y);
168                 if (ilow < 0) ilow = 0;
169 tplarson 1.1    if (ilow+1 >= nx) ilow -= 1;
170                 if (jlow < 0) jlow = 0;
171                 if (jlow+1 >= ny) jlow -= 1;
172                 p1 = x - ilow;
173                 p0 = 1.0 - p1;
174                 q1 = y - jlow;
175                 q0 = 1.0 - q1;
176                   
177                 /* Bilinear interpolation from four data points. */
178                 u = 0.0;
179                 fptr = f + jlow*nx + ilow;
180                 u += p0 * q0 * *fptr++;
181                 u += p1 * q0 * *fptr;
182                 fptr += nx - 1;
183                 u += p0 * q1 * *fptr++;
184                 u += p1 * q1 * *fptr;
185                 return u;
186              }
187              
188              
189              int SetDistort(int dist, 
190 tplarson 1.1 	       double cubic, 
191              	       double alpha, 
192              	       double beta, 
193              	       double feff,
194               	       LIBASTRO_Dist_t *dOut)
195              {
196                 int error = 0;
197              
198                 if (dOut != NULL)
199                 {
200                    int status = 0;
201              
202                    dOut->disttype = dist;
203              
204                    /* disttype == 0 is 'no distortion' and is not an error */
205                    if (dOut->disttype != 0)
206                    {
207              	 if (dOut->disttype == 1) 
208              	 { 
209              	    /* FD */
210              	    dOut->xc = 511.5;
211 tplarson 1.1 	    dOut->yc = 511.5;
212              	    dOut->scale = 1.0;
213              	 }
214              	 else if (dOut->disttype == 2)
215              	 { 
216              	    /* vw */
217              	    dOut->xc = 95.1;
218              	    dOut->yc = 95.1;
219              	    dOut->scale = 5.0;
220              	 }
221              	 else 
222              	 {
223              	    error = 1;
224              	 }
225              
226              	 if (!error)
227              	 {
228              	    dOut->feff = feff;
229              	    dOut->alpha = alpha * kRadsPerDegree;
230              	    beta *= kRadsPerDegree;
231              	    dOut->cosbeta = cos(beta);
232 tplarson 1.1 	    dOut->sinbeta = sin(beta);
233              	    dOut->cdist = cubic;
234              
235              	    if (status != CMDPARAMS_SUCCESS)
236              	    {
237              	       error = 1;
238              	    }
239              
240              	 }
241              
242                    } /* disttype != 0*/
243                 }
244                 else
245                 {
246                    error = 1;
247                 }
248              
249                 return error;
250              }
251              
252              
253 tplarson 1.1 void Distort(double x, 
254              	     double y, 
255              	     double rsun, 
256              	     int sign, 
257              	     double *xd, 
258              	     double *yd, 
259              	     const LIBASTRO_Dist_t *distpars)
260              {
261                 /*
262                   For sign=+1 transforms ideal coordinates (x,y) to distorted
263                   coordinates (xd,yd). rsun needed to correct for scale error.
264                   For sign=-1 ideal and distorted coordinates are swapped.
265                   Units are pixels from lower left corner.
266                   Adapted from /home/schou/calib/dist/distort.pro
267                   Note that rsun must be passed in FD pixels!
268                 */
269                   
270                 double xx,yy,xxl,yyl,x0,y0,x00,y00,x0l,y0l,epsx,epsy,rdist;
271                   
272                 if (distpars->disttype == 0) 
273                 {
274 tplarson 1.1       *xd=x;
275                    *yd=y;
276                    return;
277                 }
278                 /* Convert fo FD pixel scale*/
279              
280                   rsun=rsun*distpars->scale;
281              
282                 /* Don't do anyway since rsun already corrected. Don't get me started... */
283                 /* above comment no longer applies, rsun is now untouched in o2helio (v2helio) */  
284                   
285                 /* Shift zero point to nominal image center and convert to FD pixel scale */
286                   
287                 x00=distpars->scale*(x-distpars->xc);
288                 y00=distpars->scale*(y-distpars->yc);
289                   
290                 /* Radial distortion */
291                   
292                 rdist=distpars->cdist*(x00*x00+y00*y00-rsun*rsun);
293                   
294                 x0=x00+rdist*x00;
295 tplarson 1.1    y0=y00+rdist*y00;
296                   
297                 /* Distortion due to CCD tilt */
298                   
299                 /* First rotate coordinate system */
300                 x0l=x0*distpars->cosbeta+y0*distpars->sinbeta;
301                 y0l=-x0*distpars->sinbeta+y0*distpars->cosbeta;
302                   
303                 /* Apply distortion */
304                 xxl=x0l*(1.-distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
305                 yyl=y0l*(1.+distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha) - 
306                   rsun*rsun/distpars->feff*distpars->alpha;
307                   
308                 /* Rotate coordinates back */
309                 xx=xxl*distpars->cosbeta-yyl*distpars->sinbeta;
310                 yy=xxl*distpars->sinbeta+yyl*distpars->cosbeta;
311                   
312                 /* Total distortion */
313                   
314                 epsx=xx-x00;
315                 epsy=yy-y00;
316 tplarson 1.1      
317                 /* Return distorted values. Cheating with inverse transform. */
318                   
319                 *xd=x+sign*epsx/distpars->scale;
320                 *yd=y+sign*epsy/distpars->scale;
321              }
322              
323              
324              int imageinterp(
325                      float *V,  /* input velocity array    */
326                      /* assumed declared V[ypixels][xpixels]  */
327                      float *U,  /* output projected array   */
328                      
329 tplarson 1.2         int xpixels, /* x width of input array */
330                      int ypixels, /* y width of input array */
331                      double   x0, /* x pixel address of disk center		*/
332                      double   y0, /* y pixel address of disk center		*/
333                      double   P,  /* angle between CCD y-axis and solar vertical */
334                                   /* positive to the east (CCW) */
335 tplarson 1.1         double   rsun,  /* pixels */
336 tplarson 1.3         double   Rmax, /* maximum disk radius to use (e.g. 0.95)	*/
337                      int      NaN_beyond_rmax,
338 tplarson 1.1         int interpolation, /* option */
339                      int cols,  /* width of output array   */
340                      int rows,  /* height of output array   */
341                      double   vsign,
342                      const LIBASTRO_Dist_t *distpars)
343              {
344                   
345                 float u;   /* output temp    */
346                 double x, y;   /* CCD location of desired point */
347                 int row, col;   /* index into output array  */
348                 int rowoffset;
349 tplarson 1.3    double xratio, yratio, Rmax2;
350 tplarson 1.1 
351                 long i;
352                 double *vdp;
353                 float *vp;
354 tplarson 1.2    double xtmp, ytmp;
355 tplarson 1.1      
356                 double *VD;
357                 VD=(double *)(malloc(xpixels*ypixels*sizeof(double)));
358                 vdp=VD;
359                 vp=V;
360                 for (i=0;i<xpixels*ypixels;i++) *vdp++=(double)*vp++;
361 tplarson 1.3 
362                 Rmax2 = Rmax*Rmax*rsun*rsun;
363 tplarson 1.1      
364                 if (cols > kMaxCols) 
365                 {
366                    return kLIBASTRO_DimensionMismatch;
367                 }
368              
369                   
370                 if (interpolation != kLIBASTRO_InterCubic && interpolation != kLIBASTRO_InterBilinear) 
371                 {     
372                    return kLIBASTRO_UnsupportedInterp;
373                 }
374              
375                 xratio=(double)xpixels/cols;
376                 yratio=(double)ypixels/rows;
377                   
378                 for (row = 0; row < rows; row++) {
379               
380                    rowoffset = row*cols;
381                
382                    for (col = 0; col < cols; col++) {
383              
384 tplarson 1.3          xtmp = (col + 0.5)*xratio - 0.5 - x0;
385                       ytmp = (row + 0.5)*yratio - 0.5 - y0;
386                       if ((xtmp*xtmp + ytmp*ytmp) >= Rmax2)
387                       {
388                         *(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
389                         continue;
390                       }
391                       x= x0 + xtmp*cos(P) - ytmp*sin(P);
392                       y= y0 + xtmp*sin(P) + ytmp*cos(P);
393 tplarson 1.1 
394                       Distort(x, y, rsun, +1, &x, &y, distpars);
395              
396                       if (interpolation == kLIBASTRO_InterCubic) 
397                       {
398                          u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
399                       }
400                       else if (interpolation == kLIBASTRO_InterBilinear) 
401                       {
402                          u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
403                       }
404              
405                       *(U + rowoffset + col) = u;
406                    }
407                 }
408                 free(VD);
409                 return kLIBASTRO_Success;
410              }

Karen Tian
Powered by
ViewCVS 0.9.4