• Skip to content
  • Skip to link menu
KDE API Reference
  • KDE API Reference
  • edu API Reference
  • KDE Home
  • Contact Us
 

kstars

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

KDE's Doxygen guidelines are available online.

kstars

Skip menu "kstars"
  • Main Page
  • Namespace List
  • Namespace Members
  • Alphabetical List
  • Class List
  • Class Hierarchy
  • Class Members
  • File List
  • File Members
  • Modules
  • Related Pages

edu API Reference

Skip menu "edu API Reference"
  •     core
  • kstars

Search



Report problems with this website to our bug tracking system.
Contact the specific authors with questions and comments about the page contents.

KDE® and the K Desktop Environment® logo are registered trademarks of KDE e.V. | Legal