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 int xpixels, /* x width of input array */
330 int ypixels, /* y width of input array */
331 double rsun, /* pixels */
332
333 int interpolation, /* option */
334 int cols, /* width of output array */
335 int rows, /* height of output array */
336 double vsign,
337 tplarson 1.1 const LIBASTRO_Dist_t *distpars)
338 {
339
340 float u; /* output temp */
341 double x, y; /* CCD location of desired point */
342 int row, col; /* index into output array */
343 int rowoffset;
344 double xratio, yratio;
345
346 long i;
347 double *vdp;
348 float *vp;
349
350 double *VD;
351 VD=(double *)(malloc(xpixels*ypixels*sizeof(double)));
352 vdp=VD;
353 vp=V;
354 for (i=0;i<xpixels*ypixels;i++) *vdp++=(double)*vp++;
355
356 if (cols > kMaxCols)
357 {
358 tplarson 1.1 return kLIBASTRO_DimensionMismatch;
359 }
360
361
362 if (interpolation != kLIBASTRO_InterCubic && interpolation != kLIBASTRO_InterBilinear)
363 {
364 return kLIBASTRO_UnsupportedInterp;
365 }
366
367 xratio=(double)xpixels/cols;
368 yratio=(double)ypixels/rows;
369
370 for (row = 0; row < rows; row++) {
371
372 rowoffset = row*cols;
373
374 for (col = 0; col < cols; col++) {
375
376 x=col*xratio;
377 y=row*yratio;
378
379 tplarson 1.1 Distort(x, y, rsun, +1, &x, &y, distpars);
380
381 if (interpolation == kLIBASTRO_InterCubic)
382 {
383 u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
384 }
385 else if (interpolation == kLIBASTRO_InterBilinear)
386 {
387 u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
388 }
389
390 *(U + rowoffset + col) = u;
391 }
392 }
393 free(VD);
394 return kLIBASTRO_Success;
395 }
|
Return to
imageinterp.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps