Kstars

fitsdata.cpp
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#include "auxiliary/robuststatistics.h"
26
27#include <KFormat>
28#include <QApplication>
29#include <QImage>
30#include <QtConcurrent>
31#include <QImageReader>
32#include <QUrl>
33#include <QNetworkAccessManager>
34
35#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
36#include <wcshdr.h>
37#include <wcsfix.h>
38#endif
39
40#if !defined(KSTARS_LITE) && defined(HAVE_LIBRAW)
41#include <libraw/libraw.h>
42#endif
43
44#ifdef HAVE_XISF
45#include <libxisf.h>
46#endif
47
48#include <cfloat>
49#include <cmath>
50
51#include <fits_debug.h>
52
53#define ZOOM_DEFAULT 100.0
54#define ZOOM_MIN 10
55#define ZOOM_MAX 400
56#define ZOOM_LOW_INCR 10
57#define ZOOM_HIGH_INCR 50
58
59QString getTemporaryPath()
60{
61 return QDir(KSPaths::writableLocation(QStandardPaths::TempLocation) + "/" +
62 qAppName()).path();
63}
64
65const QStringList RAWFormats = { "cr2", "cr3", "crw", "nef", "raf", "dng", "arw", "orf" };
66
67bool FITSData::readableFilename(const QString &filename)
68{
69 QFileInfo info(filename);
70 QString extension = info.completeSuffix().toLower();
71 return extension.contains("fit") || extension.contains("fz") ||
72 extension.contains("xisf") ||
74 RAWFormats.contains(extension);
75}
76
77FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)
78{
79 static const QRegularExpression re("[-{}]");
80
81 qRegisterMetaType<FITSMode>("FITSMode");
82
83 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
84 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
85 debayerParams.offsetX = debayerParams.offsetY = 0;
86
87 // Reserve 3 channels
88 m_CumulativeFrequency.resize(3);
89 m_HistogramBinWidth.resize(3);
90 m_HistogramFrequency.resize(3);
91 m_HistogramIntensity.resize(3);
92
93 // Set UUID for each view
95 uuid = uuid.remove(re);
96 setObjectName(uuid);
97}
98
99FITSData::FITSData(const QSharedPointer<FITSData> &other)
100{
101 static const QRegularExpression re("[-{}]");
102 qRegisterMetaType<FITSMode>("FITSMode");
103
104 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
105 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
106 debayerParams.offsetX = debayerParams.offsetY = 0;
107
108 m_TemporaryDataFile.setFileTemplate("fits_memory_XXXXXX");
109
110 this->m_Mode = other->m_Mode;
111 this->m_Statistics.channels = other->m_Statistics.channels;
112 memcpy(&m_Statistics, &(other->m_Statistics), sizeof(m_Statistics));
113 m_ImageBuffer = new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel];
114 memcpy(m_ImageBuffer, other->m_ImageBuffer,
115 m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel);
116
117 // Set UUID for each view
119 uuid = uuid.remove(re);
120 setObjectName(uuid);
121}
122
123FITSData::~FITSData()
124{
125 int status = 0;
126
127 if (m_StarFindFuture.isRunning())
128 m_StarFindFuture.waitForFinished();
129
130 clearImageBuffers();
131
132#ifdef HAVE_WCSLIB
133 if (m_WCSHandle != nullptr)
134 {
135 wcsvfree(&m_nwcs, &m_WCSHandle);
136 m_WCSHandle = nullptr;
137 m_nwcs = 0;
138 }
139#endif
140
141 if (starCenters.count() > 0)
142 qDeleteAll(starCenters);
143 starCenters.clear();
144
145 if (m_SkyObjects.count() > 0)
146 qDeleteAll(m_SkyObjects);
147 m_SkyObjects.clear();
148
149 m_CatObjects.clear();
150
151 if (fptr != nullptr)
152 {
153 fits_flush_file(fptr, &status);
154 fits_close_file(fptr, &status);
155 free(m_PackBuffer);
156 m_PackBuffer = nullptr;
157 fptr = nullptr;
158 }
159}
160
161void FITSData::loadCommon(const QString &inFilename)
162{
163 int status = 0;
164 qDeleteAll(starCenters);
165 starCenters.clear();
166
167 if (fptr != nullptr)
168 {
169 fits_flush_file(fptr, &status);
170 fits_close_file(fptr, &status);
171 free(m_PackBuffer);
172 m_PackBuffer = nullptr;
173 fptr = nullptr;
174 }
175
176 m_Filename = inFilename;
177}
178
179bool FITSData::loadFromBuffer(const QByteArray &buffer)
180{
181 loadCommon("");
182 qCDebug(KSTARS_FITS) << "Reading file buffer (" << KFormat().formatByteSize(buffer.size()) << ")";
183 return privateLoad(buffer);
184}
185
186QFuture<bool> FITSData::loadFromFile(const QString &inFilename)
187{
188 loadCommon(inFilename);
189 QFileInfo info(m_Filename);
190 m_Extension = info.completeSuffix().toLower();
191 qCDebug(KSTARS_FITS) << "Loading file " << m_Filename;
192#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
193 return QtConcurrent::run(&FITSData::privateLoad, this, QByteArray());
194#else
195 return QtConcurrent::run(this, &FITSData::privateLoad, QByteArray());
196#endif
197}
198
199namespace
200{
201// Common code for reporting fits read errors. Always returns false.
202QString fitsErrorToString(int status)
203{
204 char error_status[512] = {0};
205 fits_report_error(stderr, status);
206 fits_get_errstatus(status, error_status);
207 QString message = error_status;
208 return message;
209}
210}
211
212bool FITSData::privateLoad(const QByteArray &buffer)
213{
214 m_isTemporary = m_Filename.startsWith(KSPaths::writableLocation(QStandardPaths::TempLocation));
215 cacheHFR = -1;
216 cacheEccentricity = -1;
217
218 if (m_Extension.contains("fit") || m_Extension.contains("fz"))
219 return loadFITSImage(buffer);
220 if (m_Extension.contains("xisf"))
221 return loadXISFImage(buffer);
222 if (QImageReader::supportedImageFormats().contains(m_Extension.toLatin1()))
223 return loadCanonicalImage(buffer);
224 else if (RAWFormats.contains(m_Extension))
225 return loadRAWImage(buffer);
226
227 return false;
228}
229
230bool FITSData::loadFITSImage(const QByteArray &buffer, const bool isCompressed)
231{
232 int status = 0, anynull = 0;
233 long naxes[3];
234
235 m_HistogramConstructed = false;
236
237 if (m_Extension.contains(".fz") || isCompressed)
238 {
239 fpstate fpvar;
240 fp_init (&fpvar);
241 bool rc = false;
242
243 if (buffer.isEmpty())
244 {
245 // Store so we don't lose.
246 m_compressedFilename = m_Filename;
247
248 QString uncompressedFile = QDir::tempPath() + QString("/%1").arg(QUuid::createUuid().toString().remove(
249 QRegularExpression("[-{}]")));
250
251 rc = fp_unpack_file_to_fits(m_Filename.toLocal8Bit().data(), &fptr, fpvar) == 0;
252 if (rc)
253 {
254 m_Filename = uncompressedFile;
255 }
256 }
257 else
258 {
259 size_t m_PackBufferSize = 100000;
260 free(m_PackBuffer);
261 m_PackBuffer = (uint8_t *)malloc(m_PackBufferSize);
262 rc = fp_unpack_data_to_data(buffer.data(), buffer.size(), &m_PackBuffer, &m_PackBufferSize, fpvar) == 0;
263
264 if (rc)
265 {
266 void *data = reinterpret_cast<void *>(m_PackBuffer);
267 if (fits_open_memfile(&fptr, m_Filename.toLocal8Bit().data(), READONLY, &data, &m_PackBufferSize, 0,
268 nullptr, &status))
269 {
270 m_LastError = i18n("Error reading fits buffer: %1.", fitsErrorToString(status));
271 return false;
272 }
273
274 m_Statistics.size = m_PackBufferSize;
275 }
276 //rc = fp_unpack_data_to_fits(buffer.data(), buffer.size(), &fptr, fpvar) == 0;
277 }
278
279 if (rc == false)
280 {
281 free(m_PackBuffer);
282 m_PackBuffer = nullptr;
283 m_LastError = i18n("Failed to unpack compressed fits");
284 qCCritical(KSTARS_FITS) << m_LastError;
285 return false;
286 }
287
288 m_isTemporary = true;
289 m_isCompressed = true;
290 m_Statistics.size = fptr->Fptr->logfilesize;
291
292 }
293 else if (buffer.isEmpty())
294 {
295 // Use open diskfile as it does not use extended file names which has problems opening
296 // files with [ ] or ( ) in their names.
297 if (fits_open_diskfile(&fptr, m_Filename.toLocal8Bit(), READONLY, &status))
298 {
299 m_LastError = i18n("Error opening fits file %1 : %2", m_Filename, fitsErrorToString(status));
300 }
301
302 m_Statistics.size = QFile(m_Filename).size();
303 }
304 else
305 {
306 // Read the FITS file from a memory buffer.
307 void *temp_buffer = const_cast<void *>(reinterpret_cast<const void *>(buffer.data()));
308 size_t temp_size = buffer.size();
309 if (fits_open_memfile(&fptr, m_Filename.toLocal8Bit().data(), READONLY,
310 &temp_buffer, &temp_size, 0, nullptr, &status))
311 {
312 m_LastError = i18n("Error reading fits buffer: %1", fitsErrorToString(status));
313 return false;
314 }
315
316 m_Statistics.size = temp_size;
317 }
318
319 if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status))
320 {
321
322 free(m_PackBuffer);
323 m_PackBuffer = nullptr;
324 m_LastError = i18n("Could not locate image HDU: %1", fitsErrorToString(status));
325 }
326
327 if (fits_get_img_param(fptr, 3, &m_FITSBITPIX, &(m_Statistics.ndim), naxes, &status))
328 {
329 free(m_PackBuffer);
330 m_PackBuffer = nullptr;
331 m_LastError = i18n("FITS file open error (fits_get_img_param): %1", fitsErrorToString(status));
332 return false;
333 }
334
335 // Reload if it is transparently compressed.
336 if ((fits_is_compressed_image(fptr, &status) || m_Statistics.ndim <= 0) && !isCompressed)
337 {
338 loadCommon(m_Filename);
339 qCDebug(KSTARS_FITS) << "Image is compressed. Reloading...";
340 return loadFITSImage(buffer, true);
341 }
342
343 if (m_Statistics.ndim < 2)
344 {
345 m_LastError = i18n("1D FITS images are not supported in KStars.");
346 qCCritical(KSTARS_FITS) << m_LastError;
347 free(m_PackBuffer);
348 m_PackBuffer = nullptr;
349 return false;
350 }
351
352 switch (m_FITSBITPIX)
353 {
354 case BYTE_IMG:
355 m_Statistics.dataType = TBYTE;
356 m_Statistics.bytesPerPixel = sizeof(uint8_t);
357 break;
358 case SHORT_IMG:
359 // Read SHORT image as USHORT
360 m_FITSBITPIX = USHORT_IMG;
361 m_Statistics.dataType = TUSHORT;
362 m_Statistics.bytesPerPixel = sizeof(int16_t);
363 break;
364 case USHORT_IMG:
365 m_Statistics.dataType = TUSHORT;
366 m_Statistics.bytesPerPixel = sizeof(uint16_t);
367 break;
368 case LONG_IMG:
369 // Read LONG image as ULONG
370 m_FITSBITPIX = ULONG_IMG;
371 m_Statistics.dataType = TULONG;
372 m_Statistics.bytesPerPixel = sizeof(int32_t);
373 break;
374 case ULONG_IMG:
375 m_Statistics.dataType = TULONG;
376 m_Statistics.bytesPerPixel = sizeof(uint32_t);
377 break;
378 case FLOAT_IMG:
379 m_Statistics.dataType = TFLOAT;
380 m_Statistics.bytesPerPixel = sizeof(float);
381 break;
382 case LONGLONG_IMG:
383 m_Statistics.dataType = TLONGLONG;
384 m_Statistics.bytesPerPixel = sizeof(int64_t);
385 break;
386 case DOUBLE_IMG:
387 m_Statistics.dataType = TDOUBLE;
388 m_Statistics.bytesPerPixel = sizeof(double);
389 break;
390 default:
391 m_LastError = i18n("Bit depth %1 is not supported.", m_FITSBITPIX);
392 qCCritical(KSTARS_FITS) << m_LastError;
393 return false;
394 }
395
396 if (m_Statistics.ndim < 3)
397 naxes[2] = 1;
398
399 if (naxes[0] == 0 || naxes[1] == 0)
400 {
401 m_LastError = i18n("Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
402 qCCritical(KSTARS_FITS) << m_LastError;
403 free(m_PackBuffer);
404 m_PackBuffer = nullptr;
405 return false;
406 }
407
408 m_Statistics.width = naxes[0];
409 m_Statistics.height = naxes[1];
410 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
411 roiCenter.setX(m_Statistics.width / 2);
412 roiCenter.setY(m_Statistics.height / 2);
413 if(m_Statistics.width % 2)
414 roiCenter.setX(roiCenter.x() + 1);
415 if(m_Statistics.height % 2)
416 roiCenter.setY(roiCenter.y() + 1);
417
418 clearImageBuffers();
419
420 m_Statistics.channels = naxes[2];
421
422 // Channels always set to #1 if we are not required to process 3D Cubes
423 // Or if mode is not FITS_NORMAL (guide, focus..etc)
424 if ( (m_Mode != FITS_NORMAL && m_Mode != FITS_CALIBRATE) || !Options::auto3DCube())
425 m_Statistics.channels = 1;
426
427 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
428 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
429 if (m_ImageBuffer == nullptr)
430 {
431 qCWarning(KSTARS_FITS) << "FITSData: Not enough memory for image_buffer channel. Requested: "
432 << m_ImageBufferSize << " bytes.";
433 clearImageBuffers();
434 free(m_PackBuffer);
435 m_PackBuffer = nullptr;
436 return false;
437 }
438
439 rotCounter = 0;
440 flipHCounter = 0;
441 flipVCounter = 0;
442 long nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
443
444 if (fits_read_img(fptr, m_Statistics.dataType, 1, nelements, nullptr, m_ImageBuffer, &anynull, &status))
445 {
446 m_LastError = i18n("Error reading image: %1", fitsErrorToString(status));
447 return false;
448 }
449
450 parseHeader();
451
452 // Get UTC date time
453 QVariant value;
454 if (getRecordValue("DATE-OBS", value) && value.isValid())
455 {
456 QDateTime ts = value.toDateTime();
457 m_DateTime = KStarsDateTime(ts.date(), ts.time());
458 }
459
460 // Only check for debayed IF the original naxes[2] is 1
461 // which is for single channels.
462 if (naxes[2] == 1 && m_Statistics.channels == 1 && Options::autoDebayer() && checkDebayer())
463 {
464 // Save bayer image on disk in case we need to save it later since debayer destorys this data
465 if (m_isTemporary && m_TemporaryDataFile.open())
466 {
467 m_TemporaryDataFile.write(buffer);
468 m_TemporaryDataFile.close();
469 m_Filename = m_TemporaryDataFile.fileName();
470 }
471
472 if (debayer())
473 calculateStats(false, false);
474 }
475 else
476 calculateStats(false, false);
477
478 if (m_Mode == FITS_NORMAL || m_Mode == FITS_ALIGN)
479 loadWCS();
480
481 starsSearched = false;
482
483 return true;
484}
485
486bool FITSData::loadXISFImage(const QByteArray &buffer)
487{
488 m_HistogramConstructed = false;
489 clearImageBuffers();
490
491#ifdef HAVE_XISF
492 try
493 {
494 LibXISF::XISFReader xisfReader;
495 if (buffer.isEmpty())
496 {
497 xisfReader.open(m_Filename.toLocal8Bit().data());
498 }
499 else
500 {
501 LibXISF::ByteArray byteArray(buffer.constData(), buffer.size());
502 xisfReader.open(byteArray);
503 }
504
505 if (xisfReader.imagesCount() == 0)
506 {
507 m_LastError = i18n("File contain no images");
508 return false;
509 }
510
511 const LibXISF::Image &image = xisfReader.getImage(0);
512
513 switch (image.sampleFormat())
514 {
515 case LibXISF::Image::UInt8:
516 m_Statistics.dataType = TBYTE;
517 m_Statistics.bytesPerPixel = sizeof(LibXISF::UInt8);
518 m_FITSBITPIX = TBYTE;
519 break;
520 case LibXISF::Image::UInt16:
521 m_Statistics.dataType = TUSHORT;
522 m_Statistics.bytesPerPixel = sizeof(LibXISF::UInt16);
523 m_FITSBITPIX = TUSHORT;
524 break;
525 case LibXISF::Image::UInt32:
526 m_Statistics.dataType = TULONG;
527 m_Statistics.bytesPerPixel = sizeof(LibXISF::UInt32);
528 m_FITSBITPIX = TULONG;
529 break;
530 case LibXISF::Image::Float32:
531 m_Statistics.dataType = TFLOAT;
532 m_Statistics.bytesPerPixel = sizeof(LibXISF::Float32);
533 m_FITSBITPIX = TFLOAT;
534 break;
535 default:
536 m_LastError = i18n("Sample format %1 is not supported.", LibXISF::Image::sampleFormatString(image.sampleFormat()).c_str());
537 qCCritical(KSTARS_FITS) << m_LastError;
538 return false;
539 }
540
541 m_Statistics.width = image.width();
542 m_Statistics.height = image.height();
543 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
544 m_Statistics.channels = image.channelCount();
545 m_Statistics.size = buffer.size();
546 roiCenter.setX(m_Statistics.width / 2);
547 roiCenter.setY(m_Statistics.height / 2);
548 if(m_Statistics.width % 2)
549 roiCenter.setX(roiCenter.x() + 1);
550 if(m_Statistics.height % 2)
551 roiCenter.setY(roiCenter.y() + 1);
552
553 m_HeaderRecords.clear();
554 auto &fitsKeywords = image.fitsKeywords();
555 for(auto &fitsKeyword : fitsKeywords)
556 {
557 m_HeaderRecords.push_back({QString::fromStdString(fitsKeyword.name), QString::fromStdString(fitsKeyword.value), QString::fromStdString(fitsKeyword.comment)});
558 }
559
560 setupWCSParams();
561
562 m_ImageBufferSize = image.imageDataSize();
563 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
564 std::memcpy(m_ImageBuffer, image.imageData(), m_ImageBufferSize);
565
566 calculateStats(false, false);
567 loadWCS();
568 }
569 catch (LibXISF::Error &error)
570 {
571 m_LastError = i18n("XISF file open error: ") + error.what();
572 return false;
573 }
574 return true;
575#else
576 Q_UNUSED(buffer)
577 return false;
578#endif
579
580}
581
582bool FITSData::saveXISFImage(const QString &newFilename)
583{
584#ifdef HAVE_XISF
585 try
586 {
587 LibXISF::XISFWriter xisfWriter;
588 LibXISF::Image image;
589 image.setGeometry(m_Statistics.width, m_Statistics.height, m_Statistics.channels);
590 if (m_Statistics.channels > 1)
591 image.setColorSpace(LibXISF::Image::RGB);
592
593 switch (m_FITSBITPIX)
594 {
595 case BYTE_IMG:
596 image.setSampleFormat(LibXISF::Image::UInt8);
597 break;
598 case USHORT_IMG:
599 image.setSampleFormat(LibXISF::Image::UInt16);
600 break;
601 case ULONG_IMG:
602 image.setSampleFormat(LibXISF::Image::UInt32);
603 break;
604 case FLOAT_IMG:
605 image.setSampleFormat(LibXISF::Image::Float32);
606 break;
607 default:
608 m_LastError = i18n("Bit depth %1 is not supported.", m_FITSBITPIX);
609 qCCritical(KSTARS_FITS) << m_LastError;
610 return false;
611 }
612
613 std::memcpy(image.imageData(), m_ImageBuffer, m_ImageBufferSize);
614 for (auto &fitsKeyword : m_HeaderRecords)
615 image.addFITSKeyword({fitsKeyword.key.toUtf8().data(), fitsKeyword.value.toString().toUtf8().data(), fitsKeyword.comment.toUtf8().data()});
616
617 xisfWriter.writeImage(image);
618 xisfWriter.save(newFilename.toLocal8Bit().data());
619 m_Filename = newFilename;
620 }
621 catch (LibXISF::Error &err)
622 {
623 m_LastError = i18n("Error saving XISF image") + err.what();
624 return false;
625 }
626 return true;
627#else
628 Q_UNUSED(newFilename)
629 return false;
630#endif
631}
632
633bool FITSData::loadCanonicalImage(const QByteArray &buffer)
634{
635 QImage imageFromFile;
636 if (!buffer.isEmpty())
637 {
638 if(!imageFromFile.loadFromData(reinterpret_cast<const uint8_t*>(buffer.data()), buffer.size()))
639 {
640 qCCritical(KSTARS_FITS) << "Failed to open image.";
641 return false;
642 }
643
644 m_Statistics.size = buffer.size();
645 }
646 else
647 {
648 if(!imageFromFile.load(m_Filename.toLocal8Bit()))
649 {
650 qCCritical(KSTARS_FITS) << "Failed to open image.";
651 return false;
652 }
653
654 m_Statistics.size = QFileInfo(m_Filename).size();
655 }
656
657 //imageFromFile = imageFromFile.convertToFormat(QImage::Format_RGB32);
658
659 // Note: This will need to be changed. I think QT only loads 8 bpp images.
660 // Also the depth method gives the total bits per pixel in the image not just the bits per
661 // pixel in each channel.
662 m_FITSBITPIX = BYTE_IMG;
663 switch (m_FITSBITPIX)
664 {
665 case BYTE_IMG:
666 m_Statistics.dataType = TBYTE;
667 m_Statistics.bytesPerPixel = sizeof(uint8_t);
668 break;
669 case SHORT_IMG:
670 // Read SHORT image as USHORT
671 m_Statistics.dataType = TUSHORT;
672 m_Statistics.bytesPerPixel = sizeof(int16_t);
673 break;
674 case USHORT_IMG:
675 m_Statistics.dataType = TUSHORT;
676 m_Statistics.bytesPerPixel = sizeof(uint16_t);
677 break;
678 case LONG_IMG:
679 // Read LONG image as ULONG
680 m_Statistics.dataType = TULONG;
681 m_Statistics.bytesPerPixel = sizeof(int32_t);
682 break;
683 case ULONG_IMG:
684 m_Statistics.dataType = TULONG;
685 m_Statistics.bytesPerPixel = sizeof(uint32_t);
686 break;
687 case FLOAT_IMG:
688 m_Statistics.dataType = TFLOAT;
689 m_Statistics.bytesPerPixel = sizeof(float);
690 break;
691 case LONGLONG_IMG:
692 m_Statistics.dataType = TLONGLONG;
693 m_Statistics.bytesPerPixel = sizeof(int64_t);
694 break;
695 case DOUBLE_IMG:
696 m_Statistics.dataType = TDOUBLE;
697 m_Statistics.bytesPerPixel = sizeof(double);
698 break;
699 default:
700 m_LastError = QString("Bit depth %1 is not supported.").arg(m_FITSBITPIX);
701 qCCritical(KSTARS_FITS) << m_LastError;
702 return false;
703 }
704
705 m_Statistics.width = static_cast<uint16_t>(imageFromFile.width());
706 m_Statistics.height = static_cast<uint16_t>(imageFromFile.height());
707 m_Statistics.channels = imageFromFile.format() == QImage::Format_Grayscale8 ? 1 : 3;
708 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
709 clearImageBuffers();
710 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * static_cast<uint16_t>
711 (m_Statistics.bytesPerPixel);
712 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
713 if (m_ImageBuffer == nullptr)
714 {
715 m_LastError = i18n("FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
716 qCCritical(KSTARS_FITS) << m_LastError;
717 clearImageBuffers();
718 return false;
719 }
720
721 if (m_Statistics.channels == 1)
722 {
723 memcpy(m_ImageBuffer, imageFromFile.bits(), m_ImageBufferSize);
724 }
725 else
726 {
727
728 auto debayered_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
729 auto * original_bayered_buffer = reinterpret_cast<uint8_t *>(imageFromFile.bits());
730
731 // Data in RGBA, with bytes in the order of B,G,R we need to copy them into 3 layers for FITS
732 uint8_t * rBuff = debayered_buffer;
733 uint8_t * gBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height);
734 uint8_t * bBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
735
736 int imax = m_Statistics.samples_per_channel * 4 - 4;
737 for (int i = 0; i <= imax; i += 4)
738 {
739 *rBuff++ = original_bayered_buffer[i + 2];
740 *gBuff++ = original_bayered_buffer[i + 1];
741 *bBuff++ = original_bayered_buffer[i + 0];
742 }
743 }
744
745 m_HeaderRecords.clear();
746 setupWCSParams();
747
748 calculateStats(false, false);
749 loadWCS();
750 return true;
751}
752
753bool FITSData::loadRAWImage(const QByteArray &buffer)
754{
755#if !defined(KSTARS_LITE) && !defined(HAVE_LIBRAW)
756 m_LastError = i18n("Unable to find dcraw and cjpeg. Please install the required tools to convert CR2/NEF to JPEG.");
757 return false;
758#else
759
760 int ret = 0;
761 // Creation of image processing object
762 LibRaw RawProcessor;
763
764 // Let us open the file/buffer
765 if (buffer.isEmpty())
766 {
767 // Open file
768 if ((ret = RawProcessor.open_file(m_Filename.toLocal8Bit().constData())) != LIBRAW_SUCCESS)
769 {
770 m_LastError = i18n("Cannot open file %1: %2", m_Filename, libraw_strerror(ret));
771 RawProcessor.recycle();
772 return false;
773 }
774
775 m_Statistics.size = QFileInfo(m_Filename).size();
776 }
777 // Open Buffer
778 else
779 {
780 if ((ret = RawProcessor.open_buffer(const_cast<void *>(reinterpret_cast<const void *>(buffer.data())),
781 buffer.size())) != LIBRAW_SUCCESS)
782 {
783 m_LastError = i18n("Cannot open buffer: %1", libraw_strerror(ret));
784 RawProcessor.recycle();
785 return false;
786 }
787
788 m_Statistics.size = buffer.size();
789 }
790
791 // Let us unpack the thumbnail
792 if ((ret = RawProcessor.unpack()) != LIBRAW_SUCCESS)
793 {
794 m_LastError = i18n("Cannot unpack_thumb: %1", libraw_strerror(ret));
795 RawProcessor.recycle();
796 return false;
797 }
798
799 if ((ret = RawProcessor.dcraw_process()) != LIBRAW_SUCCESS)
800 {
801 m_LastError = i18n("Cannot dcraw_process: %1", libraw_strerror(ret));
802 RawProcessor.recycle();
803 return false;
804 }
805
806 libraw_processed_image_t *image = RawProcessor.dcraw_make_mem_image(&ret);
807 if (ret != LIBRAW_SUCCESS)
808 {
809 m_LastError = i18n("Cannot load to memory: %1", libraw_strerror(ret));
810 RawProcessor.recycle();
811 return false;
812 }
813
814 RawProcessor.recycle();
815
816 m_Statistics.bytesPerPixel = image->bits / 8;
817 // We only support two types now
818 if (m_Statistics.bytesPerPixel == 1)
819 m_Statistics.dataType = TBYTE;
820 else
821 m_Statistics.dataType = TUSHORT;
822 m_Statistics.width = image->width;
823 m_Statistics.height = image->height;
824 m_Statistics.channels = image->colors;
825 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
826 clearImageBuffers();
827 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
828 m_ImageBuffer = new uint8_t[m_ImageBufferSize];
829 if (m_ImageBuffer == nullptr)
830 {
831 m_LastError = i18n("FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
832 qCCritical(KSTARS_FITS) << m_LastError;
833 libraw_dcraw_clear_mem(image);
834 clearImageBuffers();
835 return false;
836 }
837
838 auto destination_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
839 auto source_buffer = reinterpret_cast<uint8_t *>(image->data);
840
841 // For mono, we memcpy directly
842 if (image->colors == 1)
843 {
844 memcpy(destination_buffer, source_buffer, m_ImageBufferSize);
845 }
846 else
847 {
848 // Data in RGB24, with bytes in the order of R,G,B. We copy them copy them into 3 layers for FITS
849 uint8_t * rBuff = destination_buffer;
850 uint8_t * gBuff = destination_buffer + (m_Statistics.width * m_Statistics.height);
851 uint8_t * bBuff = destination_buffer + (m_Statistics.width * m_Statistics.height * 2);
852
853 int imax = m_Statistics.samples_per_channel * 3 - 3;
854 for (int i = 0; i <= imax; i += 3)
855 {
856 *rBuff++ = source_buffer[i + 0];
857 *gBuff++ = source_buffer[i + 1];
858 *bBuff++ = source_buffer[i + 2];
859 }
860 }
861 libraw_dcraw_clear_mem(image);
862
863 m_HeaderRecords.clear();
864 setupWCSParams();
865
866 calculateStats(false, false);
867 loadWCS();
868 return true;
869#endif
870}
871
872bool FITSData::saveImage(const QString &newFilename)
873{
874 if (newFilename == m_Filename)
875 return true;
876
877 const QString ext = QFileInfo(newFilename).suffix();
878
879 if (ext == "jpg" || ext == "png")
880 {
881 double min, max;
882 QImage fitsImage = FITSToImage(newFilename);
883 getMinMax(&min, &max);
884
885 if (min == max)
886 {
887 fitsImage.fill(Qt::white);
888 }
889 else if (channels() == 1)
890 {
891 fitsImage = QImage(width(), height(), QImage::Format_Indexed8);
892
893 fitsImage.setColorCount(256);
894 for (int i = 0; i < 256; i++)
895 fitsImage.setColor(i, qRgb(i, i, i));
896 }
897 else
898 {
899 fitsImage = QImage(width(), height(), QImage::Format_RGB32);
900 }
901
902 double dataMin = m_Statistics.mean[0] - m_Statistics.stddev[0];
903 double dataMax = m_Statistics.mean[0] + m_Statistics.stddev[0] * 3;
904
905 double bscale = 255. / (dataMax - dataMin);
906 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
907
908 // Long way to do this since we do not want to use templated functions here
909 switch (m_Statistics.dataType)
910 {
911 case TBYTE:
912 convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
913 break;
914
915 case TSHORT:
916 convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
917 break;
918
919 case TUSHORT:
920 convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
921 break;
922
923 case TLONG:
924 convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
925 break;
926
927 case TULONG:
928 convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
929 break;
930
931 case TFLOAT:
932 convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
933 break;
934
935 case TLONGLONG:
936 convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
937 break;
938
939 case TDOUBLE:
940 convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
941 break;
942
943 default:
944 break;
945 }
946
947 fitsImage.save(newFilename, ext.toLatin1().constData());
948
949 m_Filename = newFilename;
950 qCInfo(KSTARS_FITS) << "Saved image file:" << m_Filename;
951 return true;
952 }
953
954 if (m_isCompressed)
955 {
956 m_LastError = i18n("Saving compressed files is not supported.");
957 return false;
958 }
959
960 int status = 0;
961 long nelements;
962 fitsfile * new_fptr;
963
964 if (ext == "xisf")
965 {
966 if(fptr)
967 {
968 fits_close_file(fptr, &status);
969 fptr = nullptr;
970 }
971 rotCounter = flipHCounter = flipVCounter = 0;
972 return saveXISFImage(newFilename);
973 }
974
975 // if (HasDebayer && m_Filename.isEmpty() == false)
976 // {
977 // fits_flush_file(fptr, &status);
978 // /* close current file */
979 // if (fits_close_file(fptr, &status))
980 // {
981 // recordLastError(status);
982 // return status;
983 // }
984
985 // // Skip "!" in the beginning of the new file name
986 // QString finalFileName(newFilename);
987
988 // finalFileName.remove('!');
989
990 // // Remove first otherwise copy will fail below if file exists
991 // QFile::remove(finalFileName);
992
993 // if (!QFile::copy(m_Filename, finalFileName))
994 // {
995 // qCCritical(KSTARS_FITS()) << "FITS: Failed to copy " << m_Filename << " to " << finalFileName;
996 // fptr = nullptr;
997 // return false;
998 // }
999
1000 // m_Filename = finalFileName;
1001
1002 // // Use open diskfile as it does not use extended file names which has problems opening
1003 // // files with [ ] or ( ) in their names.
1004 // fits_open_diskfile(&fptr, m_Filename.toLocal8Bit(), READONLY, &status);
1005 // fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status);
1006
1007 // return true;
1008 // }
1009
1010 // Read the image back into buffer in case we debyayed
1011 nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
1012
1013 /* close current file */
1014 if (fptr && fits_close_file(fptr, &status))
1015 {
1016 m_LastError = i18n("Failed to close file: %1", fitsErrorToString(status));
1017 return false;
1018 }
1019
1020 /* Create a new File, overwriting existing*/
1021 if (fits_create_file(&new_fptr, QString("!%1").arg(newFilename).toLocal8Bit(), &status))
1022 {
1023 m_LastError = i18n("Failed to create file: %1", fitsErrorToString(status));
1024 return status;
1025 }
1026
1027 status = 0;
1028
1029 fptr = new_fptr;
1030
1031 // Create image
1032 long naxis = m_Statistics.channels == 1 ? 2 : 3;
1033 long naxes[3] = {m_Statistics.width, m_Statistics.height, naxis};
1034
1035 // JM 2020-12-28: Here we to use bitpix values
1036 if (fits_create_img(fptr, m_FITSBITPIX, naxis, naxes, &status))
1037 {
1038 m_LastError = i18n("Failed to create image: %1", fitsErrorToString(status));
1039 return false;
1040 }
1041
1042 /* Write keywords */
1043
1044 // Minimum
1045 if (fits_update_key(fptr, TDOUBLE, "DATAMIN", &(m_Statistics.min), "Minimum value", &status))
1046 {
1047 m_LastError = i18n("Failed to update key: %1", fitsErrorToString(status));
1048 return false;
1049 }
1050
1051 // Maximum
1052 if (fits_update_key(fptr, TDOUBLE, "DATAMAX", &(m_Statistics.max), "Maximum value", &status))
1053 {
1054 m_LastError = i18n("Failed to update key: %1", fitsErrorToString(status));
1055 return false;
1056 }
1057
1058 // KStars Min, for 3 channels
1059 fits_write_key(fptr, TDOUBLE, "MIN1", &m_Statistics.min[0], "Min Channel 1", &status);
1060 if (channels() > 1)
1061 {
1062 fits_write_key(fptr, TDOUBLE, "MIN2", &m_Statistics.min[1], "Min Channel 2", &status);
1063 fits_write_key(fptr, TDOUBLE, "MIN3", &m_Statistics.min[2], "Min Channel 3", &status);
1064 }
1065
1066 // KStars max, for 3 channels
1067 fits_write_key(fptr, TDOUBLE, "MAX1", &m_Statistics.max[0], "Max Channel 1", &status);
1068 if (channels() > 1)
1069 {
1070 fits_write_key(fptr, TDOUBLE, "MAX2", &m_Statistics.max[1], "Max Channel 2", &status);
1071 fits_write_key(fptr, TDOUBLE, "MAX3", &m_Statistics.max[2], "Max Channel 3", &status);
1072 }
1073
1074 // Mean
1075 if (m_Statistics.mean[0] > 0)
1076 {
1077 fits_write_key(fptr, TDOUBLE, "MEAN1", &m_Statistics.mean[0], "Mean Channel 1", &status);
1078 if (channels() > 1)
1079 {
1080 fits_write_key(fptr, TDOUBLE, "MEAN2", &m_Statistics.mean[1], "Mean Channel 2", &status);
1081 fits_write_key(fptr, TDOUBLE, "MEAN3", &m_Statistics.mean[2], "Mean Channel 3", &status);
1082 }
1083 }
1084
1085 // Median
1086 if (m_Statistics.median[0] > 0)
1087 {
1088 fits_write_key(fptr, TDOUBLE, "MEDIAN1", &m_Statistics.median[0], "Median Channel 1", &status);
1089 if (channels() > 1)
1090 {
1091 fits_write_key(fptr, TDOUBLE, "MEDIAN2", &m_Statistics.median[1], "Median Channel 2", &status);
1092 fits_write_key(fptr, TDOUBLE, "MEDIAN3", &m_Statistics.median[2], "Median Channel 3", &status);
1093 }
1094 }
1095
1096 // Standard Deviation
1097 if (m_Statistics.stddev[0] > 0)
1098 {
1099 fits_write_key(fptr, TDOUBLE, "STDDEV1", &m_Statistics.stddev[0], "Standard Deviation Channel 1", &status);
1100 if (channels() > 1)
1101 {
1102 fits_write_key(fptr, TDOUBLE, "STDDEV2", &m_Statistics.stddev[1], "Standard Deviation Channel 2", &status);
1103 fits_write_key(fptr, TDOUBLE, "STDDEV3", &m_Statistics.stddev[2], "Standard Deviation Channel 3", &status);
1104 }
1105 }
1106
1107 // Skip first 10 standard records and copy the rest.
1108 for (int i = 10; i < m_HeaderRecords.count(); i++)
1109 {
1110 QString key = m_HeaderRecords[i].key;
1111 const char *comment = m_HeaderRecords[i].comment.toLatin1().constBegin();
1112 QVariant value = m_HeaderRecords[i].value;
1113
1114 switch (value.type())
1115 {
1116 case QVariant::Int:
1117 {
1118 int number = value.toInt();
1119 fits_write_key(fptr, TINT, key.toLatin1().constData(), &number, comment, &status);
1120 }
1121 break;
1122
1123 case QVariant::Double:
1124 {
1125 double number = value.toDouble();
1126 fits_write_key(fptr, TDOUBLE, key.toLatin1().constData(), &number, comment, &status);
1127 }
1128 break;
1129
1130 case QVariant::String:
1131 default:
1132 {
1133 char valueBuffer[256] = {0};
1134 strncpy(valueBuffer, value.toString().toLatin1().constData(), 256 - 1);
1135 fits_write_key(fptr, TSTRING, key.toLatin1().constData(), valueBuffer, comment, &status);
1136 }
1137 }
1138 }
1139
1140 // ISO Date
1141 if (fits_write_date(fptr, &status))
1142 {
1143 m_LastError = i18n("Failed to update date: %1", fitsErrorToString(status));
1144 return false;
1145 }
1146
1147 QString history =
1148 QString("Modified by KStars on %1").arg(QDateTime::currentDateTime().toString("yyyy-MM-ddThh:mm:ss"));
1149 // History
1150 if (fits_write_history(fptr, history.toLatin1(), &status))
1151 {
1152 m_LastError = i18n("Failed to update history: %1", fitsErrorToString(status));
1153 return false;
1154 }
1155
1156 int rot = 0, mirror = 0;
1157 if (rotCounter > 0)
1158 {
1159 rot = (90 * rotCounter) % 360;
1160 if (rot < 0)
1161 rot += 360;
1162 }
1163 if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0)
1164 mirror = 1;
1165
1166 if ((rot != 0) || (mirror != 0))
1167 rotWCSFITS(rot, mirror);
1168
1169 rotCounter = flipHCounter = flipVCounter = 0;
1170
1171 // Here we need to use the actual data type
1172 if (fits_write_img(fptr, m_Statistics.dataType, 1, nelements, m_ImageBuffer, &status))
1173 {
1174 m_LastError = i18n("Failed to write image: %1", fitsErrorToString(status));
1175 return false;
1176 }
1177
1178 m_Filename = newFilename;
1179
1180 fits_flush_file(fptr, &status);
1181
1182 qCInfo(KSTARS_FITS) << "Saved FITS file:" << m_Filename;
1183
1184 return true;
1185}
1186
1187void FITSData::clearImageBuffers()
1188{
1189 delete[] m_ImageBuffer;
1190 m_ImageBuffer = nullptr;
1191 if(m_ImageRoiBuffer != nullptr )
1192 {
1193 delete[] m_ImageRoiBuffer;
1194 m_ImageRoiBuffer = nullptr;
1195
1196 }
1197 //m_BayerBuffer = nullptr;
1198}
1199
1200void FITSData::makeRoiBuffer(QRect roi)
1201{
1202 if (!roi.isValid())
1203 return;
1204
1205 uint32_t channelSize = roi.height() * roi.width();
1206 if(channelSize > m_Statistics.samples_per_channel || channelSize == 1)
1207 {
1208 return;
1209 }
1210 if(m_ImageRoiBuffer != nullptr )
1211 {
1212 delete[] m_ImageRoiBuffer;
1213 m_ImageRoiBuffer = nullptr;
1214 }
1215 int xoffset = roi.topLeft().x() - 1;
1216 int yoffset = roi.topLeft().y() - 1;
1217 uint32_t bpp = m_Statistics.bytesPerPixel;
1218 m_ImageRoiBuffer = new uint8_t[roi.height()*roi.width()*m_Statistics.channels * m_Statistics.bytesPerPixel]();
1219 for(int n = 0 ; n < m_Statistics.channels ; n++)
1220 {
1221 for(int i = 0; i < roi.height(); i++)
1222 {
1223 size_t i1 = n * channelSize * bpp + i * roi.width() * bpp;
1224 size_t i2 = n * m_Statistics.samples_per_channel * bpp + (yoffset + i) * width() * bpp + xoffset * bpp;
1225 memcpy(&m_ImageRoiBuffer[i1],
1226 &m_ImageBuffer[i2],
1227 roi.width() * bpp);
1228 }
1229
1230 }
1231 memcpy(&m_ROIStatistics, &m_Statistics, sizeof(FITSImage::Statistic));
1232 m_ROIStatistics.samples_per_channel = roi.height() * roi.width();
1233 m_ROIStatistics.width = roi.width();
1234 m_ROIStatistics.height = roi.height();
1235 calculateStats(false, true);
1236}
1237
1238void FITSData::setupWCSParams()
1239{
1240 FITSImage::Solution solution;
1241 if (parseSolution(solution))
1242 {
1243 const bool eastToTheRight = solution.parity == FITSImage::POSITIVE ? false : true;
1244 updateWCSHeaderData(solution.orientation, solution.ra, solution.dec, solution.pixscale, eastToTheRight);
1245
1246 QVariant value;
1247 bool validObservationDate = false;
1248
1249 if (getRecordValue("DATE-OBS", value))
1250 {
1251 QString tsString(value.toString());
1252 tsString = tsString.remove('\'').trimmed();
1253 // Add Zulu time to indicate UTC
1254 tsString += "Z";
1255
1257
1258 if (ts.isValid())
1259 {
1260 validObservationDate = true;
1261 m_DateTime = KStarsDateTime(ts.date(), ts.time());
1262 }
1263 }
1264
1265 // Set to current datetime if no valid observation datetime. Not 100% accurate but close enough in most cases
1266 if (!validObservationDate)
1268 }
1269}
1270
1271void FITSData::calculateStats(bool refresh, bool roi)
1272{
1273 // Calculate min max
1274 if(roi == false)
1275 {
1276 calculateMinMax(refresh);
1277 calculateMedian(refresh);
1278
1279 // Try to read mean/median/stddev if in file
1280 if (refresh == false && fptr)
1281 {
1282 int status = 0;
1283 int nfound = 0;
1284 if (fits_read_key_dbl(fptr, "MEAN1", &m_Statistics.mean[0], nullptr, &status) == 0)
1285 nfound++;
1286 // NB. These could fail if missing, which is OK.
1287 fits_read_key_dbl(fptr, "MEAN2", & m_Statistics.mean[1], nullptr, &status);
1288 fits_read_key_dbl(fptr, "MEAN3", &m_Statistics.mean[2], nullptr, &status);
1289
1290 status = 0;
1291 if (fits_read_key_dbl(fptr, "STDDEV1", &m_Statistics.stddev[0], nullptr, &status) == 0)
1292 nfound++;
1293 // NB. These could fail if missing, which is OK.
1294 fits_read_key_dbl(fptr, "STDDEV2", &m_Statistics.stddev[1], nullptr, &status);
1295 fits_read_key_dbl(fptr, "STDDEV3", &m_Statistics.stddev[2], nullptr, &status);
1296
1297 // If all is OK, we're done
1298 if (nfound == 2)
1299 return;
1300 }
1301
1302 // Get standard deviation and mean in one run
1303 switch (m_Statistics.dataType)
1304 {
1305 case TBYTE:
1306 calculateStdDev<uint8_t>();
1307 break;
1308
1309 case TSHORT:
1310 calculateStdDev<int16_t>();
1311 break;
1312
1313 case TUSHORT:
1314 calculateStdDev<uint16_t>();
1315 break;
1316
1317 case TLONG:
1318 calculateStdDev<int32_t>();
1319 break;
1320
1321 case TULONG:
1322 calculateStdDev<uint32_t>();
1323 break;
1324
1325 case TFLOAT:
1326 calculateStdDev<float>();
1327 break;
1328
1329 case TLONGLONG:
1330 calculateStdDev<int64_t>();
1331 break;
1332
1333 case TDOUBLE:
1334 calculateStdDev<double>();
1335 break;
1336
1337 default:
1338 return;
1339 }
1340
1341 // FIXME That's not really SNR, must implement a proper solution for this value
1342 m_Statistics.SNR = m_Statistics.mean[0] / m_Statistics.stddev[0];
1343 }
1344 else
1345 {
1346 calculateMinMax(refresh, roi);
1347 calculateMedian(refresh, roi);
1348
1349 switch (m_ROIStatistics.dataType)
1350 {
1351 case TBYTE:
1352 calculateStdDev<uint8_t>(roi);
1353 break;
1354
1355 case TSHORT:
1356 calculateStdDev<int16_t>(roi);
1357 break;
1358
1359 case TUSHORT:
1360 calculateStdDev<uint16_t>(roi);
1361 break;
1362
1363 case TLONG:
1364 calculateStdDev<int32_t>(roi);
1365 break;
1366
1367 case TULONG:
1368 calculateStdDev<uint32_t>(roi);
1369 break;
1370
1371 case TFLOAT:
1372 calculateStdDev<float>(roi);
1373 break;
1374
1375 case TLONGLONG:
1376 calculateStdDev<int64_t>(roi);
1377 break;
1378
1379 case TDOUBLE:
1380 calculateStdDev<double>(roi);
1381 break;
1382
1383 default:
1384 return;
1385 }
1386 }
1387}
1388
1389void FITSData::calculateMinMax(bool refresh, bool roi)
1390{
1391 if(roi == false)
1392 {
1393 int status, nfound = 0;
1394
1395 status = 0;
1396
1397 // Only fetch from header if we have a single channel
1398 // Otherwise, calculate manually.
1399 if (fptr != nullptr && !refresh)
1400 {
1401
1402 status = 0;
1403
1404 if (fits_read_key_dbl(fptr, "DATAMIN", &(m_Statistics.min[0]), nullptr, &status) == 0)
1405 nfound++;
1406 else if (fits_read_key_dbl(fptr, "MIN1", &(m_Statistics.min[0]), nullptr, &status) == 0)
1407 nfound++;
1408
1409 // NB. These could fail if missing, which is OK.
1410 fits_read_key_dbl(fptr, "MIN2", &m_Statistics.min[1], nullptr, &status);
1411 fits_read_key_dbl(fptr, "MIN3", &m_Statistics.min[2], nullptr, &status);
1412
1413 status = 0;
1414
1415 if (fits_read_key_dbl(fptr, "DATAMAX", &(m_Statistics.max[0]), nullptr, &status) == 0)
1416 nfound++;
1417 else if (fits_read_key_dbl(fptr, "MAX1", &(m_Statistics.max[0]), nullptr, &status) == 0)
1418 nfound++;
1419
1420 // NB. These could fail if missing, which is OK.
1421 fits_read_key_dbl(fptr, "MAX2", &m_Statistics.max[1], nullptr, &status);
1422 fits_read_key_dbl(fptr, "MAX3", &m_Statistics.max[2], nullptr, &status);
1423
1424 // If we found both keywords, no need to calculate them, unless they are both zeros
1425 if (nfound == 2 && !(m_Statistics.min[0] == 0 && m_Statistics.max[0] == 0))
1426 return;
1427 }
1428
1429 m_Statistics.min[0] = 1.0E30;
1430 m_Statistics.max[0] = -1.0E30;
1431
1432 m_Statistics.min[1] = 1.0E30;
1433 m_Statistics.max[1] = -1.0E30;
1434
1435 m_Statistics.min[2] = 1.0E30;
1436 m_Statistics.max[2] = -1.0E30;
1437
1438 switch (m_Statistics.dataType)
1439 {
1440 case TBYTE:
1441 calculateMinMax<uint8_t>();
1442 break;
1443
1444 case TSHORT:
1445 calculateMinMax<int16_t>();
1446 break;
1447
1448 case TUSHORT:
1449 calculateMinMax<uint16_t>();
1450 break;
1451
1452 case TLONG:
1453 calculateMinMax<int32_t>();
1454 break;
1455
1456 case TULONG:
1457 calculateMinMax<uint32_t>();
1458 break;
1459
1460 case TFLOAT:
1461 calculateMinMax<float>();
1462 break;
1463
1464 case TLONGLONG:
1465 calculateMinMax<int64_t>();
1466 break;
1467
1468 case TDOUBLE:
1469 calculateMinMax<double>();
1470 break;
1471
1472 default:
1473 break;
1474 }
1475 }
1476 else
1477 {
1478 m_ROIStatistics.min[0] = 1.0E30;
1479 m_ROIStatistics.max[0] = -1.0E30;
1480
1481 m_ROIStatistics.min[1] = 1.0E30;
1482 m_ROIStatistics.max[1] = -1.0E30;
1483
1484 m_ROIStatistics.min[2] = 1.0E30;
1485 m_ROIStatistics.max[2] = -1.0E30;
1486
1487 switch (m_Statistics.dataType)
1488 {
1489 case TBYTE:
1490 calculateMinMax<uint8_t>(roi);
1491 break;
1492
1493 case TSHORT:
1494 calculateMinMax<int16_t>(roi);
1495 break;
1496
1497 case TUSHORT:
1498 calculateMinMax<uint16_t>(roi);
1499 break;
1500
1501 case TLONG:
1502 calculateMinMax<int32_t>(roi);
1503 break;
1504
1505 case TULONG:
1506 calculateMinMax<uint32_t>(roi);
1507 break;
1508
1509 case TFLOAT:
1510 calculateMinMax<float>(roi);
1511 break;
1512
1513 case TLONGLONG:
1514 calculateMinMax<int64_t>(roi);
1515 break;
1516
1517 case TDOUBLE:
1518 calculateMinMax<double>(roi);
1519 break;
1520
1521 default:
1522 break;
1523 }
1524
1525 }
1526}
1527
1528void FITSData::calculateMedian(bool refresh, bool roi)
1529{
1530 if(!roi)
1531 {
1532 int status, nfound = 0;
1533
1534 status = 0;
1535
1536 // Only fetch from header if we have a single channel
1537 // Otherwise, calculate manually.
1538 if (fptr != nullptr && !refresh)
1539 {
1540 status = 0;
1541 if (fits_read_key_dbl(fptr, "MEDIAN1", &m_Statistics.median[0], nullptr, &status) == 0)
1542 nfound++;
1543
1544 // NB. These could fail if missing, which is OK.
1545 fits_read_key_dbl(fptr, "MEDIAN2", &m_Statistics.median[1], nullptr, &status);
1546 fits_read_key_dbl(fptr, "MEDIAN3", &m_Statistics.median[2], nullptr, &status);
1547
1548 if (nfound == 1)
1549 return;
1550 }
1551
1552 m_Statistics.median[RED_CHANNEL] = 0;
1553 m_Statistics.median[GREEN_CHANNEL] = 0;
1554 m_Statistics.median[BLUE_CHANNEL] = 0;
1555
1556 switch (m_Statistics.dataType)
1557 {
1558 case TBYTE:
1559 calculateMedian<uint8_t>();
1560 break;
1561
1562 case TSHORT:
1563 calculateMedian<int16_t>();
1564 break;
1565
1566 case TUSHORT:
1567 calculateMedian<uint16_t>();
1568 break;
1569
1570 case TLONG:
1571 calculateMedian<int32_t>();
1572 break;
1573
1574 case TULONG:
1575 calculateMedian<uint32_t>();
1576 break;
1577
1578 case TFLOAT:
1579 calculateMedian<float>();
1580 break;
1581
1582 case TLONGLONG:
1583 calculateMedian<int64_t>();
1584 break;
1585
1586 case TDOUBLE:
1587 calculateMedian<double>();
1588 break;
1589
1590 default:
1591 break;
1592 }
1593 }
1594 else
1595 {
1596 m_ROIStatistics.median[RED_CHANNEL] = 0;
1597 m_ROIStatistics.median[GREEN_CHANNEL] = 0;
1598 m_ROIStatistics.median[BLUE_CHANNEL] = 0;
1599
1600 switch (m_ROIStatistics.dataType)
1601 {
1602 case TBYTE:
1603 calculateMedian<uint8_t>(roi);
1604 break;
1605
1606 case TSHORT:
1607 calculateMedian<int16_t>(roi);
1608 break;
1609
1610 case TUSHORT:
1611 calculateMedian<uint16_t>(roi);
1612 break;
1613
1614 case TLONG:
1615 calculateMedian<int32_t>(roi);
1616 break;
1617
1618 case TULONG:
1619 calculateMedian<uint32_t>(roi);
1620 break;
1621
1622 case TFLOAT:
1623 calculateMedian<float>(roi);
1624 break;
1625
1626 case TLONGLONG:
1627 calculateMedian<int64_t>(roi);
1628 break;
1629
1630 case TDOUBLE:
1631 calculateMedian<double>(roi);
1632 break;
1633
1634 default:
1635 break;
1636 }
1637
1638 }
1639}
1640
1641template <typename T>
1642void FITSData::calculateMedian(bool roi)
1643{
1644 auto * buffer = reinterpret_cast<T *>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1645 const uint32_t maxMedianSize = 500000;
1646 uint32_t medianSize = roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel;
1647 uint8_t downsample = 1;
1648 if (medianSize > maxMedianSize)
1649 {
1650 downsample = (static_cast<double>(medianSize) / maxMedianSize) + 0.999;
1651 medianSize /= downsample;
1652 }
1653 // Ideally samples would be declared like this...
1654 // std::vector<T> samples;
1655 // Unfortunately this doesn't compile on Mac - see the comments in robuststatistics.cpp for more details
1656 // So for now declare samples like this...
1657 std::vector<int32_t> samples;
1658 samples.reserve(medianSize);
1659
1660 for (uint8_t n = 0; n < m_Statistics.channels; n++)
1661 {
1662 auto *oneChannel = buffer + n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1663 for (uint32_t upto = 0; upto < (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1664 upto += downsample)
1665 samples.push_back(oneChannel[upto]);
1666 auto median = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_MEDIAN, samples);
1667 roi ? m_ROIStatistics.median[n] = median : m_Statistics.median[n] = median;
1668 }
1669}
1670
1671template <typename T>
1672QPair<T, T> FITSData::getParitionMinMax(uint32_t start, uint32_t stride, bool roi)
1673{
1674 auto * buffer = reinterpret_cast<T *>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1675 T min = std::numeric_limits<T>::max();
1676 T max = std::numeric_limits<T>::min();
1677
1678 uint32_t end = start + stride;
1679
1680 for (uint32_t i = start; i < end; i++)
1681 {
1682 min = qMin(buffer[i], min);
1683 max = qMax(buffer[i], max);
1684 // if (buffer[i] < min)
1685 // min = buffer[i];
1686 // else if (buffer[i] > max)
1687 // max = buffer[i];
1688 }
1689
1690 return qMakePair(min, max);
1691}
1692
1693template <typename T>
1694void FITSData::calculateMinMax(bool roi)
1695{
1696 T min = std::numeric_limits<T>::max();
1697 T max = std::numeric_limits<T>::min();
1698
1699
1700 // Create N threads
1701 const uint8_t nThreads = 16;
1702
1703 for (int n = 0; n < m_Statistics.channels; n++)
1704 {
1705 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1706
1707 // Calculate how many elements we process per thread
1708 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1709
1710 // Calculate the final stride since we can have some left over due to division above
1711 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1712 (tStride * nThreads));
1713
1714 // Start location for inspecting elements
1715 uint32_t tStart = cStart;
1716
1717 // List of futures
1718 QList<QFuture<QPair<T, T>>> futures;
1719
1720 for (int i = 0; i < nThreads; i++)
1721 {
1722 // Run threads
1723#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
1724 futures.append(QtConcurrent::run(&FITSData::getParitionMinMax<T>, this, tStart, (i == (nThreads - 1)) ? fStride : tStride,
1725 roi));
1726#else
1727 futures.append(QtConcurrent::run(this, &FITSData::getParitionMinMax<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride,
1728 roi));
1729#endif
1730 tStart += tStride;
1731 }
1732
1733 // Now wait for results
1734 for (int i = 0; i < nThreads; i++)
1735 {
1736 QPair<T, T> result = futures[i].result();
1737 min = qMin(result.first, min);
1738 max = qMax(result.second, max);
1739 }
1740
1741 if(!roi)
1742 {
1743 m_Statistics.min[n] = min;
1744 m_Statistics.max[n] = max;
1745 }
1746 else
1747 {
1748 m_ROIStatistics.min[n] = min;
1749 m_ROIStatistics.max[n] = max;
1750 }
1751 }
1752
1753}
1754
1755// This struct is used when returning results from the threaded getSumAndSquaredSum calculations
1756// used to compute the mean and variance of the image.
1757struct SumData
1758{
1759 double sum;
1760 double squaredSum;
1761 double numSamples;
1762 SumData(double s, double sq, int n) : sum(s), squaredSum(sq), numSamples(n) {}
1763 SumData() : sum(0), squaredSum(0), numSamples(0) {}
1764};
1765
1766template <typename T>
1767SumData getSumAndSquaredSum(uint32_t start, uint32_t stride, uint8_t *buff)
1768{
1769 auto * buffer = reinterpret_cast<T *>(buff);
1770 const uint32_t end = start + stride;
1771 double sum = 0;
1772 double squaredSum = 0;
1773 for (uint32_t i = start; i < end; i++)
1774 {
1775 double sample = buffer[i];
1776 sum += sample;
1777 squaredSum += sample * sample;
1778 }
1779 const double numSamples = end - start;
1780 return SumData(sum, squaredSum, numSamples);
1781}
1782
1783template <typename T>
1784void FITSData::calculateStdDev(bool roi )
1785{
1786 // Create N threads
1787 const uint8_t nThreads = 16;
1788
1789 for (int n = 0; n < m_Statistics.channels; n++)
1790 {
1791 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1792
1793 // Calculate how many elements we process per thread
1794 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1795
1796 // Calculate the final stride since we can have some left over due to division above
1797 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1798 (tStride * nThreads));
1799
1800 // Start location for inspecting elements
1801 uint32_t tStart = cStart;
1802
1803 // List of futures
1804 QList<QFuture<SumData>> futures;
1805
1806 for (int i = 0; i < nThreads; i++)
1807 {
1808 // Run threads
1809 uint8_t *buff = roi ? m_ImageRoiBuffer : m_ImageBuffer;
1810 futures.append(QtConcurrent::run(&getSumAndSquaredSum<T>, tStart,
1811 (i == (nThreads - 1)) ? fStride : tStride, buff));
1812 tStart += tStride;
1813 }
1814
1815 // Now wait for results
1816 double sum = 0, squared_sum = 0;
1817 double numSamples = 0;
1818 for (int i = 0; i < nThreads; i++)
1819 {
1820 const SumData result = futures[i].result();
1821 sum += result.sum;
1822 squared_sum += result.squaredSum;
1823 numSamples += result.numSamples;
1824
1825 }
1826 if (numSamples <= 0) continue;
1827 const double mean = sum / numSamples;
1828 const double variance = squared_sum / numSamples - mean * mean;
1829 if(!roi)
1830 {
1831 m_Statistics.mean[n] = mean;
1832 m_Statistics.stddev[n] = sqrt(variance);
1833 }
1834 else
1835 {
1836 m_ROIStatistics.mean[n] = mean;
1837 m_ROIStatistics.stddev[n] = sqrt(variance);
1838 }
1839 }
1840}
1841
1842QVector<double> FITSData::createGaussianKernel(int size, double sigma)
1843{
1844 QVector<double> kernel(size * size);
1845 kernel.fill(0.0, size * size);
1846
1847 double kernelSum = 0.0;
1848 int fOff = (size - 1) / 2;
1849 double normal = 1.0 / (2.0 * M_PI * sigma * sigma);
1850 for (int y = -fOff; y <= fOff; y++)
1851 {
1852 for (int x = -fOff; x <= fOff; x++)
1853 {
1854 double distance = ((y * y) + (x * x)) / (2.0 * sigma * sigma);
1855 int index = (y + fOff) * size + (x + fOff);
1856 kernel[index] = normal * qExp(-distance);
1857 kernelSum += kernel.at(index);
1858 }
1859 }
1860 for (int y = 0; y < size; y++)
1861 {
1862 for (int x = 0; x < size; x++)
1863 {
1864 int index = y * size + x;
1865 kernel[index] = kernel.at(index) * 1.0 / kernelSum;
1866 }
1867 }
1868
1869 return kernel;
1870}
1871
1872template <typename T>
1873void FITSData::convolutionFilter(const QVector<double> &kernel, int kernelSize)
1874{
1875 T * imagePtr = reinterpret_cast<T *>(m_ImageBuffer);
1876
1877 // Create variable for pixel data for each kernel
1878 T gt = 0;
1879
1880 // This is how much your center pixel is offset from the border of your kernel
1881 int fOff = (kernelSize - 1) / 2;
1882
1883 // Start with the pixel that is offset fOff from top and fOff from the left side
1884 // this is so entire kernel is on your image
1885 for (int offsetY = 0; offsetY < m_Statistics.height; offsetY++)
1886 {
1887 for (int offsetX = 0; offsetX < m_Statistics.width; offsetX++)
1888 {
1889 // reset gray value to 0
1890 gt = 0;
1891 // position of the kernel center pixel
1892 int byteOffset = offsetY * m_Statistics.width + offsetX;
1893
1894 // kernel calculations
1895 for (int filterY = -fOff; filterY <= fOff; filterY++)
1896 {
1897 for (int filterX = -fOff; filterX <= fOff; filterX++)
1898 {
1899 if ((offsetY + filterY) >= 0 && (offsetY + filterY) < m_Statistics.height
1900 && ((offsetX + filterX) >= 0 && (offsetX + filterX) < m_Statistics.width ))
1901 {
1902
1903 int calcOffset = byteOffset + filterX + filterY * m_Statistics.width;
1904 int index = (filterY + fOff) * kernelSize + (filterX + fOff);
1905 double kernelValue = kernel.at(index);
1906 gt += (imagePtr[calcOffset]) * kernelValue;
1907 }
1908 }
1909 }
1910
1911 // set new data in the other byte array for your image data
1912 imagePtr[byteOffset] = gt;
1913 }
1914 }
1915}
1916
1917template <typename T>
1918void FITSData::gaussianBlur(int kernelSize, double sigma)
1919{
1920 // Size must be an odd number!
1921 if (kernelSize % 2 == 0)
1922 {
1923 kernelSize--;
1924 qCInfo(KSTARS_FITS) << "Warning, size must be an odd number, correcting size to " << kernelSize;
1925 }
1926 // Edge must be a positive number!
1927 if (kernelSize < 1)
1928 {
1929 kernelSize = 1;
1930 }
1931
1932 QVector<double> gaussianKernel = createGaussianKernel(kernelSize, sigma);
1933 convolutionFilter<T>(gaussianKernel, kernelSize);
1934}
1935
1936void FITSData::setMinMax(double newMin, double newMax, uint8_t channel)
1937{
1938 m_Statistics.min[channel] = newMin;
1939 m_Statistics.max[channel] = newMax;
1940}
1941
1942bool FITSData::parseHeader()
1943{
1944 char * header = nullptr;
1945 int status = 0, nkeys = 0;
1946
1947 if (fits_hdr2str(fptr, 0, nullptr, 0, &header, &nkeys, &status))
1948 {
1949 fits_report_error(stderr, status);
1950 free(header);
1951 return false;
1952 }
1953
1954 m_HeaderRecords.clear();
1955 QString recordList = QString(header);
1956
1957 for (int i = 0; i < nkeys; i++)
1958 {
1959 Record oneRecord;
1960 // Quotes cause issues for simplified below so we're removing them.
1961 QString record = recordList.mid(i * 80, 80).remove("'");
1963 // If it is only a comment
1964 if (properties.size() == 1)
1965 {
1966 oneRecord.key = properties[0].mid(0, 7);
1967 oneRecord.comment = properties[0].mid(8).simplified();
1968 }
1969 else
1970 {
1971 oneRecord.key = properties[0].simplified();
1972 oneRecord.value = properties[1].simplified();
1973 if (properties.size() > 2)
1974 oneRecord.comment = properties[2].simplified();
1975
1976 // Try to guess the value.
1977 // Test for integer & double. If neither, then leave it as "string".
1978 bool ok = false;
1979
1980 // Is it Integer?
1981 oneRecord.value.toInt(&ok);
1982 if (ok)
1983 oneRecord.value.convert(QMetaType::Int);
1984 else
1985 {
1986 // Is it double?
1987 oneRecord.value.toDouble(&ok);
1988 if (ok)
1989 oneRecord.value.convert(QMetaType::Double);
1990 }
1991 }
1992
1993 m_HeaderRecords.append(oneRecord);
1994 }
1995
1996 free(header);
1997
1998 return true;
1999}
2000
2001bool FITSData::getRecordValue(const QString &key, QVariant &value) const
2002{
2003 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](const Record & oneRecord)
2004 {
2005 return (oneRecord.key == key && oneRecord.value.isValid());
2006 });
2007
2008 if (result != m_HeaderRecords.end())
2009 {
2010 value = (*result).value;
2011 return true;
2012 }
2013 return false;
2014}
2015
2016void FITSData::updateRecordValue(const QString &key, QVariant value, const QString &comment)
2017{
2018 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](const Record & oneRecord)
2019 {
2020 return (oneRecord.key == key && oneRecord.value.isValid());
2021 });
2022
2023 if (result != m_HeaderRecords.end())
2024 {
2025 (*result).value = value;
2026 }
2027 else
2028 {
2029 // Add item as penultimate entry (END is usually the last one).
2030 FITSData::Record record = {key, value.toString(), comment};
2031#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
2032 m_HeaderRecords.insert(std::max(0LL, m_HeaderRecords.size() - 1), record);
2033#else
2034 m_HeaderRecords.insert(std::max(0, m_HeaderRecords.size() - 1), record);
2035#endif
2036 }
2037}
2038
2039bool FITSData::parseSolution(FITSImage::Solution &solution) const
2040{
2041 dms angleValue;
2042 bool raOK = false, deOK = false, coordOK = false, scaleOK = false;
2043 QVariant value;
2044
2045 // Reset all
2046 solution.fieldWidth = solution.fieldHeight = solution.pixscale = solution.ra = solution.dec = -1;
2047
2048 // RA
2049 if (getRecordValue("OBJCTRA", value))
2050 {
2051 angleValue = dms::fromString(value.toString(), false);
2052 solution.ra = angleValue.Degrees();
2053 raOK = true;
2054 }
2055 else if (getRecordValue("RA", value))
2056 {
2057 solution.ra = value.toDouble(&raOK);
2058 }
2059
2060 // DE
2061 if (getRecordValue("OBJCTDEC", value))
2062 {
2063 angleValue = dms::fromString(value.toString(), true);
2064 solution.dec = angleValue.Degrees();
2065 deOK = true;
2066 }
2067 else if (getRecordValue("DEC", value))
2068 {
2069 solution.dec = value.toDouble(&deOK);
2070 }
2071
2072 coordOK = raOK && deOK;
2073
2074 // PixScale
2075 double scale = -1;
2076 if (getRecordValue("SCALE", value))
2077 {
2078 scale = value.toDouble();
2079 }
2080
2081 double focal_length = -1;
2082 if (getRecordValue("FOCALLEN", value))
2083 {
2084 focal_length = value.toDouble();
2085 }
2086
2087 double pixsize1 = -1, pixsize2 = -1;
2088 // Pixel Size 1
2089 if (getRecordValue("PIXSIZE1", value))
2090 {
2091 pixsize1 = value.toDouble();
2092 }
2093 else if (getRecordValue("XPIXSZ", value))
2094 {
2095 pixsize1 = value.toDouble();
2096 }
2097 // Pixel Size 2
2098 if (getRecordValue("PIXSIZE2", value))
2099 {
2100 pixsize2 = value.toDouble();
2101 }
2102 else if (getRecordValue("YPIXSZ", value))
2103 {
2104 pixsize2 = value.toDouble();
2105 }
2106
2107 int binx = 1, biny = 1;
2108 // Binning X
2109 if (getRecordValue("XBINNING", value))
2110 {
2111 binx = value.toDouble();
2112 }
2113 // Binning Y
2114 if (getRecordValue("YBINNING", value))
2115 {
2116 biny = value.toDouble();
2117 }
2118
2119 if (pixsize1 > 0 && pixsize2 > 0)
2120 {
2121 // If we have scale, then that's it
2122 if (scale > 0)
2123 {
2124 // Arcsecs per pixel
2125 solution.pixscale = scale;
2126 // Arcmins
2127 solution.fieldWidth = (m_Statistics.width * scale) / 60.0;
2128 // Arcmins, and account for pixel ratio if non-squared.
2129 solution.fieldHeight = (m_Statistics.height * scale * (pixsize2 / pixsize1)) / 60.0;
2130 }
2131 else if (focal_length > 0)
2132 {
2133 // Arcmins
2134 solution.fieldWidth = ((206264.8062470963552 * m_Statistics.width * (pixsize1 / 1000.0)) / (focal_length * binx)) / 60.0;
2135 // Arsmins
2136 solution.fieldHeight = ((206264.8062470963552 * m_Statistics.height * (pixsize2 / 1000.0)) / (focal_length * biny)) / 60.0;
2137 // Arcsecs per pixel
2138 solution.pixscale = (206264.8062470963552 * (pixsize1 / 1000.0)) / (focal_length * binx);
2139 }
2140
2141 scaleOK = true;
2142 }
2143
2144 return (coordOK || scaleOK);
2145}
2146
2147QFuture<bool> FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox)
2148{
2149 if (m_StarFindFuture.isRunning())
2150 m_StarFindFuture.waitForFinished();
2151
2152 starAlgorithm = algorithm;
2153 qDeleteAll(starCenters);
2154 starCenters.clear();
2155 starsSearched = true;
2156
2157 switch (algorithm)
2158 {
2159 case ALGORITHM_SEP:
2160 {
2161 m_StarDetector.reset(new FITSSEPDetector(this));
2162 m_StarDetector->setSettings(m_SourceExtractorSettings);
2163 if (m_Mode == FITS_NORMAL && trackingBox.isNull())
2164 {
2165 if (Options::quickHFR())
2166 {
2167 //Just finds stars in the center 25% of the image.
2168 const int w = getStatistics().width;
2169 const int h = getStatistics().height;
2170 QRect middle(static_cast<int>(w * 0.25), static_cast<int>(h * 0.25), w / 2, h / 2);
2171 m_StarFindFuture = m_StarDetector->findSources(middle);
2172 return m_StarFindFuture;
2173 }
2174 }
2175 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2176 return m_StarFindFuture;
2177 }
2178
2179 case ALGORITHM_GRADIENT:
2180 default:
2181 {
2182 m_StarDetector.reset(new FITSGradientDetector(this));
2183 m_StarDetector->setSettings(m_SourceExtractorSettings);
2184 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2185 return m_StarFindFuture;
2186 }
2187
2188 case ALGORITHM_CENTROID:
2189 {
2190#ifndef KSTARS_LITE
2191 m_StarDetector.reset(new FITSCentroidDetector(this));
2192 m_StarDetector->setSettings(m_SourceExtractorSettings);
2193 // We need JMIndex calculated from histogram
2194 if (!isHistogramConstructed())
2195 constructHistogram();
2196 m_StarDetector->configure("JMINDEX", m_JMIndex);
2197 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2198 return m_StarFindFuture;
2199 }
2200#else
2201 {
2202 m_StarDetector.reset(new FITSCentroidDetector(this));
2203 m_StarFindFuture = starDetector->findSources(trackingBox);
2204 return m_StarFindFuture;
2205 }
2206#endif
2207
2208 case ALGORITHM_THRESHOLD:
2209 {
2210 m_StarDetector.reset(new FITSThresholdDetector(this));
2211 m_StarDetector->setSettings(m_SourceExtractorSettings);
2212 m_StarDetector->configure("THRESHOLD_PERCENTAGE", Options::focusThreshold());
2213 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2214 return m_StarFindFuture;
2215 }
2216
2217 case ALGORITHM_BAHTINOV:
2218 {
2219 m_StarDetector.reset(new FITSBahtinovDetector(this));
2220 m_StarDetector->setSettings(m_SourceExtractorSettings);
2221 m_StarDetector->configure("NUMBER_OF_AVERAGE_ROWS", Options::focusMultiRowAverage());
2222 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2223 return m_StarFindFuture;
2224 }
2225 }
2226}
2227
2228int FITSData::filterStars(QSharedPointer<ImageMask> mask)
2229{
2230 if (mask.isNull() == false)
2231 {
2232 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(),
2233 [&](Edge * edge)
2234 {
2235 return (mask->isVisible(edge->x, edge->y) == false);
2236 }), starCenters.end());
2237 }
2238
2239 return starCenters.count();
2240
2241}
2242
2243double FITSData::getHFR(HFRType type)
2244{
2245 if (starCenters.empty())
2246 return -1;
2247
2248 if (cacheHFR >= 0 && cacheHFRType == type)
2249 return cacheHFR;
2250
2251 m_SelectedHFRStar.invalidate();
2252
2253 if (type == HFR_MAX)
2254 {
2255
2256 int maxVal = 0;
2257 int maxIndex = 0;
2258 for (int i = 0; i < starCenters.count(); i++)
2259 {
2260 if (starCenters[i]->HFR > maxVal)
2261 {
2262 maxIndex = i;
2263 maxVal = starCenters[i]->HFR;
2264 }
2265 }
2266
2267 m_SelectedHFRStar = *starCenters[maxIndex];
2268 cacheHFR = starCenters[maxIndex]->HFR;
2269 cacheHFRType = type;
2270 return cacheHFR;
2271 }
2272 else if (type == HFR_HIGH)
2273 {
2274 // Reject all stars within 10% of border
2275 int minX = width() / 10;
2276 int minY = height() / 10;
2277 int maxX = width() - minX;
2278 int maxY = height() - minY;
2279 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(), [minX, minY, maxX, maxY](Edge * oneStar)
2280 {
2281 return (oneStar->x < minX || oneStar->x > maxX || oneStar->y < minY || oneStar->y > maxY);
2282 }), starCenters.end());
2283 // Top 5%
2284 if (starCenters.empty())
2285 return -1;
2286
2287 m_SelectedHFRStar = *starCenters[static_cast<int>(starCenters.size() * 0.05)];
2288 cacheHFR = m_SelectedHFRStar.HFR;
2289 cacheHFRType = type;
2290 return cacheHFR;
2291 }
2292 else if (type == HFR_MEDIAN)
2293 {
2294 std::nth_element(starCenters.begin(), starCenters.begin() + starCenters.size() / 2, starCenters.end());
2295 m_SelectedHFRStar = *starCenters[starCenters.size() / 2];
2296
2297 cacheHFR = m_SelectedHFRStar.HFR;
2298 cacheHFRType = type;
2299 return cacheHFR;
2300 }
2301
2302 // We may remove saturated stars from the HFR calculation, if we have enough stars.
2303 // No real way to tell the scale, so only remove saturated stars with range 0 -> 2**16
2304 // for non-byte types. Unsigned types and floating types, or really any pixels whose
2305 // range isn't 0-64 (or 0-255 for bytes) won't have their saturated stars removed.
2306 const int saturationValue = m_Statistics.dataType == TBYTE ? 250 : 50000;
2307 int numSaturated = 0;
2308 for (auto center : starCenters)
2309 if (center->val > saturationValue)
2310 numSaturated++;
2311 const bool removeSaturatedStars = starCenters.size() - numSaturated > 20;
2312 if (removeSaturatedStars && numSaturated > 0)
2313 qCDebug(KSTARS_FITS) << "Removing " << numSaturated << " stars from HFR calculation";
2314
2315 std::vector<double> HFRs;
2316
2317 for (auto center : starCenters)
2318 {
2319 if (removeSaturatedStars && center->val > saturationValue) continue;
2320
2321 if (type == HFR_AVERAGE)
2322 HFRs.push_back(center->HFR);
2323 else
2324 {
2325 // HFR_ADJ_AVERAGE - so adjust the HFR based on the stars brightness
2326 // HFRadj = HFR.erf(sqrt(ln(peak/background)))/(1 - background/peak)
2327 // Sanity check inputs to prevent equation blowing up
2328 if (m_SkyBackground.mean <= 0.0 || center->val < m_SkyBackground.mean)
2329 {
2330 HFRs.push_back(center->HFR);
2331 qCDebug(KSTARS_FITS) << "HFR Adj, sky background " << m_SkyBackground.mean << " star peak " << center->val <<
2332 " not adjusting";
2333 }
2334 else
2335 {
2336 const double a_div_b = center->val / m_SkyBackground.mean;
2337 const double factor = erf(sqrt(log(a_div_b))) / (1 - (1 / a_div_b));
2338 HFRs.push_back(center->HFR * factor);
2339 qCDebug(KSTARS_FITS) << "HFR Adj, brightness adjusted from " << center->HFR << " to " << center->HFR * factor;
2340 }
2341 }
2342 }
2343
2344 auto m = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_SIGMACLIPPING,
2345 HFRs, 2);
2346
2347 cacheHFR = m;
2348 cacheHFRType = HFR_AVERAGE;
2349 return cacheHFR;
2350}
2351
2352double FITSData::getHFR(int x, int y, double scale)
2353{
2354 if (starCenters.empty())
2355 return -1;
2356
2357 for (int i = 0; i < starCenters.count(); i++)
2358 {
2359 const int maxDist = std::max(2, static_cast<int>(0.5 + 2 * starCenters[i]->width / scale));
2360 const int dx = std::fabs(starCenters[i]->x - x);
2361 const int dy = std::fabs(starCenters[i]->y - y);
2362 if (dx <= maxDist && dy <= maxDist)
2363 {
2364 m_SelectedHFRStar = *starCenters[i];
2365 return starCenters[i]->HFR;
2366 }
2367 }
2368
2369 return -1;
2370}
2371
2372double FITSData::getEccentricity()
2373{
2374 if (starCenters.empty())
2375 return -1;
2376 if (cacheEccentricity >= 0)
2377 return cacheEccentricity;
2378 std::vector<float> eccs;
2379 for (const auto &s : starCenters)
2380 eccs.push_back(s->ellipticity);
2381 int middle = eccs.size() / 2;
2382 std::nth_element(eccs.begin(), eccs.begin() + middle, eccs.end());
2383 float medianEllipticity = eccs[middle];
2384
2385 // SEP gives us the ellipticity (flattening) directly.
2386 // To get the eccentricity, use this formula:
2387 // e = sqrt(ellipticity * (2 - ellipticity))
2388 // https://en.wikipedia.org/wiki/Eccentricity_(mathematics)#Ellipses
2389 const float eccentricity = sqrt(medianEllipticity * (2 - medianEllipticity));
2390 cacheEccentricity = eccentricity;
2391 return eccentricity;
2392}
2393
2394void FITSData::applyFilter(FITSScale type, uint8_t * image, QVector<double> * min, QVector<double> * max)
2395{
2396 if (type == FITS_NONE)
2397 return;
2398
2399 QVector<double> dataMin(3);
2400 QVector<double> dataMax(3);
2401
2402 if (min)
2403 dataMin = *min;
2404 if (max)
2405 dataMax = *max;
2406
2407 switch (type)
2408 {
2409 case FITS_AUTO_STRETCH:
2410 {
2411 for (int i = 0; i < 3; i++)
2412 {
2413 dataMin[i] = m_Statistics.mean[i] - m_Statistics.stddev[i];
2414 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2415 }
2416 }
2417 break;
2418
2419 case FITS_HIGH_CONTRAST:
2420 {
2421 for (int i = 0; i < 3; i++)
2422 {
2423 dataMin[i] = m_Statistics.mean[i] + m_Statistics.stddev[i];
2424 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2425 }
2426 }
2427 break;
2428
2429 case FITS_HIGH_PASS:
2430 {
2431 for (int i = 0; i < 3; i++)
2432 {
2433 dataMin[i] = m_Statistics.mean[i];
2434 }
2435 }
2436 break;
2437
2438 default:
2439 break;
2440 }
2441
2442 switch (m_Statistics.dataType)
2443 {
2444 case TBYTE:
2445 {
2446 for (int i = 0; i < 3; i++)
2447 {
2448 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2449 dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i];
2450 }
2451 applyFilter<uint8_t>(type, image, &dataMin, &dataMax);
2452 }
2453 break;
2454
2455 case TSHORT:
2456 {
2457 for (int i = 0; i < 3; i++)
2458 {
2459 dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i];
2460 dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i];
2461 }
2462 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2463 }
2464
2465 break;
2466
2467 case TUSHORT:
2468 {
2469 for (int i = 0; i < 3; i++)
2470 {
2471 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2472 dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i];
2473 }
2474 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2475 }
2476 break;
2477
2478 case TLONG:
2479 {
2480 for (int i = 0; i < 3; i++)
2481 {
2482 dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i];
2483 dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i];
2484 }
2485 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2486 }
2487 break;
2488
2489 case TULONG:
2490 {
2491 for (int i = 0; i < 3; i++)
2492 {
2493 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2494 dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i];
2495 }
2496 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2497 }
2498 break;
2499
2500 case TFLOAT:
2501 {
2502 for (int i = 0; i < 3; i++)
2503 {
2504 dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i];
2505 dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i];
2506 }
2507 applyFilter<float>(type, image, &dataMin, &dataMax);
2508 }
2509 break;
2510
2511 case TLONGLONG:
2512 {
2513 for (int i = 0; i < 3; i++)
2514 {
2515 dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i];
2516 dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i];
2517 }
2518
2519 applyFilter<long>(type, image, &dataMin, &dataMax);
2520 }
2521 break;
2522
2523 case TDOUBLE:
2524 {
2525 for (int i = 0; i < 3; i++)
2526 {
2527 dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i];
2528 dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i];
2529 }
2530 applyFilter<double>(type, image, &dataMin, &dataMax);
2531 }
2532
2533 break;
2534
2535 default:
2536 return;
2537 }
2538
2539 if (min != nullptr)
2540 *min = dataMin;
2541 if (max != nullptr)
2542 *max = dataMax;
2543
2544 emit dataChanged();
2545}
2546
2547template <typename T>
2548void FITSData::applyFilter(FITSScale type, uint8_t * targetImage, QVector<double> * targetMin, QVector<double> * targetMax)
2549{
2550 bool calcStats = false;
2551 T * image = nullptr;
2552
2553 if (targetImage)
2554 image = reinterpret_cast<T *>(targetImage);
2555 else
2556 {
2557 image = reinterpret_cast<T *>(m_ImageBuffer);
2558 calcStats = true;
2559 }
2560
2561 T min[3], max[3];
2562 for (int i = 0; i < 3; i++)
2563 {
2564 min[i] = (*targetMin)[i] < std::numeric_limits<T>::min() ? std::numeric_limits<T>::min() : (*targetMin)[i];
2565 max[i] = (*targetMax)[i] > std::numeric_limits<T>::max() ? std::numeric_limits<T>::max() : (*targetMax)[i];
2566 }
2567
2568
2569 // Create N threads
2570 const uint8_t nThreads = 16;
2571
2572 uint32_t width = m_Statistics.width;
2573 uint32_t height = m_Statistics.height;
2574
2575 //QTime timer;
2576 //timer.start();
2577 switch (type)
2578 {
2579 case FITS_AUTO:
2580 case FITS_LINEAR:
2581 case FITS_AUTO_STRETCH:
2582 case FITS_HIGH_CONTRAST:
2583 case FITS_LOG:
2584 case FITS_SQRT:
2585 case FITS_HIGH_PASS:
2586 {
2587 // List of futures
2588 QList<QFuture<void>> futures;
2589 QVector<double> coeff(3);
2590
2591 if (type == FITS_LOG)
2592 {
2593 for (int i = 0; i < 3; i++)
2594 coeff[i] = max[i] / std::log(1 + max[i]);
2595 }
2596 else if (type == FITS_SQRT)
2597 {
2598 for (int i = 0; i < 3; i++)
2599 coeff[i] = max[i] / sqrt(max[i]);
2600 }
2601
2602 for (int n = 0; n < m_Statistics.channels; n++)
2603 {
2604 if (type == FITS_HIGH_PASS)
2605 min[n] = m_Statistics.mean[n];
2606
2607 uint32_t cStart = n * m_Statistics.samples_per_channel;
2608
2609 // Calculate how many elements we process per thread
2610 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
2611
2612 // Calculate the final stride since we can have some left over due to division above
2613 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
2614
2615 T * runningBuffer = image + cStart;
2616
2617 if (type == FITS_LOG)
2618 {
2619 for (int i = 0; i < nThreads; i++)
2620 {
2621 // Run threads
2622 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2623 coeff, n](T & a)
2624 {
2625 a = qBound(min[n], static_cast<T>(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]);
2626 }));
2627
2628 runningBuffer += tStride;
2629 }
2630 }
2631 else if (type == FITS_SQRT)
2632 {
2633 for (int i = 0; i < nThreads; i++)
2634 {
2635 // Run threads
2636 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2637 coeff, n](T & a)
2638 {
2639 a = qBound(min[n], static_cast<T>(round(coeff[n] * a)), max[n]);
2640 }));
2641 }
2642
2643 runningBuffer += tStride;
2644 }
2645 else
2646 {
2647 for (int i = 0; i < nThreads; i++)
2648 {
2649 // Run threads
2650 futures.append(QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2651 n](T & a)
2652 {
2653 a = qBound(min[n], a, max[n]);
2654 }));
2655 runningBuffer += tStride;
2656 }
2657 }
2658 }
2659
2660 for (int i = 0; i < nThreads * m_Statistics.channels; i++)
2661 futures[i].waitForFinished();
2662
2663 if (calcStats)
2664 {
2665 for (int i = 0; i < 3; i++)
2666 {
2667 m_Statistics.min[i] = min[i];
2668 m_Statistics.max[i] = max[i];
2669 }
2670 calculateStdDev<T>();
2671 }
2672 }
2673 break;
2674
2675 case FITS_EQUALIZE:
2676 {
2677#ifndef KSTARS_LITE
2678 if (!isHistogramConstructed())
2679 constructHistogram();
2680
2681 T bufferVal = 0;
2682
2683 double coeff = 255.0 / (height * width);
2684 uint32_t row = 0;
2685 uint32_t index = 0;
2686
2687 for (int i = 0; i < m_Statistics.channels; i++)
2688 {
2689 uint32_t offset = i * m_Statistics.samples_per_channel;
2690 for (uint32_t j = 0; j < height; j++)
2691 {
2692 row = offset + j * width;
2693 for (uint32_t k = 0; k < width; k++)
2694 {
2695 index = k + row;
2696 bufferVal = (image[index] - min[i]) / m_HistogramBinWidth[i];
2697
2698 if (bufferVal >= m_CumulativeFrequency[i].size())
2699 bufferVal = m_CumulativeFrequency[i].size() - 1;
2700
2701 image[index] = qBound(min[i], static_cast<T>(round(coeff * m_CumulativeFrequency[i][bufferVal])), max[i]);
2702 }
2703 }
2704 }
2705#endif
2706 }
2707 if (calcStats)
2708 calculateStats(true, false);
2709 break;
2710
2711 // Based on http://www.librow.com/articles/article-1
2712 case FITS_MEDIAN:
2713 {
2714 uint8_t BBP = m_Statistics.bytesPerPixel;
2715 auto * extension = new T[(width + 2) * (height + 2)];
2716 // Check memory allocation
2717 if (!extension)
2718 return;
2719 // Create image extension
2720 for (uint32_t ch = 0; ch < m_Statistics.channels; ch++)
2721 {
2722 uint32_t offset = ch * m_Statistics.samples_per_channel;
2723 uint32_t N = width, M = height;
2724
2725 for (uint32_t i = 0; i < M; ++i)
2726 {
2727 memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP);
2728 extension[(N + 2) * (i + 1)] = image[N * i + offset];
2729 extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset];
2730 }
2731 // Fill first line of image extension
2732 memcpy(extension, extension + N + 2, (N + 2) * BBP);
2733 // Fill last line of image extension
2734 memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP);
2735 // Call median filter implementation
2736
2737 N = width + 2;
2738 M = height + 2;
2739
2740 // Move window through all elements of the image
2741 for (uint32_t m = 1; m < M - 1; ++m)
2742 for (uint32_t n = 1; n < N - 1; ++n)
2743 {
2744 // Pick up window elements
2745 int k = 0;
2746 float window[9];
2747
2748 memset(&window[0], 0, 9 * sizeof(float));
2749 for (uint32_t j = m - 1; j < m + 2; ++j)
2750 for (uint32_t i = n - 1; i < n + 2; ++i)
2751 window[k++] = extension[j * N + i];
2752 // Order elements (only half of them)
2753 for (uint32_t j = 0; j < 5; ++j)
2754 {
2755 // Find position of minimum element
2756 int mine = j;
2757 for (uint32_t l = j + 1; l < 9; ++l)
2758 if (window[l] < window[mine])
2759 mine = l;
2760 // Put found minimum element in its place
2761 const float temp = window[j];
2762 window[j] = window[mine];
2763 window[mine] = temp;
2764 }
2765 // Get result - the middle element
2766 image[(m - 1) * (N - 2) + n - 1 + offset] = window[4];
2767 }
2768 }
2769
2770 // Free memory
2771 delete[] extension;
2772
2773 if (calcStats)
2774 calculateStdDev<T>();
2775 }
2776 break;
2777
2778 case FITS_GAUSSIAN:
2779 gaussianBlur<T>(Options::focusGaussianKernelSize(), Options::focusGaussianSigma());
2780 if (calcStats)
2781 calculateStats(true, false);
2782 break;
2783
2784 case FITS_ROTATE_CW:
2785 rotFITS<T>(90, 0);
2786 rotCounter++;
2787 break;
2788
2789 case FITS_ROTATE_CCW:
2790 rotFITS<T>(270, 0);
2791 rotCounter--;
2792 break;
2793
2794 case FITS_MOUNT_FLIP_H:
2795 rotFITS<T>(0, 1);
2796 flipHCounter++;
2797 break;
2798
2799 case FITS_MOUNT_FLIP_V:
2800 rotFITS<T>(0, 2);
2801 flipVCounter++;
2802 break;
2803
2804 default:
2805 break;
2806 }
2807}
2808
2809QList<Edge *> FITSData::getStarCentersInSubFrame(QRect subFrame) const
2810{
2811 QList<Edge *> starCentersInSubFrame;
2812 for (int i = 0; i < starCenters.count(); i++)
2813 {
2814 int x = static_cast<int>(starCenters[i]->x);
2815 int y = static_cast<int>(starCenters[i]->y);
2816 if(subFrame.contains(x, y))
2817 {
2818 starCentersInSubFrame.append(starCenters[i]);
2819 }
2820 }
2821 return starCentersInSubFrame;
2822}
2823
2824bool FITSData::loadWCS()
2825{
2826#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2827
2828 if (m_WCSState == Success)
2829 {
2830 qCWarning(KSTARS_FITS) << "WCS data already loaded";
2831 return true;
2832 }
2833
2834 if (m_WCSHandle != nullptr)
2835 {
2836 wcsvfree(&m_nwcs, &m_WCSHandle);
2837 m_nwcs = 0;
2838 m_WCSHandle = nullptr;
2839 }
2840
2841 qCDebug(KSTARS_FITS) << "Started WCS Data Processing...";
2842
2843 QByteArray header_str;
2844 int status = 0;
2845 int nkeyrec = 0, nreject = 0;
2846 if (fptr)
2847 {
2848 char *header = nullptr;
2849 if (fits_hdr2str(fptr, 1, nullptr, 0, &header, &nkeyrec, &status))
2850 {
2851 char errmsg[512];
2852 fits_get_errstatus(status, errmsg);
2853 m_LastError = errmsg;
2854 m_WCSState = Failure;
2855 return false;
2856 }
2857 header_str = QByteArray(header);
2858 fits_free_memory(header, &status);
2859 }
2860 else
2861 {
2862 nkeyrec = 1;
2863 for(auto &fitsKeyword : m_HeaderRecords)
2864 {
2865 QByteArray rec;
2866 rec.append(fitsKeyword.key.leftJustified(8, ' ').toLatin1());
2867 rec.append("= ");
2868 rec.append(fitsKeyword.value.toByteArray());
2869 rec.append(" / ");
2870 rec.append(fitsKeyword.comment.toLatin1());
2871 header_str.append(rec.leftJustified(80, ' ', true));
2872 nkeyrec++;
2873 }
2874 header_str.append(QByteArray("END").leftJustified(80));
2875 }
2876
2877 if ((status = wcspih(header_str.data(), nkeyrec, WCSHDR_all, 0, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2878 {
2879 wcsvfree(&m_nwcs, &m_WCSHandle);
2880 m_WCSHandle = nullptr;
2881 m_nwcs = 0;
2882 m_LastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]);
2883 m_WCSState = Failure;
2884 return false;
2885 }
2886
2887 if (m_WCSHandle == nullptr)
2888 {
2889 m_LastError = i18n("No world coordinate systems found.");
2890 m_WCSState = Failure;
2891 return false;
2892 }
2893
2894 // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now.
2895 if (m_WCSHandle->crpix[0] == 0)
2896 {
2897 wcsvfree(&m_nwcs, &m_WCSHandle);
2898 m_WCSHandle = nullptr;
2899 m_nwcs = 0;
2900 m_LastError = i18n("No world coordinate systems found.");
2901 m_WCSState = Failure;
2902 return false;
2903 }
2904
2905 cdfix(m_WCSHandle);
2906 if ((status = wcsset(m_WCSHandle)) != 0)
2907 {
2908 wcsvfree(&m_nwcs, &m_WCSHandle);
2909 m_WCSHandle = nullptr;
2910 m_nwcs = 0;
2911 m_LastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2912 m_WCSState = Failure;
2913 return false;
2914 }
2915
2916 m_ObjectsSearched = false;
2917 m_CatObjectsSearched = false;
2918 m_WCSState = Success;
2919 HasWCS = true;
2920
2921 qCDebug(KSTARS_FITS) << "Finished WCS Data processing...";
2922
2923 return true;
2924#else
2925 return false;
2926#endif
2927}
2928
2929bool FITSData::wcsToPixel(const SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint)
2930{
2931#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2932 int status, stat[NWCSFIX];
2933 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2934
2935 if (m_WCSHandle == nullptr)
2936 {
2937 m_LastError = i18n("No world coordinate systems found.");
2938 return false;
2939 }
2940
2941 worldcrd[0] = wcsCoord.ra0().Degrees();
2942 worldcrd[1] = wcsCoord.dec0().Degrees();
2943
2944 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2945 {
2946 m_LastError = QString("wcss2p error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2947 return false;
2948 }
2949
2950 wcsImagePoint.setX(imgcrd[0]);
2951 wcsImagePoint.setY(imgcrd[1]);
2952
2953 wcsPixelPoint.setX(pixcrd[0]);
2954 wcsPixelPoint.setY(pixcrd[1]);
2955
2956 return true;
2957#else
2958 Q_UNUSED(wcsCoord);
2959 Q_UNUSED(wcsPixelPoint);
2960 Q_UNUSED(wcsImagePoint);
2961 return false;
2962#endif
2963}
2964
2965bool FITSData::pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord)
2966{
2967#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2968 int status, stat[NWCSFIX];
2969 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2970
2971 if (m_WCSHandle == nullptr)
2972 {
2973 m_LastError = i18n("No world coordinate systems found.");
2974 return false;
2975 }
2976
2977 pixcrd[0] = wcsPixelPoint.x();
2978 pixcrd[1] = wcsPixelPoint.y();
2979
2980 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2981 {
2982 m_LastError = QString("wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2983 return false;
2984 }
2985 else
2986 {
2987 wcsCoord.setRA0(world[0] / 15.0);
2988 wcsCoord.setDec0(world[1]);
2989 }
2990
2991 return true;
2992#else
2993 Q_UNUSED(wcsPixelPoint);
2994 Q_UNUSED(wcsCoord);
2995 return false;
2996#endif
2997}
2998
2999#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3000bool FITSData::searchObjects()
3001{
3002 return (Options::fitsCatalog() == CAT_SKYMAP) ? searchSkyMapObjects() : searchCatObjects();
3003}
3004
3005bool FITSData::searchSkyMapObjects()
3006{
3007 if (m_ObjectsSearched)
3008 return true;
3009
3010 m_ObjectsSearched = true;
3011
3012 SkyPoint startPoint;
3013 SkyPoint endPoint;
3014
3015 pixelToWCS(QPointF(0, 0), startPoint);
3016 pixelToWCS(QPointF(width() - 1, height() - 1), endPoint);
3017
3018 return findObjectsInImage(startPoint, endPoint);
3019}
3020
3021bool FITSData::searchCatObjects()
3022{
3023 if (m_CatObjectsSearched)
3024 return true;
3025
3026 m_CatObjectsSearched = true;
3027
3028 SkyPoint searchCenter;
3029 double radius;
3030 QPoint pt;
3031 bool ok = true;
3032 if (catROIRadius() > 0)
3033 {
3034 // A ROI has been set so this is a request where the user set a ROI
3035 pt = catROIPt();
3036 QPoint edgePt = QPoint(pt.x() + catROIRadius(), pt.y());
3037 SkyPoint searchEdge;
3038 ok = pixelToWCS(pt, searchCenter);
3039 if (ok)
3040 ok = pixelToWCS(edgePt, searchEdge);
3041 if (ok)
3042 {
3043 QVariant date;
3044 KSNumbers * num = nullptr;
3045
3046 if (getRecordValue("DATE-OBS", date))
3047 {
3048 QString tsString(date.toString());
3049 tsString = tsString.remove('\'').trimmed();
3050 // Add Zulu time to indicate UTC
3051 tsString += "Z";
3052
3054
3055 if (ts.isValid())
3056 num = new KSNumbers(KStarsDateTime(ts).djd());
3057 }
3058
3059 //Set to current time if the above does not work.
3060 if (num == nullptr)
3061 num = new KSNumbers(KStarsData::Instance()->ut().djd());
3062
3063 searchCenter.updateCoordsNow(num);
3064 searchEdge.updateCoordsNow(num);
3065 radius = searchCenter.angularDistanceTo(&searchEdge).Degrees() * 60; // Arc minutes
3066
3067 delete num;
3068 }
3069 }
3070 else
3071 {
3072 // No ROI setup so this is a first call to use a circle of 0.5 arcmins in the center of the image
3073 // Lets calculate the number of pixels that correspond to the 0.5 arcmin radius
3074 pt = QPoint((width() / 2) - 1, (height() / 2) - 1);
3075 ok = pixelToWCS(pt, searchCenter);
3076 if (ok)
3077 {
3078 radius = 1.0; // Use 1.0 arcmins as a starting point... not too big to avoid getting swamped with objects
3079 double raEdge = searchCenter.ra0().Degrees() + (radius / 60.0);
3080 double decEdge = searchCenter.dec0().Degrees();
3081 SkyPoint searchEdge(raEdge / 15.0, decEdge);
3082 QPointF edgePoint, pEdge;
3083 ok = wcsToPixel(searchEdge, pEdge, edgePoint);
3084 if (ok)
3085 {
3086 // Set the ROI which will be drawn
3087 const double radiusPix = std::hypot((pt.x() - pEdge.x()), pt.y() - pEdge.y());
3088 setCatSearchROI(pt, radiusPix);
3089 }
3090 }
3091 }
3092 if (!ok)
3093 {
3094 qCDebug(KSTARS_FITS) << "Unable to process Catalog Object request...";
3095 return false;
3096 }
3097 if (Options::fitsCatalog() == CAT_SIMBAD)
3098 return findSimbadObjectsInImage(searchCenter, radius);
3099 else
3100 return false;
3101}
3102
3103bool FITSData::findWCSBounds(double &minRA, double &maxRA, double &minDec, double &maxDec)
3104{
3105 if (m_WCSHandle == nullptr)
3106 {
3107 m_LastError = i18n("No world coordinate systems found.");
3108 return false;
3109 }
3110
3111 maxRA = -1000;
3112 minRA = 1000;
3113 maxDec = -1000;
3114 minDec = 1000;
3115
3116 auto updateMinMax = [&](int x, int y)
3117 {
3118 int stat[2];
3119 double imgcrd[2], phi, pixcrd[2], theta, world[2];
3120
3121 pixcrd[0] = x;
3122 pixcrd[1] = y;
3123
3124 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
3125 return;
3126
3127 minRA = std::min(minRA, world[0]);
3128 maxRA = std::max(maxRA, world[0]);
3129 minDec = std::min(minDec, world[1]);
3130 maxDec = std::max(maxDec, world[1]);
3131 };
3132
3133 // Find min and max values from edges
3134 for (int y = 0; y < height(); y++)
3135 {
3136 updateMinMax(0, y);
3137 updateMinMax(width() - 1, y);
3138 }
3139
3140 for (int x = 1; x < width() - 1; x++)
3141 {
3142 updateMinMax(x, 0);
3143 updateMinMax(x, height() - 1);
3144 }
3145
3146 // Check if either pole is in the image
3147 SkyPoint NCP(0, 90);
3148 SkyPoint SCP(0, -90);
3149 QPointF imagePoint, pPoint;
3150 if (wcsToPixel(NCP, pPoint, imagePoint))
3151 {
3152 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
3153 {
3154 maxDec = 90;
3155 }
3156 }
3157 if (wcsToPixel(SCP, pPoint, imagePoint))
3158 {
3159 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
3160 {
3161 minDec = -90;
3162 }
3163 }
3164 return true;
3165}
3166#endif
3167
3168#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3169bool FITSData::findObjectsInImage(SkyPoint startPoint, SkyPoint endPoint)
3170{
3171 if (KStarsData::Instance() == nullptr)
3172 return false;
3173
3174 int w = width();
3175 int h = height();
3176 QVariant date;
3177 KSNumbers * num = nullptr;
3178
3179 if (getRecordValue("DATE-OBS", date))
3180 {
3181 QString tsString(date.toString());
3182 tsString = tsString.remove('\'').trimmed();
3183 // Add Zulu time to indicate UTC
3184 tsString += "Z";
3185
3187
3188 if (ts.isValid())
3189 num = new KSNumbers(KStarsDateTime(ts).djd());
3190 }
3191
3192 //Set to current time if the above does not work.
3193 if (num == nullptr)
3194 num = new KSNumbers(KStarsData::Instance()->ut().djd());
3195
3196 startPoint.updateCoordsNow(num);
3197 endPoint.updateCoordsNow(num);
3198
3199 m_SkyObjects.clear();
3200
3201 QList<SkyObject *> list = KStarsData::Instance()->skyComposite()->findObjectsInArea(startPoint, endPoint);
3202 list.erase(std::remove_if(list.begin(), list.end(), [](SkyObject * oneObject)
3203 {
3204 int type = oneObject->type();
3205 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
3206 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
3207 type == SkyObject::SATELLITE);
3208 }), list.end());
3209
3210 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3211 int stat[2];
3212 for (auto &object : list)
3213 {
3214 world[0] = object->ra0().Degrees();
3215 world[1] = object->dec0().Degrees();
3216
3217 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
3218 {
3219 //The X and Y are set to the found position if it does work.
3220 int x = pixcrd[0];
3221 int y = pixcrd[1];
3222 if (x > 0 && y > 0 && x < w && y < h)
3223 m_SkyObjects.append(new FITSSkyObject(object, x, y));
3224 }
3225 }
3226
3227 delete (num);
3228 return true;
3229}
3230
3231// See https://simbad.u-strasbg.fr/Pages/guide/sim-q.htx and
3232// https://simbad.u-strasbg.fr/Pages/guide/sim-url.htx
3233// for details of how to query Simbad and the various options
3234//
3235// Set the epoch from the datetime of the image and ensure that the reply is in this same epoch
3236bool FITSData::findSimbadObjectsInImage(SkyPoint searchCenter, double radius)
3237{
3238 m_CatObjects.clear();
3239 m_CatUpdateTable = false;
3240
3241 // Query Simbad
3242 QUrl simbadURL = QUrl("https://simbad.cds.unistra.fr/simbad/sim-coo");
3243 QUrlQuery simbadQuery(simbadURL.query());
3244
3245 QString coord = QString("%1 %2").arg(searchCenter.ra0().toHMSString(true, true))
3246 .arg(searchCenter.dec0().toDMSString(true, true, true));
3247
3248 QString radiusStr = QString("%1").arg(radius * 60.0, 0, 'f', 5);
3249 m_CatObjQuery = QString("%1 %2").arg(coord).arg(radiusStr);
3250
3251 coord.replace("+", "%2B"); // Need to replace the + otherwise it gets lost by URL processing
3252
3253 QString epoch = QString("J%1").arg(m_DateTime.epoch());
3254
3255 simbadQuery.addQueryItem("Coord", coord);
3256 simbadQuery.addQueryItem("Radius", radiusStr);
3257 simbadQuery.addQueryItem("Radius.unit", "arcsec");
3258 simbadQuery.addQueryItem("CooFrame", "ICRS");
3259 simbadQuery.addQueryItem("CooEpoch", epoch);
3260 simbadQuery.addQueryItem("output.format", "ASCII");
3261 simbadQuery.addQueryItem("output.max", QString("%1").arg(10000));
3262 simbadQuery.addQueryItem("list.otypesel", "on");
3263 simbadQuery.addQueryItem("otypedisp", "3"); // Use Simbad's 3 char object type
3264 simbadQuery.addQueryItem("list.coo1", "on"); // Display coord 1
3265 simbadQuery.addQueryItem("frame1", "ICRS");
3266 simbadQuery.addQueryItem("epoch1", "J2000");
3267 simbadQuery.addQueryItem("equi1", epoch);
3268 simbadQuery.addQueryItem("list.spsel", "off"); // Spectral Type
3269 simbadQuery.addQueryItem("list.sizesel", "on"); // Angular size
3270 simbadQuery.addQueryItem("list.bibsel", "off"); // Bibliography
3271
3272 simbadURL.setQuery(simbadQuery);
3273 qCDebug(KSTARS_FITS) << "Simbad query:" << simbadURL;
3274
3275 m_NetworkAccessManager.reset(new QNetworkAccessManager(this));
3276 connect(m_NetworkAccessManager.get(), &QNetworkAccessManager::finished, this, &FITSData::simbadResponseReady);
3277 QNetworkReply *response = m_NetworkAccessManager->get(QNetworkRequest(simbadURL));
3278 if (response->error() != QNetworkReply::NoError)
3279 {
3280 qCDebug(KSTARS_FITS) << "Error:" << response->errorString() << " occured querying SIMBAD with" << simbadURL;
3281 m_CatQueryInProgress = false;
3282 emit catalogQueryFailed(i18n("Error querying Simbad"));
3283 return false;
3284 }
3285
3286 m_CatQueryInProgress = true;
3287 m_CatQueryTimer.setSingleShot(true);
3288 connect(&m_CatQueryTimer, &QTimer::timeout, this, &FITSData::catTimeout);
3289 m_CatQueryTimer.start(60 * 1000);
3290 emit loadingCatalogData();
3291 return true;
3292}
3293#endif
3294
3295#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3296void FITSData::catTimeout()
3297{
3298 m_CatQueryInProgress = false;
3299 QString text = i18n("Simbad query timed out");
3300 qCDebug(KSTARS_FITS) << text;
3301 emit catalogQueryFailed(text);
3302}
3303#endif
3304
3305// There are 3 types of responses... zero objects, 1 object and >1 objects. The formats are as follows:
3306// (Note that Line xx isn't part of the reply, these are added for convenience)
3307//
3308// Zero Objects
3309// Line 0 "!! No astronomical object found : "
3310//
3311// 1 Object
3312// Line 0 "C.D.S. - SIMBAD4 rel 1.8 - 2024.09.30CEST14:07:24"
3313// Line 1 ""
3314// Line 2 "coord 21:36:49 +57:31:23 (ICRS, J2024.7, 2000.0), radius: 5.71083 arcsec"
3315// Line 3 "------------------------------------------------------------------------"
3316// Line 4 ""
3317// Line 5 "Object EM* LkHA 349C --- RS* --- OID=@157406 (@@122902,11) --- coobox=594"
3318// Line 6 ""
3319// Line 7 "Coordinates(ICRS,ep=J2000,eq=J2024.7): 21 36 49.4008740691 +57 31 21.834541600 (Opt ) A [0.0658 0.0673 90] 2020yCat.1350....0G"
3320// Line 8 "Coordinates(FK4,ep=B1950,eq=1950): 21 35 17.0484337158 +57 17 51.458438447"
3321// Line 9 "Coordinates(Gal,ep=J2000,eq=2000): 099.0986092951232 +03.9548388170815"
3322// Line 10 "hierarchy counts: #parents=2, #children=0, #siblings=0"
3323// Line 11 "Proper motions: -2.232 -5.540 [0.087 0.086 90] A 2020yCat.1350....0G"
3324// Line 12 "Parallax: 1.0825 [0.0718] A 2020yCat.1350....0G"
3325// Line 13 "Radial Velocity: ~ [~ ~ ] ~ ~ "
3326// Line 14 "Redshift: ~ [~ ~ ] ~ ~ "
3327// Line 15 "cz: ~ [~ ~ ] ~ ~ "
3328// Line 16 "Flux U : 19.219 [0.083] C 2010ApJ...710..597S"
3329// Line 17 "Flux V : 15.81 [0.01] B 2012MNRAS.426.2917G"
3330// Line 18 "Flux G : 14.948605 [0.007309] C 2020yCat.1350....0G"
3331// Line 19 "Flux R : 14.508 [0.030] C 2010ApJ...710..597S"
3332// Line 20 "Flux I : 13.094 [0.003] C 2010ApJ...710..597S"
3333// Line 21 "Flux J : 11.918 [0.022] C 2003yCat.2246....0C"
3334// Line 22 "Flux H : 10.914 [0.029] C 2003yCat.2246....0C"
3335// Line 23 "Flux K : 10.355 [0.021] C 2003yCat.2246....0C"
3336// Line 24 "Flux g : 16.702 [~] D 2020ApJS..249...18C"
3337// Line 25 "Flux r : 15.057 [~] D 2020ApJS..249...18C"
3338// Line 26 "Flux i : 14.84 [~] D 2012AJ....143...61N"
3339// Line 27 "Spectral type: K6 D 2006AJ....132.2135S"
3340// Line 28 "Morphological type: ~ ~ ~"
3341// Line 29 "Angular size: ~ ~ ~ (~) ~ ~"
3342// Line 30 ""
3343// Line 31 "Identifiers (12):"
3344// Line 32 " ZTF J213649.40+573121.9 Gaia DR3 2178442322028512512 TIC 469533750 "
3345// Line 33 " CoKu LkHA 349 c EM* LkHA 349C 2MASS J21364941+5731220 "
3346// Line 34 " [SHB2004] Trumpler 37 14-141 [MSR2009] IC1396A-29 [NSW2012] 125 "
3347// Line 35 " CXOIC1396A J213649.39+573122.1 [GFS2012] 173 Gaia DR2 2178442322028512512 "
3348// Line 36 ""
3349// Line 37 "Bibcodes 1850-2024 () (20):"
3350// Line 38 " 2024MNRAS.532.2108S 2020ApJS..249...18C 2019ApJ...878....7M 2018MNRAS.478.5091F"
3351// Line 39 " 2012AJ....143...61N 2012MNRAS.426.2917G 2011ApJ...742...39S 2010ApJ...710..597S"
3352// Line 40 " 2009ApJ...690..683R 2009ApJ...702.1507M 2006AJ....132.2135S 2006ApJ...638..897S"
3353// Line 41 " 2005A&A...443..535V 2005AJ....130..188S 2004AJ....128..805S 2004ApJS..154..385R"
3354// Line 42 " 1997A&A...325.1001S 1996ApJ...463L.105M 1995A&A...299..464H 1979ApJS...41..743C"
3355// Line 43 ""
3356// Line 44 "Measures (distance:1 PLX:2 PM:2 ROT:1 V*:3 velocities:1 ):"
3357// Line 45 "distance:1PLX:2PM:2ROT:1V*:3velocities:1"
3358// Line 46 ""
3359// Line 47 "Notes (0) :"
3360// Line 48 ""
3361// Line 49 "================================================================================"
3362// Line 50 ""
3363// Line 51 ""
3364//
3365// >1 Objects
3366// Line 0 "C.D.S. - SIMBAD4 rel 1.8 - 2024.09.30CEST14:13:23"
3367// Line 1 ""
3368// Line 2 "coord 21:36:51 +57:31:13 (ICRS, J2024.7, 2000.0), radius: 11.42166 arcsec"
3369// Line 3 "-------------------------------------------------------------------------"
3370// Line 4 ""
3371// Line 5 "Number of objects : 2"
3372// Line 6 ""
3373// Line 7 "#|dist(asec)| identifier |typ| coord1 (ICRS,J2000/J2024.7) |Mag U | Mag B | Mag V | Mag R | Mag I | ang. size |#not"
3374// Line 8 "-|----------|-----------------------------------|---|---------------------------------------|------|---------|---------|---------|---------|---------------|----"
3375// Line 9 "1| 2.62|DOBASHI 3187 |DNe|21 36 51.3 +57 31 12 | ~| ~ | ~ | ~ | ~ | ~ ~ | 1"
3376// Line 10 "2| 3.19|EM* LkHA 349 |Or*|21 36 50.7149106104 +57 31 10.580260473| ~|15.13 |13.45 |12.88 |10.91 | ~ ~ | 1"
3377// Line 11 "================================================================================"
3378// Line 12 ""
3379// Line 13 ""
3380
3381#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3382void FITSData::simbadResponseReady(QNetworkReply * simbadReply)
3383{
3384 m_CatQueryTimer.stop();
3385 m_CatQueryInProgress = false;
3386 if (simbadReply->error() != QNetworkReply::NoError)
3387 {
3388 qCDebug(KSTARS_FITS) << "Error:" << simbadReply->errorString() << " occured in SIMBAD query reply";
3389 m_CatQueryInProgress = false;
3390 emit catalogQueryFailed(i18n("Error querying Simbad"));
3391 return;
3392 }
3393
3394 auto data = simbadReply->readAll();
3395 QString dataStr = QString(data);
3396
3397 // Break the line into comma-separated components
3398 int numObjs = 0;
3399 bool ok;
3400 QString name, type, coord, sizeStr;
3401 double dist = 0, magnitude = 0;
3402 int num = 0;
3403 QStringList replyLines = dataStr.split(QLatin1Char('\n'));
3404 int dataRow = -1;
3405 for (int i = 0; i < replyLines.size(); i++)
3406 {
3407 // JEE qCDebug(KSTARS_FITS) << "Line" << i << replyLines[i];
3408 if (i == 0 && replyLines[i].contains("No astronomical object found : "))
3409 // No objects in the reply so we can just ignore the rest of the data
3410 break;
3411 if (i == 2)
3412 {
3413 // Check that the query coord & radius match the reply - otherwise ignore...
3414 QString replyStr;
3415 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3416 if (replyData.size() == 9 && replyData[0] == "coord" && replyData[6] == "radius:")
3417 replyStr = QString("%1 %2 %3").arg(replyData[1]).arg(replyData[2]).arg(replyData[7]);
3418 if (replyStr != m_CatObjQuery)
3419 {
3420 qCDebug(KSTARS_FITS) << "Simbad query:" << m_CatObjQuery << "Reply:" << replyStr << ". Ignoring...";
3421 break;
3422 }
3423 }
3424 else if (i == 5)
3425 {
3426 if (replyLines[i].contains("Number of objects : "))
3427 numObjs = QString(replyLines[i].mid(20)).toInt(&ok);
3428 else if (replyLines[i].startsWith("Object "))
3429 {
3430 // Dealing with a single object
3431 num = numObjs = 1;
3432 const int firstBreak = replyLines[i].indexOf("---");
3433 const int secondBreak = replyLines[i].indexOf("---", firstBreak + 3);
3434 name = replyLines[i].mid(7, firstBreak - 7).trimmed();
3435 type = replyLines[i].mid(firstBreak + 3, secondBreak - firstBreak - 3).trimmed();
3436 }
3437 else
3438 {
3439 qCDebug(KSTARS_FITS) << "Bad Simbad Reply, Ignoring...";
3440 break;
3441 }
3442 }
3443 else if (numObjs == 1 && i >= 7)
3444 {
3445 if (replyLines[i].startsWith("Coordinates(ICRS"))
3446 {
3447 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3448 if (replyData.size() >= 8)
3449 coord = replyData[1] + " " + replyData[2] + " " + replyData[3] + " " +
3450 replyData[5] + " " + replyData[6] + " " + replyData[7];
3451 }
3452 else if (replyLines[i].startsWith("Flux V :"))
3453 {
3454 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3455 if (replyData.size() >= 4)
3456 magnitude = QString(replyData[3]).toDouble(&ok);
3457 }
3458 else if (replyLines[i].startsWith("Angular size:"))
3459 sizeStr = replyLines[i].mid(13, replyLines[i].indexOf("~") - 13).trimmed();
3460 else if (replyLines[i].startsWith("=========================================="))
3461 {
3462 // End of data so add the object
3463 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3464 dataRow++;
3465 }
3466 }
3467 else if (numObjs > 1 && i >= 9)
3468 {
3469 // These are the data lines - 1 per object
3470 QStringList replyData = replyLines[i].split(QLatin1Char('|'));
3471 num = 0;
3472 name = type = coord = sizeStr = "";
3473 magnitude = dist = 0.0;
3474 for (int j = 0; j < replyData.size(); j++)
3475 {
3476 if (j == 0)
3477 num = QString(replyData[j]).toInt(&ok);
3478 else if (j == 1)
3479 dist = QString(replyData[j]).toDouble(&ok) / 60.0; // Convert from arcsecs to arcmins
3480 else if (j == 2)
3481 name = (QString(replyData[j]).trimmed());
3482 else if (j == 3)
3483 type = QString(replyData[j]).trimmed();
3484 else if (j == 4)
3485 coord = QString(replyData[j]).trimmed();
3486 else if (j == 7) // Mag V
3487 magnitude = QString(replyData[j]).toDouble(&ok);
3488 else if (j == 10) // Angular size
3489 sizeStr = QString(replyData[j]).trimmed();
3490 }
3491
3492 if (num == 0 || name.isEmpty() || type.isEmpty() || coord.isEmpty())
3493 continue;
3494
3495 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3496 dataRow++;
3497 }
3498 }
3499 if (numObjs != ++dataRow)
3500 qCDebug(KSTARS_FITS) << "Simbad Header rows:" << numObjs << ". Data rows:" << dataRow;
3501 // Signal fitsview to update based on new data. Fitsview will signal fitstab
3502 m_CatObjectsSearched = false;
3503 // Signal FITSView that there is new data to process
3504 emit dataChanged();
3505 // Signal FITSTab that there is new data to process
3506 emit loadedCatalogData();
3507 simbadReply->deleteLater();
3508}
3509#endif
3510
3511bool FITSData::addCatObject(const int num, const QString name, const QString type, const QString coord, const double dist,
3512 const double magnitude, const QString sizeStr)
3513{
3514 bool ok;
3515 dms r(0.0), d(0.0);
3516
3517 QStringList split = coord.split(QLatin1Char(' '));
3518 if (split.size() != 6)
3519 {
3520 qCDebug(KSTARS_FITS) << "Coordinates for " << name << "invalid: " << coord;
3521 return false;
3522 }
3523
3524 QString raStr = QString("%1 %2 %3").arg(split[0]).arg(split[1]).arg(split[2]);
3525 QString decStr = QString("%1 %2 %3").arg(split[3]).arg(split[4]).arg(split[5]);
3526 ok = r.setFromString(raStr, false); // angle in hours
3527 if (ok)
3528 ok = d.setFromString(decStr, true); // angle in degrees
3529 if (!ok)
3530 {
3531 qCDebug(KSTARS_FITS) << "Coordinates for " << name << "invalid: " << coord;
3532 return false;
3533 }
3534
3535 double size = 0.0;
3536 if (!sizeStr.contains("~"))
3537 {
3538 QStringList sizeSplit = sizeStr.split(QLatin1Char(' '));
3539 if (sizeSplit.size() >= 1)
3540 // Take the first (major axis) size. Note sometimes there's a double space between axes so
3541 // number of elements can be 3.
3542 size = QString(sizeSplit[0]).toDouble(&ok);
3543 else
3544 qCDebug(KSTARS_FITS) << "Angular size for " << name << "invalid: " << sizeStr;
3545 }
3546 CatObject catObject;
3547 catObject.num = num;
3548 catObject.catType = CAT_SIMBAD;
3549 catObject.name = name;
3550 catObject.typeCode = type;
3551 catObject.typeLabel = getCatObjectLabel(type);
3552 catObject.r = r;
3553 catObject.d = d;
3554 catObject.dist = dist;
3555 catObject.magnitude = magnitude;
3556 catObject.size = size;
3557 catObject.highlight = false;
3558 catObject.show = getCatObjectFilter(type);
3559
3560 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3561 int stat[2];
3562
3563 world[0] = r.Degrees();
3564 world[1] = d.Degrees();
3565
3566 int status = wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]);
3567 if (status != 0)
3568 {
3569 qCDebug(KSTARS_FITS) << QString("wcss2p error processing object %1 %2: %3.").arg(name).arg(status).arg(wcs_errmsg[status]);
3570 return false;
3571 }
3572
3573 //X and Y are set to the found position if in image
3574 double x = pixcrd[0];
3575 double y = pixcrd[1];
3576 if (x > 0 && y > 0 && x < width() && y < height())
3577 {
3578 catObject.x = x;
3579 catObject.y = y;
3580 m_CatObjects.append(catObject);
3581 }
3582 return true;
3583}
3584
3585void FITSData::setCatObjectsFilters(const QStringList filters)
3586{
3587 m_CatObjectsFilters = filters;
3588}
3589
3590#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3591void FITSData::setCatSearchROI(const QPoint searchCenter, const int radius)
3592{
3593 m_CatROIPt = searchCenter;
3594 m_CatROIRadius = radius;
3595 Q_UNUSED(searchCenter);
3596 Q_UNUSED(radius);
3597}
3598#endif
3599
3600#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3601QString FITSData::getCatObjectLabel(const QString code) const
3602{
3603 QString label;
3604 if (code != "")
3605 {
3606 for (int i = 0; i < MAX_CAT_OBJ_TYPES; i++)
3607 {
3608 if (code == catObjTypes[i].code)
3609 {
3610 label = catObjTypes[i].label;
3611 break;
3612 }
3613 else if (code == catObjTypes[i].candidateCode)
3614 {
3615 label = catObjTypes[i].label + "_Cand";
3616 break;
3617 }
3618 }
3619 }
3620 return label;
3621}
3622#endif
3623
3624#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3625bool FITSData::getCatObjectFilter(const QString type) const
3626{
3627 if (m_CatObjectsFilters.isEmpty() || m_CatObjectsFilters.contains(type))
3628 return true;
3629 return false;
3630}
3631#endif
3632
3633void FITSData::filterCatObjects()
3634{
3635#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3636 bool changed = false;
3637 for (int i = 0; i < m_CatObjects.size(); i++)
3638 {
3639 bool showState = getCatObjectFilter(m_CatObjects[i].typeCode);
3640 if (m_CatObjects[i].show != showState)
3641 {
3642 changed = true;
3643 m_CatObjects[i].show = showState;
3644 }
3645 }
3646 if (changed)
3647 emit dataChanged();
3648#endif
3649}
3650
3651bool FITSData::highlightCatObject(const int hilite, const int lolite)
3652{
3653#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3654 if (hilite < 0 || hilite >= m_CatObjects.size())
3655 return false;
3656
3657 // Lowlight previous highlighted catalog object if passed
3658 if (lolite >= 0 && lolite < m_CatObjects.size())
3659 m_CatObjects[lolite].highlight = false;
3660 else
3661 {
3662 // Loop through data to find current highlighted item
3663 for (int i = 0; i < m_CatObjects.size(); i++)
3664 {
3665 if (m_CatObjects[i].highlight)
3666 {
3667 m_CatObjects[i].highlight = false;
3668 break;
3669 }
3670 }
3671 }
3672
3673 // Highlight the catalog object
3674 m_CatObjects[hilite].highlight = true;
3675#endif
3676 return true;
3677}
3678
3679int FITSData::getFlipVCounter() const
3680{
3681 return flipVCounter;
3682}
3683
3684void FITSData::setFlipVCounter(int value)
3685{
3686 flipVCounter = value;
3687}
3688
3689int FITSData::getFlipHCounter() const
3690{
3691 return flipHCounter;
3692}
3693
3694void FITSData::setFlipHCounter(int value)
3695{
3696 flipHCounter = value;
3697}
3698
3699int FITSData::getRotCounter() const
3700{
3701 return rotCounter;
3702}
3703
3704void FITSData::setRotCounter(int value)
3705{
3706 rotCounter = value;
3707}
3708
3709/* Rotate an image by 90, 180, or 270 degrees, with an optional
3710 * reflection across the vertical or horizontal axis.
3711 * verbose generates extra info on stdout.
3712 * return nullptr if successful or rotated image.
3713 */
3714template <typename T>
3715bool FITSData::rotFITS(int rotate, int mirror)
3716{
3717 int ny, nx;
3718 int x1, y1, x2, y2;
3719 uint8_t * rotimage = nullptr;
3720 int offset = 0;
3721
3722 if (rotate == 1)
3723 rotate = 90;
3724 else if (rotate == 2)
3725 rotate = 180;
3726 else if (rotate == 3)
3727 rotate = 270;
3728 else if (rotate < 0)
3729 rotate = rotate + 360;
3730
3731 nx = m_Statistics.width;
3732 ny = m_Statistics.height;
3733
3734 int BBP = m_Statistics.bytesPerPixel;
3735
3736 /* Allocate buffer for rotated image */
3737 rotimage = new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
3738
3739 if (rotimage == nullptr)
3740 {
3741 qWarning() << "Unable to allocate memory for rotated image buffer!";
3742 return false;
3743 }
3744
3745 auto * rotBuffer = reinterpret_cast<T *>(rotimage);
3746 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
3747
3748 /* Mirror image without rotation */
3749 if (rotate < 45 && rotate > -45)
3750 {
3751 if (mirror == 1)
3752 {
3753 for (int i = 0; i < m_Statistics.channels; i++)
3754 {
3755 offset = m_Statistics.samples_per_channel * i;
3756 for (x1 = 0; x1 < nx; x1++)
3757 {
3758 x2 = nx - x1 - 1;
3759 for (y1 = 0; y1 < ny; y1++)
3760 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3761 }
3762 }
3763 }
3764 else if (mirror == 2)
3765 {
3766 for (int i = 0; i < m_Statistics.channels; i++)
3767 {
3768 offset = m_Statistics.samples_per_channel * i;
3769 for (y1 = 0; y1 < ny; y1++)
3770 {
3771 y2 = ny - y1 - 1;
3772 for (x1 = 0; x1 < nx; x1++)
3773 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3774 }
3775 }
3776 }
3777 else
3778 {
3779 for (int i = 0; i < m_Statistics.channels; i++)
3780 {
3781 offset = m_Statistics.samples_per_channel * i;
3782 for (y1 = 0; y1 < ny; y1++)
3783 {
3784 for (x1 = 0; x1 < nx; x1++)
3785 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3786 }
3787 }
3788 }
3789 }
3790
3791 /* Rotate by 90 degrees */
3792 else if (rotate >= 45 && rotate < 135)
3793 {
3794 if (mirror == 1)
3795 {
3796 for (int i = 0; i < m_Statistics.channels; i++)
3797 {
3798 offset = m_Statistics.samples_per_channel * i;
3799 for (y1 = 0; y1 < ny; y1++)
3800 {
3801 x2 = ny - y1 - 1;
3802 for (x1 = 0; x1 < nx; x1++)
3803 {
3804 y2 = nx - x1 - 1;
3805 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3806 }
3807 }
3808 }
3809 }
3810 else if (mirror == 2)
3811 {
3812 for (int i = 0; i < m_Statistics.channels; i++)
3813 {
3814 offset = m_Statistics.samples_per_channel * i;
3815 for (y1 = 0; y1 < ny; y1++)
3816 {
3817 for (x1 = 0; x1 < nx; x1++)
3818 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3819 }
3820 }
3821 }
3822 else
3823 {
3824 for (int i = 0; i < m_Statistics.channels; i++)
3825 {
3826 offset = m_Statistics.samples_per_channel * i;
3827 for (y1 = 0; y1 < ny; y1++)
3828 {
3829 x2 = ny - y1 - 1;
3830 for (x1 = 0; x1 < nx; x1++)
3831 {
3832 y2 = x1;
3833 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3834 }
3835 }
3836 }
3837 }
3838
3839 m_Statistics.width = ny;
3840 m_Statistics.height = nx;
3841 }
3842
3843 /* Rotate by 180 degrees */
3844 else if (rotate >= 135 && rotate < 225)
3845 {
3846 if (mirror == 1)
3847 {
3848 for (int i = 0; i < m_Statistics.channels; i++)
3849 {
3850 offset = m_Statistics.samples_per_channel * i;
3851 for (y1 = 0; y1 < ny; y1++)
3852 {
3853 y2 = ny - y1 - 1;
3854 for (x1 = 0; x1 < nx; x1++)
3855 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3856 }
3857 }
3858 }
3859 else if (mirror == 2)
3860 {
3861 for (int i = 0; i < m_Statistics.channels; i++)
3862 {
3863 offset = m_Statistics.samples_per_channel * i;
3864 for (x1 = 0; x1 < nx; x1++)
3865 {
3866 x2 = nx - x1 - 1;
3867 for (y1 = 0; y1 < ny; y1++)
3868 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3869 }
3870 }
3871 }
3872 else
3873 {
3874 for (int i = 0; i < m_Statistics.channels; i++)
3875 {
3876 offset = m_Statistics.samples_per_channel * i;
3877 for (y1 = 0; y1 < ny; y1++)
3878 {
3879 y2 = ny - y1 - 1;
3880 for (x1 = 0; x1 < nx; x1++)
3881 {
3882 x2 = nx - x1 - 1;
3883 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3884 }
3885 }
3886 }
3887 }
3888 }
3889
3890 /* Rotate by 270 degrees */
3891 else if (rotate >= 225 && rotate < 315)
3892 {
3893 if (mirror == 1)
3894 {
3895 for (int i = 0; i < m_Statistics.channels; i++)
3896 {
3897 offset = m_Statistics.samples_per_channel * i;
3898 for (y1 = 0; y1 < ny; y1++)
3899 {
3900 for (x1 = 0; x1 < nx; x1++)
3901 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3902 }
3903 }
3904 }
3905 else if (mirror == 2)
3906 {
3907 for (int i = 0; i < m_Statistics.channels; i++)
3908 {
3909 offset = m_Statistics.samples_per_channel * i;
3910 for (y1 = 0; y1 < ny; y1++)
3911 {
3912 x2 = ny - y1 - 1;
3913 for (x1 = 0; x1 < nx; x1++)
3914 {
3915 y2 = nx - x1 - 1;
3916 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3917 }
3918 }
3919 }
3920 }
3921 else
3922 {
3923 for (int i = 0; i < m_Statistics.channels; i++)
3924 {
3925 offset = m_Statistics.samples_per_channel * i;
3926 for (y1 = 0; y1 < ny; y1++)
3927 {
3928 x2 = y1;
3929 for (x1 = 0; x1 < nx; x1++)
3930 {
3931 y2 = nx - x1 - 1;
3932 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3933 }
3934 }
3935 }
3936 }
3937
3938 m_Statistics.width = ny;
3939 m_Statistics.height = nx;
3940 }
3941
3942 /* If rotating by more than 315 degrees, assume top-bottom reflection */
3943 else if (rotate >= 315 && mirror)
3944 {
3945 for (int i = 0; i < m_Statistics.channels; i++)
3946 {
3947 offset = m_Statistics.samples_per_channel * i;
3948 for (y1 = 0; y1 < ny; y1++)
3949 {
3950 for (x1 = 0; x1 < nx; x1++)
3951 {
3952 x2 = y1;
3953 y2 = x1;
3954 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3955 }
3956 }
3957 }
3958 }
3959
3960 delete[] m_ImageBuffer;
3961 m_ImageBuffer = rotimage;
3962
3963 return true;
3964}
3965
3966void FITSData::rotWCSFITS(int angle, int mirror)
3967{
3968 int status = 0;
3969 char comment[100];
3970 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
3971 int WCS_DECIMALS = 6;
3972
3973 naxis1 = m_Statistics.width;
3974 naxis2 = m_Statistics.height;
3975
3976 if (fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
3977 {
3978 // No WCS keywords
3979 return;
3980 }
3981
3982 /* Reset CROTAn and CD matrix if axes have been exchanged */
3983 if (angle == 90)
3984 {
3985 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
3986 fits_update_key_dbl(fptr, "CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3987
3988 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
3989 fits_update_key_dbl(fptr, "CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3990 }
3991
3992 status = 0;
3993
3994 /* Negate rotation angle if mirrored */
3995 if (mirror != 0)
3996 {
3997 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
3998 fits_update_key_dbl(fptr, "CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
3999
4000 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
4001 fits_update_key_dbl(fptr, "CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
4002
4003 status = 0;
4004
4005 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
4006 fits_update_key_dbl(fptr, "LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4007
4008 status = 0;
4009
4010 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
4011 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4012
4013 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp1, comment, &status))
4014 fits_update_key_dbl(fptr, "CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
4015
4016 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp1, comment, &status))
4017 fits_update_key_dbl(fptr, "CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
4018 }
4019
4020 status = 0;
4021
4022 /* Unbin CRPIX and CD matrix */
4023 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
4024 {
4025 if (ctemp1 != 1.0)
4026 {
4027 if (!fits_read_key_dbl(fptr, "LTM2_2", &ctemp2, comment, &status))
4028 if (ctemp1 == ctemp2)
4029 {
4030 double ltv1 = 0.0;
4031 double ltv2 = 0.0;
4032 status = 0;
4033 if (!fits_read_key_dbl(fptr, "LTV1", &ltv1, comment, &status))
4034 fits_delete_key(fptr, "LTV1", &status);
4035 if (!fits_read_key_dbl(fptr, "LTV2", &ltv2, comment, &status))
4036 fits_delete_key(fptr, "LTV2", &status);
4037
4038 status = 0;
4039
4040 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp3, comment, &status))
4041 fits_update_key_dbl(fptr, "CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
4042
4043 if (!fits_read_key_dbl(fptr, "CRPIX2", &ctemp3, comment, &status))
4044 fits_update_key_dbl(fptr, "CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
4045
4046 status = 0;
4047
4048 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp3, comment, &status))
4049 fits_update_key_dbl(fptr, "CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4050
4051 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp3, comment, &status))
4052 fits_update_key_dbl(fptr, "CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4053
4054 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status))
4055 fits_update_key_dbl(fptr, "CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4056
4057 if (!fits_read_key_dbl(fptr, "CD2_2", &ctemp3, comment, &status))
4058 fits_update_key_dbl(fptr, "CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4059
4060 status = 0;
4061
4062 fits_delete_key(fptr, "LTM1_1", &status);
4063 fits_delete_key(fptr, "LTM1_2", &status);
4064 }
4065 }
4066 }
4067
4068 status = 0;
4069
4070 /* Reset CRPIXn */
4071 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp1, comment, &status) &&
4072 !fits_read_key_dbl(fptr, "CRPIX2", &ctemp2, comment, &status))
4073 {
4074 if (mirror != 0)
4075 {
4076 if (angle == 0)
4077 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4078 else if (angle == 90)
4079 {
4080 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4081 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4082 }
4083 else if (angle == 180)
4084 {
4085 fits_update_key_dbl(fptr, "CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
4086 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4087 }
4088 else if (angle == 270)
4089 {
4090 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4091 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4092 }
4093 }
4094 else
4095 {
4096 if (angle == 90)
4097 {
4098 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4099 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4100 }
4101 else if (angle == 180)
4102 {
4103 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4104 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4105 }
4106 else if (angle == 270)
4107 {
4108 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4109 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4110 }
4111 }
4112 }
4113
4114 status = 0;
4115
4116 /* Reset CDELTn (degrees per pixel) */
4117 if (!fits_read_key_dbl(fptr, "CDELT1", &ctemp1, comment, &status) &&
4118 !fits_read_key_dbl(fptr, "CDELT2", &ctemp2, comment, &status))
4119 {
4120 if (mirror != 0)
4121 {
4122 if (angle == 0)
4123 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4124 else if (angle == 90)
4125 {
4126 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4127 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4128 }
4129 else if (angle == 180)
4130 {
4131 fits_update_key_dbl(fptr, "CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
4132 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4133 }
4134 else if (angle == 270)
4135 {
4136 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4137 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4138 }
4139 }
4140 else
4141 {
4142 if (angle == 90)
4143 {
4144 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4145 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4146 }
4147 else if (angle == 180)
4148 {
4149 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4150 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4151 }
4152 else if (angle == 270)
4153 {
4154 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4155 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4156 }
4157 }
4158 }
4159
4160 /* Reset CD matrix, if present */
4161 ctemp1 = 0.0;
4162 ctemp2 = 0.0;
4163 ctemp3 = 0.0;
4164 ctemp4 = 0.0;
4165 status = 0;
4166 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
4167 {
4168 fits_read_key_dbl(fptr, "CD1_2", &ctemp2, comment, &status);
4169 fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status);
4170 fits_read_key_dbl(fptr, "CD2_2", &ctemp4, comment, &status);
4171 status = 0;
4172 if (mirror != 0)
4173 {
4174 if (angle == 0)
4175 {
4176 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4177 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4178 }
4179 else if (angle == 90)
4180 {
4181 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4182 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4183 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4184 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4185 }
4186 else if (angle == 180)
4187 {
4188 fits_update_key_dbl(fptr, "CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
4189 fits_update_key_dbl(fptr, "CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
4190 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4191 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4192 }
4193 else if (angle == 270)
4194 {
4195 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4196 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4197 fits_update_key_dbl(fptr, "CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
4198 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4199 }
4200 }
4201 else
4202 {
4203 if (angle == 90)
4204 {
4205 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4206 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4207 fits_update_key_dbl(fptr, "CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
4208 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4209 }
4210 else if (angle == 180)
4211 {
4212 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4213 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4214 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4215 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4216 }
4217 else if (angle == 270)
4218 {
4219 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4220 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4221 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4222 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4223 }
4224 }
4225 }
4226
4227 /* Delete any polynomial solution */
4228 /* (These could maybe be switched, but I don't want to work them out yet */
4229 status = 0;
4230 if (!fits_read_key_dbl(fptr, "CO1_1", &ctemp1, comment, &status))
4231 {
4232 int i;
4233 char keyword[16];
4234
4235 for (i = 1; i < 13; i++)
4236 {
4237 sprintf(keyword, "CO1_%d", i);
4238 fits_delete_key(fptr, keyword, &status);
4239 }
4240 for (i = 1; i < 13; i++)
4241 {
4242 sprintf(keyword, "CO2_%d", i);
4243 fits_delete_key(fptr, keyword, &status);
4244 }
4245 }
4246
4247}
4248
4249uint8_t * FITSData::getWritableImageBuffer()
4250{
4251 return m_ImageBuffer;
4252}
4253
4254uint8_t const * FITSData::getImageBuffer() const
4255{
4256 return m_ImageBuffer;
4257}
4258
4259void FITSData::setImageBuffer(uint8_t * buffer)
4260{
4261 delete[] m_ImageBuffer;
4262 m_ImageBuffer = buffer;
4263}
4264
4265bool FITSData::checkDebayer()
4266{
4267 int status = 0;
4268 char bayerPattern[64], roworder[64];
4269
4270 // Let's search for BAYERPAT keyword, if it's not found we return as there is no bayer pattern in this image
4271 if (fits_read_keyword(fptr, "BAYERPAT", bayerPattern, nullptr, &status))
4272 return false;
4273
4274 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
4275 {
4276 m_LastError = i18n("Only 8 and 16 bits bayered images supported.");
4277 return false;
4278 }
4279 QString pattern(bayerPattern);
4280 pattern = pattern.remove('\'').trimmed();
4281
4282 QString order(roworder);
4283 order = order.remove('\'').trimmed();
4284
4285 if (order == "BOTTOM-UP" && !(m_Statistics.height % 2))
4286 {
4287 if (pattern == "RGGB")
4288 pattern = "GBRG";
4289 else if (pattern == "GBRG")
4290 pattern = "RGGB";
4291 else if (pattern == "GRBG")
4292 pattern = "BGGR";
4293 else if (pattern == "BGGR")
4294 pattern = "GRBG";
4295 else return false;
4296 }
4297
4298 if (pattern == "RGGB")
4299 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4300 else if (pattern == "GBRG")
4301 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4302 else if (pattern == "GRBG")
4303 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4304 else if (pattern == "BGGR")
4305 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4306 // We return unless we find a valid pattern
4307 else
4308 {
4309 m_LastError = i18n("Unsupported bayer pattern %1.", pattern);
4310 return false;
4311 }
4312
4313 fits_read_key(fptr, TINT, "XBAYROFF", &debayerParams.offsetX, nullptr, &status);
4314 fits_read_key(fptr, TINT, "YBAYROFF", &debayerParams.offsetY, nullptr, &status);
4315
4316 if (debayerParams.offsetX == 1)
4317 {
4318 // This may leave odd values in the 0th column if the color filter is not there
4319 // in the sensor, but otherwise should process the offset correctly.
4320 // Only offsets of 0 or 1 are implemented in debayer_8bit() and debayer_16bit().
4321 switch (debayerParams.filter)
4322 {
4323 case DC1394_COLOR_FILTER_RGGB:
4324 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4325 break;
4326 case DC1394_COLOR_FILTER_GBRG:
4327 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4328 break;
4329 case DC1394_COLOR_FILTER_GRBG:
4330 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4331 break;
4332 case DC1394_COLOR_FILTER_BGGR:
4333 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4334 break;
4335 }
4336 debayerParams.offsetX = 0;
4337 }
4338 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
4339 {
4340 m_LastError = i18n("Unsupported bayer offsets %1 %2.", debayerParams.offsetX, debayerParams.offsetY);
4341 return false;
4342 }
4343
4344 HasDebayer = true;
4345
4346 return true;
4347}
4348
4349void FITSData::getBayerParams(BayerParams * param)
4350{
4351 param->method = debayerParams.method;
4352 param->filter = debayerParams.filter;
4353 param->offsetX = debayerParams.offsetX;
4354 param->offsetY = debayerParams.offsetY;
4355}
4356
4357void FITSData::setBayerParams(BayerParams * param)
4358{
4359 debayerParams.method = param->method;
4360 debayerParams.filter = param->filter;
4361 debayerParams.offsetX = param->offsetX;
4362 debayerParams.offsetY = param->offsetY;
4363}
4364
4365bool FITSData::debayer(bool reload)
4366{
4367 if (reload)
4368 {
4369 int anynull = 0, status = 0;
4370
4371 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel, nullptr, m_ImageBuffer,
4372 &anynull, &status))
4373 {
4374 // char errmsg[512];
4375 // fits_get_errstatus(status, errmsg);
4376 // KSNotification::error(i18n("Error reading image: %1", QString(errmsg)), i18n("Debayer error"));
4377 return false;
4378 }
4379 }
4380
4381 switch (m_Statistics.dataType)
4382 {
4383 case TBYTE:
4384 return debayer_8bit();
4385
4386 case TUSHORT:
4387 return debayer_16bit();
4388
4389 default:
4390 return false;
4391 }
4392}
4393
4394bool FITSData::debayer_8bit()
4395{
4396 dc1394error_t error_code;
4397
4398 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4399 uint8_t * destinationBuffer = nullptr;
4400
4401 try
4402 {
4403 destinationBuffer = new uint8_t[rgb_size];
4404 }
4405 catch (const std::bad_alloc &e)
4406 {
4407 logOOMError(rgb_size);
4408 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4409 return false;
4410 }
4411
4412 auto * bayer_source_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
4413 auto * bayer_destination_buffer = reinterpret_cast<uint8_t *>(destinationBuffer);
4414
4415 if (bayer_destination_buffer == nullptr)
4416 {
4417 logOOMError(rgb_size);
4418 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer.");
4419 return false;
4420 }
4421
4422 int ds1394_height = m_Statistics.height;
4423 auto dc1394_source = bayer_source_buffer;
4424
4425 if (debayerParams.offsetY == 1)
4426 {
4427 dc1394_source += m_Statistics.width;
4428 ds1394_height--;
4429 }
4430 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
4431
4432 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4433 debayerParams.filter,
4434 debayerParams.method);
4435
4436 if (error_code != DC1394_SUCCESS)
4437 {
4438 m_LastError = i18n("Debayer failed (%1)", error_code);
4439 m_Statistics.channels = 1;
4440 delete[] destinationBuffer;
4441 return false;
4442 }
4443
4444 if (m_ImageBufferSize != rgb_size)
4445 {
4446 delete[] m_ImageBuffer;
4447 try
4448 {
4449 m_ImageBuffer = new uint8_t[rgb_size];
4450 }
4451 catch (const std::bad_alloc &e)
4452 {
4453 delete[] destinationBuffer;
4454 logOOMError(rgb_size);
4455 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4456 return false;
4457 }
4458
4459 m_ImageBufferSize = rgb_size;
4460 }
4461
4462 auto bayered_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
4463
4464 // Data in R1G1B1, we need to copy them into 3 layers for FITS
4465
4466 uint8_t * rBuff = bayered_buffer;
4467 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4468 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4469
4470 int imax = m_Statistics.samples_per_channel * 3 - 3;
4471 for (int i = 0; i <= imax; i += 3)
4472 {
4473 *rBuff++ = bayer_destination_buffer[i];
4474 *gBuff++ = bayer_destination_buffer[i + 1];
4475 *bBuff++ = bayer_destination_buffer[i + 2];
4476 }
4477
4478 // TODO Maybe all should be treated the same
4479 // Doing single channel saves lots of memory though for non-essential
4480 // frames
4481 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4482 m_Statistics.dataType = TBYTE;
4483 delete[] destinationBuffer;
4484 return true;
4485}
4486
4487bool FITSData::debayer_16bit()
4488{
4489 dc1394error_t error_code;
4490
4491 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4492 uint8_t *destinationBuffer = nullptr;
4493 try
4494 {
4495 destinationBuffer = new uint8_t[rgb_size];
4496 }
4497 catch (const std::bad_alloc &e)
4498 {
4499 logOOMError(rgb_size);
4500 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4501 return false;
4502 }
4503
4504 auto * bayer_source_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
4505 auto * bayer_destination_buffer = reinterpret_cast<uint16_t *>(destinationBuffer);
4506
4507 if (bayer_destination_buffer == nullptr)
4508 {
4509 logOOMError(rgb_size);
4510 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer.");
4511 return false;
4512 }
4513
4514 int ds1394_height = m_Statistics.height;
4515 auto dc1394_source = bayer_source_buffer;
4516
4517 if (debayerParams.offsetY == 1)
4518 {
4519 dc1394_source += m_Statistics.width;
4520 ds1394_height--;
4521 }
4522 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
4523
4524 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4525 debayerParams.filter,
4526 debayerParams.method, 16);
4527
4528 if (error_code != DC1394_SUCCESS)
4529 {
4530 m_LastError = i18n("Debayer failed (%1)");
4531 m_Statistics.channels = 1;
4532 delete[] destinationBuffer;
4533 return false;
4534 }
4535
4536 if (m_ImageBufferSize != rgb_size)
4537 {
4538 delete[] m_ImageBuffer;
4539 try
4540 {
4541 m_ImageBuffer = new uint8_t[rgb_size];
4542 }
4543 catch (const std::bad_alloc &e)
4544 {
4545 logOOMError(rgb_size);
4546 delete[] destinationBuffer;
4547 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4548 return false;
4549 }
4550
4551 m_ImageBufferSize = rgb_size;
4552 }
4553
4554 auto bayered_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
4555
4556 // Data in R1G1B1, we need to copy them into 3 layers for FITS
4557
4558 uint16_t * rBuff = bayered_buffer;
4559 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4560 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4561
4562 int imax = m_Statistics.samples_per_channel * 3 - 3;
4563 for (int i = 0; i <= imax; i += 3)
4564 {
4565 *rBuff++ = bayer_destination_buffer[i];
4566 *gBuff++ = bayer_destination_buffer[i + 1];
4567 *bBuff++ = bayer_destination_buffer[i + 2];
4568 }
4569
4570 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4571 m_Statistics.dataType = TUSHORT;
4572 delete[] destinationBuffer;
4573 return true;
4574}
4575
4576void FITSData::logOOMError(uint32_t requiredMemory)
4577{
4578 qCCritical(KSTARS_FITS) << "Debayed memory allocation failure. Required Memory:" << KFormat().formatByteSize(requiredMemory)
4579 << "Available system memory:" << KSUtils::getAvailableRAM();
4580}
4581
4582double FITSData::getADU() const
4583{
4584 double adu = 0;
4585 for (int i = 0; i < m_Statistics.channels; i++)
4586 adu += m_Statistics.mean[i];
4587
4588 return (adu / static_cast<double>(m_Statistics.channels));
4589}
4590
4591QString FITSData::getLastError() const
4592{
4593 return m_LastError;
4594}
4595
4596template <typename T>
4597void FITSData::convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image)
4598{
4599 const auto * buffer = reinterpret_cast<const T *>(getImageBuffer());
4600 const T limit = std::numeric_limits<T>::max();
4601 T bMin = dataMin < 0 ? 0 : dataMin;
4602 T bMax = dataMax > limit ? limit : dataMax;
4603 uint16_t w = width();
4604 uint16_t h = height();
4605 uint32_t size = w * h;
4606 double val;
4607
4608 if (channels() == 1)
4609 {
4610 /* Fill in pixel values using indexed map, linear scale */
4611 for (int j = 0; j < h; j++)
4612 {
4613 unsigned char * scanLine = image.scanLine(j);
4614
4615 for (int i = 0; i < w; i++)
4616 {
4617 val = qBound(bMin, buffer[j * w + i], bMax);
4618 val = val * scale + zero;
4619 scanLine[i] = qBound<unsigned char>(0, static_cast<uint8_t>(val), 255);
4620 }
4621 }
4622 }
4623 else
4624 {
4625 double rval = 0, gval = 0, bval = 0;
4626 QRgb value;
4627 /* Fill in pixel values using indexed map, linear scale */
4628 for (int j = 0; j < h; j++)
4629 {
4630 auto * scanLine = reinterpret_cast<QRgb *>((image.scanLine(j)));
4631
4632 for (int i = 0; i < w; i++)
4633 {
4634 rval = qBound(bMin, buffer[j * w + i], bMax);
4635 gval = qBound(bMin, buffer[j * w + i + size], bMax);
4636 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
4637
4638 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
4639
4640 scanLine[i] = value;
4641 }
4642 }
4643 }
4644}
4645
4646QImage FITSData::FITSToImage(const QString &filename)
4647{
4648 QImage fitsImage;
4649 double min, max;
4650
4651 FITSData data;
4652
4653 QFuture<bool> future = data.loadFromFile(filename);
4654
4655 // Wait synchronously
4656 future.waitForFinished();
4657 if (future.result() == false)
4658 return fitsImage;
4659
4660 data.getMinMax(&min, &max);
4661
4662 if (min == max)
4663 {
4664 fitsImage.fill(Qt::white);
4665 return fitsImage;
4666 }
4667
4668 if (data.channels() == 1)
4669 {
4670 fitsImage = QImage(data.width(), data.height(), QImage::Format_Indexed8);
4671
4672 fitsImage.setColorCount(256);
4673 for (int i = 0; i < 256; i++)
4674 fitsImage.setColor(i, qRgb(i, i, i));
4675 }
4676 else
4677 {
4678 fitsImage = QImage(data.width(), data.height(), QImage::Format_RGB32);
4679 }
4680
4681 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
4682 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
4683
4684 double bscale = 255. / (dataMax - dataMin);
4685 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
4686
4687 // Long way to do this since we do not want to use templated functions here
4688 switch (data.m_Statistics.dataType)
4689 {
4690 case TBYTE:
4691 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4692 break;
4693
4694 case TSHORT:
4695 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4696 break;
4697
4698 case TUSHORT:
4699 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4700 break;
4701
4702 case TLONG:
4703 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4704 break;
4705
4706 case TULONG:
4707 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4708 break;
4709
4710 case TFLOAT:
4711 data.convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
4712 break;
4713
4714 case TLONGLONG:
4715 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4716 break;
4717
4718 case TDOUBLE:
4719 data.convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
4720 break;
4721
4722 default:
4723 break;
4724 }
4725
4726 return fitsImage;
4727}
4728
4729bool FITSData::ImageToFITS(const QString &filename, const QString &format, QString &output)
4730{
4731 if (QImageReader::supportedImageFormats().contains(format.toLatin1()) == false)
4732 {
4733 qCCritical(KSTARS_FITS) << "Failed to convert" << filename << "to FITS since format" << format << "is not supported in Qt";
4734 return false;
4735 }
4736
4737 QImage input;
4738
4739 if (input.load(filename, format.toLatin1()) == false)
4740 {
4741 qCCritical(KSTARS_FITS) << "Failed to open image" << filename;
4742 return false;
4743 }
4744
4745 output = QDir(getTemporaryPath()).filePath(filename + ".fits");
4746
4747 //This section sets up the FITS File
4748 fitsfile *fptr = nullptr;
4749 int status = 0;
4750 long fpixel = 1, naxis = input.allGray() ? 2 : 3, nelements, exposure;
4751 long naxes[3] = { input.width(), input.height(), naxis == 3 ? 3 : 1 };
4752 char error_status[512] = {0};
4753
4754 if (fits_create_file(&fptr, QString('!' + output).toLocal8Bit().data(), &status))
4755 {
4756 fits_get_errstatus(status, error_status);
4757 qCCritical(KSTARS_FITS) << "Failed to create FITS file. Error:" << error_status;
4758 return false;
4759 }
4760
4761 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
4762 {
4763 qCWarning(KSTARS_FITS) << "fits_create_img failed:" << error_status;
4764 status = 0;
4765 fits_flush_file(fptr, &status);
4766 fits_close_file(fptr, &status);
4767 return false;
4768 }
4769
4770 exposure = 1;
4771 fits_update_key(fptr, TLONG, "EXPOSURE", &exposure, "Total Exposure Time", &status);
4772
4773 // Gray image
4774 if (naxis == 2)
4775 {
4776 nelements = naxes[0] * naxes[1];
4777 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.bits(), &status))
4778 {
4779 fits_get_errstatus(status, error_status);
4780 qCWarning(KSTARS_FITS) << "fits_write_img GRAY failed:" << error_status;
4781 status = 0;
4782 fits_flush_file(fptr, &status);
4783 fits_close_file(fptr, &status);
4784 return false;
4785 }
4786 }
4787 // RGB image, we have to convert from ARGB format to R G B for each plane
4788 else
4789 {
4790 nelements = naxes[0] * naxes[1] * 3;
4791
4792 uint8_t *srcBuffer = input.bits();
4793 // ARGB
4794 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
4795
4796 uint8_t *rgbBuffer = new uint8_t[nelements];
4797 if (rgbBuffer == nullptr)
4798 {
4799 qCWarning(KSTARS_FITS) << "Not enough memory for RGB buffer";
4800 fits_flush_file(fptr, &status);
4801 fits_close_file(fptr, &status);
4802 return false;
4803 }
4804
4805 uint8_t *subR = rgbBuffer;
4806 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
4807 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
4808 for (uint32_t i = 0; i < srcBytes; i += 4)
4809 {
4810 *subB++ = srcBuffer[i];
4811 *subG++ = srcBuffer[i + 1];
4812 *subR++ = srcBuffer[i + 2];
4813 }
4814
4815 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
4816 {
4817 fits_get_errstatus(status, error_status);
4818 qCWarning(KSTARS_FITS) << "fits_write_img RGB failed:" << error_status;
4819 status = 0;
4820 fits_flush_file(fptr, &status);
4821 fits_close_file(fptr, &status);
4822 delete [] rgbBuffer;
4823 return false;
4824 }
4825
4826 delete [] rgbBuffer;
4827 }
4828
4829 if (fits_flush_file(fptr, &status))
4830 {
4831 fits_get_errstatus(status, error_status);
4832 qCWarning(KSTARS_FITS) << "fits_flush_file failed:" << error_status;
4833 status = 0;
4834 fits_close_file(fptr, &status);
4835 return false;
4836 }
4837
4838 if (fits_close_file(fptr, &status))
4839 {
4840 fits_get_errstatus(status, error_status);
4841 qCWarning(KSTARS_FITS) << "fits_close_file failed:" << error_status;
4842 return false;
4843 }
4844
4845 return true;
4846}
4847
4848void FITSData::injectWCS(double orientation, double ra, double dec, double pixscale, bool eastToTheRight)
4849{
4850 int status = 0;
4851
4852 updateWCSHeaderData(orientation, ra, dec, pixscale, eastToTheRight);
4853
4854 if (fptr)
4855 {
4856 fits_update_key(fptr, TDOUBLE, "OBJCTRA", &ra, "Object RA", &status);
4857 fits_update_key(fptr, TDOUBLE, "OBJCTDEC", &dec, "Object DEC", &status);
4858
4859 int epoch = 2000;
4860
4861 fits_update_key(fptr, TINT, "EQUINOX", &epoch, "Equinox", &status);
4862
4863 fits_update_key(fptr, TDOUBLE, "CRVAL1", &ra, "CRVAL1", &status);
4864 fits_update_key(fptr, TDOUBLE, "CRVAL2", &dec, "CRVAL1", &status);
4865
4866 char radecsys[8] = "FK5";
4867 char ctype1[16] = "RA---TAN";
4868 char ctype2[16] = "DEC--TAN";
4869
4870 fits_update_key(fptr, TSTRING, "RADECSYS", radecsys, "RADECSYS", &status);
4871 fits_update_key(fptr, TSTRING, "CTYPE1", ctype1, "CTYPE1", &status);
4872 fits_update_key(fptr, TSTRING, "CTYPE2", ctype2, "CTYPE2", &status);
4873
4874 double crpix1 = width() / 2.0;
4875 double crpix2 = height() / 2.0;
4876
4877 fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix1, "CRPIX1", &status);
4878 fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix2, "CRPIX2", &status);
4879
4880 // Arcsecs per Pixel
4881 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4882 double secpix2 = pixscale;
4883
4884 fits_update_key(fptr, TDOUBLE, "SECPIX1", &secpix1, "SECPIX1", &status);
4885 fits_update_key(fptr, TDOUBLE, "SECPIX2", &secpix2, "SECPIX2", &status);
4886
4887 double degpix1 = secpix1 / 3600.0;
4888 double degpix2 = secpix2 / 3600.0;
4889
4890 fits_update_key(fptr, TDOUBLE, "CDELT1", &degpix1, "CDELT1", &status);
4891 fits_update_key(fptr, TDOUBLE, "CDELT2", &degpix2, "CDELT2", &status);
4892
4893 // Rotation is CW, we need to convert it to CCW per CROTA1 definition
4894 double rotation = 360 - orientation;
4895 if (rotation > 360)
4896 rotation -= 360;
4897
4898 fits_update_key(fptr, TDOUBLE, "CROTA1", &rotation, "CROTA1", &status);
4899 fits_update_key(fptr, TDOUBLE, "CROTA2", &rotation, "CROTA2", &status);
4900 }
4901
4902 emit headerChanged();
4903 m_WCSState = Idle;
4904}
4905
4906// Update header records based on plate solver solution
4907void FITSData::updateWCSHeaderData(const double orientation, const double ra, const double dec, const double pixscale,
4908 const bool eastToTheRight)
4909{
4910 updateRecordValue("OBJCTRA", ra, "Object RA");
4911 updateRecordValue("OBJCTDEC", dec, "Object DEC");
4912 updateRecordValue("EQUINOX", 2000, "Equinox");
4913 updateRecordValue("CRVAL1", ra, "CRVAL1");
4914 updateRecordValue("CRVAL2", dec, "CRVAL2");
4915
4916 updateRecordValue("RADECSYS", "'FK5'", "RADECSYS");
4917 updateRecordValue("CTYPE1", "'RA---TAN'", "CTYPE1");
4918 updateRecordValue("CTYPE2", "'DEC--TAN'", "CTYPE2");
4919
4920 updateRecordValue("CRPIX1", m_Statistics.width / 2.0, "CRPIX1");
4921 updateRecordValue("CRPIX2", m_Statistics.height / 2.0, "CRPIX2");
4922
4923 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4924 double secpix2 = pixscale;
4925
4926 double degpix1 = secpix1 / 3600.0;
4927 double degpix2 = secpix2 / 3600.0;
4928 updateRecordValue("CDELT1", degpix1, "CDELT1");
4929 updateRecordValue("CDELT2", degpix2, "CDELT2");
4930
4931 // Rotation is CW, we need to convert it to CCW per CROTA1 definition
4932 double rotation = 360 - orientation;
4933 if (rotation > 360)
4934 rotation -= 360;
4935 updateRecordValue("CROTA1", rotation, "CROTA1");
4936 updateRecordValue("CROTA2", rotation, "CROTA2");
4937}
4938
4939bool FITSData::contains(const QPointF &point) const
4940{
4941 return (point.x() >= 0 && point.y() >= 0 && point.x() <= m_Statistics.width && point.y() <= m_Statistics.height);
4942}
4943
4944void FITSData::saveStatistics(FITSImage::Statistic &other)
4945{
4946 other = m_Statistics;
4947}
4948
4949void FITSData::restoreStatistics(FITSImage::Statistic &other)
4950{
4951 m_Statistics = other;
4952
4953 emit dataChanged();
4954}
4955
4956void FITSData::constructHistogram()
4957{
4958 switch (m_Statistics.dataType)
4959 {
4960 case TBYTE:
4961 constructHistogramInternal<uint8_t>();
4962 break;
4963
4964 case TSHORT:
4965 constructHistogramInternal<int16_t>();
4966 break;
4967
4968 case TUSHORT:
4969 constructHistogramInternal<uint16_t>();
4970 break;
4971
4972 case TLONG:
4973 constructHistogramInternal<int32_t>();
4974 break;
4975
4976 case TULONG:
4977 constructHistogramInternal<uint32_t>();
4978 break;
4979
4980 case TFLOAT:
4981 constructHistogramInternal<float>();
4982 break;
4983
4984 case TLONGLONG:
4985 constructHistogramInternal<int64_t>();
4986 break;
4987
4988 case TDOUBLE:
4989 constructHistogramInternal<double>();
4990 break;
4991
4992 default:
4993 break;
4994 }
4995}
4996
4997template <typename T> int32_t FITSData::histogramBinInternal(T value, int channel) const
4998{
4999 return qMax(static_cast<T>(0), qMin(static_cast<T>(m_HistogramBinCount),
5000 static_cast<T>(rint((value - m_Statistics.min[channel]) / m_HistogramBinWidth[channel]))));
5001}
5002
5003template <typename T> int32_t FITSData::histogramBinInternal(int x, int y, int channel) const
5004{
5005 if (!m_ImageBuffer || !isHistogramConstructed())
5006 return 0;
5007 uint32_t samples = m_Statistics.width * m_Statistics.height;
5008 uint32_t offset = channel * samples;
5009 auto * const buffer = reinterpret_cast<T const *>(m_ImageBuffer);
5010 int index = y * m_Statistics.width + x;
5011 const T &sample = buffer[index + offset];
5012 return histogramBinInternal(sample, channel);
5013}
5014
5015int32_t FITSData::histogramBin(int x, int y, int channel) const
5016{
5017 switch (m_Statistics.dataType)
5018 {
5019 case TBYTE:
5020 return histogramBinInternal<uint8_t>(x, y, channel);
5021 break;
5022
5023 case TSHORT:
5024 return histogramBinInternal<int16_t>(x, y, channel);
5025 break;
5026
5027 case TUSHORT:
5028 return histogramBinInternal<uint16_t>(x, y, channel);
5029 break;
5030
5031 case TLONG:
5032 return histogramBinInternal<int32_t>(x, y, channel);
5033 break;
5034
5035 case TULONG:
5036 return histogramBinInternal<uint32_t>(x, y, channel);
5037 break;
5038
5039 case TFLOAT:
5040 return histogramBinInternal<float>(x, y, channel);
5041 break;
5042
5043 case TLONGLONG:
5044 return histogramBinInternal<int64_t>(x, y, channel);
5045 break;
5046
5047 case TDOUBLE:
5048 return histogramBinInternal<double>(x, y, channel);
5049 break;
5050
5051 default:
5052 return 0;
5053 break;
5054 }
5055}
5056
5057template <typename T> void FITSData::constructHistogramInternal()
5058{
5059 auto * const buffer = reinterpret_cast<T const *>(m_ImageBuffer);
5060 uint32_t samples = m_Statistics.width * m_Statistics.height;
5061 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
5062
5063 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
5064 if (m_HistogramBinCount <= 0)
5065 m_HistogramBinCount = 256;
5066
5067 for (int n = 0; n < m_Statistics.channels; n++)
5068 {
5069 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
5070 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
5071 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
5072 // Distinguish between 0-1.0 ranges and ranges with integer values.
5073 const double minBinSize = (m_Statistics.max[n] > 1.1) ? 1.0 : .0001;
5074 m_HistogramBinWidth[n] = qMax(minBinSize, (m_Statistics.max[n] - m_Statistics.min[n]) / m_HistogramBinCount);
5075 }
5076
5077 QVector<QFuture<void>> futures;
5078
5079 for (int n = 0; n < m_Statistics.channels; n++)
5080 {
5081 futures.append(QtConcurrent::run([ = ]()
5082 {
5083 for (int i = 0; i < m_HistogramBinCount; i++)
5084 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
5085 }));
5086 }
5087
5088 for (int n = 0; n < m_Statistics.channels; n++)
5089 {
5090 futures.append(QtConcurrent::run([ = ]()
5091 {
5092 uint32_t offset = n * samples;
5093
5094 for (uint32_t i = 0; i < samples; i += sampleBy)
5095 {
5096 int32_t id = histogramBinInternal<T>(buffer[i + offset], n);
5097 m_HistogramFrequency[n][id] += sampleBy;
5098 }
5099 }));
5100 }
5101
5102 for (QFuture<void> future : futures)
5103 future.waitForFinished();
5104
5105 futures.clear();
5106
5107 for (int n = 0; n < m_Statistics.channels; n++)
5108 {
5109 futures.append(QtConcurrent::run([ = ]()
5110 {
5111 uint32_t accumulator = 0;
5112 for (int i = 0; i < m_HistogramBinCount; i++)
5113 {
5114 accumulator += m_HistogramFrequency[n][i];
5115 m_CumulativeFrequency[n].replace(i, accumulator);
5116 }
5117 }));
5118 }
5119
5120 for (QFuture<void> future : futures)
5121 future.waitForFinished();
5122
5123 futures.clear();
5124
5125 // Custom index to indicate the overall contrast of the image
5126 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
5127 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] / static_cast<double>
5128 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
5129 4]);
5130 else
5131 m_JMIndex = 1;
5132
5133 qCDebug(KSTARS_FITS) << "FITHistogram: JMIndex " << m_JMIndex;
5134
5135 m_HistogramConstructed = true;
5136 emit histogramReady();
5137}
5138
5139void FITSData::recordLastError(int errorCode)
5140{
5141 char fitsErrorMessage[512] = {0};
5142 fits_get_errstatus(errorCode, fitsErrorMessage);
5143 m_LastError = fitsErrorMessage;
5144}
5145
5146double FITSData::getAverageMean(bool roi) const
5147{
5148 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5149 if (ptr->channels == 1)
5150 return ptr->mean[0];
5151 else
5152 return (ptr->mean[0] + ptr->mean[1] + ptr->mean[2]) / 3.0;
5153}
5154
5155double FITSData::getAverageMedian(bool roi) const
5156{
5157 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5158 if (ptr->channels == 1)
5159 return ptr->median[0];
5160 else
5161 return (ptr->median[0] + ptr->median[1] + ptr->median[2]) / 3.0;
5162}
5163
5164double FITSData::getAverageStdDev(bool roi) const
5165{
5166 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5167 if (ptr->channels == 1)
5168 return ptr->stddev[0];
5169 else
5170 return (ptr->stddev[0] + ptr->stddev[1] + ptr->stddev[2]) / 3.0;
5171}
QString formatByteSize(double size, int precision=1, KFormat::BinaryUnitDialect dialect=KFormat::DefaultBinaryDialect, KFormat::BinarySizeUnits units=KFormat::DefaultBinaryUnits) const
There are several time-dependent values used in position calculations, that are not specific to an ob...
Definition ksnumbers.h:43
SkyMapComposite * skyComposite()
Definition kstarsdata.h:168
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day,...
double epoch() const
This is (approximately) the year expressed as a floating-point value.
static KStarsDateTime currentDateTimeUtc()
QList< SkyObject * > findObjectsInArea(const SkyPoint &p1, const SkyPoint &p2)
Provides all necessary information about an object in the sky: its coordinates, name(s),...
Definition skyobject.h:42
The sky coordinates of a point in the sky.
Definition skypoint.h:45
const CachingDms & ra0() const
Definition skypoint.h:251
virtual void updateCoordsNow(const KSNumbers *num)
updateCoordsNow Shortcut for updateCoords( const KSNumbers *num, false, nullptr, nullptr,...
Definition skypoint.h:391
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
Definition skypoint.cpp:899
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition skypoint.h:94
const CachingDms & dec0() const
Definition skypoint.h:257
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition skypoint.h:119
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
Definition dms.cpp:429
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:287
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:378
const double & Degrees() const
Definition dms.h:141
QString i18n(const char *text, const TYPE &arg...)
Type type(const QSqlDatabase &db)
bool remove(const QString &column, const QVariant &value)
char * toString(const EngineQuery &query)
KIOCORE_EXPORT StatJob * stat(const QUrl &url, JobFlags flags=DefaultFlags)
KIOCORE_EXPORT QString number(KIO::filesize_t size)
KGUIADDONS_EXPORT QWindow * window(QObject *job)
void error(QWidget *parent, const QString &text, const QString &title, const KGuiItem &buttonOk, Options options=Notify)
KIOCORE_EXPORT QStringList list(const QString &fileClass)
QString name(StandardAction id)
KGuiItem properties()
QString label(StandardShortcut id)
const QList< QKeySequence > & end()
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
QByteArray & append(QByteArrayView data)
const_iterator constBegin() const const
const char * constData() const const
char * data()
bool isEmpty() const const
QByteArray leftJustified(qsizetype width, char fill, bool truncate) const const
void push_back(QByteArrayView str)
qsizetype size() const const
QDateTime currentDateTime()
QDate date() const const
QDateTime fromString(QStringView string, QStringView format, QCalendar cal)
bool isValid() const const
QTime time() const const
QString filePath(const QString &fileName) const const
QString path() const const
QString tempPath()
virtual qint64 size() const const override
virtual void close() override
qint64 size() const const
QString suffix() const const
bool isRunning() const const
T result() const const
void waitForFinished()
Format_Grayscale8
bool allGray() const const
uchar * bits()
void fill(Qt::GlobalColor color)
Format format() const const
int height() const const
bool load(QIODevice *device, const char *format)
bool loadFromData(QByteArrayView data, const char *format)
bool save(QIODevice *device, const char *format, int quality) const const
uchar * scanLine(int i)
void setColor(int index, QRgb colorValue)
void setColorCount(int colorCount)
int width() const const
QList< QByteArray > supportedImageFormats()
QString errorString() const const
QByteArray readAll()
qint64 write(const QByteArray &data)
void append(QList< T > &&value)
const_reference at(qsizetype i) const const
iterator begin()
void clear()
bool contains(const AT &value) const const
qsizetype count() const const
bool empty() const const
iterator end()
iterator erase(const_iterator begin, const_iterator end)
bool isEmpty() const const
QList< T > mid(qsizetype pos, qsizetype length) const const
qsizetype size() const const
void finished(QNetworkReply *reply)
NetworkError error() const const
QMetaObject::Connection connect(const QObject *sender, PointerToMemberFunction signal, Functor functor)
void deleteLater()
void setObjectName(QAnyStringView name)
void setX(int x)
void setY(int y)
int x() const const
int y() const const
void setX(qreal x)
void setY(qreal y)
qreal x() const const
qreal y() const const
bool contains(const QPoint &point, bool proper) const const
int height() const const
bool isNull() const const
bool isValid() const const
QPoint topLeft() const const
int width() const const
void reset(T *other)
T * get() const const
bool isNull() const const
QString arg(Args &&... args) const const
bool contains(QChar ch, Qt::CaseSensitivity cs) const const
QString fromStdString(const std::string &str)
QString & insert(qsizetype position, QChar ch)
bool isEmpty() const const
QString mid(qsizetype position, qsizetype n) const const
QString & remove(QChar ch, Qt::CaseSensitivity cs)
QString & replace(QChar before, QChar after, Qt::CaseSensitivity cs)
QStringList split(QChar sep, Qt::SplitBehavior behavior, Qt::CaseSensitivity cs) const const
bool startsWith(QChar c, Qt::CaseSensitivity cs) const const
double toDouble(bool *ok) const const
int toInt(bool *ok, int base) const const
QByteArray toLatin1() const const
QByteArray toLocal8Bit() const const
QString trimmed() const const
bool contains(QLatin1StringView str, Qt::CaseSensitivity cs) const const
qsizetype indexOf(const QRegularExpression &re, qsizetype from) const const
QTextStream & center(QTextStream &stream)
QFuture< void > map(Iterator begin, Iterator end, MapFunctor &&function)
QFuture< T > run(Function function,...)
virtual QString fileName() const const override
void setFileTemplate(const QString &name)
void setSingleShot(bool singleShot)
void start()
void stop()
void timeout()
QString query(ComponentFormattingOptions options) const const
void setQuery(const QString &query, ParsingMode mode)
QUuid createUuid()
QString toString(StringFormat mode) const const
Type type() const const
bool isValid() const const
QDateTime toDateTime() const const
double toDouble(bool *ok) const const
int toInt(bool *ok) const const
QString toString() const const
T value() const const
Object to hold FITS Header records.
Definition fitsdata.h:90
This file is part of the KDE documentation.
Documentation copyright © 1996-2025 The KDE developers.
Generated on Fri Jan 3 2025 11:47:15 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.