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 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.5*s3 - 2.5*s2 + 1.0;
67 u[2] = -1.5*s3 + 2.0*s2 + 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 tplarson 1.1 * 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 tplarson 1.1 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 tplarson 1.1 * 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 /* 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1
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*(x00*x00+y00*y00-rsun*rsun);
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=x0*distpars->cosbeta+y0*distpars->sinbeta;
271 y0l=-x0*distpars->sinbeta+y0*distpars->cosbeta;
272 tplarson 1.1
273 /* Apply distortion */
274 xxl=x0l*(1.-distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
275 yyl=y0l*(1.+distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha) -
276 rsun*rsun/distpars->feff*distpars->alpha;
277
278 /* Rotate coordinates back */
279 xx=xxl*distpars->cosbeta-yyl*distpars->sinbeta;
280 yy=xxl*distpars->sinbeta+yyl*distpars->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+sign*epsx/distpars->scale;
290 *yd=y+sign*epsy/distpars->scale;
291 }
292
293 tplarson 1.1
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(xpixels*ypixels*sizeof(double)));
328 vdp=VD;
329 vp=V;
330 for (i=0;i<xpixels*ypixels;i++) *vdp++=(double)*vp++;
|
331 tplarson 1.3
332 Rmax2 = Rmax*Rmax*rsun*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 ((xtmp*xtmp + ytmp*ytmp) >= Rmax2)
357 {
358 *(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
359 continue;
360 }
361 x= x0 + xtmp*cos(P) - ytmp*sin(P);
362 y= y0 + xtmp*sin(P) + ytmp*cos(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 }
|
Return to
imageinterp.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps