21 #define FREE_ARG char *
23 #define TWOPI 6.28318530717959
28 float ***
f3tensorSP(
long nrl,
long nrh,
long ncl,
long nch,
long ndl,
long ndh);
29 void free_f3tensorSP(
float ***t,
long nrl,
long nrh,
long ncl,
long nch,
long ndl,
long ndh);
31 float **
matrixSP(
long nrl,
long nrh,
long ncl,
long nch);
33 void free_matrixSP(
float **m,
long nrl,
long nrh,
long ncl,
long nch);
36 void fournNR(
float data[],
unsigned long nn[],
long ndim,
long isign);
38 void rlft3NR(
float ***data,
float **speq,
unsigned long nn1,
unsigned long nn2,
unsigned long nn3,
long isign);
40 void ShiftEST(
float ***testimage,
float ***refimage,
int n,
float *xshift,
float *yshift,
int k);
42 namespace ImageAutoGuiding
46 float ***RefImage, ***TestImage;
59 for (j = 1; j <= n; ++j)
61 for (i = 1; i <= n; ++i)
63 RefImage[1][j][i] = ref[k];
64 TestImage[1][j][i] = im[k];
72 ShiftEST(TestImage, RefImage, n, &x, &y, 1);
86 void ShiftEST(
float ***testimage,
float ***refimage,
int n,
float *xshift,
float *yshift,
int k)
88 int ix, iy, nh, nhplusone;
89 double deltax, deltay, fx2sum, fy2sum, phifxsum, phifysum, fxfysum;
90 double fx, fy, ff, fn, re, im, testre, testim, rev, imv, phi;
91 double power, dem, f2, f2limit;
106 for (ix = 1; ix <= 2 * n; ++ix)
111 rlft3NR(refimage, speq, 1, n, n, 1);
115 for (ix = 1; ix <= 2 * n; ++ix)
120 rlft3NR(testimage, speq, 1, n, n, 1);
132 for (ix = 1; ix <= n; ++ix)
136 fx = ff * ((float)(ix - 1));
140 fx = -ff * ((float)(n + 1 - ix));
143 for (iy = 1; iy <= nh; ++iy)
145 fy = ff * ((float)(iy - 1));
147 f2 = fx * fx + fy * fy;
155 re = refimage[k][ix][2 * iy - 1];
159 im = refimage[k][ix][2 * iy];
161 power = re * re + im * im;
165 testre = testimage[k][ix][2 * iy - 1];
169 testim = testimage[k][ix][2 * iy];
171 rev = re * testre + im * testim;
172 imv = re * testim - im * testre;
176 phi = atan2(imv, rev);
178 fx2sum += power * fx * fx;
179 fy2sum += power * fy * fy;
181 phifxsum += power * fx * phi;
182 phifysum += power * fy * phi;
184 fxfysum += power * fx * fy;
191 dem = fx2sum * fy2sum - fxfysum * fxfysum;
193 deltax = (phifxsum * fy2sum - fxfysum * phifysum) / (dem *
TWOPI);
194 deltay = (phifysum * fx2sum - fxfysum * phifxsum) / (dem * TWOPI);
209 void rlft3NR(
float ***data,
float **speq,
unsigned long nn1,
unsigned long nn2,
unsigned long nn3,
long isign)
211 void fournNR(
float data[],
unsigned long nn[],
long ndim,
long isign);
213 unsigned long i1, i2, i3, j1, j2, j3, nn[4], ii3;
214 double theta, wpi, wpr, wtemp;
215 float c1, c2, h1r, h1i, h2r, h2i;
218 if ((
unsigned long)(1 + &data[nn1][nn2][nn3] - &data[1][1][1]) != nn1 * nn2 * nn3)
222 theta = isign * (6.28318530717959 / nn3);
223 wtemp = sin(0.5 * theta);
224 wpr = -2.0 * wtemp * wtemp;
231 fournNR(&data[1][1][1] - 1, nn, 3, isign);
232 for (i1 = 1; i1 <= nn1; i1++)
233 for (i2 = 1, j2 = 0; i2 <= nn2; i2++)
235 speq[i1][++j2] = data[i1][i2][1];
236 speq[i1][++j2] = data[i1][i2][2];
239 for (i1 = 1; i1 <= nn1; i1++)
241 j1 = (i1 != 1 ? nn1 - i1 + 2 : 1);
244 for (ii3 = 1, i3 = 1; i3 <= (nn3 >> 2) + 1; i3++, ii3 += 2)
246 for (i2 = 1; i2 <= nn2; i2++)
250 j2 = (i2 != 1 ? ((nn2 - i2) << 1) + 3 : 1);
251 h1r = c1 * (data[i1][i2][1] + speq[j1][j2]);
252 h1i = c1 * (data[i1][i2][2] - speq[j1][j2 + 1]);
253 h2i = c2 * (data[i1][i2][1] - speq[j1][j2]);
254 h2r = -c2 * (data[i1][i2][2] + speq[j1][j2 + 1]);
255 data[i1][i2][1] = h1r + h2r;
256 data[i1][i2][2] = h1i + h2i;
257 speq[j1][j2] = h1r - h2r;
258 speq[j1][j2 + 1] = h2i - h1i;
262 j2 = (i2 != 1 ? nn2 - i2 + 2 : 1);
263 j3 = nn3 + 3 - (i3 << 1);
264 h1r = c1 * (data[i1][i2][ii3] + data[j1][j2][j3]);
265 h1i = c1 * (data[i1][i2][ii3 + 1] - data[j1][j2][j3 + 1]);
266 h2i = c2 * (data[i1][i2][ii3] - data[j1][j2][j3]);
267 h2r = -c2 * (data[i1][i2][ii3 + 1] + data[j1][j2][j3 + 1]);
268 data[i1][i2][ii3] = h1r + wr * h2r - wi * h2i;
269 data[i1][i2][ii3 + 1] = h1i + wr * h2i + wi * h2r;
270 data[j1][j2][j3] = h1r - wr * h2r + wi * h2i;
271 data[j1][j2][j3 + 1] = -h1i + wr * h2i + wi * h2r;
274 wr = (wtemp = wr) * wpr - wi * wpi + wr;
275 wi = wi * wpr + wtemp * wpi + wi;
279 fournNR(&data[1][1][1] - 1, nn, 3, isign);
288 void fournNR(
float data[],
unsigned long nn[],
long ndim,
long isign)
291 unsigned long i1, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
292 unsigned long i2, i3;
293 unsigned long ibit, k1 = 0, k2, n, nprev, nrem, ntot;
294 float wi, wr, tempi, tempr;
295 double theta, wpi, wpr, wtemp;
297 for (ntot = 1, idim = 1; idim <= ndim; idim++)
300 for (idim = ndim; idim >= 1; idim--)
303 nrem = ntot / (n * nprev);
308 for (i2 = 1; i2 <= ip2; i2 += ip1)
312 for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2)
314 for (i3 = i1; i3 <= ip3; i3 += ip2)
316 i3rev = i2rev + i3 - i2;
317 SWAP(data[i3], data[i3rev]);
318 SWAP(data[i3 + 1], data[i3rev + 1]);
323 while (ibit >= ip1 && i2rev > ibit)
334 theta = isign * 6.28318530717959 / (ifp2 / ip1);
335 wtemp = sin(0.5 * theta);
336 wpr = -2.0 * wtemp * wtemp;
340 for (i3 = 1; i3 <= ifp1; i3 += ip1)
342 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
344 for (i2 = i1; i2 <= ip3; i2 += ifp2)
348 tempr = (double)wr * data[k2] - (
double)wi * data[k2 + 1];
349 tempi = (double)wr * data[k2 + 1] + (
double)wi * data[k2];
350 data[k2] = data[k1] - tempr;
351 data[k2 + 1] = data[k1 + 1] - tempi;
353 data[k1 + 1] += tempi;
356 wr = (wtemp = wr) * wpr - wi * wpi + wr;
357 wi = wi * wpr + wtemp * wpi + wi;
372 float **
matrixSP(
long nrl,
long nrh,
long ncl,
long nch)
375 long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
379 m = (
float **)malloc((
size_t)((nrow +
NR_END) *
sizeof(
float *)));
386 m[nrl] = (
float *)malloc((
size_t)((nrow * ncol +
NR_END) *
sizeof(
float)));
392 for (i = nrl + 1; i <= nrh; i++)
393 m[i] = m[i - 1] + ncol;
399 float ***
f3tensorSP(
long nrl,
long nrh,
long ncl,
long nch,
long ndl,
long ndh)
402 long i, j, nrow = nrh - nrl + 1, ncol = nch - ncl + 1, ndep = ndh - ndl + 1;
406 t = (
float ***)malloc((
size_t)((nrow +
NR_END) *
sizeof(
float **)));
413 t[nrl] = (
float **)malloc((
size_t)((nrow * ncol +
NR_END) *
sizeof(
float *)));
420 t[nrl][ncl] = (
float *)malloc((
size_t)((nrow * ncol * ndep +
NR_END) *
sizeof(
float)));
426 for (j = ncl + 1; j <= nch; j++)
427 t[nrl][j] = t[nrl][j - 1] + ndep;
428 for (i = nrl + 1; i <= nrh; i++)
430 t[i] = t[i - 1] + ncol;
431 t[i][ncl] = t[i - 1][ncl] + ncol * ndep;
432 for (j = ncl + 1; j <= nch; j++)
433 t[i][j] = t[i][j - 1] + ndep;
449 void free_f3tensorSP(
float ***t,
long nrl,
long nrh,
long ncl,
long nch,
long ndl,
long ndh)
void free_f3tensorSP(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, long isign)
void free_matrixSP(float **m, long nrl, long nrh, long ncl, long nch)
float ** matrixSP(long nrl, long nrh, long ncl, long nch)
void fournNR(float data[], unsigned long nn[], long ndim, long isign)
void ShiftEST(float ***testimage, float ***refimage, int n, float *xshift, float *yshift, int k)
float *** f3tensorSP(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
void ImageAutoGuiding1(float *ref, float *im, int n, float *xshift, float *yshift)