Kstars

imageautoguiding.cpp
1 /*
2  SPDX-FileCopyrightText: 2017 Bob Majewski <[email protected]>
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 
25 float ***f3tensorSP(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
26 void free_f3tensorSP(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
27 
28 float **matrixSP(long nrl, long nrh, long ncl, long nch);
29 
30 void free_matrixSP(float **m, long nrl, long nrh, long ncl, long nch);
31 void nrerrorNR(void);
32 
33 void fournNR(float data[], unsigned long nn[], long ndim, long isign);
34 
35 void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, long isign);
36 
37 void ShiftEST(float ***testimage, float ***refimage, int n, float *xshift, float *yshift, int k);
38 
39 namespace ImageAutoGuiding
40 {
41 void 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 
83 void ShiftEST(float ***testimage, float ***refimage, int n, float *xshift, float *yshift, int k)
84 {
85  int ix, iy, nh, nhplusone;
86  double deltax, deltay, fx2sum, fy2sum, phifxsum, phifysum, fxfysum;
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 
188  dem = fx2sum * fy2sum - fxfysum * fxfysum;
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 
206 void 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 
285 void 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 
362 void nrerrorNR(void)
363 {
364 }
365 
366 // Numerical Recipes memory allocation routines
367 // One based arrays
368 
369 float **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 
396 float ***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 
437 void 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 
446 void 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 }
void ref()
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:54 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.