1 /*
2 SPDX-FileCopyrightText: 2004 Jasem Mutlaq <mutlaqja@ikarustech.com>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5
6 Some code fragments were adapted from Peter Kirchgessner's FITS plugin
7 SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
8 */
9
10 #include "fitsdata.h"
11 #include "fitsbahtinovdetector.h"
12 #include "fitsthresholddetector.h"
13 #include "fitsgradientdetector.h"
14 #include "fitscentroiddetector.h"
15 #include "fitssepdetector.h"
16
17 #include "fpack.h"
18
19 #include "kstarsdata.h"
20 #include "ksutils.h"
21 #include "kspaths.h"
22 #include "Options.h"
23 #include "skymapcomposite.h"
24 #include "auxiliary/ksnotification.h"
25
26 #include <KFormat>
27 #include <QApplication>
28 #include <QImage>
29 #include <QtConcurrent>
30 #include <QImageReader>
31
32 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
33 #include <wcshdr.h>
34 #include <wcsfix.h>
35 #endif
36
37 #if !defined(KSTARS_LITE) && defined(HAVE_LIBRAW)
38 #include <libraw/libraw.h>
39 #endif
40
41 #include <cfloat>
42 #include <cmath>
43
44 #include <fits_debug.h>
45
46 #define ZOOM_DEFAULT 100.0
47 #define ZOOM_MIN 10
48 #define ZOOM_MAX 400
49 #define ZOOM_LOW_INCR 10
50 #define ZOOM_HIGH_INCR 50
51
getTemporaryPath()52 QString getTemporaryPath()
53 {
54 return QDir(KSPaths::writableLocation(QStandardPaths::TempLocation) + "/" +
55 qAppName()).path();
56 }
57
58 const QStringList RAWFormats = { "cr2", "cr3", "crw", "nef", "raf", "dng", "arw", "orf" };
59
FITSData(FITSMode fitsMode)60 FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)
61 {
62 qRegisterMetaType<FITSMode>("FITSMode");
63
64 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
65 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
66 debayerParams.offsetX = debayerParams.offsetY = 0;
67
68 // Reserve 3 channels
69 m_CumulativeFrequency.resize(3);
70 m_HistogramBinWidth.resize(3);
71 m_HistogramFrequency.resize(3);
72 m_HistogramIntensity.resize(3);
73 }
74
FITSData(const QSharedPointer<FITSData> & other)75 FITSData::FITSData(const QSharedPointer<FITSData> &other)
76 {
77 qRegisterMetaType<FITSMode>("FITSMode");
78
79 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
80 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
81 debayerParams.offsetX = debayerParams.offsetY = 0;
82
83 m_TemporaryDataFile.setFileTemplate("fits_memory_XXXXXX");
84
85 this->m_Mode = other->m_Mode;
86 this->m_Statistics.channels = other->m_Statistics.channels;
87 memcpy(&m_Statistics, &(other->m_Statistics), sizeof(m_Statistics));
88 m_ImageBuffer = new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel];
89 memcpy(m_ImageBuffer, other->m_ImageBuffer,
90 m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel);
91 }
92
~FITSData()93 FITSData::~FITSData()
94 {
95 int status = 0;
96
97 if (m_StarFindFuture.isRunning())
98 m_StarFindFuture.waitForFinished();
99
100 clearImageBuffers();
101
102 #ifdef HAVE_WCSLIB
103 if (m_WCSHandle != nullptr)
104 wcsvfree(&m_nwcs, &m_WCSHandle);
105 #endif
106
107 if (starCenters.count() > 0)
108 qDeleteAll(starCenters);
109
110 if (m_SkyObjects.count() > 0)
111 qDeleteAll(m_SkyObjects);
112
113 if (fptr != nullptr)
114 {
115 fits_flush_file(fptr, &status);
116 fits_close_file(fptr, &status);
117 fptr = nullptr;
118 }
119 }
120
loadCommon(const QString & inFilename)121 void FITSData::loadCommon(const QString &inFilename)
122 {
123 int status = 0;
124 qDeleteAll(starCenters);
125 starCenters.clear();
126
127 if (fptr != nullptr)
128 {
129 fits_flush_file(fptr, &status);
130 fits_close_file(fptr, &status);
131 fptr = nullptr;
132 }
133
134 m_Filename = inFilename;
135 }
136
loadFromBuffer(const QByteArray & buffer,const QString & extension,const QString & inFilename,bool silent)137 bool FITSData::loadFromBuffer(const QByteArray &buffer, const QString &extension, const QString &inFilename, bool silent)
138 {
139 loadCommon(inFilename);
140 qCDebug(KSTARS_FITS) << "Reading file buffer (" << KFormat().formatByteSize(buffer.size()) << ")";
141 return privateLoad(buffer, extension, silent);
142 }
143
loadFromFile(const QString & inFilename,bool silent)144 QFuture<bool> FITSData::loadFromFile(const QString &inFilename, bool silent)
145 {
146 loadCommon(inFilename);
147 QFileInfo info(m_Filename);
148 QString extension = info.completeSuffix().toLower();
149 qCInfo(KSTARS_FITS) << "Loading file " << m_Filename;
150 return QtConcurrent::run(this, &FITSData::privateLoad, QByteArray(), extension, silent);
151 }
152
153 namespace
154 {
155 // Common code for reporting fits read errors. Always returns false.
fitsOpenError(int status,const QString & message,bool silent)156 bool fitsOpenError(int status, const QString &message, bool silent)
157 {
158 char error_status[512];
159 fits_report_error(stderr, status);
160 fits_get_errstatus(status, error_status);
161 QString errMessage = message;
162 errMessage.append(i18n(" Error: %1", QString::fromUtf8(error_status)));
163 if (!silent)
164 KSNotification::error(errMessage, i18n("FITS Open"));
165 qCCritical(KSTARS_FITS) << errMessage;
166 return false;
167 }
168 }
169
privateLoad(const QByteArray & buffer,const QString & extension,bool silent)170 bool FITSData::privateLoad(const QByteArray &buffer, const QString &extension, bool silent)
171 {
172 m_isTemporary = m_Filename.startsWith(KSPaths::writableLocation(QStandardPaths::TempLocation));
173 cacheHFR = -1;
174 cacheEccentricity = -1;
175
176 if (extension.contains("fit"))
177 return loadFITSImage(buffer, extension, silent);
178 if (QImageReader::supportedImageFormats().contains(extension.toLatin1()))
179 return loadCanonicalImage(buffer, extension, silent);
180 else if (RAWFormats.contains(extension))
181 return loadRAWImage(buffer, extension, silent);
182
183 return false;
184 }
185
loadFITSImage(const QByteArray & buffer,const QString & extension,bool silent)186 bool FITSData::loadFITSImage(const QByteArray &buffer, const QString &extension, bool silent)
187 {
188 int status = 0, anynull = 0;
189 long naxes[3];
190 QString errMessage;
191
192 m_HistogramConstructed = false;
193
194 if (buffer.isEmpty() && extension.contains(".fz"))
195 {
196 // Store so we don't lose.
197 m_compressedFilename = m_Filename;
198
199 QString uncompressedFile = QDir::tempPath() + QString("/%1").arg(QUuid::createUuid().toString().remove(
200 QRegularExpression("[-{}]")));
201 fpstate fpvar;
202 fp_init (&fpvar);
203 if (fp_unpack(m_Filename.toLocal8Bit().data(), uncompressedFile.toLocal8Bit().data(), fpvar) < 0)
204 {
205 errMessage = i18n("Failed to unpack compressed fits");
206 if (!silent)
207 KSNotification::error(errMessage, i18n("FITS Open"));
208 qCCritical(KSTARS_FITS) << errMessage;
209 return false;
210 }
211
212 m_Filename = uncompressedFile;
213 m_isTemporary = true;
214 m_isCompressed = true;
215 }
216
217 if (buffer.isEmpty())
218 {
219 // Use open diskfile as it does not use extended file names which has problems opening
220 // files with [ ] or ( ) in their names.
221 if (fits_open_diskfile(&fptr, m_Filename.toLocal8Bit(), READONLY, &status))
222 {
223 recordLastError(status);
224 return fitsOpenError(status, i18n("Error opening fits file %1", m_Filename), silent);
225 }
226
227 m_Statistics.size = QFile(m_Filename).size();
228 }
229 else
230 {
231 // Read the FITS file from a memory buffer.
232 void *temp_buffer = const_cast<void *>(reinterpret_cast<const void *>(buffer.data()));
233 size_t temp_size = buffer.size();
234 if (fits_open_memfile(&fptr, m_Filename.toLocal8Bit().data(), READONLY,
235 &temp_buffer, &temp_size, 0, nullptr, &status))
236 {
237 recordLastError(status);
238 return fitsOpenError(status, i18n("Error reading fits buffer."), silent);
239 }
240
241 m_Statistics.size = temp_size;
242 }
243
244 if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status))
245 {
246 recordLastError(status);
247 return fitsOpenError(status, i18n("Could not locate image HDU."), silent);
248 }
249
250 if (fits_get_img_param(fptr, 3, &m_FITSBITPIX, &(m_Statistics.ndim), naxes, &status))
251 {
252 recordLastError(status);
253 return fitsOpenError(status, i18n("FITS file open error (fits_get_img_param)."), silent);
254 }
255
256 if (m_Statistics.ndim < 2)
257 {
258 m_LastError = i18n("1D FITS images are not supported in KStars.");
259 if (!silent)
260 KSNotification::error(m_LastError, i18n("FITS Open"));
261 qCCritical(KSTARS_FITS) << m_LastError;
262 return false;
263 }
264
265 switch (m_FITSBITPIX)
266 {
267 case BYTE_IMG:
268 m_Statistics.dataType = TBYTE;
269 m_Statistics.bytesPerPixel = sizeof(uint8_t);
270 break;
271 case SHORT_IMG:
272 // Read SHORT image as USHORT
273 m_FITSBITPIX = USHORT_IMG;
274 m_Statistics.dataType = TUSHORT;
275 m_Statistics.bytesPerPixel = sizeof(int16_t);
276 break;
277 case USHORT_IMG:
278 m_Statistics.dataType = TUSHORT;
279 m_Statistics.bytesPerPixel = sizeof(uint16_t);
280 break;
281 case LONG_IMG:
282 // Read LONG image as ULONG
283 m_FITSBITPIX = ULONG_IMG;
284 m_Statistics.dataType = TULONG;
285 m_Statistics.bytesPerPixel = sizeof(int32_t);
286 break;
287 case ULONG_IMG:
288 m_Statistics.dataType = TULONG;
289 m_Statistics.bytesPerPixel = sizeof(uint32_t);
290 break;
291 case FLOAT_IMG:
292 m_Statistics.dataType = TFLOAT;
293 m_Statistics.bytesPerPixel = sizeof(float);
294 break;
295 case LONGLONG_IMG:
296 m_Statistics.dataType = TLONGLONG;
297 m_Statistics.bytesPerPixel = sizeof(int64_t);
298 break;
299 case DOUBLE_IMG:
300 m_Statistics.dataType = TDOUBLE;
301 m_Statistics.bytesPerPixel = sizeof(double);
302 break;
303 default:
304 m_LastError = i18n("Bit depth %1 is not supported.", m_FITSBITPIX);
305 if (!silent)
306 KSNotification::error(m_LastError, i18n("FITS Open"));
307 qCCritical(KSTARS_FITS) << m_LastError;
308 return false;
309 }
310
311 if (m_Statistics.ndim < 3)
312 naxes[2] = 1;
313
314 if (naxes[0] == 0 || naxes[1] == 0)
315 {
316 m_LastError = i18n("Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
317 if (!silent)
318 KSNotification::error(m_LastError, i18n("FITS Open"));
319 qCCritical(KSTARS_FITS) << m_LastError;
320 return false;
321 }
322
323 m_Statistics.width = naxes[0];
324 m_Statistics.height = naxes[1];
325 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
326
327 clearImageBuffers();
328
329 m_Statistics.channels = naxes[2];
330
331 // Channels always set to #1 if we are not required to process 3D Cubes
332 // Or if mode is not FITS_NORMAL (guide, focus..etc)
333 if (m_Mode != FITS_NORMAL || !Options::auto3DCube())
334 m_Statistics.channels = 1;
335
336 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
337 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
338 if (m_ImageBuffer == nullptr)
339 {
340 qCWarning(KSTARS_FITS) << "FITSData: Not enough memory for image_buffer channel. Requested: "
341 << m_ImageBufferSize << " bytes.";
342 clearImageBuffers();
343 return false;
344 }
345
346 rotCounter = 0;
347 flipHCounter = 0;
348 flipVCounter = 0;
349 long nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
350
351 if (fits_read_img(fptr, m_Statistics.dataType, 1, nelements, nullptr, m_ImageBuffer, &anynull, &status))
352 {
353 recordLastError(status);
354 return fitsOpenError(status, i18n("Error reading image."), silent);
355 }
356
357 parseHeader();
358
359 // Get UTC date time
360 QVariant value;
361 if (getRecordValue("DATE-OBS", value) && value.isValid())
362 {
363 QDateTime ts = value.toDateTime();
364 m_DateTime = KStarsDateTime(ts.date(), ts.time());
365 }
366
367 // Only check for debayed IF the original naxes[2] is 1
368 // which is for single channels.
369 if (naxes[2] == 1 && m_Statistics.channels == 1 && Options::autoDebayer() && checkDebayer())
370 {
371 // Save bayer image on disk in case we need to save it later since debayer destorys this data
372 if (m_isTemporary && m_TemporaryDataFile.open())
373 {
374 m_TemporaryDataFile.write(buffer);
375 m_TemporaryDataFile.close();
376 m_Filename = m_TemporaryDataFile.fileName();
377 }
378
379 if (debayer())
380 calculateStats();
381 }
382 else
383 calculateStats();
384
385 if (m_Mode == FITS_NORMAL || m_Mode == FITS_ALIGN)
386 checkForWCS();
387
388 starsSearched = false;
389
390 return true;
391 }
392
loadCanonicalImage(const QByteArray & buffer,const QString & extension,bool silent)393 bool FITSData::loadCanonicalImage(const QByteArray &buffer, const QString &extension, bool silent)
394 {
395 // TODO need to add error popups as well later on
396 Q_UNUSED(silent);
397 QString errMessage;
398 QImage imageFromFile;
399 if (!buffer.isEmpty())
400 {
401 if(!imageFromFile.loadFromData(reinterpret_cast<const uint8_t*>(buffer.data()), buffer.size(),
402 extension.toLatin1().constData()))
403 {
404 qCCritical(KSTARS_FITS) << "Failed to open image.";
405 return false;
406 }
407
408 m_Statistics.size = buffer.size();
409 }
410 else
411 {
412 if(!imageFromFile.load(m_Filename.toLocal8Bit()))
413 {
414 qCCritical(KSTARS_FITS) << "Failed to open image.";
415 return false;
416 }
417
418 m_Statistics.size = QFileInfo(m_Filename).size();
419 }
420
421 //imageFromFile = imageFromFile.convertToFormat(QImage::Format_RGB32);
422
423 // Note: This will need to be changed. I think QT only loads 8 bpp images.
424 // Also the depth method gives the total bits per pixel in the image not just the bits per
425 // pixel in each channel.
426 const int m_FITSBITPIX = 8;
427 switch (m_FITSBITPIX)
428 {
429 case BYTE_IMG:
430 m_Statistics.dataType = TBYTE;
431 m_Statistics.bytesPerPixel = sizeof(uint8_t);
432 break;
433 case SHORT_IMG:
434 // Read SHORT image as USHORT
435 m_Statistics.dataType = TUSHORT;
436 m_Statistics.bytesPerPixel = sizeof(int16_t);
437 break;
438 case USHORT_IMG:
439 m_Statistics.dataType = TUSHORT;
440 m_Statistics.bytesPerPixel = sizeof(uint16_t);
441 break;
442 case LONG_IMG:
443 // Read LONG image as ULONG
444 m_Statistics.dataType = TULONG;
445 m_Statistics.bytesPerPixel = sizeof(int32_t);
446 break;
447 case ULONG_IMG:
448 m_Statistics.dataType = TULONG;
449 m_Statistics.bytesPerPixel = sizeof(uint32_t);
450 break;
451 case FLOAT_IMG:
452 m_Statistics.dataType = TFLOAT;
453 m_Statistics.bytesPerPixel = sizeof(float);
454 break;
455 case LONGLONG_IMG:
456 m_Statistics.dataType = TLONGLONG;
457 m_Statistics.bytesPerPixel = sizeof(int64_t);
458 break;
459 case DOUBLE_IMG:
460 m_Statistics.dataType = TDOUBLE;
461 m_Statistics.bytesPerPixel = sizeof(double);
462 break;
463 default:
464 errMessage = QString("Bit depth %1 is not supported.").arg(m_FITSBITPIX);
465 QMessageBox::critical(nullptr, "Message", errMessage);
466 qCCritical(KSTARS_FITS) << errMessage;
467 return false;
468 }
469
470 m_Statistics.width = static_cast<uint16_t>(imageFromFile.width());
471 m_Statistics.height = static_cast<uint16_t>(imageFromFile.height());
472 m_Statistics.channels = imageFromFile.isGrayscale() ? 1 : 3;
473 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
474 clearImageBuffers();
475 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * static_cast<uint16_t>
476 (m_Statistics.bytesPerPixel);
477 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
478 if (m_ImageBuffer == nullptr)
479 {
480 qCCritical(KSTARS_FITS) << QString("FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ").arg(
481 m_ImageBufferSize);
482 clearImageBuffers();
483 return false;
484 }
485
486 if (m_Statistics.channels == 1)
487 {
488 memcpy(m_ImageBuffer, imageFromFile.bits(), m_ImageBufferSize);
489 }
490 else
491 {
492
493 auto debayered_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
494 auto * original_bayered_buffer = reinterpret_cast<uint8_t *>(imageFromFile.bits());
495
496 // Data in RGBA, with bytes in the order of B,G,R we need to copy them into 3 layers for FITS
497 uint8_t * rBuff = debayered_buffer;
498 uint8_t * gBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height);
499 uint8_t * bBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
500
501 int imax = m_Statistics.samples_per_channel * 4 - 4;
502 for (int i = 0; i <= imax; i += 4)
503 {
504 *rBuff++ = original_bayered_buffer[i + 2];
505 *gBuff++ = original_bayered_buffer[i + 1];
506 *bBuff++ = original_bayered_buffer[i + 0];
507 }
508 }
509
510 calculateStats();
511 return true;
512 }
513
loadRAWImage(const QByteArray & buffer,const QString & extension,bool silent)514 bool FITSData::loadRAWImage(const QByteArray &buffer, const QString &extension, bool silent)
515 {
516 // TODO need to add error popups as well later on
517 Q_UNUSED(silent);
518 Q_UNUSED(extension);
519
520 #if !defined(KSTARS_LITE) && !defined(HAVE_LIBRAW)
521 m_LastError = i18n("Unable to find dcraw and cjpeg. Please install the required tools to convert CR2/NEF to JPEG.");
522 return false;
523 #else
524
525 int ret = 0;
526 // Creation of image processing object
527 LibRaw RawProcessor;
528
529 // Let us open the file/buffer
530 if (buffer.isEmpty())
531 {
532 // Open file
533 if ((ret = RawProcessor.open_file(m_Filename.toLocal8Bit().constData())) != LIBRAW_SUCCESS)
534 {
535 m_LastError = i18n("Cannot open file %1: %2", m_Filename, libraw_strerror(ret));
536 RawProcessor.recycle();
537 return false;
538 }
539
540 m_Statistics.size = QFileInfo(m_Filename).size();
541 }
542 // Open Buffer
543 else
544 {
545 if ((ret = RawProcessor.open_buffer(const_cast<void *>(reinterpret_cast<const void *>(buffer.data())),
546 buffer.size())) != LIBRAW_SUCCESS)
547 {
548 m_LastError = i18n("Cannot open buffer: %1", libraw_strerror(ret));
549 RawProcessor.recycle();
550 return false;
551 }
552
553 m_Statistics.size = buffer.size();
554 }
555
556 // Let us unpack the thumbnail
557 if ((ret = RawProcessor.unpack()) != LIBRAW_SUCCESS)
558 {
559 m_LastError = i18n("Cannot unpack_thumb: %1", libraw_strerror(ret));
560 RawProcessor.recycle();
561 return false;
562 }
563
564 if ((ret = RawProcessor.dcraw_process()) != LIBRAW_SUCCESS)
565 {
566 m_LastError = i18n("Cannot dcraw_process: %1", libraw_strerror(ret));
567 RawProcessor.recycle();
568 return false;
569 }
570
571 libraw_processed_image_t *image = RawProcessor.dcraw_make_mem_image(&ret);
572 if (ret != LIBRAW_SUCCESS)
573 {
574 m_LastError = i18n("Cannot load to memory: %1", libraw_strerror(ret));
575 RawProcessor.recycle();
576 return false;
577 }
578
579 RawProcessor.recycle();
580
581 m_Statistics.bytesPerPixel = image->bits / 8;
582 // We only support two types now
583 if (m_Statistics.bytesPerPixel == 1)
584 m_Statistics.dataType = TBYTE;
585 else
586 m_Statistics.dataType = TUSHORT;
587 m_Statistics.width = image->width;
588 m_Statistics.height = image->height;
589 m_Statistics.channels = image->colors;
590 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
591 clearImageBuffers();
592 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
593 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
594 if (m_ImageBuffer == nullptr)
595 {
596 qCCritical(KSTARS_FITS) << QString("FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ").arg(
597 m_ImageBufferSize);
598 libraw_dcraw_clear_mem(image);
599 clearImageBuffers();
600 return false;
601 }
602
603 auto destination_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
604 auto source_buffer = reinterpret_cast<uint8_t *>(image->data);
605
606 // For mono, we memcpy directly
607 if (image->colors == 1)
608 {
609 memcpy(destination_buffer, source_buffer, m_ImageBufferSize);
610 }
611 else
612 {
613 // Data in RGB24, with bytes in the order of R,G,B. We copy them copy them into 3 layers for FITS
614 uint8_t * rBuff = destination_buffer;
615 uint8_t * gBuff = destination_buffer + (m_Statistics.width * m_Statistics.height);
616 uint8_t * bBuff = destination_buffer + (m_Statistics.width * m_Statistics.height * 2);
617
618 int imax = m_Statistics.samples_per_channel * 3 - 3;
619 for (int i = 0; i <= imax; i += 3)
620 {
621 *rBuff++ = source_buffer[i + 0];
622 *gBuff++ = source_buffer[i + 1];
623 *bBuff++ = source_buffer[i + 2];
624 }
625 }
626 libraw_dcraw_clear_mem(image);
627
628 calculateStats();
629 return true;
630 #endif
631 }
632
saveImage(const QString & newFilename)633 bool FITSData::saveImage(const QString &newFilename)
634 {
635 if (newFilename == m_Filename)
636 return true;
637
638 const QString ext = QFileInfo(newFilename).suffix();
639
640 if (ext == "jpg" || ext == "png")
641 {
642 double min, max;
643 QImage fitsImage = FITSToImage(newFilename);
644 getMinMax(&min, &max);
645
646 if (min == max)
647 {
648 fitsImage.fill(Qt::white);
649 }
650 else if (channels() == 1)
651 {
652 fitsImage = QImage(width(), height(), QImage::Format_Indexed8);
653
654 fitsImage.setColorCount(256);
655 for (int i = 0; i < 256; i++)
656 fitsImage.setColor(i, qRgb(i, i, i));
657 }
658 else
659 {
660 fitsImage = QImage(width(), height(), QImage::Format_RGB32);
661 }
662
663 double dataMin = m_Statistics.mean[0] - m_Statistics.stddev[0];
664 double dataMax = m_Statistics.mean[0] + m_Statistics.stddev[0] * 3;
665
666 double bscale = 255. / (dataMax - dataMin);
667 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
668
669 // Long way to do this since we do not want to use templated functions here
670 switch (m_Statistics.dataType)
671 {
672 case TBYTE:
673 convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
674 break;
675
676 case TSHORT:
677 convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
678 break;
679
680 case TUSHORT:
681 convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
682 break;
683
684 case TLONG:
685 convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
686 break;
687
688 case TULONG:
689 convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
690 break;
691
692 case TFLOAT:
693 convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
694 break;
695
696 case TLONGLONG:
697 convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
698 break;
699
700 case TDOUBLE:
701 convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
702 break;
703
704 default:
705 break;
706 }
707
708 fitsImage.save(newFilename, ext.toLatin1().constData());
709
710 m_Filename = newFilename;
711 qCInfo(KSTARS_FITS) << "Saved image file:" << m_Filename;
712 return true;
713 }
714
715 if (m_isCompressed)
716 {
717 KSNotification::error(i18n("Saving compressed files is not supported."));
718 return false;
719 }
720
721 int status = 0;
722 long nelements;
723 fitsfile * new_fptr;
724
725 // if (HasDebayer && m_Filename.isEmpty() == false)
726 // {
727 // fits_flush_file(fptr, &status);
728 // /* close current file */
729 // if (fits_close_file(fptr, &status))
730 // {
731 // recordLastError(status);
732 // return status;
733 // }
734
735 // // Skip "!" in the beginning of the new file name
736 // QString finalFileName(newFilename);
737
738 // finalFileName.remove('!');
739
740 // // Remove first otherwise copy will fail below if file exists
741 // QFile::remove(finalFileName);
742
743 // if (!QFile::copy(m_Filename, finalFileName))
744 // {
745 // qCCritical(KSTARS_FITS()) << "FITS: Failed to copy " << m_Filename << " to " << finalFileName;
746 // fptr = nullptr;
747 // return false;
748 // }
749
750 // m_Filename = finalFileName;
751
752 // // Use open diskfile as it does not use extended file names which has problems opening
753 // // files with [ ] or ( ) in their names.
754 // fits_open_diskfile(&fptr, m_Filename.toLocal8Bit(), READONLY, &status);
755 // fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status);
756
757 // return true;
758 // }
759
760 // Read the image back into buffer in case we debyayed
761 nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
762
763 /* close current file */
764 if (fptr && fits_close_file(fptr, &status))
765 {
766 recordLastError(status);
767 return false;
768 }
769
770 /* Create a new File, overwriting existing*/
771 if (fits_create_file(&new_fptr, QString("!%1").arg(newFilename).toLocal8Bit(), &status))
772 {
773 recordLastError(status);
774 return status;
775 }
776
777 status = 0;
778
779 fptr = new_fptr;
780
781 // Create image
782 long naxis = m_Statistics.channels == 1 ? 2 : 3;
783 long naxes[3] = {m_Statistics.width, m_Statistics.height, naxis};
784
785 // JM 2020-12-28: Here we to use bitpix values
786 if (fits_create_img(fptr, m_FITSBITPIX, naxis, naxes, &status))
787 {
788 recordLastError(status);
789 return false;
790 }
791
792 /* Write keywords */
793
794 // Minimum
795 if (fits_update_key(fptr, TDOUBLE, "DATAMIN", &(m_Statistics.min), "Minimum value", &status))
796 {
797 recordLastError(status);
798 return false;
799 }
800
801 // Maximum
802 if (fits_update_key(fptr, TDOUBLE, "DATAMAX", &(m_Statistics.max), "Maximum value", &status))
803 {
804 recordLastError(status);
805 return false;
806 }
807
808 // KStars Min, for 3 channels
809 fits_write_key(fptr, TDOUBLE, "MIN1", &m_Statistics.min[0], "Min Channel 1", &status);
810 if (channels() > 1)
811 {
812 fits_write_key(fptr, TDOUBLE, "MIN2", &m_Statistics.min[1], "Min Channel 2", &status);
813 fits_write_key(fptr, TDOUBLE, "MIN3", &m_Statistics.min[2], "Min Channel 3", &status);
814 }
815
816 // KStars max, for 3 channels
817 fits_write_key(fptr, TDOUBLE, "MAX1", &m_Statistics.max[0], "Max Channel 1", &status);
818 if (channels() > 1)
819 {
820 fits_write_key(fptr, TDOUBLE, "MAX2", &m_Statistics.max[1], "Max Channel 2", &status);
821 fits_write_key(fptr, TDOUBLE, "MAX3", &m_Statistics.max[2], "Max Channel 3", &status);
822 }
823
824 // Mean
825 if (m_Statistics.mean[0] > 0)
826 {
827 fits_write_key(fptr, TDOUBLE, "MEAN1", &m_Statistics.mean[0], "Mean Channel 1", &status);
828 if (channels() > 1)
829 {
830 fits_write_key(fptr, TDOUBLE, "MEAN2", &m_Statistics.mean[1], "Mean Channel 2", &status);
831 fits_write_key(fptr, TDOUBLE, "MEAN3", &m_Statistics.mean[2], "Mean Channel 3", &status);
832 }
833 }
834
835 // Median
836 if (m_Statistics.median[0] > 0)
837 {
838 fits_write_key(fptr, TDOUBLE, "MEDIAN1", &m_Statistics.median[0], "Median Channel 1", &status);
839 if (channels() > 1)
840 {
841 fits_write_key(fptr, TDOUBLE, "MEDIAN2", &m_Statistics.median[1], "Median Channel 2", &status);
842 fits_write_key(fptr, TDOUBLE, "MEDIAN3", &m_Statistics.median[2], "Median Channel 3", &status);
843 }
844 }
845
846 // Standard Deviation
847 if (m_Statistics.stddev[0] > 0)
848 {
849 fits_write_key(fptr, TDOUBLE, "STDDEV1", &m_Statistics.stddev[0], "Standard Deviation Channel 1", &status);
850 if (channels() > 1)
851 {
852 fits_write_key(fptr, TDOUBLE, "STDDEV2", &m_Statistics.stddev[1], "Standard Deviation Channel 2", &status);
853 fits_write_key(fptr, TDOUBLE, "STDDEV3", &m_Statistics.stddev[2], "Standard Deviation Channel 3", &status);
854 }
855 }
856
857 // Skip first 10 standard records and copy the rest.
858 for (int i = 10; i < m_HeaderRecords.count(); i++)
859 {
860 QString key = m_HeaderRecords[i].key;
861 const char *comment = m_HeaderRecords[i].comment.toLatin1().constBegin();
862 QVariant value = m_HeaderRecords[i].value;
863
864 switch (value.type())
865 {
866 case QVariant::Int:
867 {
868 int number = value.toInt();
869 fits_write_key(fptr, TINT, key.toLatin1().constData(), &number, comment, &status);
870 }
871 break;
872
873 case QVariant::Double:
874 {
875 double number = value.toDouble();
876 fits_write_key(fptr, TDOUBLE, key.toLatin1().constData(), &number, comment, &status);
877 }
878 break;
879
880 case QVariant::String:
881 default:
882 {
883 char valueBuffer[256] = {0};
884 strncpy(valueBuffer, value.toString().toLatin1().constData(), 256 - 1);
885 fits_write_key(fptr, TSTRING, key.toLatin1().constData(), valueBuffer, comment, &status);
886 }
887 }
888 }
889
890 // ISO Date
891 if (fits_write_date(fptr, &status))
892 {
893 recordLastError(status);
894 return false;
895 }
896
897 QString history =
898 QString("Modified by KStars on %1").arg(QDateTime::currentDateTime().toString("yyyy-MM-ddThh:mm:ss"));
899 // History
900 if (fits_write_history(fptr, history.toLatin1(), &status))
901 {
902 recordLastError(status);
903 return false;
904 }
905
906 int rot = 0, mirror = 0;
907 if (rotCounter > 0)
908 {
909 rot = (90 * rotCounter) % 360;
910 if (rot < 0)
911 rot += 360;
912 }
913 if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0)
914 mirror = 1;
915
916 if ((rot != 0) || (mirror != 0))
917 rotWCSFITS(rot, mirror);
918
919 rotCounter = flipHCounter = flipVCounter = 0;
920
921 // Here we need to use the actual data type
922 if (fits_write_img(fptr, m_Statistics.dataType, 1, nelements, m_ImageBuffer, &status))
923 {
924 recordLastError(status);
925 return false;
926 }
927
928 m_Filename = newFilename;
929
930 fits_flush_file(fptr, &status);
931
932 qCInfo(KSTARS_FITS) << "Saved FITS file:" << m_Filename;
933
934 return true;
935 }
936
clearImageBuffers()937 void FITSData::clearImageBuffers()
938 {
939 delete[] m_ImageBuffer;
940 m_ImageBuffer = nullptr;
941 //m_BayerBuffer = nullptr;
942 }
943
calculateStats(bool refresh)944 void FITSData::calculateStats(bool refresh)
945 {
946 // Calculate min max
947 calculateMinMax(refresh);
948 calculateMedian(refresh);
949
950 // Try to read mean/median/stddev if in file
951 if (refresh == false && fptr)
952 {
953 int status = 0;
954 int nfound = 0;
955 if (fits_read_key_dbl(fptr, "MEAN1", &m_Statistics.mean[0], nullptr, &status) == 0)
956 nfound++;
957 // NB. These could fail if missing, which is OK.
958 fits_read_key_dbl(fptr, "MEAN2", & m_Statistics.mean[1], nullptr, &status);
959 fits_read_key_dbl(fptr, "MEAN3", &m_Statistics.mean[2], nullptr, &status);
960
961 status = 0;
962 if (fits_read_key_dbl(fptr, "STDDEV1", &m_Statistics.stddev[0], nullptr, &status) == 0)
963 nfound++;
964 // NB. These could fail if missing, which is OK.
965 fits_read_key_dbl(fptr, "STDDEV2", &m_Statistics.stddev[1], nullptr, &status);
966 fits_read_key_dbl(fptr, "STDDEV3", &m_Statistics.stddev[2], nullptr, &status);
967
968 // If all is OK, we're done
969 if (nfound == 2)
970 return;
971 }
972
973 // Get standard deviation and mean in one run
974 switch (m_Statistics.dataType)
975 {
976 case TBYTE:
977 runningAverageStdDev<uint8_t>();
978 break;
979
980 case TSHORT:
981 runningAverageStdDev<int16_t>();
982 break;
983
984 case TUSHORT:
985 runningAverageStdDev<uint16_t>();
986 break;
987
988 case TLONG:
989 runningAverageStdDev<int32_t>();
990 break;
991
992 case TULONG:
993 runningAverageStdDev<uint32_t>();
994 break;
995
996 case TFLOAT:
997 runningAverageStdDev<float>();
998 break;
999
1000 case TLONGLONG:
1001 runningAverageStdDev<int64_t>();
1002 break;
1003
1004 case TDOUBLE:
1005 runningAverageStdDev<double>();
1006 break;
1007
1008 default:
1009 return;
1010 }
1011
1012 // FIXME That's not really SNR, must implement a proper solution for this value
1013 m_Statistics.SNR = m_Statistics.mean[0] / m_Statistics.stddev[0];
1014 }
1015
calculateMinMax(bool refresh)1016 void FITSData::calculateMinMax(bool refresh)
1017 {
1018 int status, nfound = 0;
1019
1020 status = 0;
1021
1022 // Only fetch from header if we have a single channel
1023 // Otherwise, calculate manually.
1024 if (fptr != nullptr && !refresh)
1025 {
1026
1027 status = 0;
1028
1029 if (fits_read_key_dbl(fptr, "DATAMIN", &(m_Statistics.min[0]), nullptr, &status) == 0)
1030 nfound++;
1031 else if (fits_read_key_dbl(fptr, "MIN1", &(m_Statistics.min[0]), nullptr, &status) == 0)
1032 nfound++;
1033
1034 // NB. These could fail if missing, which is OK.
1035 fits_read_key_dbl(fptr, "MIN2", &m_Statistics.min[1], nullptr, &status);
1036 fits_read_key_dbl(fptr, "MIN3", &m_Statistics.min[2], nullptr, &status);
1037
1038 status = 0;
1039
1040 if (fits_read_key_dbl(fptr, "DATAMAX", &(m_Statistics.max[0]), nullptr, &status) == 0)
1041 nfound++;
1042 else if (fits_read_key_dbl(fptr, "MAX1", &(m_Statistics.max[0]), nullptr, &status) == 0)
1043 nfound++;
1044
1045 // NB. These could fail if missing, which is OK.
1046 fits_read_key_dbl(fptr, "MAX2", &m_Statistics.max[1], nullptr, &status);
1047 fits_read_key_dbl(fptr, "MAX3", &m_Statistics.max[2], nullptr, &status);
1048
1049 // If we found both keywords, no need to calculate them, unless they are both zeros
1050 if (nfound == 2 && !(m_Statistics.min[0] == 0 && m_Statistics.max[0] == 0))
1051 return;
1052 }
1053
1054 m_Statistics.min[0] = 1.0E30;
1055 m_Statistics.max[0] = -1.0E30;
1056
1057 m_Statistics.min[1] = 1.0E30;
1058 m_Statistics.max[1] = -1.0E30;
1059
1060 m_Statistics.min[2] = 1.0E30;
1061 m_Statistics.max[2] = -1.0E30;
1062
1063 switch (m_Statistics.dataType)
1064 {
1065 case TBYTE:
1066 calculateMinMax<uint8_t>();
1067 break;
1068
1069 case TSHORT:
1070 calculateMinMax<int16_t>();
1071 break;
1072
1073 case TUSHORT:
1074 calculateMinMax<uint16_t>();
1075 break;
1076
1077 case TLONG:
1078 calculateMinMax<int32_t>();
1079 break;
1080
1081 case TULONG:
1082 calculateMinMax<uint32_t>();
1083 break;
1084
1085 case TFLOAT:
1086 calculateMinMax<float>();
1087 break;
1088
1089 case TLONGLONG:
1090 calculateMinMax<int64_t>();
1091 break;
1092
1093 case TDOUBLE:
1094 calculateMinMax<double>();
1095 break;
1096
1097 default:
1098 break;
1099 }
1100 }
1101
calculateMedian(bool refresh)1102 void FITSData::calculateMedian(bool refresh)
1103 {
1104 int status, nfound = 0;
1105
1106 status = 0;
1107
1108 // Only fetch from header if we have a single channel
1109 // Otherwise, calculate manually.
1110 if (fptr != nullptr && !refresh)
1111 {
1112 status = 0;
1113 if (fits_read_key_dbl(fptr, "MEDIAN1", &m_Statistics.median[0], nullptr, &status) == 0)
1114 nfound++;
1115
1116 // NB. These could fail if missing, which is OK.
1117 fits_read_key_dbl(fptr, "MEDIAN2", &m_Statistics.median[1], nullptr, &status);
1118 fits_read_key_dbl(fptr, "MEDIAN3", &m_Statistics.median[2], nullptr, &status);
1119
1120 if (nfound == 1)
1121 return;
1122 }
1123
1124 m_Statistics.median[RED_CHANNEL] = 0;
1125 m_Statistics.median[GREEN_CHANNEL] = 0;
1126 m_Statistics.median[BLUE_CHANNEL] = 0;
1127
1128 switch (m_Statistics.dataType)
1129 {
1130 case TBYTE:
1131 calculateMedian<uint8_t>();
1132 break;
1133
1134 case TSHORT:
1135 calculateMedian<int16_t>();
1136 break;
1137
1138 case TUSHORT:
1139 calculateMedian<uint16_t>();
1140 break;
1141
1142 case TLONG:
1143 calculateMedian<int32_t>();
1144 break;
1145
1146 case TULONG:
1147 calculateMedian<uint32_t>();
1148 break;
1149
1150 case TFLOAT:
1151 calculateMedian<float>();
1152 break;
1153
1154 case TLONGLONG:
1155 calculateMedian<int64_t>();
1156 break;
1157
1158 case TDOUBLE:
1159 calculateMedian<double>();
1160 break;
1161
1162 default:
1163 break;
1164 }
1165 }
1166
1167 template <typename T>
calculateMedian()1168 void FITSData::calculateMedian()
1169 {
1170 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
1171 const uint32_t maxMedianSize = 500000;
1172 uint32_t medianSize = m_Statistics.samples_per_channel;
1173 uint8_t downsample = 1;
1174 if (medianSize > maxMedianSize)
1175 {
1176 downsample = (static_cast<double>(medianSize) / maxMedianSize) + 0.999;
1177 medianSize /= downsample;
1178 }
1179 std::vector<T> samples;
1180 samples.reserve(medianSize);
1181
1182 for (uint8_t n = 0; n < m_Statistics.channels; n++)
1183 {
1184 auto *oneChannel = buffer + n * m_Statistics.samples_per_channel;
1185 for (uint32_t upto = 0; upto < m_Statistics.samples_per_channel; upto += downsample)
1186 samples.push_back(oneChannel[upto]);
1187 const uint32_t middle = samples.size() / 2;
1188 std::nth_element(samples.begin(), samples.begin() + middle, samples.end());
1189 m_Statistics.median[n] = samples[middle];
1190 }
1191 }
1192
1193 template <typename T>
getParitionMinMax(uint32_t start,uint32_t stride)1194 QPair<T, T> FITSData::getParitionMinMax(uint32_t start, uint32_t stride)
1195 {
1196 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
1197 T min = std::numeric_limits<T>::max();
1198 T max = std::numeric_limits<T>::min();
1199
1200 uint32_t end = start + stride;
1201
1202 for (uint32_t i = start; i < end; i++)
1203 {
1204 min = qMin(buffer[i], min);
1205 max = qMax(buffer[i], max);
1206 // if (buffer[i] < min)
1207 // min = buffer[i];
1208 // else if (buffer[i] > max)
1209 // max = buffer[i];
1210 }
1211
1212 return qMakePair(min, max);
1213 }
1214
1215 template <typename T>
calculateMinMax()1216 void FITSData::calculateMinMax()
1217 {
1218 T min = std::numeric_limits<T>::max();
1219 T max = std::numeric_limits<T>::min();
1220
1221
1222 // Create N threads
1223 const uint8_t nThreads = 16;
1224
1225 for (int n = 0; n < m_Statistics.channels; n++)
1226 {
1227 uint32_t cStart = n * m_Statistics.samples_per_channel;
1228
1229 // Calculate how many elements we process per thread
1230 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
1231
1232 // Calculate the final stride since we can have some left over due to division above
1233 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
1234
1235 // Start location for inspecting elements
1236 uint32_t tStart = cStart;
1237
1238 // List of futures
1239 QList<QFuture<QPair<T, T>>> futures;
1240
1241 for (int i = 0; i < nThreads; i++)
1242 {
1243 // Run threads
1244 futures.append(QtConcurrent::run(this, &FITSData::getParitionMinMax<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride));
1245 tStart += tStride;
1246 }
1247
1248 // Now wait for results
1249 for (int i = 0; i < nThreads; i++)
1250 {
1251 QPair<T, T> result = futures[i].result();
1252 min = qMin(result.first, min);
1253 max = qMax(result.second, max);
1254 }
1255
1256 m_Statistics.min[n] = min;
1257 m_Statistics.max[n] = max;
1258 }
1259 }
1260
1261 template <typename T>
getSquaredSumAndMean(uint32_t start,uint32_t stride)1262 QPair<double, double> FITSData::getSquaredSumAndMean(uint32_t start, uint32_t stride)
1263 {
1264 uint32_t m_n = 2;
1265 double m_oldM = 0, m_newM = 0, m_oldS = 0, m_newS = 0;
1266
1267 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
1268 uint32_t end = start + stride;
1269
1270 for (uint32_t i = start; i < end; i++)
1271 {
1272 m_newM = m_oldM + (buffer[i] - m_oldM) / m_n;
1273 m_newS = m_oldS + (buffer[i] - m_oldM) * (buffer[i] - m_newM);
1274
1275 m_oldM = m_newM;
1276 m_oldS = m_newS;
1277 m_n++;
1278 }
1279
1280 return qMakePair<double, double>(m_newM, m_newS);
1281 }
1282
1283 template <typename T>
runningAverageStdDev()1284 void FITSData::runningAverageStdDev()
1285 {
1286 // Create N threads
1287 const uint8_t nThreads = 16;
1288
1289 for (int n = 0; n < m_Statistics.channels; n++)
1290 {
1291 uint32_t cStart = n * m_Statistics.samples_per_channel;
1292
1293 // Calculate how many elements we process per thread
1294 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
1295
1296 // Calculate the final stride since we can have some left over due to division above
1297 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
1298
1299 // Start location for inspecting elements
1300 uint32_t tStart = cStart;
1301
1302 // List of futures
1303 QList<QFuture<QPair<double, double>>> futures;
1304
1305 for (int i = 0; i < nThreads; i++)
1306 {
1307 // Run threads
1308 futures.append(QtConcurrent::run(this, &FITSData::getSquaredSumAndMean<T>, tStart,
1309 (i == (nThreads - 1)) ? fStride : tStride));
1310 tStart += tStride;
1311 }
1312
1313 double mean = 0, squared_sum = 0;
1314
1315 // Now wait for results
1316 for (int i = 0; i < nThreads; i++)
1317 {
1318 QPair<double, double> result = futures[i].result();
1319 mean += result.first;
1320 squared_sum += result.second;
1321 }
1322
1323 double variance = squared_sum / m_Statistics.samples_per_channel;
1324
1325 m_Statistics.mean[n] = mean / nThreads;
1326 m_Statistics.stddev[n] = sqrt(variance);
1327 }
1328 }
1329
createGaussianKernel(int size,double sigma)1330 QVector<double> FITSData::createGaussianKernel(int size, double sigma)
1331 {
1332 QVector<double> kernel(size * size);
1333 kernel.fill(0.0, size * size);
1334
1335 double kernelSum = 0.0;
1336 int fOff = (size - 1) / 2;
1337 double normal = 1.0 / (2.0 * M_PI * sigma * sigma);
1338 for (int y = -fOff; y <= fOff; y++)
1339 {
1340 for (int x = -fOff; x <= fOff; x++)
1341 {
1342 double distance = ((y * y) + (x * x)) / (2.0 * sigma * sigma);
1343 int index = (y + fOff) * size + (x + fOff);
1344 kernel[index] = normal * qExp(-distance);
1345 kernelSum += kernel.at(index);
1346 }
1347 }
1348 for (int y = 0; y < size; y++)
1349 {
1350 for (int x = 0; x < size; x++)
1351 {
1352 int index = y * size + x;
1353 kernel[index] = kernel.at(index) * 1.0 / kernelSum;
1354 }
1355 }
1356
1357 return kernel;
1358 }
1359
1360 template <typename T>
convolutionFilter(const QVector<double> & kernel,int kernelSize)1361 void FITSData::convolutionFilter(const QVector<double> &kernel, int kernelSize)
1362 {
1363 T * imagePtr = reinterpret_cast<T *>(m_ImageBuffer);
1364
1365 // Create variable for pixel data for each kernel
1366 T gt = 0;
1367
1368 // This is how much your center pixel is offset from the border of your kernel
1369 int fOff = (kernelSize - 1) / 2;
1370
1371 // Start with the pixel that is offset fOff from top and fOff from the left side
1372 // this is so entire kernel is on your image
1373 for (int offsetY = 0; offsetY < m_Statistics.height; offsetY++)
1374 {
1375 for (int offsetX = 0; offsetX < m_Statistics.width; offsetX++)
1376 {
1377 // reset gray value to 0
1378 gt = 0;
1379 // position of the kernel center pixel
1380 int byteOffset = offsetY * m_Statistics.width + offsetX;
1381
1382 // kernel calculations
1383 for (int filterY = -fOff; filterY <= fOff; filterY++)
1384 {
1385 for (int filterX = -fOff; filterX <= fOff; filterX++)
1386 {
1387 if ((offsetY + filterY) >= 0 && (offsetY + filterY) < m_Statistics.height
1388 && ((offsetX + filterX) >= 0 && (offsetX + filterX) < m_Statistics.width ))
1389 {
1390
1391 int calcOffset = byteOffset + filterX + filterY * m_Statistics.width;
1392 int index = (filterY + fOff) * kernelSize + (filterX + fOff);
1393 double kernelValue = kernel.at(index);
1394 gt += (imagePtr[calcOffset]) * kernelValue;
1395 }
1396 }
1397 }
1398
1399 // set new data in the other byte array for your image data
1400 imagePtr[byteOffset] = gt;
1401 }
1402 }
1403 }
1404
1405 template <typename T>
gaussianBlur(int kernelSize,double sigma)1406 void FITSData::gaussianBlur(int kernelSize, double sigma)
1407 {
1408 // Size must be an odd number!
1409 if (kernelSize % 2 == 0)
1410 {
1411 kernelSize--;
1412 qCInfo(KSTARS_FITS) << "Warning, size must be an odd number, correcting size to " << kernelSize;
1413 }
1414 // Edge must be a positive number!
1415 if (kernelSize < 1)
1416 {
1417 kernelSize = 1;
1418 }
1419
1420 QVector<double> gaussianKernel = createGaussianKernel(kernelSize, sigma);
1421 convolutionFilter<T>(gaussianKernel, kernelSize);
1422 }
1423
setMinMax(double newMin,double newMax,uint8_t channel)1424 void FITSData::setMinMax(double newMin, double newMax, uint8_t channel)
1425 {
1426 m_Statistics.min[channel] = newMin;
1427 m_Statistics.max[channel] = newMax;
1428 }
1429
parseHeader()1430 bool FITSData::parseHeader()
1431 {
1432 char * header = nullptr;
1433 int status = 0, nkeys = 0;
1434
1435 if (fits_hdr2str(fptr, 0, nullptr, 0, &header, &nkeys, &status))
1436 {
1437 fits_report_error(stderr, status);
1438 free(header);
1439 return false;
1440 }
1441
1442 m_HeaderRecords.clear();
1443 QString recordList = QString(header);
1444
1445 for (int i = 0; i < nkeys; i++)
1446 {
1447 Record oneRecord;
1448 // Quotes cause issues for simplified below so we're removing them.
1449 QString record = recordList.mid(i * 80, 80).remove("'");
1450 QStringList properties = record.split(QRegExp("[=/]"));
1451 // If it is only a comment
1452 if (properties.size() == 1)
1453 {
1454 oneRecord.key = properties[0].mid(0, 7);
1455 oneRecord.comment = properties[0].mid(8).simplified();
1456 }
1457 else
1458 {
1459 oneRecord.key = properties[0].simplified();
1460 oneRecord.value = properties[1].simplified();
1461 if (properties.size() > 2)
1462 oneRecord.comment = properties[2].simplified();
1463
1464 // Try to guess the value.
1465 // Test for integer & double. If neither, then leave it as "string".
1466 bool ok = false;
1467
1468 // Is it Integer?
1469 oneRecord.value.toInt(&ok);
1470 if (ok)
1471 oneRecord.value.convert(QMetaType::Int);
1472 else
1473 {
1474 // Is it double?
1475 oneRecord.value.toDouble(&ok);
1476 if (ok)
1477 oneRecord.value.convert(QMetaType::Double);
1478 }
1479 }
1480
1481 m_HeaderRecords.append(oneRecord);
1482 }
1483
1484 free(header);
1485
1486 return true;
1487 }
1488
getRecordValue(const QString & key,QVariant & value) const1489 bool FITSData::getRecordValue(const QString &key, QVariant &value) const
1490 {
1491 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](const Record & oneRecord)
1492 {
1493 return (oneRecord.key == key && oneRecord.value.isValid());
1494 });
1495
1496 if (result != m_HeaderRecords.end())
1497 {
1498 value = (*result).value;
1499 return true;
1500 }
1501 return false;
1502 }
1503
parseSolution(FITSImage::Solution & solution) const1504 bool FITSData::parseSolution(FITSImage::Solution &solution) const
1505 {
1506 dms angleValue;
1507 bool raOK = false, deOK = false, coordOK = false, scaleOK = false;
1508 QVariant value;
1509
1510 // Reset all
1511 solution.fieldWidth = solution.fieldHeight = solution.pixscale = solution.ra = solution.dec = -1;
1512
1513 // RA
1514 if (getRecordValue("OBJCTRA", value))
1515 {
1516 angleValue = dms::fromString(value.toString(), false);
1517 solution.ra = angleValue.Degrees();
1518 raOK = true;
1519 }
1520 else if (getRecordValue("RA", value))
1521 {
1522 solution.ra = value.toDouble(&raOK);
1523 }
1524
1525 // DE
1526 if (getRecordValue("OBJCTDEC", value))
1527 {
1528 angleValue = dms::fromString(value.toString(), true);
1529 solution.dec = angleValue.Degrees();
1530 deOK = true;
1531 }
1532 else if (getRecordValue("DEC", value))
1533 {
1534 solution.dec = value.toDouble(&deOK);
1535 }
1536
1537 coordOK = raOK && deOK;
1538
1539 // PixScale
1540 double scale = -1;
1541 if (getRecordValue("SCALE", value))
1542 {
1543 scale = value.toDouble();
1544 }
1545
1546 double focal_length = -1;
1547 if (getRecordValue("FOCALLEN", value))
1548 {
1549 focal_length = value.toDouble();
1550 }
1551
1552 double pixsize1 = -1, pixsize2 = -1;
1553 // Pixel Size 1
1554 if (getRecordValue("PIXSIZE1", value))
1555 {
1556 pixsize1 = value.toDouble();
1557 }
1558 // Pixel Size 2
1559 if (getRecordValue("PIXSIZE2", value))
1560 {
1561 pixsize2 = value.toDouble();
1562 }
1563
1564 int binx = 1, biny = 1;
1565 // Binning X
1566 if (getRecordValue("XBINNING", value))
1567 {
1568 binx = value.toDouble();
1569 }
1570 // Binning Y
1571 if (getRecordValue("YBINNING", value))
1572 {
1573 biny = value.toDouble();
1574 }
1575
1576 if (pixsize1 > 0 && pixsize2 > 0)
1577 {
1578 // If we have scale, then that's it
1579 if (scale > 0)
1580 {
1581 // Arcsecs per pixel
1582 solution.pixscale = scale;
1583 // Arcmins
1584 solution.fieldWidth = (m_Statistics.width * scale) / 60.0;
1585 // Arcmins, and account for pixel ratio if non-squared.
1586 solution.fieldHeight = (m_Statistics.height * scale * (pixsize2 / pixsize1)) / 60.0;
1587 }
1588 else if (focal_length > 0)
1589 {
1590 // Arcmins
1591 solution.fieldWidth = ((206264.8062470963552 * m_Statistics.width * (pixsize1 / 1000.0)) / (focal_length * binx)) / 60.0;
1592 // Arsmins
1593 solution.fieldHeight = ((206264.8062470963552 * m_Statistics.height * (pixsize2 / 1000.0)) / (focal_length * biny)) / 60.0;
1594 // Arcsecs per pixel
1595 solution.pixscale = (206264.8062470963552 * (pixsize1 / 1000.0)) / (focal_length * binx);
1596 }
1597
1598 scaleOK = true;
1599 }
1600
1601 return (coordOK || scaleOK);
1602 }
1603
findStars(StarAlgorithm algorithm,const QRect & trackingBox)1604 QFuture<bool> FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox)
1605 {
1606 if (m_StarFindFuture.isRunning())
1607 m_StarFindFuture.waitForFinished();
1608
1609 starAlgorithm = algorithm;
1610 qDeleteAll(starCenters);
1611 starCenters.clear();
1612 starsSearched = true;
1613
1614 switch (algorithm)
1615 {
1616 case ALGORITHM_SEP:
1617 {
1618 QPointer<FITSSEPDetector> detector = new FITSSEPDetector(this);
1619 detector->setSettings(m_SourceExtractorSettings);
1620 if (m_Mode == FITS_NORMAL && trackingBox.isNull())
1621 {
1622 if (Options::quickHFR())
1623 {
1624 //Just finds stars in the center 25% of the image.
1625 const int w = getStatistics().width;
1626 const int h = getStatistics().height;
1627 QRect middle(static_cast<int>(w * 0.25), static_cast<int>(h * 0.25), w / 2, h / 2);
1628 return detector->findSources(middle);
1629 }
1630 }
1631 m_StarFindFuture = detector->findSources(trackingBox);
1632 return m_StarFindFuture;
1633 }
1634
1635 case ALGORITHM_GRADIENT:
1636 default:
1637 {
1638 QPointer<FITSGradientDetector> detector = new FITSGradientDetector(this);
1639 detector->setSettings(m_SourceExtractorSettings);
1640 m_StarFindFuture = detector->findSources(trackingBox);
1641 return m_StarFindFuture;
1642 }
1643
1644 case ALGORITHM_CENTROID:
1645 {
1646 #ifndef KSTARS_LITE
1647 QPointer<FITSCentroidDetector> detector = new FITSCentroidDetector(this);
1648 detector->setSettings(m_SourceExtractorSettings);
1649 // We need JMIndex calculated from histogram
1650 if (!isHistogramConstructed())
1651 constructHistogram();
1652 detector->configure("JMINDEX", m_JMIndex);
1653 m_StarFindFuture = detector->findSources(trackingBox);
1654 return m_StarFindFuture;
1655 }
1656 #else
1657 {
1658 QPointer<FITSCentroidDetector> detector = new FITSCentroidDetector(this);
1659 return detector->findSources(starCenters, trackingBox);
1660 }
1661 #endif
1662
1663 case ALGORITHM_THRESHOLD:
1664 {
1665 QPointer<FITSThresholdDetector> detector = new FITSThresholdDetector(this);
1666 detector->setSettings(m_SourceExtractorSettings);
1667 detector->configure("THRESHOLD_PERCENTAGE", Options::focusThreshold());
1668 m_StarFindFuture = detector->findSources(trackingBox);
1669 return m_StarFindFuture;
1670 }
1671
1672 case ALGORITHM_BAHTINOV:
1673 {
1674 QPointer<FITSBahtinovDetector> detector = new FITSBahtinovDetector(this);
1675 detector->setSettings(m_SourceExtractorSettings);
1676 detector->configure("NUMBER_OF_AVERAGE_ROWS", Options::focusMultiRowAverage());
1677 m_StarFindFuture = detector->findSources(trackingBox);
1678 return m_StarFindFuture;
1679 }
1680 }
1681 }
1682
filterStars(const float innerRadius,const float outerRadius)1683 int FITSData::filterStars(const float innerRadius, const float outerRadius)
1684 {
1685 long const sqDiagonal = (long) this->width() * (long) this->width() / 4 + (long) this->height() * (long) this->height() / 4;
1686 long const sqInnerRadius = std::lround(sqDiagonal * innerRadius * innerRadius);
1687 long const sqOuterRadius = std::lround(sqDiagonal * outerRadius * outerRadius);
1688
1689 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(),
1690 [&](Edge * edge)
1691 {
1692 long const x = edge->x - this->width() / 2;
1693 long const y = edge->y - this->height() / 2;
1694 long const sqRadius = x * x + y * y;
1695 return sqRadius < sqInnerRadius || sqOuterRadius < sqRadius;
1696 }), starCenters.end());
1697
1698 return starCenters.count();
1699 }
1700
getHFR(HFRType type)1701 double FITSData::getHFR(HFRType type)
1702 {
1703 if (starCenters.empty())
1704 return -1;
1705
1706 if (cacheHFR >= 0 && cacheHFRType == type)
1707 return cacheHFR;
1708
1709 m_SelectedHFRStar = nullptr;
1710
1711 if (type == HFR_MAX)
1712 {
1713
1714 int maxVal = 0;
1715 int maxIndex = 0;
1716 for (int i = 0; i < starCenters.count(); i++)
1717 {
1718 if (starCenters[i]->HFR > maxVal)
1719 {
1720 maxIndex = i;
1721 maxVal = starCenters[i]->HFR;
1722 }
1723 }
1724
1725 m_SelectedHFRStar = starCenters[maxIndex];
1726 cacheHFR = starCenters[maxIndex]->HFR;
1727 cacheHFRType = type;
1728 return cacheHFR;
1729 }
1730 else if (type == HFR_HIGH)
1731 {
1732 // Reject all stars within 10% of border
1733 int minX = width() / 10;
1734 int minY = height() / 10;
1735 int maxX = width() - minX;
1736 int maxY = height() - minY;
1737 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(), [minX, minY, maxX, maxY](Edge * oneStar)
1738 {
1739 return (oneStar->x < minX || oneStar->x > maxX || oneStar->y < minY || oneStar->y > maxY);
1740 }), starCenters.end());
1741 // Top 5%
1742 if (starCenters.empty())
1743 return -1;
1744
1745 m_SelectedHFRStar = starCenters[static_cast<int>(starCenters.size() * 0.05)];
1746 cacheHFR = m_SelectedHFRStar->HFR;
1747 cacheHFRType = type;
1748 return cacheHFR;
1749 }
1750 else if (type == HFR_MEDIAN)
1751 {
1752 std::nth_element(starCenters.begin(), starCenters.begin() + starCenters.size() / 2, starCenters.end());
1753 m_SelectedHFRStar = starCenters[starCenters.size() / 2];
1754
1755 cacheHFR = m_SelectedHFRStar->HFR;
1756 cacheHFRType = type;
1757 return cacheHFR;
1758 }
1759
1760 // We may remove saturated stars from the HFR calculation, if we have enough stars.
1761 // No real way to tell the scale, so only remove saturated stars with range 0 -> 2**16
1762 // for non-byte types. Unsigned types and floating types, or really any pixels whose
1763 // range isn't 0-64 (or 0-255 for bytes) won't have their saturated stars removed.
1764 const int saturationValue = m_Statistics.dataType == TBYTE ? 250 : 50000;
1765 int numSaturated = 0;
1766 for (auto center : starCenters)
1767 if (center->val > saturationValue)
1768 numSaturated++;
1769 const bool removeSaturatedStars = starCenters.size() - numSaturated > 20;
1770 if (removeSaturatedStars && numSaturated > 0)
1771 qCDebug(KSTARS_FITS) << "Removing " << numSaturated << " stars from HFR calculation";
1772
1773 QVector<double> HFRs;
1774 for (auto center : starCenters)
1775 {
1776 if (removeSaturatedStars && center->val > saturationValue) continue;
1777 HFRs << center->HFR;
1778 }
1779 std::sort(HFRs.begin(), HFRs.end());
1780
1781 double sum = std::accumulate(HFRs.begin(), HFRs.end(), 0.0);
1782 double m = sum / HFRs.size();
1783
1784 if (HFRs.size() > 3)
1785 {
1786 double accum = 0.0;
1787 std::for_each (HFRs.begin(), HFRs.end(), [&](const double d)
1788 {
1789 accum += (d - m) * (d - m);
1790 });
1791 double stddev = sqrt(accum / (HFRs.size() - 1));
1792
1793 // Remove stars over 2 standard deviations away.
1794 auto end1 = std::remove_if(HFRs.begin(), HFRs.end(), [m, stddev](const double hfr)
1795 {
1796 return hfr > (m + stddev * 2);
1797 });
1798 auto end2 = std::remove_if(HFRs.begin(), end1, [m, stddev](const double hfr)
1799 {
1800 return hfr < (m - stddev * 2);
1801 });
1802
1803 // New mean
1804 sum = std::accumulate(HFRs.begin(), end2, 0.0);
1805 const int num_remaining = std::distance(HFRs.begin(), end2);
1806 if (num_remaining > 0) m = sum / num_remaining;
1807 }
1808
1809 cacheHFR = m;
1810 cacheHFRType = HFR_AVERAGE;
1811 return cacheHFR;
1812 }
1813
getHFR(int x,int y)1814 double FITSData::getHFR(int x, int y)
1815 {
1816 if (starCenters.empty())
1817 return -1;
1818
1819 for (int i = 0; i < starCenters.count(); i++)
1820 {
1821 if (std::fabs(starCenters[i]->x - x) <= starCenters[i]->width / 2 &&
1822 std::fabs(starCenters[i]->y - y) <= starCenters[i]->width / 2)
1823 {
1824 m_SelectedHFRStar = starCenters[i];
1825 return starCenters[i]->HFR;
1826 }
1827 }
1828
1829 return -1;
1830 }
1831
getEccentricity()1832 double FITSData::getEccentricity()
1833 {
1834 if (starCenters.empty())
1835 return -1;
1836 if (cacheEccentricity >= 0)
1837 return cacheEccentricity;
1838 std::vector<float> eccs;
1839 for (const auto &s : starCenters)
1840 eccs.push_back(s->ellipticity);
1841 int middle = eccs.size() / 2;
1842 std::nth_element(eccs.begin(), eccs.begin() + middle, eccs.end());
1843 float medianEllipticity = eccs[middle];
1844
1845 // SEP gives us the ellipticity (flattening) directly.
1846 // To get the eccentricity, use this formula:
1847 // e = sqrt(ellipticity * (2 - ellipticity))
1848 // https://en.wikipedia.org/wiki/Eccentricity_(mathematics)#Ellipses
1849 const float eccentricity = sqrt(medianEllipticity * (2 - medianEllipticity));
1850 cacheEccentricity = eccentricity;
1851 return eccentricity;
1852 }
1853
applyFilter(FITSScale type,uint8_t * image,QVector<double> * min,QVector<double> * max)1854 void FITSData::applyFilter(FITSScale type, uint8_t * image, QVector<double> * min, QVector<double> * max)
1855 {
1856 if (type == FITS_NONE)
1857 return;
1858
1859 QVector<double> dataMin(3);
1860 QVector<double> dataMax(3);
1861
1862 if (min)
1863 dataMin = *min;
1864 if (max)
1865 dataMax = *max;
1866
1867 switch (type)
1868 {
1869 case FITS_AUTO_STRETCH:
1870 {
1871 for (int i = 0; i < 3; i++)
1872 {
1873 dataMin[i] = m_Statistics.mean[i] - m_Statistics.stddev[i];
1874 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
1875 }
1876 }
1877 break;
1878
1879 case FITS_HIGH_CONTRAST:
1880 {
1881 for (int i = 0; i < 3; i++)
1882 {
1883 dataMin[i] = m_Statistics.mean[i] + m_Statistics.stddev[i];
1884 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
1885 }
1886 }
1887 break;
1888
1889 case FITS_HIGH_PASS:
1890 {
1891 for (int i = 0; i < 3; i++)
1892 {
1893 dataMin[i] = m_Statistics.mean[i];
1894 }
1895 }
1896 break;
1897
1898 default:
1899 break;
1900 }
1901
1902 switch (m_Statistics.dataType)
1903 {
1904 case TBYTE:
1905 {
1906 for (int i = 0; i < 3; i++)
1907 {
1908 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1909 dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i];
1910 }
1911 applyFilter<uint8_t>(type, image, &dataMin, &dataMax);
1912 }
1913 break;
1914
1915 case TSHORT:
1916 {
1917 for (int i = 0; i < 3; i++)
1918 {
1919 dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i];
1920 dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i];
1921 }
1922 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1923 }
1924
1925 break;
1926
1927 case TUSHORT:
1928 {
1929 for (int i = 0; i < 3; i++)
1930 {
1931 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1932 dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i];
1933 }
1934 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1935 }
1936 break;
1937
1938 case TLONG:
1939 {
1940 for (int i = 0; i < 3; i++)
1941 {
1942 dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i];
1943 dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i];
1944 }
1945 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1946 }
1947 break;
1948
1949 case TULONG:
1950 {
1951 for (int i = 0; i < 3; i++)
1952 {
1953 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1954 dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i];
1955 }
1956 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1957 }
1958 break;
1959
1960 case TFLOAT:
1961 {
1962 for (int i = 0; i < 3; i++)
1963 {
1964 dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i];
1965 dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i];
1966 }
1967 applyFilter<float>(type, image, &dataMin, &dataMax);
1968 }
1969 break;
1970
1971 case TLONGLONG:
1972 {
1973 for (int i = 0; i < 3; i++)
1974 {
1975 dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i];
1976 dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i];
1977 }
1978
1979 applyFilter<long>(type, image, &dataMin, &dataMax);
1980 }
1981 break;
1982
1983 case TDOUBLE:
1984 {
1985 for (int i = 0; i < 3; i++)
1986 {
1987 dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i];
1988 dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i];
1989 }
1990 applyFilter<double>(type, image, &dataMin, &dataMax);
1991 }
1992
1993 break;
1994
1995 default:
1996 return;
1997 }
1998
1999 if (min != nullptr)
2000 *min = dataMin;
2001 if (max != nullptr)
2002 *max = dataMax;
2003
2004 emit dataChanged();
2005 }
2006
2007 template <typename T>
applyFilter(FITSScale type,uint8_t * targetImage,QVector<double> * targetMin,QVector<double> * targetMax)2008 void FITSData::applyFilter(FITSScale type, uint8_t * targetImage, QVector<double> * targetMin, QVector<double> * targetMax)
2009 {
2010 bool calcStats = false;
2011 T * image = nullptr;
2012
2013 if (targetImage)
2014 image = reinterpret_cast<T *>(targetImage);
2015 else
2016 {
2017 image = reinterpret_cast<T *>(m_ImageBuffer);
2018 calcStats = true;
2019 }
2020
2021 T min[3], max[3];
2022 for (int i = 0; i < 3; i++)
2023 {
2024 min[i] = (*targetMin)[i] < std::numeric_limits<T>::min() ? std::numeric_limits<T>::min() : (*targetMin)[i];
2025 max[i] = (*targetMax)[i] > std::numeric_limits<T>::max() ? std::numeric_limits<T>::max() : (*targetMax)[i];
2026 }
2027
2028
2029 // Create N threads
2030 const uint8_t nThreads = 16;
2031
2032 uint32_t width = m_Statistics.width;
2033 uint32_t height = m_Statistics.height;
2034
2035 //QTime timer;
2036 //timer.start();
2037 switch (type)
2038 {
2039 case FITS_AUTO:
2040 case FITS_LINEAR:
2041 case FITS_AUTO_STRETCH:
2042 case FITS_HIGH_CONTRAST:
2043 case FITS_LOG:
2044 case FITS_SQRT:
2045 case FITS_HIGH_PASS:
2046 {
2047 // List of futures
2048 QList<QFuture<void>> futures;
2049 QVector<double> coeff(3);
2050
2051 if (type == FITS_LOG)
2052 {
2053 for (int i = 0; i < 3; i++)
2054 coeff[i] = max[i] / std::log(1 + max[i]);
2055 }
2056 else if (type == FITS_SQRT)
2057 {
2058 for (int i = 0; i < 3; i++)
2059 coeff[i] = max[i] / sqrt(max[i]);
2060 }
2061
2062 for (int n = 0; n < m_Statistics.channels; n++)
2063 {
2064 if (type == FITS_HIGH_PASS)
2065 min[n] = m_Statistics.mean[n];
2066
2067 uint32_t cStart = n * m_Statistics.samples_per_channel;
2068
2069 // Calculate how many elements we process per thread
2070 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
2071
2072 // Calculate the final stride since we can have some left over due to division above
2073 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
2074
2075 T * runningBuffer = image + cStart;
2076
2077 if (type == FITS_LOG)
2078 {
2079 for (int i = 0; i < nThreads; i++)
2080 {
2081 // Run threads
2082 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2083 coeff, n](T & a)
2084 {
2085 a = qBound(min[n], static_cast<T>(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]);
2086 }));
2087
2088 runningBuffer += tStride;
2089 }
2090 }
2091 else if (type == FITS_SQRT)
2092 {
2093 for (int i = 0; i < nThreads; i++)
2094 {
2095 // Run threads
2096 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2097 coeff, n](T & a)
2098 {
2099 a = qBound(min[n], static_cast<T>(round(coeff[n] * a)), max[n]);
2100 }));
2101 }
2102
2103 runningBuffer += tStride;
2104 }
2105 else
2106 {
2107 for (int i = 0; i < nThreads; i++)
2108 {
2109 // Run threads
2110 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2111 n](T & a)
2112 {
2113 a = qBound(min[n], a, max[n]);
2114 }));
2115 runningBuffer += tStride;
2116 }
2117 }
2118 }
2119
2120 for (int i = 0; i < nThreads * m_Statistics.channels; i++)
2121 futures[i].waitForFinished();
2122
2123 if (calcStats)
2124 {
2125 for (int i = 0; i < 3; i++)
2126 {
2127 m_Statistics.min[i] = min[i];
2128 m_Statistics.max[i] = max[i];
2129 }
2130 //if (type != FITS_AUTO && type != FITS_LINEAR)
2131 runningAverageStdDev<T>();
2132 //QtConcurrent::run(this, &FITSData::runningAverageStdDev<T>);
2133 }
2134 }
2135 break;
2136
2137 case FITS_EQUALIZE:
2138 {
2139 #ifndef KSTARS_LITE
2140 if (!isHistogramConstructed())
2141 constructHistogram();
2142
2143 T bufferVal = 0;
2144
2145 double coeff = 255.0 / (height * width);
2146 uint32_t row = 0;
2147 uint32_t index = 0;
2148
2149 for (int i = 0; i < m_Statistics.channels; i++)
2150 {
2151 uint32_t offset = i * m_Statistics.samples_per_channel;
2152 for (uint32_t j = 0; j < height; j++)
2153 {
2154 row = offset + j * width;
2155 for (uint32_t k = 0; k < width; k++)
2156 {
2157 index = k + row;
2158 bufferVal = (image[index] - min[i]) / m_HistogramBinWidth[i];
2159
2160 if (bufferVal >= m_CumulativeFrequency[i].size())
2161 bufferVal = m_CumulativeFrequency[i].size() - 1;
2162
2163 image[index] = qBound(min[i], static_cast<T>(round(coeff * m_CumulativeFrequency[i][bufferVal])), max[i]);
2164 }
2165 }
2166 }
2167 #endif
2168 }
2169 if (calcStats)
2170 calculateStats(true);
2171 break;
2172
2173 // Based on http://www.librow.com/articles/article-1
2174 case FITS_MEDIAN:
2175 {
2176 uint8_t BBP = m_Statistics.bytesPerPixel;
2177 auto * extension = new T[(width + 2) * (height + 2)];
2178 // Check memory allocation
2179 if (!extension)
2180 return;
2181 // Create image extension
2182 for (uint32_t ch = 0; ch < m_Statistics.channels; ch++)
2183 {
2184 uint32_t offset = ch * m_Statistics.samples_per_channel;
2185 uint32_t N = width, M = height;
2186
2187 for (uint32_t i = 0; i < M; ++i)
2188 {
2189 memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP);
2190 extension[(N + 2) * (i + 1)] = image[N * i + offset];
2191 extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset];
2192 }
2193 // Fill first line of image extension
2194 memcpy(extension, extension + N + 2, (N + 2) * BBP);
2195 // Fill last line of image extension
2196 memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP);
2197 // Call median filter implementation
2198
2199 N = width + 2;
2200 M = height + 2;
2201
2202 // Move window through all elements of the image
2203 for (uint32_t m = 1; m < M - 1; ++m)
2204 for (uint32_t n = 1; n < N - 1; ++n)
2205 {
2206 // Pick up window elements
2207 int k = 0;
2208 float window[9];
2209
2210 memset(&window[0], 0, 9 * sizeof(float));
2211 for (uint32_t j = m - 1; j < m + 2; ++j)
2212 for (uint32_t i = n - 1; i < n + 2; ++i)
2213 window[k++] = extension[j * N + i];
2214 // Order elements (only half of them)
2215 for (uint32_t j = 0; j < 5; ++j)
2216 {
2217 // Find position of minimum element
2218 int mine = j;
2219 for (uint32_t l = j + 1; l < 9; ++l)
2220 if (window[l] < window[mine])
2221 mine = l;
2222 // Put found minimum element in its place
2223 const float temp = window[j];
2224 window[j] = window[mine];
2225 window[mine] = temp;
2226 }
2227 // Get result - the middle element
2228 image[(m - 1) * (N - 2) + n - 1 + offset] = window[4];
2229 }
2230 }
2231
2232 // Free memory
2233 delete[] extension;
2234
2235 if (calcStats)
2236 runningAverageStdDev<T>();
2237 }
2238 break;
2239
2240 case FITS_GAUSSIAN:
2241 gaussianBlur<T>(Options::focusGaussianKernelSize(), Options::focusGaussianSigma());
2242 if (calcStats)
2243 calculateStats(true);
2244 break;
2245
2246 case FITS_ROTATE_CW:
2247 rotFITS<T>(90, 0);
2248 rotCounter++;
2249 break;
2250
2251 case FITS_ROTATE_CCW:
2252 rotFITS<T>(270, 0);
2253 rotCounter--;
2254 break;
2255
2256 case FITS_FLIP_H:
2257 rotFITS<T>(0, 1);
2258 flipHCounter++;
2259 break;
2260
2261 case FITS_FLIP_V:
2262 rotFITS<T>(0, 2);
2263 flipVCounter++;
2264 break;
2265
2266 default:
2267 break;
2268 }
2269 }
2270
getStarCentersInSubFrame(QRect subFrame) const2271 QList<Edge *> FITSData::getStarCentersInSubFrame(QRect subFrame) const
2272 {
2273 QList<Edge *> starCentersInSubFrame;
2274 for (int i = 0; i < starCenters.count(); i++)
2275 {
2276 int x = static_cast<int>(starCenters[i]->x);
2277 int y = static_cast<int>(starCenters[i]->y);
2278 if(subFrame.contains(x, y))
2279 {
2280 starCentersInSubFrame.append(starCenters[i]);
2281 }
2282 }
2283 return starCentersInSubFrame;
2284 }
2285
checkForWCS()2286 bool FITSData::checkForWCS()
2287 {
2288 #ifndef KSTARS_LITE
2289 #ifdef HAVE_WCSLIB
2290
2291 int status = 0;
2292 char * header;
2293 int nkeyrec, nreject;
2294
2295 // Free wcs before re-use
2296 if (m_WCSHandle != nullptr)
2297 {
2298 wcsvfree(&m_nwcs, &m_WCSHandle);
2299 m_WCSHandle = nullptr;
2300 }
2301
2302 if (fits_hdr2str(fptr, 1, nullptr, 0, &header, &nkeyrec, &status))
2303 {
2304 char errmsg[512];
2305 fits_get_errstatus(status, errmsg);
2306 m_LastError = errmsg;
2307 return false;
2308 }
2309
2310 if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2311 {
2312 free(header);
2313 wcsvfree(&m_nwcs, &m_WCSHandle);
2314 m_LastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]);
2315 return false;
2316 }
2317
2318 free(header);
2319
2320 if (m_WCSHandle == nullptr)
2321 {
2322 m_LastError = i18n("No world coordinate systems found.");
2323 return false;
2324 }
2325
2326 // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now.
2327 if (m_WCSHandle->crpix[0] == 0)
2328 {
2329 wcsvfree(&m_nwcs, &m_WCSHandle);
2330 m_WCSHandle = nullptr;
2331 m_LastError = i18n("No world coordinate systems found.");
2332 return false;
2333 }
2334
2335 cdfix(m_WCSHandle);
2336 if ((status = wcsset(m_WCSHandle)) != 0)
2337 {
2338 wcsvfree(&m_nwcs, &m_WCSHandle);
2339 m_WCSHandle = nullptr;
2340 m_LastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2341 return false;
2342 }
2343
2344 HasWCS = true;
2345 #endif
2346 #endif
2347 return HasWCS;
2348 }
2349
loadWCS(bool extras)2350 bool FITSData::loadWCS(bool extras)
2351 {
2352 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2353
2354 if (m_WCSState == Success)
2355 {
2356 qCWarning(KSTARS_FITS) << "WCS data already loaded";
2357 return true;
2358 }
2359
2360 if (m_WCSHandle != nullptr)
2361 {
2362 wcsvfree(&m_nwcs, &m_WCSHandle);
2363 m_WCSHandle = nullptr;
2364 }
2365
2366 qCDebug(KSTARS_FITS) << "Started WCS Data Processing...";
2367
2368 int status = 0;
2369 char * header;
2370 int nkeyrec, nreject, nwcs;
2371
2372 if (fits_hdr2str(fptr, 1, nullptr, 0, &header, &nkeyrec, &status))
2373 {
2374 char errmsg[512];
2375 fits_get_errstatus(status, errmsg);
2376 m_LastError = errmsg;
2377 m_WCSState = Failure;
2378 return false;
2379 }
2380
2381 if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &m_WCSHandle)) != 0)
2382 {
2383 free(header);
2384 wcsvfree(&m_nwcs, &m_WCSHandle);
2385 m_WCSHandle = nullptr;
2386 m_LastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]);
2387 m_WCSState = Failure;
2388 return false;
2389 }
2390
2391 free(header);
2392
2393 if (m_WCSHandle == nullptr)
2394 {
2395 m_LastError = i18n("No world coordinate systems found.");
2396 m_WCSState = Failure;
2397 return false;
2398 }
2399
2400 // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now.
2401 if (m_WCSHandle->crpix[0] == 0)
2402 {
2403 wcsvfree(&m_nwcs, &m_WCSHandle);
2404 m_WCSHandle = nullptr;
2405 m_LastError = i18n("No world coordinate systems found.");
2406 m_WCSState = Failure;
2407 return false;
2408 }
2409
2410 cdfix(m_WCSHandle);
2411 if ((status = wcsset(m_WCSHandle)) != 0)
2412 {
2413 wcsvfree(&m_nwcs, &m_WCSHandle);
2414 m_WCSHandle = nullptr;
2415 m_LastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2416 m_WCSState = Failure;
2417 return false;
2418 }
2419
2420 m_ObjectsSearched = false;
2421 m_WCSState = Success;
2422 FullWCS = extras;
2423 HasWCS = true;
2424
2425 qCDebug(KSTARS_FITS) << "Finished WCS Data processing...";
2426
2427 return true;
2428 #else
2429 return false;
2430 #endif
2431 }
2432
wcsToPixel(const SkyPoint & wcsCoord,QPointF & wcsPixelPoint,QPointF & wcsImagePoint)2433 bool FITSData::wcsToPixel(const SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint)
2434 {
2435 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2436 int status, stat[NWCSFIX];
2437 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2438
2439 if (m_WCSHandle == nullptr)
2440 {
2441 m_LastError = i18n("No world coordinate systems found.");
2442 return false;
2443 }
2444
2445 worldcrd[0] = wcsCoord.ra0().Degrees();
2446 worldcrd[1] = wcsCoord.dec0().Degrees();
2447
2448 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2449 {
2450 m_LastError = QString("wcss2p error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2451 return false;
2452 }
2453
2454 wcsImagePoint.setX(imgcrd[0]);
2455 wcsImagePoint.setY(imgcrd[1]);
2456
2457 wcsPixelPoint.setX(pixcrd[0]);
2458 wcsPixelPoint.setY(pixcrd[1]);
2459
2460 return true;
2461 #else
2462 Q_UNUSED(wcsCoord);
2463 Q_UNUSED(wcsPixelPoint);
2464 Q_UNUSED(wcsImagePoint);
2465 return false;
2466 #endif
2467 }
2468
pixelToWCS(const QPointF & wcsPixelPoint,SkyPoint & wcsCoord)2469 bool FITSData::pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord)
2470 {
2471 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2472 int status, stat[NWCSFIX];
2473 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2474
2475 if (m_WCSHandle == nullptr)
2476 {
2477 m_LastError = i18n("No world coordinate systems found.");
2478 return false;
2479 }
2480
2481 pixcrd[0] = wcsPixelPoint.x();
2482 pixcrd[1] = wcsPixelPoint.y();
2483
2484 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2485 {
2486 m_LastError = QString("wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2487 return false;
2488 }
2489 else
2490 {
2491 wcsCoord.setRA0(world[0] / 15.0);
2492 wcsCoord.setDec0(world[1]);
2493 }
2494
2495 return true;
2496 #else
2497 Q_UNUSED(wcsPixelPoint);
2498 Q_UNUSED(wcsCoord);
2499 return false;
2500 #endif
2501 }
2502
2503 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
searchObjects()2504 bool FITSData::searchObjects()
2505 {
2506 if (m_ObjectsSearched)
2507 return true;
2508
2509 m_ObjectsSearched = true;
2510
2511 SkyPoint startPoint;
2512 SkyPoint endPoint;
2513
2514 pixelToWCS(QPointF(0, 0), startPoint);
2515 pixelToWCS(QPointF(width() - 1, height() - 1), endPoint);
2516
2517 return findObjectsInImage(startPoint, endPoint);
2518 }
2519
findWCSBounds(double & minRA,double & maxRA,double & minDec,double & maxDec)2520 bool FITSData::findWCSBounds(double &minRA, double &maxRA, double &minDec, double &maxDec)
2521 {
2522 if (m_WCSHandle == nullptr)
2523 {
2524 m_LastError = i18n("No world coordinate systems found.");
2525 return false;
2526 }
2527
2528 maxRA = -1000;
2529 minRA = 1000;
2530 maxDec = -1000;
2531 minDec = 1000;
2532
2533 auto updateMinMax = [&](int x, int y)
2534 {
2535 int stat[2];
2536 double imgcrd[2], phi, pixcrd[2], theta, world[2];
2537
2538 pixcrd[0] = x;
2539 pixcrd[1] = y;
2540
2541 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
2542 return;
2543
2544 minRA = std::min(minRA, world[0]);
2545 maxRA = std::max(maxRA, world[0]);
2546 minDec = std::min(minDec, world[1]);
2547 maxDec = std::max(maxDec, world[1]);
2548 };
2549
2550 // Find min and max values from edges
2551 for (int y = 0; y < height(); y++)
2552 {
2553 updateMinMax(0, y);
2554 updateMinMax(width() - 1, y);
2555 }
2556
2557 for (int x = 1; x < width() - 1; x++)
2558 {
2559 updateMinMax(x, 0);
2560 updateMinMax(x, height() - 1);
2561 }
2562
2563 // Check if either pole is in the image
2564 SkyPoint NCP(0, 90);
2565 SkyPoint SCP(0, -90);
2566 QPointF pixelPoint, imagePoint, pPoint;
2567 if (wcsToPixel(NCP, pPoint, imagePoint))
2568 {
2569 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
2570 {
2571 maxDec = 90;
2572 }
2573 }
2574 if (wcsToPixel(SCP, pPoint, imagePoint))
2575 {
2576 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
2577 {
2578 minDec = -90;
2579 }
2580 }
2581 return true;
2582 }
2583 #endif
2584
2585 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
findObjectsInImage(SkyPoint startPoint,SkyPoint endPoint)2586 bool FITSData::findObjectsInImage(SkyPoint startPoint, SkyPoint endPoint)
2587 {
2588 if (KStarsData::Instance() == nullptr)
2589 return false;
2590
2591 int w = width();
2592 int h = height();
2593 QVariant date;
2594 KSNumbers * num = nullptr;
2595
2596 if (getRecordValue("DATE-OBS", date))
2597 {
2598 QString tsString(date.toString());
2599 tsString = tsString.remove('\'').trimmed();
2600 // Add Zulu time to indicate UTC
2601 tsString += "Z";
2602
2603 QDateTime ts = QDateTime::fromString(tsString, Qt::ISODate);
2604
2605 if (ts.isValid())
2606 num = new KSNumbers(KStarsDateTime(ts).djd());
2607 }
2608
2609 //Set to current time if the above does not work.
2610 if (num == nullptr)
2611 num = new KSNumbers(KStarsData::Instance()->ut().djd());
2612
2613 startPoint.updateCoordsNow(num);
2614 endPoint.updateCoordsNow(num);
2615
2616 m_SkyObjects.clear();
2617
2618 QList<SkyObject *> list = KStarsData::Instance()->skyComposite()->findObjectsInArea(startPoint, endPoint);
2619 list.erase(std::remove_if(list.begin(), list.end(), [](SkyObject * oneObject)
2620 {
2621 int type = oneObject->type();
2622 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
2623 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
2624 type == SkyObject::SATELLITE);
2625 }), list.end());
2626
2627 double world[2], phi, theta, imgcrd[2], pixcrd[2];
2628 int stat[2];
2629 for (auto &object : list)
2630 {
2631 world[0] = object->ra0().Degrees();
2632 world[1] = object->dec0().Degrees();
2633
2634 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
2635 {
2636 //The X and Y are set to the found position if it does work.
2637 int x = pixcrd[0];
2638 int y = pixcrd[1];
2639 if (x > 0 && y > 0 && x < w && y < h)
2640 m_SkyObjects.append(new FITSSkyObject(object, x, y));
2641 }
2642 }
2643
2644 delete (num);
2645 return true;
2646 }
2647 #endif
2648
getFlipVCounter() const2649 int FITSData::getFlipVCounter() const
2650 {
2651 return flipVCounter;
2652 }
2653
setFlipVCounter(int value)2654 void FITSData::setFlipVCounter(int value)
2655 {
2656 flipVCounter = value;
2657 }
2658
getFlipHCounter() const2659 int FITSData::getFlipHCounter() const
2660 {
2661 return flipHCounter;
2662 }
2663
setFlipHCounter(int value)2664 void FITSData::setFlipHCounter(int value)
2665 {
2666 flipHCounter = value;
2667 }
2668
getRotCounter() const2669 int FITSData::getRotCounter() const
2670 {
2671 return rotCounter;
2672 }
2673
setRotCounter(int value)2674 void FITSData::setRotCounter(int value)
2675 {
2676 rotCounter = value;
2677 }
2678
2679 /* Rotate an image by 90, 180, or 270 degrees, with an optional
2680 * reflection across the vertical or horizontal axis.
2681 * verbose generates extra info on stdout.
2682 * return nullptr if successful or rotated image.
2683 */
2684 template <typename T>
rotFITS(int rotate,int mirror)2685 bool FITSData::rotFITS(int rotate, int mirror)
2686 {
2687 int ny, nx;
2688 int x1, y1, x2, y2;
2689 uint8_t * rotimage = nullptr;
2690 int offset = 0;
2691
2692 if (rotate == 1)
2693 rotate = 90;
2694 else if (rotate == 2)
2695 rotate = 180;
2696 else if (rotate == 3)
2697 rotate = 270;
2698 else if (rotate < 0)
2699 rotate = rotate + 360;
2700
2701 nx = m_Statistics.width;
2702 ny = m_Statistics.height;
2703
2704 int BBP = m_Statistics.bytesPerPixel;
2705
2706 /* Allocate buffer for rotated image */
2707 rotimage = new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
2708
2709 if (rotimage == nullptr)
2710 {
2711 qWarning() << "Unable to allocate memory for rotated image buffer!";
2712 return false;
2713 }
2714
2715 auto * rotBuffer = reinterpret_cast<T *>(rotimage);
2716 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
2717
2718 /* Mirror image without rotation */
2719 if (rotate < 45 && rotate > -45)
2720 {
2721 if (mirror == 1)
2722 {
2723 for (int i = 0; i < m_Statistics.channels; i++)
2724 {
2725 offset = m_Statistics.samples_per_channel * i;
2726 for (x1 = 0; x1 < nx; x1++)
2727 {
2728 x2 = nx - x1 - 1;
2729 for (y1 = 0; y1 < ny; y1++)
2730 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2731 }
2732 }
2733 }
2734 else if (mirror == 2)
2735 {
2736 for (int i = 0; i < m_Statistics.channels; i++)
2737 {
2738 offset = m_Statistics.samples_per_channel * i;
2739 for (y1 = 0; y1 < ny; y1++)
2740 {
2741 y2 = ny - y1 - 1;
2742 for (x1 = 0; x1 < nx; x1++)
2743 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2744 }
2745 }
2746 }
2747 else
2748 {
2749 for (int i = 0; i < m_Statistics.channels; i++)
2750 {
2751 offset = m_Statistics.samples_per_channel * i;
2752 for (y1 = 0; y1 < ny; y1++)
2753 {
2754 for (x1 = 0; x1 < nx; x1++)
2755 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2756 }
2757 }
2758 }
2759 }
2760
2761 /* Rotate by 90 degrees */
2762 else if (rotate >= 45 && rotate < 135)
2763 {
2764 if (mirror == 1)
2765 {
2766 for (int i = 0; i < m_Statistics.channels; i++)
2767 {
2768 offset = m_Statistics.samples_per_channel * i;
2769 for (y1 = 0; y1 < ny; y1++)
2770 {
2771 x2 = ny - y1 - 1;
2772 for (x1 = 0; x1 < nx; x1++)
2773 {
2774 y2 = nx - x1 - 1;
2775 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2776 }
2777 }
2778 }
2779 }
2780 else if (mirror == 2)
2781 {
2782 for (int i = 0; i < m_Statistics.channels; i++)
2783 {
2784 offset = m_Statistics.samples_per_channel * i;
2785 for (y1 = 0; y1 < ny; y1++)
2786 {
2787 for (x1 = 0; x1 < nx; x1++)
2788 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
2789 }
2790 }
2791 }
2792 else
2793 {
2794 for (int i = 0; i < m_Statistics.channels; i++)
2795 {
2796 offset = m_Statistics.samples_per_channel * i;
2797 for (y1 = 0; y1 < ny; y1++)
2798 {
2799 x2 = ny - y1 - 1;
2800 for (x1 = 0; x1 < nx; x1++)
2801 {
2802 y2 = x1;
2803 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2804 }
2805 }
2806 }
2807 }
2808
2809 m_Statistics.width = ny;
2810 m_Statistics.height = nx;
2811 }
2812
2813 /* Rotate by 180 degrees */
2814 else if (rotate >= 135 && rotate < 225)
2815 {
2816 if (mirror == 1)
2817 {
2818 for (int i = 0; i < m_Statistics.channels; i++)
2819 {
2820 offset = m_Statistics.samples_per_channel * i;
2821 for (y1 = 0; y1 < ny; y1++)
2822 {
2823 y2 = ny - y1 - 1;
2824 for (x1 = 0; x1 < nx; x1++)
2825 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2826 }
2827 }
2828 }
2829 else if (mirror == 2)
2830 {
2831 for (int i = 0; i < m_Statistics.channels; i++)
2832 {
2833 offset = m_Statistics.samples_per_channel * i;
2834 for (x1 = 0; x1 < nx; x1++)
2835 {
2836 x2 = nx - x1 - 1;
2837 for (y1 = 0; y1 < ny; y1++)
2838 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2839 }
2840 }
2841 }
2842 else
2843 {
2844 for (int i = 0; i < m_Statistics.channels; i++)
2845 {
2846 offset = m_Statistics.samples_per_channel * i;
2847 for (y1 = 0; y1 < ny; y1++)
2848 {
2849 y2 = ny - y1 - 1;
2850 for (x1 = 0; x1 < nx; x1++)
2851 {
2852 x2 = nx - x1 - 1;
2853 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2854 }
2855 }
2856 }
2857 }
2858 }
2859
2860 /* Rotate by 270 degrees */
2861 else if (rotate >= 225 && rotate < 315)
2862 {
2863 if (mirror == 1)
2864 {
2865 for (int i = 0; i < m_Statistics.channels; i++)
2866 {
2867 offset = m_Statistics.samples_per_channel * i;
2868 for (y1 = 0; y1 < ny; y1++)
2869 {
2870 for (x1 = 0; x1 < nx; x1++)
2871 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
2872 }
2873 }
2874 }
2875 else if (mirror == 2)
2876 {
2877 for (int i = 0; i < m_Statistics.channels; i++)
2878 {
2879 offset = m_Statistics.samples_per_channel * i;
2880 for (y1 = 0; y1 < ny; y1++)
2881 {
2882 x2 = ny - y1 - 1;
2883 for (x1 = 0; x1 < nx; x1++)
2884 {
2885 y2 = nx - x1 - 1;
2886 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2887 }
2888 }
2889 }
2890 }
2891 else
2892 {
2893 for (int i = 0; i < m_Statistics.channels; i++)
2894 {
2895 offset = m_Statistics.samples_per_channel * i;
2896 for (y1 = 0; y1 < ny; y1++)
2897 {
2898 x2 = y1;
2899 for (x1 = 0; x1 < nx; x1++)
2900 {
2901 y2 = nx - x1 - 1;
2902 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2903 }
2904 }
2905 }
2906 }
2907
2908 m_Statistics.width = ny;
2909 m_Statistics.height = nx;
2910 }
2911
2912 /* If rotating by more than 315 degrees, assume top-bottom reflection */
2913 else if (rotate >= 315 && mirror)
2914 {
2915 for (int i = 0; i < m_Statistics.channels; i++)
2916 {
2917 offset = m_Statistics.samples_per_channel * i;
2918 for (y1 = 0; y1 < ny; y1++)
2919 {
2920 for (x1 = 0; x1 < nx; x1++)
2921 {
2922 x2 = y1;
2923 y2 = x1;
2924 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2925 }
2926 }
2927 }
2928 }
2929
2930 delete[] m_ImageBuffer;
2931 m_ImageBuffer = rotimage;
2932
2933 return true;
2934 }
2935
rotWCSFITS(int angle,int mirror)2936 void FITSData::rotWCSFITS(int angle, int mirror)
2937 {
2938 int status = 0;
2939 char comment[100];
2940 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
2941 int WCS_DECIMALS = 6;
2942
2943 naxis1 = m_Statistics.width;
2944 naxis2 = m_Statistics.height;
2945
2946 if (fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
2947 {
2948 // No WCS keywords
2949 return;
2950 }
2951
2952 /* Reset CROTAn and CD matrix if axes have been exchanged */
2953 if (angle == 90)
2954 {
2955 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
2956 fits_update_key_dbl(fptr, "CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
2957
2958 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
2959 fits_update_key_dbl(fptr, "CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
2960 }
2961
2962 status = 0;
2963
2964 /* Negate rotation angle if mirrored */
2965 if (mirror != 0)
2966 {
2967 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
2968 fits_update_key_dbl(fptr, "CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
2969
2970 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
2971 fits_update_key_dbl(fptr, "CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
2972
2973 status = 0;
2974
2975 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
2976 fits_update_key_dbl(fptr, "LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
2977
2978 status = 0;
2979
2980 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
2981 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
2982
2983 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp1, comment, &status))
2984 fits_update_key_dbl(fptr, "CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
2985
2986 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp1, comment, &status))
2987 fits_update_key_dbl(fptr, "CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
2988 }
2989
2990 status = 0;
2991
2992 /* Unbin CRPIX and CD matrix */
2993 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
2994 {
2995 if (ctemp1 != 1.0)
2996 {
2997 if (!fits_read_key_dbl(fptr, "LTM2_2", &ctemp2, comment, &status))
2998 if (ctemp1 == ctemp2)
2999 {
3000 double ltv1 = 0.0;
3001 double ltv2 = 0.0;
3002 status = 0;
3003 if (!fits_read_key_dbl(fptr, "LTV1", <v1, comment, &status))
3004 fits_delete_key(fptr, "LTV1", &status);
3005 if (!fits_read_key_dbl(fptr, "LTV2", <v2, comment, &status))
3006 fits_delete_key(fptr, "LTV2", &status);
3007
3008 status = 0;
3009
3010 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp3, comment, &status))
3011 fits_update_key_dbl(fptr, "CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
3012
3013 if (!fits_read_key_dbl(fptr, "CRPIX2", &ctemp3, comment, &status))
3014 fits_update_key_dbl(fptr, "CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
3015
3016 status = 0;
3017
3018 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp3, comment, &status))
3019 fits_update_key_dbl(fptr, "CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3020
3021 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp3, comment, &status))
3022 fits_update_key_dbl(fptr, "CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3023
3024 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status))
3025 fits_update_key_dbl(fptr, "CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3026
3027 if (!fits_read_key_dbl(fptr, "CD2_2", &ctemp3, comment, &status))
3028 fits_update_key_dbl(fptr, "CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3029
3030 status = 0;
3031
3032 fits_delete_key(fptr, "LTM1_1", &status);
3033 fits_delete_key(fptr, "LTM1_2", &status);
3034 }
3035 }
3036 }
3037
3038 status = 0;
3039
3040 /* Reset CRPIXn */
3041 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp1, comment, &status) &&
3042 !fits_read_key_dbl(fptr, "CRPIX2", &ctemp2, comment, &status))
3043 {
3044 if (mirror != 0)
3045 {
3046 if (angle == 0)
3047 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3048 else if (angle == 90)
3049 {
3050 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3051 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3052 }
3053 else if (angle == 180)
3054 {
3055 fits_update_key_dbl(fptr, "CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
3056 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3057 }
3058 else if (angle == 270)
3059 {
3060 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
3061 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
3062 }
3063 }
3064 else
3065 {
3066 if (angle == 90)
3067 {
3068 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3069 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
3070 }
3071 else if (angle == 180)
3072 {
3073 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3074 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3075 }
3076 else if (angle == 270)
3077 {
3078 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
3079 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3080 }
3081 }
3082 }
3083
3084 status = 0;
3085
3086 /* Reset CDELTn (degrees per pixel) */
3087 if (!fits_read_key_dbl(fptr, "CDELT1", &ctemp1, comment, &status) &&
3088 !fits_read_key_dbl(fptr, "CDELT2", &ctemp2, comment, &status))
3089 {
3090 if (mirror != 0)
3091 {
3092 if (angle == 0)
3093 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3094 else if (angle == 90)
3095 {
3096 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3097 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3098 }
3099 else if (angle == 180)
3100 {
3101 fits_update_key_dbl(fptr, "CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
3102 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3103 }
3104 else if (angle == 270)
3105 {
3106 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3107 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3108 }
3109 }
3110 else
3111 {
3112 if (angle == 90)
3113 {
3114 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3115 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3116 }
3117 else if (angle == 180)
3118 {
3119 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3120 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3121 }
3122 else if (angle == 270)
3123 {
3124 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3125 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3126 }
3127 }
3128 }
3129
3130 /* Reset CD matrix, if present */
3131 ctemp1 = 0.0;
3132 ctemp2 = 0.0;
3133 ctemp3 = 0.0;
3134 ctemp4 = 0.0;
3135 status = 0;
3136 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
3137 {
3138 fits_read_key_dbl(fptr, "CD1_2", &ctemp2, comment, &status);
3139 fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status);
3140 fits_read_key_dbl(fptr, "CD2_2", &ctemp4, comment, &status);
3141 status = 0;
3142 if (mirror != 0)
3143 {
3144 if (angle == 0)
3145 {
3146 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
3147 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3148 }
3149 else if (angle == 90)
3150 {
3151 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
3152 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
3153 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
3154 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
3155 }
3156 else if (angle == 180)
3157 {
3158 fits_update_key_dbl(fptr, "CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
3159 fits_update_key_dbl(fptr, "CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
3160 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3161 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
3162 }
3163 else if (angle == 270)
3164 {
3165 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
3166 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
3167 fits_update_key_dbl(fptr, "CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
3168 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
3169 }
3170 }
3171 else
3172 {
3173 if (angle == 90)
3174 {
3175 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
3176 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
3177 fits_update_key_dbl(fptr, "CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
3178 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
3179 }
3180 else if (angle == 180)
3181 {
3182 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
3183 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
3184 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3185 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
3186 }
3187 else if (angle == 270)
3188 {
3189 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
3190 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
3191 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
3192 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
3193 }
3194 }
3195 }
3196
3197 /* Delete any polynomial solution */
3198 /* (These could maybe be switched, but I don't want to work them out yet */
3199 status = 0;
3200 if (!fits_read_key_dbl(fptr, "CO1_1", &ctemp1, comment, &status))
3201 {
3202 int i;
3203 char keyword[16];
3204
3205 for (i = 1; i < 13; i++)
3206 {
3207 sprintf(keyword, "CO1_%d", i);
3208 fits_delete_key(fptr, keyword, &status);
3209 }
3210 for (i = 1; i < 13; i++)
3211 {
3212 sprintf(keyword, "CO2_%d", i);
3213 fits_delete_key(fptr, keyword, &status);
3214 }
3215 }
3216
3217 }
3218
getWritableImageBuffer()3219 uint8_t * FITSData::getWritableImageBuffer()
3220 {
3221 return m_ImageBuffer;
3222 }
3223
getImageBuffer() const3224 uint8_t const * FITSData::getImageBuffer() const
3225 {
3226 return m_ImageBuffer;
3227 }
3228
setImageBuffer(uint8_t * buffer)3229 void FITSData::setImageBuffer(uint8_t * buffer)
3230 {
3231 delete[] m_ImageBuffer;
3232 m_ImageBuffer = buffer;
3233 }
3234
checkDebayer()3235 bool FITSData::checkDebayer()
3236 {
3237 int status = 0;
3238 char bayerPattern[64], roworder[64];
3239
3240 // Let's search for BAYERPAT keyword, if it's not found we return as there is no bayer pattern in this image
3241 if (fits_read_keyword(fptr, "BAYERPAT", bayerPattern, nullptr, &status))
3242 return false;
3243
3244 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
3245 {
3246 KSNotification::error(i18n("Only 8 and 16 bits bayered images supported."), i18n("Debayer error"));
3247 return false;
3248 }
3249 QString pattern(bayerPattern);
3250 pattern = pattern.remove('\'').trimmed();
3251
3252 QString order(roworder);
3253 order = order.remove('\'').trimmed();
3254
3255 if (order == "BOTTOM-UP" && !(m_Statistics.height % 2))
3256 {
3257 if (pattern == "RGGB")
3258 pattern = "GBRG";
3259 else if (pattern == "GBRG")
3260 pattern = "RGGB";
3261 else if (pattern == "GRBG")
3262 pattern = "BGGR";
3263 else if (pattern == "BGGR")
3264 pattern = "GRBG";
3265 else return false;
3266 }
3267
3268 if (pattern == "RGGB")
3269 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
3270 else if (pattern == "GBRG")
3271 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
3272 else if (pattern == "GRBG")
3273 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
3274 else if (pattern == "BGGR")
3275 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
3276 // We return unless we find a valid pattern
3277 else
3278 {
3279 KSNotification::error(i18n("Unsupported bayer pattern %1.", pattern), i18n("Debayer error"));
3280 return false;
3281 }
3282
3283 fits_read_key(fptr, TINT, "XBAYROFF", &debayerParams.offsetX, nullptr, &status);
3284 fits_read_key(fptr, TINT, "YBAYROFF", &debayerParams.offsetY, nullptr, &status);
3285
3286 if (debayerParams.offsetX == 1)
3287 {
3288 // This may leave odd values in the 0th column if the color filter is not there
3289 // in the sensor, but otherwise should process the offset correctly.
3290 // Only offsets of 0 or 1 are implemented in debayer_8bit() and debayer_16bit().
3291 switch (debayerParams.filter)
3292 {
3293 case DC1394_COLOR_FILTER_RGGB:
3294 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
3295 break;
3296 case DC1394_COLOR_FILTER_GBRG:
3297 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
3298 break;
3299 case DC1394_COLOR_FILTER_GRBG:
3300 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
3301 break;
3302 case DC1394_COLOR_FILTER_BGGR:
3303 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
3304 break;
3305 }
3306 debayerParams.offsetX = 0;
3307 }
3308 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
3309 {
3310 KSNotification::error(i18n("Unsupported bayer offsets %1 %2.",
3311 debayerParams.offsetX, debayerParams.offsetY), i18n("Debayer error"));
3312 return false;
3313 }
3314
3315 HasDebayer = true;
3316
3317 return true;
3318 }
3319
getBayerParams(BayerParams * param)3320 void FITSData::getBayerParams(BayerParams * param)
3321 {
3322 param->method = debayerParams.method;
3323 param->filter = debayerParams.filter;
3324 param->offsetX = debayerParams.offsetX;
3325 param->offsetY = debayerParams.offsetY;
3326 }
3327
setBayerParams(BayerParams * param)3328 void FITSData::setBayerParams(BayerParams * param)
3329 {
3330 debayerParams.method = param->method;
3331 debayerParams.filter = param->filter;
3332 debayerParams.offsetX = param->offsetX;
3333 debayerParams.offsetY = param->offsetY;
3334 }
3335
debayer(bool reload)3336 bool FITSData::debayer(bool reload)
3337 {
3338 if (reload)
3339 {
3340 int anynull = 0, status = 0;
3341
3342 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel, nullptr, m_ImageBuffer,
3343 &anynull, &status))
3344 {
3345 // char errmsg[512];
3346 // fits_get_errstatus(status, errmsg);
3347 // KSNotification::error(i18n("Error reading image: %1", QString(errmsg)), i18n("Debayer error"));
3348 return false;
3349 }
3350 }
3351
3352 switch (m_Statistics.dataType)
3353 {
3354 case TBYTE:
3355 return debayer_8bit();
3356
3357 case TUSHORT:
3358 return debayer_16bit();
3359
3360 default:
3361 return false;
3362 }
3363 }
3364
debayer_8bit()3365 bool FITSData::debayer_8bit()
3366 {
3367 dc1394error_t error_code;
3368
3369 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
3370 uint8_t * destinationBuffer = nullptr;
3371
3372 try
3373 {
3374 destinationBuffer = new uint8_t[rgb_size];
3375 }
3376 catch (const std::bad_alloc &e)
3377 {
3378 logOOMError(rgb_size);
3379 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what()), i18n("Debayer error"), 10);
3380 return false;
3381 }
3382
3383 auto * bayer_source_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
3384 auto * bayer_destination_buffer = reinterpret_cast<uint8_t *>(destinationBuffer);
3385
3386 if (bayer_destination_buffer == nullptr)
3387 {
3388 logOOMError(rgb_size);
3389 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error"), 10);
3390 return false;
3391 }
3392
3393 int ds1394_height = m_Statistics.height;
3394 auto dc1394_source = bayer_source_buffer;
3395
3396 if (debayerParams.offsetY == 1)
3397 {
3398 dc1394_source += m_Statistics.width;
3399 ds1394_height--;
3400 }
3401 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
3402
3403 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
3404 debayerParams.filter,
3405 debayerParams.method);
3406
3407 if (error_code != DC1394_SUCCESS)
3408 {
3409 KSNotification::error(i18n("Debayer failed (%1)", error_code), i18n("Debayer error"), 10);
3410 m_Statistics.channels = 1;
3411 delete[] destinationBuffer;
3412 return false;
3413 }
3414
3415 if (m_ImageBufferSize != rgb_size)
3416 {
3417 delete[] m_ImageBuffer;
3418 try
3419 {
3420 m_ImageBuffer = new uint8_t[rgb_size];
3421 }
3422 catch (const std::bad_alloc &e)
3423 {
3424 delete[] destinationBuffer;
3425 logOOMError(rgb_size);
3426 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what()), i18n("Debayer error"), 10);
3427 return false;
3428 }
3429
3430 m_ImageBufferSize = rgb_size;
3431 }
3432
3433 auto bayered_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
3434
3435 // Data in R1G1B1, we need to copy them into 3 layers for FITS
3436
3437 uint8_t * rBuff = bayered_buffer;
3438 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
3439 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
3440
3441 int imax = m_Statistics.samples_per_channel * 3 - 3;
3442 for (int i = 0; i <= imax; i += 3)
3443 {
3444 *rBuff++ = bayer_destination_buffer[i];
3445 *gBuff++ = bayer_destination_buffer[i + 1];
3446 *bBuff++ = bayer_destination_buffer[i + 2];
3447 }
3448
3449 m_Statistics.channels = (m_Mode == FITS_NORMAL) ? 3 : 1;
3450 m_Statistics.dataType = TBYTE;
3451 delete[] destinationBuffer;
3452 return true;
3453 }
3454
debayer_16bit()3455 bool FITSData::debayer_16bit()
3456 {
3457 dc1394error_t error_code;
3458
3459 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
3460 uint8_t *destinationBuffer = nullptr;
3461 try
3462 {
3463 destinationBuffer = new uint8_t[rgb_size];
3464 }
3465 catch (const std::bad_alloc &e)
3466 {
3467 logOOMError(rgb_size);
3468 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what()), i18n("Debayer error"), 10);
3469 return false;
3470 }
3471
3472 auto * bayer_source_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
3473 auto * bayer_destination_buffer = reinterpret_cast<uint16_t *>(destinationBuffer);
3474
3475 if (bayer_destination_buffer == nullptr)
3476 {
3477 logOOMError(rgb_size);
3478 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer."), i18n("Debayer error"), 10);
3479 return false;
3480 }
3481
3482 int ds1394_height = m_Statistics.height;
3483 auto dc1394_source = bayer_source_buffer;
3484
3485 if (debayerParams.offsetY == 1)
3486 {
3487 dc1394_source += m_Statistics.width;
3488 ds1394_height--;
3489 }
3490 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
3491
3492 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
3493 debayerParams.filter,
3494 debayerParams.method, 16);
3495
3496 if (error_code != DC1394_SUCCESS)
3497 {
3498 KSNotification::error(i18n("Debayer failed (%1)", error_code), i18n("Debayer error"));
3499 m_Statistics.channels = 1;
3500 delete[] destinationBuffer;
3501 return false;
3502 }
3503
3504 if (m_ImageBufferSize != rgb_size)
3505 {
3506 delete[] m_ImageBuffer;
3507 try
3508 {
3509 m_ImageBuffer = new uint8_t[rgb_size];
3510 }
3511 catch (const std::bad_alloc &e)
3512 {
3513 logOOMError(rgb_size);
3514 delete[] destinationBuffer;
3515 KSNotification::error(i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what()), i18n("Debayer error"), 10);
3516 return false;
3517 }
3518
3519 m_ImageBufferSize = rgb_size;
3520 }
3521
3522 auto bayered_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
3523
3524 // Data in R1G1B1, we need to copy them into 3 layers for FITS
3525
3526 uint16_t * rBuff = bayered_buffer;
3527 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
3528 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
3529
3530 int imax = m_Statistics.samples_per_channel * 3 - 3;
3531 for (int i = 0; i <= imax; i += 3)
3532 {
3533 *rBuff++ = bayer_destination_buffer[i];
3534 *gBuff++ = bayer_destination_buffer[i + 1];
3535 *bBuff++ = bayer_destination_buffer[i + 2];
3536 }
3537
3538 m_Statistics.channels = (m_Mode == FITS_NORMAL) ? 3 : 1;
3539 m_Statistics.dataType = TUSHORT;
3540 delete[] destinationBuffer;
3541 return true;
3542 }
3543
logOOMError(uint32_t requiredMemory)3544 void FITSData::logOOMError(uint32_t requiredMemory)
3545 {
3546 qCCritical(KSTARS_FITS) << "Debayed memory allocation failure. Required Memory:" << KFormat().formatByteSize(requiredMemory)
3547 << "Available system memory:" << KSUtils::getAvailableRAM();
3548 }
3549
getADU() const3550 double FITSData::getADU() const
3551 {
3552 double adu = 0;
3553 for (int i = 0; i < m_Statistics.channels; i++)
3554 adu += m_Statistics.mean[i];
3555
3556 return (adu / static_cast<double>(m_Statistics.channels));
3557 }
3558
getLastError() const3559 QString FITSData::getLastError() const
3560 {
3561 return m_LastError;
3562 }
3563
3564 template <typename T>
convertToQImage(double dataMin,double dataMax,double scale,double zero,QImage & image)3565 void FITSData::convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image)
3566 {
3567 const auto * buffer = reinterpret_cast<const T *>(getImageBuffer());
3568 const T limit = std::numeric_limits<T>::max();
3569 T bMin = dataMin < 0 ? 0 : dataMin;
3570 T bMax = dataMax > limit ? limit : dataMax;
3571 uint16_t w = width();
3572 uint16_t h = height();
3573 uint32_t size = w * h;
3574 double val;
3575
3576 if (channels() == 1)
3577 {
3578 /* Fill in pixel values using indexed map, linear scale */
3579 for (int j = 0; j < h; j++)
3580 {
3581 unsigned char * scanLine = image.scanLine(j);
3582
3583 for (int i = 0; i < w; i++)
3584 {
3585 val = qBound(bMin, buffer[j * w + i], bMax);
3586 val = val * scale + zero;
3587 scanLine[i] = qBound<unsigned char>(0, static_cast<uint8_t>(val), 255);
3588 }
3589 }
3590 }
3591 else
3592 {
3593 double rval = 0, gval = 0, bval = 0;
3594 QRgb value;
3595 /* Fill in pixel values using indexed map, linear scale */
3596 for (int j = 0; j < h; j++)
3597 {
3598 auto * scanLine = reinterpret_cast<QRgb *>((image.scanLine(j)));
3599
3600 for (int i = 0; i < w; i++)
3601 {
3602 rval = qBound(bMin, buffer[j * w + i], bMax);
3603 gval = qBound(bMin, buffer[j * w + i + size], bMax);
3604 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
3605
3606 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
3607
3608 scanLine[i] = value;
3609 }
3610 }
3611 }
3612 }
3613
FITSToImage(const QString & filename)3614 QImage FITSData::FITSToImage(const QString &filename)
3615 {
3616 QImage fitsImage;
3617 double min, max;
3618
3619 FITSData data;
3620
3621 QFuture<bool> future = data.loadFromFile(filename);
3622
3623 // Wait synchronously
3624 future.waitForFinished();
3625 if (future.result() == false)
3626 return fitsImage;
3627
3628 data.getMinMax(&min, &max);
3629
3630 if (min == max)
3631 {
3632 fitsImage.fill(Qt::white);
3633 return fitsImage;
3634 }
3635
3636 if (data.channels() == 1)
3637 {
3638 fitsImage = QImage(data.width(), data.height(), QImage::Format_Indexed8);
3639
3640 fitsImage.setColorCount(256);
3641 for (int i = 0; i < 256; i++)
3642 fitsImage.setColor(i, qRgb(i, i, i));
3643 }
3644 else
3645 {
3646 fitsImage = QImage(data.width(), data.height(), QImage::Format_RGB32);
3647 }
3648
3649 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
3650 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
3651
3652 double bscale = 255. / (dataMax - dataMin);
3653 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
3654
3655 // Long way to do this since we do not want to use templated functions here
3656 switch (data.m_Statistics.dataType)
3657 {
3658 case TBYTE:
3659 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3660 break;
3661
3662 case TSHORT:
3663 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3664 break;
3665
3666 case TUSHORT:
3667 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3668 break;
3669
3670 case TLONG:
3671 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3672 break;
3673
3674 case TULONG:
3675 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3676 break;
3677
3678 case TFLOAT:
3679 data.convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
3680 break;
3681
3682 case TLONGLONG:
3683 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3684 break;
3685
3686 case TDOUBLE:
3687 data.convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
3688 break;
3689
3690 default:
3691 break;
3692 }
3693
3694 return fitsImage;
3695 }
3696
ImageToFITS(const QString & filename,const QString & format,QString & output)3697 bool FITSData::ImageToFITS(const QString &filename, const QString &format, QString &output)
3698 {
3699 if (QImageReader::supportedImageFormats().contains(format.toLatin1()) == false)
3700 {
3701 qCCritical(KSTARS_FITS) << "Failed to convert" << filename << "to FITS since format" << format << "is not supported in Qt";
3702 return false;
3703 }
3704
3705 QImage input;
3706
3707 if (input.load(filename, format.toLatin1()) == false)
3708 {
3709 qCCritical(KSTARS_FITS) << "Failed to open image" << filename;
3710 return false;
3711 }
3712
3713 output = QDir(getTemporaryPath()).filePath(filename + ".fits");
3714
3715 //This section sets up the FITS File
3716 fitsfile *fptr = nullptr;
3717 int status = 0;
3718 long fpixel = 1, naxis = input.allGray() ? 2 : 3, nelements, exposure;
3719 long naxes[3] = { input.width(), input.height(), naxis == 3 ? 3 : 1 };
3720 char error_status[512] = {0};
3721
3722 if (fits_create_file(&fptr, QString('!' + output).toLocal8Bit().data(), &status))
3723 {
3724 fits_get_errstatus(status, error_status);
3725 qCCritical(KSTARS_FITS) << "Failed to create FITS file. Error:" << error_status;
3726 return false;
3727 }
3728
3729 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
3730 {
3731 qCWarning(KSTARS_FITS) << "fits_create_img failed:" << error_status;
3732 status = 0;
3733 fits_flush_file(fptr, &status);
3734 fits_close_file(fptr, &status);
3735 return false;
3736 }
3737
3738 exposure = 1;
3739 fits_update_key(fptr, TLONG, "EXPOSURE", &exposure, "Total Exposure Time", &status);
3740
3741 // Gray image
3742 if (naxis == 2)
3743 {
3744 nelements = naxes[0] * naxes[1];
3745 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.bits(), &status))
3746 {
3747 fits_get_errstatus(status, error_status);
3748 qCWarning(KSTARS_FITS) << "fits_write_img GRAY failed:" << error_status;
3749 status = 0;
3750 fits_flush_file(fptr, &status);
3751 fits_close_file(fptr, &status);
3752 return false;
3753 }
3754 }
3755 // RGB image, we have to convert from ARGB format to R G B for each plane
3756 else
3757 {
3758 nelements = naxes[0] * naxes[1] * 3;
3759
3760 uint8_t *srcBuffer = input.bits();
3761 // ARGB
3762 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
3763
3764 uint8_t *rgbBuffer = new uint8_t[nelements];
3765 if (rgbBuffer == nullptr)
3766 {
3767 qCWarning(KSTARS_FITS) << "Not enough memory for RGB buffer";
3768 fits_flush_file(fptr, &status);
3769 fits_close_file(fptr, &status);
3770 return false;
3771 }
3772
3773 uint8_t *subR = rgbBuffer;
3774 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
3775 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
3776 for (uint32_t i = 0; i < srcBytes; i += 4)
3777 {
3778 *subB++ = srcBuffer[i];
3779 *subG++ = srcBuffer[i + 1];
3780 *subR++ = srcBuffer[i + 2];
3781 }
3782
3783 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
3784 {
3785 fits_get_errstatus(status, error_status);
3786 qCWarning(KSTARS_FITS) << "fits_write_img RGB failed:" << error_status;
3787 status = 0;
3788 fits_flush_file(fptr, &status);
3789 fits_close_file(fptr, &status);
3790 delete [] rgbBuffer;
3791 return false;
3792 }
3793
3794 delete [] rgbBuffer;
3795 }
3796
3797 if (fits_flush_file(fptr, &status))
3798 {
3799 fits_get_errstatus(status, error_status);
3800 qCWarning(KSTARS_FITS) << "fits_flush_file failed:" << error_status;
3801 status = 0;
3802 fits_close_file(fptr, &status);
3803 return false;
3804 }
3805
3806 if (fits_close_file(fptr, &status))
3807 {
3808 fits_get_errstatus(status, error_status);
3809 qCWarning(KSTARS_FITS) << "fits_close_file failed:" << error_status;
3810 return false;
3811 }
3812
3813 return true;
3814 }
3815
injectWCS(double orientation,double ra,double dec,double pixscale,bool eastToTheRight)3816 bool FITSData::injectWCS(double orientation, double ra, double dec, double pixscale, bool eastToTheRight)
3817 {
3818 int status = 0;
3819
3820 fits_update_key(fptr, TDOUBLE, "OBJCTRA", &ra, "Object RA", &status);
3821 fits_update_key(fptr, TDOUBLE, "OBJCTDEC", &dec, "Object DEC", &status);
3822
3823 int epoch = 2000;
3824
3825 fits_update_key(fptr, TINT, "EQUINOX", &epoch, "Equinox", &status);
3826
3827 fits_update_key(fptr, TDOUBLE, "CRVAL1", &ra, "CRVAL1", &status);
3828 fits_update_key(fptr, TDOUBLE, "CRVAL2", &dec, "CRVAL1", &status);
3829
3830 char radecsys[8] = "FK5";
3831 char ctype1[16] = "RA---TAN";
3832 char ctype2[16] = "DEC--TAN";
3833
3834 fits_update_key(fptr, TSTRING, "RADECSYS", radecsys, "RADECSYS", &status);
3835 fits_update_key(fptr, TSTRING, "CTYPE1", ctype1, "CTYPE1", &status);
3836 fits_update_key(fptr, TSTRING, "CTYPE2", ctype2, "CTYPE2", &status);
3837
3838 double crpix1 = width() / 2.0;
3839 double crpix2 = height() / 2.0;
3840
3841 fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix1, "CRPIX1", &status);
3842 fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix2, "CRPIX2", &status);
3843
3844 // Arcsecs per Pixel
3845 double secpix1 = eastToTheRight ? pixscale : -pixscale;
3846 double secpix2 = pixscale;
3847
3848 fits_update_key(fptr, TDOUBLE, "SECPIX1", &secpix1, "SECPIX1", &status);
3849 fits_update_key(fptr, TDOUBLE, "SECPIX2", &secpix2, "SECPIX2", &status);
3850
3851 double degpix1 = secpix1 / 3600.0;
3852 double degpix2 = secpix2 / 3600.0;
3853
3854 fits_update_key(fptr, TDOUBLE, "CDELT1", °pix1, "CDELT1", &status);
3855 fits_update_key(fptr, TDOUBLE, "CDELT2", °pix2, "CDELT2", &status);
3856
3857 // Rotation is CW, we need to convert it to CCW per CROTA1 definition
3858 double rotation = 360 - orientation;
3859 if (rotation > 360)
3860 rotation -= 360;
3861
3862 fits_update_key(fptr, TDOUBLE, "CROTA1", &rotation, "CROTA1", &status);
3863 fits_update_key(fptr, TDOUBLE, "CROTA2", &rotation, "CROTA2", &status);
3864
3865 m_WCSState = Idle;
3866
3867 qCDebug(KSTARS_FITS) << "Finished update WCS info.";
3868
3869 return true;
3870 }
3871
contains(const QPointF & point) const3872 bool FITSData::contains(const QPointF &point) const
3873 {
3874 return (point.x() >= 0 && point.y() >= 0 && point.x() <= m_Statistics.width && point.y() <= m_Statistics.height);
3875 }
3876
saveStatistics(FITSImage::Statistic & other)3877 void FITSData::saveStatistics(FITSImage::Statistic &other)
3878 {
3879 other = m_Statistics;
3880 }
3881
restoreStatistics(FITSImage::Statistic & other)3882 void FITSData::restoreStatistics(FITSImage::Statistic &other)
3883 {
3884 m_Statistics = other;
3885
3886 emit dataChanged();
3887 }
3888
constructHistogram()3889 void FITSData::constructHistogram()
3890 {
3891 switch (m_Statistics.dataType)
3892 {
3893 case TBYTE:
3894 constructHistogramInternal<uint8_t>();
3895 break;
3896
3897 case TSHORT:
3898 constructHistogramInternal<int16_t>();
3899 break;
3900
3901 case TUSHORT:
3902 constructHistogramInternal<uint16_t>();
3903 break;
3904
3905 case TLONG:
3906 constructHistogramInternal<int32_t>();
3907 break;
3908
3909 case TULONG:
3910 constructHistogramInternal<uint32_t>();
3911 break;
3912
3913 case TFLOAT:
3914 constructHistogramInternal<float>();
3915 break;
3916
3917 case TLONGLONG:
3918 constructHistogramInternal<int64_t>();
3919 break;
3920
3921 case TDOUBLE:
3922 constructHistogramInternal<double>();
3923 break;
3924
3925 default:
3926 break;
3927 }
3928 }
3929
constructHistogramInternal()3930 template <typename T> void FITSData::constructHistogramInternal()
3931 {
3932 auto * const buffer = reinterpret_cast<T const *>(m_ImageBuffer);
3933 uint32_t samples = m_Statistics.width * m_Statistics.height;
3934 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
3935
3936 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
3937 if (m_HistogramBinCount <= 0)
3938 m_HistogramBinCount = 256;
3939
3940 for (int n = 0; n < m_Statistics.channels; n++)
3941 {
3942 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
3943 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
3944 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
3945 m_HistogramBinWidth[n] = qMax(1.0, (m_Statistics.max[n] - m_Statistics.min[n]) / (m_HistogramBinCount - 1));
3946 }
3947
3948 QVector<QFuture<void>> futures;
3949
3950 for (int n = 0; n < m_Statistics.channels; n++)
3951 {
3952 futures.append(QtConcurrent::run([ = ]()
3953 {
3954 for (int i = 0; i < m_HistogramBinCount; i++)
3955 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
3956 }));
3957 }
3958
3959 for (int n = 0; n < m_Statistics.channels; n++)
3960 {
3961 futures.append(QtConcurrent::run([ = ]()
3962 {
3963 uint32_t offset = n * samples;
3964
3965 for (uint32_t i = 0; i < samples; i += sampleBy)
3966 {
3967 int32_t id = qMax(static_cast<T>(0), qMin(static_cast<T>(m_HistogramBinCount),
3968 static_cast<T>(rint((buffer[i + offset] - m_Statistics.min[n]) / m_HistogramBinWidth[n]))));
3969 m_HistogramFrequency[n][id] += sampleBy;
3970 }
3971 }));
3972 }
3973
3974 for (QFuture<void> future : futures)
3975 future.waitForFinished();
3976
3977 futures.clear();
3978
3979 for (int n = 0; n < m_Statistics.channels; n++)
3980 {
3981 futures.append(QtConcurrent::run([ = ]()
3982 {
3983 uint32_t accumulator = 0;
3984 for (int i = 0; i < m_HistogramBinCount; i++)
3985 {
3986 accumulator += m_HistogramFrequency[n][i];
3987 m_CumulativeFrequency[n].replace(i, accumulator);
3988 }
3989 }));
3990 }
3991
3992 for (QFuture<void> future : futures)
3993 future.waitForFinished();
3994
3995 futures.clear();
3996
3997 #if 0
3998 for (int n = 0; n < m_Statistics.channels; n++)
3999 {
4000 futures.append(QtConcurrent::run([ = ]()
4001 {
4002 double median[3] = {0};
4003 //const bool cutoffSpikes = ui->hideSaturated->isChecked();
4004 const uint32_t halfCumulative = m_CumulativeFrequency[n][m_HistogramBinCount - 1] / 2;
4005
4006 // Find which bin contains the median.
4007 int median_bin = -1;
4008 for (int i = 0; i < m_HistogramBinCount; i++)
4009 {
4010 if (m_CumulativeFrequency[n][i] > halfCumulative)
4011 {
4012 median_bin = i;
4013 break;
4014 }
4015 }
4016
4017 if (median_bin >= 0)
4018 {
4019 // The number of items in the median bin
4020 const uint32_t median_bin_size = m_HistogramFrequency[n][median_bin] / sampleBy;
4021 if (median_bin_size > 0)
4022 {
4023 // The median is this element inside the sorted median_bin;
4024 const uint32_t samples_before_median_bin = median_bin == 0 ? 0 : m_CumulativeFrequency[n][median_bin - 1];
4025 uint32_t median_position = (halfCumulative - samples_before_median_bin) / sampleBy;
4026
4027 if (median_position >= median_bin_size)
4028 median_position = median_bin_size - 1;
4029 if (median_position >= 0 && median_position < median_bin_size)
4030 {
4031 // Fill a vector with the values in the median bin (sampled by sampleBy).
4032 std::vector<T> median_bin_samples(median_bin_size);
4033 uint32_t upto = 0;
4034 const uint32_t offset = n * samples;
4035 for (uint32_t i = 0; i < samples; i += sampleBy)
4036 {
4037 if (upto >= median_bin_size) break;
4038 const int32_t id = rint((buffer[i + offset] - m_Statistics.min[n]) / m_HistogramBinWidth[n]);
4039 if (id == median_bin)
4040 median_bin_samples[upto++] = buffer[i + offset];
4041 }
4042 // Get the Nth value using N = the median position.
4043 if (upto > 0)
4044 {
4045 if (median_position >= upto) median_position = upto - 1;
4046 std::nth_element(median_bin_samples.begin(), median_bin_samples.begin() + median_position,
4047 median_bin_samples.begin() + upto);
4048 median[n] = median_bin_samples[median_position];
4049 }
4050 }
4051 }
4052 }
4053
4054 setMedian(median[n], n);
4055
4056 }));
4057 }
4058
4059 for (QFuture<void> future : futures)
4060 future.waitForFinished();
4061 #endif
4062
4063 // Custom index to indicate the overall contrast of the image
4064 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
4065 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] / static_cast<double>
4066 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
4067 4]);
4068 else
4069 m_JMIndex = 1;
4070
4071 qCDebug(KSTARS_FITS) << "FITHistogram: JMIndex " << m_JMIndex;
4072
4073 m_HistogramConstructed = true;
4074 emit histogramReady();
4075 }
4076
recordLastError(int errorCode)4077 void FITSData::recordLastError(int errorCode)
4078 {
4079 char fitsErrorMessage[512] = {0};
4080 fits_get_errstatus(errorCode, fitsErrorMessage);
4081 m_LastError = fitsErrorMessage;
4082 }
4083
getAverageMean() const4084 double FITSData::getAverageMean() const
4085 {
4086 if (m_Statistics.channels == 1)
4087 return m_Statistics.mean[0];
4088 else
4089 return (m_Statistics.mean[0] + m_Statistics.mean[1] + m_Statistics.mean[2]) / 3.0;
4090 }
4091
getAverageMedian() const4092 double FITSData::getAverageMedian() const
4093 {
4094 if (m_Statistics.channels == 1)
4095 return m_Statistics.median[0];
4096 else
4097 return (m_Statistics.median[0] + m_Statistics.median[1] + m_Statistics.median[2]) / 3.0;
4098 }
4099
getAverageStdDev() const4100 double FITSData::getAverageStdDev() const
4101 {
4102 if (m_Statistics.channels == 1)
4103 return m_Statistics.stddev[0];
4104 else
4105 return (m_Statistics.stddev[0] + m_Statistics.stddev[1] + m_Statistics.stddev[2]) / 3.0;
4106 }
4107