Kstars

focusfourierpower.h
1/*
2 SPDX-FileCopyrightText: 2023
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5*/
6
7#pragma once
8
9#include <QList>
10#include "../fitsviewer/fitsstardetector.h"
11#include "fitsviewer/fitsview.h"
12#include "fitsviewer/fitsdata.h"
13#include "auxiliary/imagemask.h"
14#include "aberrationinspectorutils.h"
15#include "curvefit.h"
16#include "../ekos.h"
17#include <ekos_focus_debug.h>
18#include <gsl/gsl_fft_complex.h>
19
20namespace Ekos
21{
22
23// This class performs Fourier Transform processing on the passed in FITS image.
24// The approach is based on this paper: https://doi.org/10.1093/mnras/stac189
25// (preprint, see https://arxiv.org/abs/2201.12466 )
26// The idea is as follows: a focus frame consists of stars and noise. As long as the
27// exposure is long enough (a few seconds) then the star can be treated as gaussians.
28// A fourier transform (FFT) of a gaussian is another gaussian. A narrow gaussian in the space
29// domain transforms to a wider gaussian in frequency domain. So the nearer to focus we
30// get, the sharper the stars (smaller HFR) become, and thereform the wider their FFTs.
31// So the nearer to focus we get the bigger the frequency domain's power spectrum becomes
32// which means the power reaches a maximum in the frequency domain.
33//
34// Random noise transforms to random noise in the frequency domain. The power spectrum of
35// the FFT of noise typically swamps the power spectrum of the stars' signal so it needs to
36// be dealt with. Tan and Schulz suggest treating noise as "white noise" which means it adds
37// the same power component to all frequencies. This means that be averaging the power
38// high frequency where the star contribution is zero, gives the noise component. This can
39// then to subtracted from each frequency to leave the power contribution of just the stars.
40//
41// I have found better results by background subtracting the frame before fourier transforming
42// which is the same thing as all thats left to FFT is an image of the stars. Currently
43// fitsviewer will have procvessed the image prior to this routine being called so background
44// information is available here.
45//
46// So we need to perform a 2D FFT on the image. GSL only performs basic 1D FFTs so this
47// routine performs a FFT on each row and then uses the results to perform a FFT on each
48// column. An optimisation would be to use a more sophisticated FFT routine. FFTW3 could,
49// for example, be used.
50//
51// Currently just the first channel (if there is more than 1) is used by this routine. It would
52// be possible to use all channels or offer the user a choice of which channel(s) to use. If
53// this were to be done it would make sense to synchonise this functionality with fitsviwer
54// functionality on HFR, FWHM, etc.
55
56class FocusFourierPower
57{
58 public:
59
60 FocusFourierPower(Mathematics::RobustStatistics::ScaleCalculation scaleCalc);
61 ~FocusFourierPower();
62
63 template <typename T>
64 void processFourierPower(const T &imageBuffer, const QSharedPointer<FITSData> &imageData,
65 const QSharedPointer<ImageMask> &mask, const int &tile, double *fourierPower, double *weight)
66 {
67 // Initialise outputs
68 *fourierPower = INVALID_STAR_MEASURE;
69 *weight = 1.0;
70
71 // Check mask & tile inputs
72 ImageMosaicMask *mosaicMask = nullptr;
73 if (tile >= 0)
74 {
75 if (mask)
76 {
77 if (tile >= NUM_TILES)
78 {
79 qCDebug(KSTARS_EKOS_FOCUS) << QString("%1 Fourier Transform called with invalid mosaic tile %2").arg(__FUNCTION__)
80 .arg(tile);
81 return;
82 }
83
84 mosaicMask = dynamic_cast<ImageMosaicMask *>(mask.get());
85 if (!mosaicMask)
86 {
87 qCDebug(KSTARS_EKOS_FOCUS) << QString("%1 Fourier Transform called with invalid 2 mosaic tile %2").arg(__FUNCTION__)
88 .arg(tile);
89 return;
90 }
91 }
92 else
93 {
94 qCDebug(KSTARS_EKOS_FOCUS) << QString("%1 Fourier Transform called for mosaic tile but no mask").arg(__FUNCTION__);
95 return;
96 }
97 }
98
99 // Dimensions on area to perform Fourier Transform on; whole sensor or just tile
100 auto stats = imageData->getStatistics();
101 unsigned int width, height;
102 if (tile < 0)
103 {
104 width = stats.width;
105 height = stats.height;
106 }
107 else
108 {
109 // Calculating for a single (square) mosaic tile
110 width = mosaicMask->tiles()[tile].width();
111 height = width;
112 }
113
114 // Allocate memory for the Fourier Transform
115 unsigned long N = width * height;
116 double *image = new(std::nothrow) double[2 * N];
117 if (!image)
118 {
119 qCDebug(KSTARS_EKOS_FOCUS) << QString("%1 Unable to allocate memory to perform Fourier Transforms").arg(__FUNCTION__);
120 return;
121 }
122
123 /* Allocate memory for GSL complex wavetables, and workspace */
124 gsl_fft_complex_wavetable *rowWT = gsl_fft_complex_wavetable_alloc(width);
125 gsl_fft_complex_workspace *rowWS = gsl_fft_complex_workspace_alloc(width);
126 gsl_fft_complex_wavetable *colWT = gsl_fft_complex_wavetable_alloc(height);
127 gsl_fft_complex_workspace *colWS = gsl_fft_complex_workspace_alloc(height);
128 if (!rowWT || !rowWS || !colWT || !colWS)
129 {
130 qCDebug(KSTARS_EKOS_FOCUS) << QString("%1 Unable to allocate memory2 to perform Fourier Transforms").arg(__FUNCTION__);
131 delete[] image;
132 return;
133 }
134
135 // Set the gsl error handler off as it aborts the program on error.
136 auto const oldErrorHandler = gsl_set_error_handler_off();
137
138 // Convert the image to complex double datatype as required by GSL FFT
139 // The real part is just the background subtracted pixel value clipped to zero
140 // The imaginary part is zero.
141 auto skyBackground = imageData->getSkyBackground();
142 auto bg = skyBackground.mean + 3.0 * skyBackground.sigma;
143
144 // Load the "image" buffer from the passed in image data. As the loop is quite large I've created 3 loops
145 // each with minimal work inside. This makes the overall code larger but avoids repeated tests within the loop
146 if (tile < 0)
147 {
148 // Setup image for whole sensor
149 if (mask.isNull() || mask->active() == false)
150 {
151 // No active mask
152 for (unsigned long i = 0; i < N; i++)
153 {
154 image[i * 2] = std::max(0.0, (double) imageBuffer[i] - bg);
155 image[i * 2 + 1] = 0.0;
156 }
157 }
158 else
159 {
160 // There is an active mask on the sensor so honour these settings
161 unsigned int posX = 0, posY = 0;
162 for (unsigned long i = 0; i < N; i++)
163 {
164 if (mask->isVisible(posX, posY))
165 image[i * 2] = std::max(0.0, (double) imageBuffer[i] - bg);
166 else
167 image[i * 2] = 0.0;
168
169 image[i * 2 + 1] = 0.0;
170
171 if (++posX == stats.width)
172 {
173 posX = 0;
174 posY++;
175 }
176 }
177 }
178 }
179 else
180 {
181 // A mosaic tile has been specified so we know we are dealing with a mosaic mask
182 unsigned int posX = mosaicMask->tiles()[tile].topLeft().x();
183 unsigned int posY = mosaicMask->tiles()[tile].topLeft().y();
184
185 unsigned long offset = posY * stats.width + posX;
186 const unsigned int widthLimit = posX + width;
187
188 // Perform calc for a specific tile of a mosaic mask
189 for (unsigned long i = 0; i < N; i++)
190 {
191 image[i * 2] = std::max(0.0, (double) imageBuffer[offset + i] - bg);
192 image[i * 2 + 1] = 0.0;
193
194 if (++posX == widthLimit)
195 {
196 posX = mosaicMask->tiles()[tile].topLeft().x();
197 offset += stats.width - width;
198 }
199 }
200 }
201
202 // Perform FFT on all the rows
203 int status = 0;
204 for (unsigned int j = 0; j < height; j++)
205 {
206 status = gsl_fft_complex_forward(&image[j * 2 * width], 1, width, rowWT, rowWS);
207 if (status != 0)
208 {
209 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error %1 [%2] calculating FFT on row %3").arg(status).arg(gsl_strerror(status)).
210 arg(j);
211 break;
212 }
213 }
214
215 if (status == 0)
216 {
217 // Perform FFT on all the cols
218 for (unsigned int i = 0; i < width; i++)
219 {
220 status = gsl_fft_complex_forward(&image[2 * i], width, height, colWT, colWS);
221 if (status != 0)
222 {
223 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error %1 [%2] calculating FFT on col %3").arg(status).arg(gsl_strerror(status)).arg(
224 i);
225 break;
226 }
227 }
228 }
229
230 if (status == 0)
231 {
232 double power = 0.0;
233 for (unsigned long i = 0; i < N; i++)
234 power += pow(image[i * 2], 2.0) + pow(image[i * 2 + 1], 2.0);
235
236 power /= pow(N, 2.0);
237
238 if (tile < 0)
239 qCDebug(KSTARS_EKOS_FOCUS) << QString("FFT power sensor %1x%2 = %3").arg(stats.width).arg(stats.height).arg(power);
240 else
241 qCDebug(KSTARS_EKOS_FOCUS) << QString("FFT power tile %1 %2x%3 = %4").arg(tile).arg(stats.width).arg(stats.height).arg(
242 power);
243
244 *fourierPower = power;
245 }
246
247 // free memory
248 delete[] image;
249 gsl_fft_complex_wavetable_free(rowWT);
250 gsl_fft_complex_workspace_free(rowWS);
251 gsl_fft_complex_wavetable_free(colWT);
252 gsl_fft_complex_workspace_free(colWS);
253
254 // Restore old GSL error handler
255 gsl_set_error_handler(oldErrorHandler);
256 }
257
258 static double constexpr INVALID_STAR_MEASURE = -1.0;
259
260 private:
261
262 Mathematics::RobustStatistics::ScaleCalculation m_ScaleCalc;
263};
264}
Ekos is an advanced Astrophotography tool for Linux.
Definition align.cpp:83
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
T * get() const const
bool isNull() const const
QString arg(Args &&... args) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.