Kstars

imageautoguiding.cpp
1/*
2 SPDX-FileCopyrightText: 2017 Bob Majewski <cygnus257@yahoo.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5*/
6
7#include "imageautoguiding.h"
8
9#include <qglobal.h>
10
11#include <cmath>
12
13#define SWAP(a, b) \
14 tempr = (a); \
15 (a) = (b); \
16 (b) = tempr
17#define NR_END 1
18#define FREE_ARG char *
19
20#define TWOPI 6.28318530717959
21#define FFITMAX 0.05
22
23//void ImageAutoGuiding1(float *ref,float *im,int n,float *xshift,float *yshift);
24
25float ***f3tensorSP(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
26void free_f3tensorSP(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
27
28float **matrixSP(long nrl, long nrh, long ncl, long nch);
29
30void free_matrixSP(float **m, long nrl, long nrh, long ncl, long nch);
31void nrerrorNR(void);
32
33void fournNR(float data[], unsigned long nn[], long ndim, long isign);
34
35void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, long isign);
36
37void ShiftEST(float ***testimage, float ***refimage, int n, float *xshift, float *yshift, int k);
38
39namespace ImageAutoGuiding
40{
41void ImageAutoGuiding1(float *ref, float *im, int n, float *xshift, float *yshift)
42{
43 float ***RefImage, ***TestImage;
44 int i, j, k;
45 float x, y;
46
47 /* Allocate memory */
48
49 RefImage = f3tensorSP(1, 1, 1, n, 1, n);
50 TestImage = f3tensorSP(1, 1, 1, n, 1, n);
51
52 /* Load Data */
53
54 k = 0;
55
56 for (j = 1; j <= n; ++j)
57 {
58 for (i = 1; i <= n; ++i)
59 {
60 RefImage[1][j][i] = ref[k];
61 TestImage[1][j][i] = im[k];
62
63 k += 1;
64 }
65 }
66
67 /* Calculate Image Shifts */
68
69 ShiftEST(TestImage, RefImage, n, &x, &y, 1);
70
71 *xshift = x;
72 *yshift = y;
73
74 /* free memory */
75
76 free_f3tensorSP(RefImage, 1, 1, 1, n, 1, n);
77 free_f3tensorSP(TestImage, 1, 1, 1, n, 1, n);
78}
79}
80
81// Calculates Image Shifts
82
83void ShiftEST(float ***testimage, float ***refimage, int n, float *xshift, float *yshift, int k)
84{
85 int ix, iy, nh, nhplusone;
87 double fx, fy, ff, fn, re, im, testre, testim, rev, imv, phi;
88 double power, dem, f2, f2limit;
89 float **speq;
90
91 f2limit = FFITMAX * FFITMAX;
92
93 speq = matrixSP(1, 1, 1, 2 * n);
94
95 nh = n / 2;
96 nhplusone = nh + 1;
97
98 fn = ((float)n);
99 ff = 1.0 / fn;
100
101 /* FFT of Reference */
102
103 for (ix = 1; ix <= 2 * n; ++ix)
104 {
105 speq[1][ix] = 0.0;
106 }
107
108 rlft3NR(refimage, speq, 1, n, n, 1);
109
110 /* FFT of Test Image */
111
112 for (ix = 1; ix <= 2 * n; ++ix)
113 {
114 speq[1][ix] = 0.0;
115 }
116
117 rlft3NR(testimage, speq, 1, n, n, 1);
118
119 /* Solving for slopes */
120
121 fx2sum = 0.0;
122 fy2sum = 0.0;
123
124 phifxsum = 0.0;
125 phifysum = 0.0;
126
127 fxfysum = 0.0;
128
129 for (ix = 1; ix <= n; ++ix)
130 {
131 if (ix <= nhplusone)
132 {
133 fx = ff * ((float)(ix - 1));
134 }
135 else
136 {
137 fx = -ff * ((float)(n + 1 - ix));
138 }
139
140 for (iy = 1; iy <= nh; ++iy)
141 {
142 fy = ff * ((float)(iy - 1));
143
144 f2 = fx * fx + fy * fy;
145
146 /* Limit to Low Spatial Frequencies */
147
148 if (f2 < f2limit)
149 {
150 /* Real part reference image */
151
152 re = refimage[k][ix][2 * iy - 1];
153
154 /* Imaginary part reference image */
155
156 im = refimage[k][ix][2 * iy];
157
158 power = re * re + im * im;
159
160 /* Real part test image */
161
162 testre = testimage[k][ix][2 * iy - 1];
163
164 /* Imaginary part image */
165
166 testim = testimage[k][ix][2 * iy];
167
168 rev = re * testre + im * testim;
169 imv = re * testim - im * testre;
170
171 /* Find Phase */
172
173 phi = atan2(imv, rev);
174
175 fx2sum += power * fx * fx;
176 fy2sum += power * fy * fy;
177
178 phifxsum += power * fx * phi;
179 phifysum += power * fy * phi;
180
181 fxfysum += power * fx * fy;
182 }
183 }
184 }
185
186 /* calculate subpixel shift */
187
189
190 deltax = (phifxsum * fy2sum - fxfysum * phifysum) / (dem * TWOPI);
191 deltay = (phifysum * fx2sum - fxfysum * phifxsum) / (dem * TWOPI);
192
193 free_matrixSP(speq, 1, 1, 1, 2 * n);
194
195 /* You can change the shift mapping here */
196
197 *xshift = deltax;
198 *yshift = deltay;
199}
200
201// 2 and 3 Dimensional FFT Routine for Real Data
202// Numerical Recipes in C Second Edition
203// The Art of Scientific Computing
204// 1999
205
206void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, long isign)
207{
208 void fournNR(float data[], unsigned long nn[], long ndim, long isign);
209 void nrerror();
210 unsigned long i1, i2, i3, j1, j2, j3, nn[4], ii3;
211 double theta, wpi, wpr, wtemp;
212 float c1, c2, h1r, h1i, h2r, h2i;
213 float wi, wr;
214
215 if ((unsigned long)(1 + &data[nn1][nn2][nn3] - &data[1][1][1]) != nn1 * nn2 * nn3)
216 nrerrorNR();
217 c1 = 0.5;
218 c2 = -0.5 * isign;
219 theta = isign * (6.28318530717959 / nn3);
220 wtemp = sin(0.5 * theta);
221 wpr = -2.0 * wtemp * wtemp;
222 wpi = sin(theta);
223 nn[1] = nn1;
224 nn[2] = nn2;
225 nn[3] = nn3 >> 1;
226 if (isign == 1)
227 {
228 fournNR(&data[1][1][1] - 1, nn, 3, isign);
229 for (i1 = 1; i1 <= nn1; i1++)
230 for (i2 = 1, j2 = 0; i2 <= nn2; i2++)
231 {
232 speq[i1][++j2] = data[i1][i2][1];
233 speq[i1][++j2] = data[i1][i2][2];
234 }
235 }
236 for (i1 = 1; i1 <= nn1; i1++)
237 {
238 j1 = (i1 != 1 ? nn1 - i1 + 2 : 1);
239 wr = 1.0;
240 wi = 0.0;
241 for (ii3 = 1, i3 = 1; i3 <= (nn3 >> 2) + 1; i3++, ii3 += 2)
242 {
243 for (i2 = 1; i2 <= nn2; i2++)
244 {
245 if (i3 == 1)
246 {
247 j2 = (i2 != 1 ? ((nn2 - i2) << 1) + 3 : 1);
248 h1r = c1 * (data[i1][i2][1] + speq[j1][j2]);
249 h1i = c1 * (data[i1][i2][2] - speq[j1][j2 + 1]);
250 h2i = c2 * (data[i1][i2][1] - speq[j1][j2]);
251 h2r = -c2 * (data[i1][i2][2] + speq[j1][j2 + 1]);
252 data[i1][i2][1] = h1r + h2r;
253 data[i1][i2][2] = h1i + h2i;
254 speq[j1][j2] = h1r - h2r;
255 speq[j1][j2 + 1] = h2i - h1i;
256 }
257 else
258 {
259 j2 = (i2 != 1 ? nn2 - i2 + 2 : 1);
260 j3 = nn3 + 3 - (i3 << 1);
261 h1r = c1 * (data[i1][i2][ii3] + data[j1][j2][j3]);
262 h1i = c1 * (data[i1][i2][ii3 + 1] - data[j1][j2][j3 + 1]);
263 h2i = c2 * (data[i1][i2][ii3] - data[j1][j2][j3]);
264 h2r = -c2 * (data[i1][i2][ii3 + 1] + data[j1][j2][j3 + 1]);
265 data[i1][i2][ii3] = h1r + wr * h2r - wi * h2i;
266 data[i1][i2][ii3 + 1] = h1i + wr * h2i + wi * h2r;
267 data[j1][j2][j3] = h1r - wr * h2r + wi * h2i;
268 data[j1][j2][j3 + 1] = -h1i + wr * h2i + wi * h2r;
269 }
270 }
271 wr = (wtemp = wr) * wpr - wi * wpi + wr;
272 wi = wi * wpr + wtemp * wpi + wi;
273 }
274 }
275 if (isign == -1)
276 fournNR(&data[1][1][1] - 1, nn, 3, isign);
277}
278
279// Multi-Dimensional FFT Routine for Complex Data
280// Numerical Recipes in C Second Edition
281// The Art of Scientific Computing
282// 1999
283// Used in rlft3
284
285void fournNR(float data[], unsigned long nn[], long ndim, long isign)
286{
287 long idim;
288 unsigned long i1, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
289 unsigned long i2, i3;
290 unsigned long ibit, k1 = 0, k2, n, nprev, nrem, ntot;
291 float wi, wr, tempi, tempr;
292 double theta, wpi, wpr, wtemp;
293
294 for (ntot = 1, idim = 1; idim <= ndim; idim++)
295 ntot *= nn[idim];
296 nprev = 1;
297 for (idim = ndim; idim >= 1; idim--)
298 {
299 n = nn[idim];
300 nrem = ntot / (n * nprev);
301 ip1 = nprev << 1;
302 ip2 = ip1 * n;
303 ip3 = ip2 * nrem;
304 i2rev = 1;
305 for (i2 = 1; i2 <= ip2; i2 += ip1)
306 {
307 if (i2 < i2rev)
308 {
309 for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2)
310 {
311 for (i3 = i1; i3 <= ip3; i3 += ip2)
312 {
313 i3rev = i2rev + i3 - i2;
314 SWAP(data[i3], data[i3rev]);
315 SWAP(data[i3 + 1], data[i3rev + 1]);
316 }
317 }
318 }
319 ibit = ip2 >> 1;
320 while (ibit >= ip1 && i2rev > ibit)
321 {
322 i2rev -= ibit;
323 ibit >>= 1;
324 }
325 i2rev += ibit;
326 }
327 ifp1 = ip1;
328 while (ifp1 < ip2)
329 {
330 ifp2 = ifp1 << 1;
331 theta = isign * 6.28318530717959 / (ifp2 / ip1);
332 wtemp = sin(0.5 * theta);
333 wpr = -2.0 * wtemp * wtemp;
334 wpi = sin(theta);
335 wr = 1.0;
336 wi = 0.0;
337 for (i3 = 1; i3 <= ifp1; i3 += ip1)
338 {
339 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
340 {
341 for (i2 = i1; i2 <= ip3; i2 += ifp2)
342 {
343 k1 = i2;
344 k2 = k1 + ifp1;
345 tempr = (double)wr * data[k2] - (double)wi * data[k2 + 1];
346 tempi = (double)wr * data[k2 + 1] + (double)wi * data[k2];
347 data[k2] = data[k1] - tempr;
348 data[k2 + 1] = data[k1 + 1] - tempi;
349 data[k1] += tempr;
350 data[k1 + 1] += tempi;
351 }
352 }
353 wr = (wtemp = wr) * wpr - wi * wpi + wr;
354 wi = wi * wpr + wtemp * wpi + wi;
355 }
356 ifp1 = ifp2;
357 }
358 nprev *= n;
359 }
360}
361
362void nrerrorNR(void)
363{
364}
365
366// Numerical Recipes memory allocation routines
367// One based arrays
368
369float **matrixSP(long nrl, long nrh, long ncl, long nch)
370/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
371{
372 long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
373 float **m;
374
375 /* allocate pointers to rows */
376 m = (float **)malloc((size_t)((nrow + NR_END) * sizeof(float *)));
377 if (!m)
378 nrerrorNR();
379 m += NR_END;
380 m -= nrl;
381
382 /* allocate rows and set pointers to them */
383 m[nrl] = (float *)malloc((size_t)((nrow * ncol + NR_END) * sizeof(float)));
384 if (!m[nrl])
385 nrerrorNR();
386 m[nrl] += NR_END;
387 m[nrl] -= ncl;
388
389 for (i = nrl + 1; i <= nrh; i++)
390 m[i] = m[i - 1] + ncol;
391
392 /* return pointer to array of pointers to rows */
393 return m;
394}
395
396float ***f3tensorSP(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
397/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
398{
399 long i, j, nrow = nrh - nrl + 1, ncol = nch - ncl + 1, ndep = ndh - ndl + 1;
400 float ***t;
401
402 /* allocate pointers to pointers to rows */
403 t = (float ***)malloc((size_t)((nrow + NR_END) * sizeof(float **)));
404 if (!t)
405 nrerrorNR();
406 t += NR_END;
407 t -= nrl;
408
409 /* allocate pointers to rows and set pointers to them */
410 t[nrl] = (float **)malloc((size_t)((nrow * ncol + NR_END) * sizeof(float *)));
411 if (!t[nrl])
412 nrerrorNR();
413 t[nrl] += NR_END;
414 t[nrl] -= ncl;
415
416 /* allocate rows and set pointers to them */
417 t[nrl][ncl] = (float *)malloc((size_t)((nrow * ncol * ndep + NR_END) * sizeof(float)));
418 if (!t[nrl][ncl])
419 nrerrorNR();
420 t[nrl][ncl] += NR_END;
421 t[nrl][ncl] -= ndl;
422
423 for (j = ncl + 1; j <= nch; j++)
424 t[nrl][j] = t[nrl][j - 1] + ndep;
425 for (i = nrl + 1; i <= nrh; i++)
426 {
427 t[i] = t[i - 1] + ncol;
428 t[i][ncl] = t[i - 1][ncl] + ncol * ndep;
429 for (j = ncl + 1; j <= nch; j++)
430 t[i][j] = t[i][j - 1] + ndep;
431 }
432
433 /* return pointer to array of pointers to rows */
434 return t;
435}
436
437void free_matrixSP(float **m, long nrl, long nrh, long ncl, long nch)
438/* free a float matrix allocated by matrix() */
439{
440 Q_UNUSED(nrh);
441 Q_UNUSED(nch);
442 free((FREE_ARG)(m[nrl] + ncl - NR_END));
443 free((FREE_ARG)(m + nrl - NR_END));
444}
445
446void free_f3tensorSP(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
447/* free a float f3tensor allocated by f3tensor() */
448{
449 Q_UNUSED(nrh);
450 Q_UNUSED(nch);
451 Q_UNUSED(ndh);
452 free((FREE_ARG)(t[nrl][ncl] + ndl - NR_END));
453 free((FREE_ARG)(t[nrl] + ncl - NR_END));
454 free((FREE_ARG)(t + nrl - NR_END));
455}
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri May 17 2024 11:48:26 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.