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 = 1, nreject = 0;
2846 for(auto &fitsKeyword : m_HeaderRecords)
2847 {
2848 // Remove keywords PC1_1... PC2_2 and CD1_1 ... CD2_2 as these are rotational / scaling matrices used by other applications
2849 // (other than KStars) that cause WCS issues with the KStars use of CROTA1/2
2850 // Also PV1_n and PS1_n
2851 // See https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf WCS section for more details.
2852 if (fitsKeyword.key.startsWith("PC1_") || fitsKeyword.key.startsWith("PC2_") ||
2853 fitsKeyword.key.startsWith("CD1_") || fitsKeyword.key.startsWith("CD2_") ||
2854 fitsKeyword.key.startsWith("PV1_") || fitsKeyword.key.startsWith("PV2_") ||
2855 fitsKeyword.key.startsWith("PS1_") || fitsKeyword.key.startsWith("PS2_"))
2856 continue;
2857 QByteArray rec;
2858 rec.append(fitsKeyword.key.leftJustified(8, ' ').toLatin1());
2859 rec.append("= ");
2860 rec.append(fitsKeyword.value.toByteArray());
2861 rec.append(" / ");
2862 rec.append(fitsKeyword.comment.toLatin1());
2863 header_str.append(rec.leftJustified(80, ' ', true));
2864 nkeyrec++;
2865 }
2866 header_str.append(QByteArray("END").leftJustified(80));
2867
2868 if ((status = wcspih(header_str.data(), nkeyrec, WCSHDR_all, 0, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2869 {
2870 wcsvfree(&m_nwcs, &m_WCSHandle);
2871 m_WCSHandle = nullptr;
2872 m_nwcs = 0;
2873 m_LastError = QString("wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]);
2874 m_WCSState = Failure;
2875 return false;
2876 }
2877
2878 if (m_WCSHandle == nullptr)
2879 {
2880 m_LastError = i18n("No world coordinate systems found.");
2881 m_WCSState = Failure;
2882 return false;
2883 }
2884
2885 // FIXME: Call above goes through EVEN if no WCS is present, so we're adding this to return for now.
2886 if (m_WCSHandle->crpix[0] == 0)
2887 {
2888 wcsvfree(&m_nwcs, &m_WCSHandle);
2889 m_WCSHandle = nullptr;
2890 m_nwcs = 0;
2891 m_LastError = i18n("No world coordinate systems found.");
2892 m_WCSState = Failure;
2893 return false;
2894 }
2895
2896 cdfix(m_WCSHandle);
2897 if ((status = wcsset(m_WCSHandle)) != 0)
2898 {
2899 wcsvfree(&m_nwcs, &m_WCSHandle);
2900 m_WCSHandle = nullptr;
2901 m_nwcs = 0;
2902 m_LastError = QString("wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2903 m_WCSState = Failure;
2904 return false;
2905 }
2906
2907 m_ObjectsSearched = false;
2908 m_CatObjectsSearched = false;
2909 m_WCSState = Success;
2910 HasWCS = true;
2911
2912 qCDebug(KSTARS_FITS) << "Finished WCS Data processing...";
2913
2914 return true;
2915#else
2916 return false;
2917#endif
2918}
2919
2920bool FITSData::wcsToPixel(const SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint)
2921{
2922#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2923 int status, stat[NWCSFIX];
2924 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2925
2926 if (m_WCSHandle == nullptr)
2927 {
2928 m_LastError = i18n("No world coordinate systems found.");
2929 return false;
2930 }
2931
2932 worldcrd[0] = wcsCoord.ra0().Degrees();
2933 worldcrd[1] = wcsCoord.dec0().Degrees();
2934
2935 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2936 {
2937 m_LastError = QString("wcss2p error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2938 return false;
2939 }
2940
2941 wcsImagePoint.setX(imgcrd[0]);
2942 wcsImagePoint.setY(imgcrd[1]);
2943
2944 wcsPixelPoint.setX(pixcrd[0]);
2945 wcsPixelPoint.setY(pixcrd[1]);
2946
2947 return true;
2948#else
2949 Q_UNUSED(wcsCoord);
2950 Q_UNUSED(wcsPixelPoint);
2951 Q_UNUSED(wcsImagePoint);
2952 return false;
2953#endif
2954}
2955
2956bool FITSData::pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord)
2957{
2958#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2959 int status, stat[NWCSFIX];
2960 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2961
2962 if (m_WCSHandle == nullptr)
2963 {
2964 m_LastError = i18n("No world coordinate systems found.");
2965 return false;
2966 }
2967
2968 pixcrd[0] = wcsPixelPoint.x();
2969 pixcrd[1] = wcsPixelPoint.y();
2970
2971 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2972 {
2973 m_LastError = QString("wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2974 return false;
2975 }
2976 else
2977 {
2978 wcsCoord.setRA0(world[0] / 15.0);
2979 wcsCoord.setDec0(world[1]);
2980 }
2981
2982 return true;
2983#else
2984 Q_UNUSED(wcsPixelPoint);
2985 Q_UNUSED(wcsCoord);
2986 return false;
2987#endif
2988}
2989
2990#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2991bool FITSData::searchObjects()
2992{
2993 return (Options::fitsCatalog() == CAT_SKYMAP) ? searchSkyMapObjects() : searchCatObjects();
2994}
2995
2996bool FITSData::searchSkyMapObjects()
2997{
2998 if (m_ObjectsSearched)
2999 return true;
3000
3001 m_ObjectsSearched = true;
3002
3003 SkyPoint startPoint;
3004 SkyPoint endPoint;
3005
3006 pixelToWCS(QPointF(0, 0), startPoint);
3007 pixelToWCS(QPointF(width() - 1, height() - 1), endPoint);
3008
3009 return findObjectsInImage(startPoint, endPoint);
3010}
3011
3012bool FITSData::searchCatObjects()
3013{
3014 if (m_CatObjectsSearched)
3015 return true;
3016
3017 m_CatObjectsSearched = true;
3018
3019 SkyPoint searchCenter;
3020 double radius;
3021 QPoint pt;
3022 bool ok = true;
3023 if (catROIRadius() > 0)
3024 {
3025 // A ROI has been set so this is a request where the user set a ROI
3026 pt = catROIPt();
3027 QPoint edgePt = QPoint(pt.x() + catROIRadius(), pt.y());
3028 SkyPoint searchEdge;
3029 ok = pixelToWCS(pt, searchCenter);
3030 if (ok)
3031 ok = pixelToWCS(edgePt, searchEdge);
3032 if (ok)
3033 {
3034 QVariant date;
3035 KSNumbers * num = nullptr;
3036
3037 if (getRecordValue("DATE-OBS", date))
3038 {
3039 QString tsString(date.toString());
3040 tsString = tsString.remove('\'').trimmed();
3041 // Add Zulu time to indicate UTC
3042 tsString += "Z";
3043
3045
3046 if (ts.isValid())
3047 num = new KSNumbers(KStarsDateTime(ts).djd());
3048 }
3049
3050 //Set to current time if the above does not work.
3051 if (num == nullptr)
3052 num = new KSNumbers(KStarsData::Instance()->ut().djd());
3053
3054 searchCenter.updateCoordsNow(num);
3055 searchEdge.updateCoordsNow(num);
3056 radius = searchCenter.angularDistanceTo(&searchEdge).Degrees() * 60; // Arc minutes
3057
3058 delete num;
3059 }
3060 }
3061 else
3062 {
3063 // No ROI setup so this is a first call to use a circle of 0.5 arcmins in the center of the image
3064 // Lets calculate the number of pixels that correspond to the 0.5 arcmin radius
3065 pt = QPoint((width() / 2) - 1, (height() / 2) - 1);
3066 ok = pixelToWCS(pt, searchCenter);
3067 if (ok)
3068 {
3069 radius = 1.0; // Use 1.0 arcmins as a starting point... not too big to avoid getting swamped with objects
3070 double raEdge = searchCenter.ra0().Degrees() + (radius / 60.0);
3071 double decEdge = searchCenter.dec0().Degrees();
3072 SkyPoint searchEdge(raEdge / 15.0, decEdge);
3073 QPointF edgePoint, pEdge;
3074 ok = wcsToPixel(searchEdge, pEdge, edgePoint);
3075 if (ok)
3076 {
3077 // Set the ROI which will be drawn
3078 const double radiusPix = std::hypot((pt.x() - pEdge.x()), pt.y() - pEdge.y());
3079 setCatSearchROI(pt, radiusPix);
3080 }
3081 }
3082 }
3083 if (!ok)
3084 {
3085 qCDebug(KSTARS_FITS) << "Unable to process Catalog Object request...";
3086 return false;
3087 }
3088 if (Options::fitsCatalog() == CAT_SIMBAD)
3089 return findSimbadObjectsInImage(searchCenter, radius);
3090 else
3091 return false;
3092}
3093
3094bool FITSData::findWCSBounds(double &minRA, double &maxRA, double &minDec, double &maxDec)
3095{
3096 if (m_WCSHandle == nullptr)
3097 {
3098 m_LastError = i18n("No world coordinate systems found.");
3099 return false;
3100 }
3101
3102 maxRA = -1000;
3103 minRA = 1000;
3104 maxDec = -1000;
3105 minDec = 1000;
3106
3107 auto updateMinMax = [&](int x, int y)
3108 {
3109 int stat[2];
3110 double imgcrd[2], phi, pixcrd[2], theta, world[2];
3111
3112 pixcrd[0] = x;
3113 pixcrd[1] = y;
3114
3115 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
3116 return;
3117
3118 minRA = std::min(minRA, world[0]);
3119 maxRA = std::max(maxRA, world[0]);
3120 minDec = std::min(minDec, world[1]);
3121 maxDec = std::max(maxDec, world[1]);
3122 };
3123
3124 // Find min and max values from edges
3125 for (int y = 0; y < height(); y++)
3126 {
3127 updateMinMax(0, y);
3128 updateMinMax(width() - 1, y);
3129 }
3130
3131 for (int x = 1; x < width() - 1; x++)
3132 {
3133 updateMinMax(x, 0);
3134 updateMinMax(x, height() - 1);
3135 }
3136
3137 // Check if either pole is in the image
3138 SkyPoint NCP(0, 90);
3139 SkyPoint SCP(0, -90);
3140 QPointF imagePoint, pPoint;
3141 if (wcsToPixel(NCP, pPoint, imagePoint))
3142 {
3143 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
3144 {
3145 maxDec = 90;
3146 }
3147 }
3148 if (wcsToPixel(SCP, pPoint, imagePoint))
3149 {
3150 if (pPoint.x() > 0 && pPoint.x() < width() && pPoint.y() > 0 && pPoint.y() < height())
3151 {
3152 minDec = -90;
3153 }
3154 }
3155 return true;
3156}
3157#endif
3158
3159#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3160bool FITSData::findObjectsInImage(SkyPoint startPoint, SkyPoint endPoint)
3161{
3162 if (KStarsData::Instance() == nullptr)
3163 return false;
3164
3165 int w = width();
3166 int h = height();
3167 QVariant date;
3168 KSNumbers * num = nullptr;
3169
3170 if (getRecordValue("DATE-OBS", date))
3171 {
3172 QString tsString(date.toString());
3173 tsString = tsString.remove('\'').trimmed();
3174 // Add Zulu time to indicate UTC
3175 tsString += "Z";
3176
3178
3179 if (ts.isValid())
3180 num = new KSNumbers(KStarsDateTime(ts).djd());
3181 }
3182
3183 //Set to current time if the above does not work.
3184 if (num == nullptr)
3185 num = new KSNumbers(KStarsData::Instance()->ut().djd());
3186
3187 startPoint.updateCoordsNow(num);
3188 endPoint.updateCoordsNow(num);
3189
3190 m_SkyObjects.clear();
3191
3192 QList<SkyObject *> list = KStarsData::Instance()->skyComposite()->findObjectsInArea(startPoint, endPoint);
3193 list.erase(std::remove_if(list.begin(), list.end(), [](SkyObject * oneObject)
3194 {
3195 int type = oneObject->type();
3196 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
3197 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
3198 type == SkyObject::SATELLITE);
3199 }), list.end());
3200
3201 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3202 int stat[2];
3203 for (auto &object : list)
3204 {
3205 world[0] = object->ra0().Degrees();
3206 world[1] = object->dec0().Degrees();
3207
3208 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
3209 {
3210 //The X and Y are set to the found position if it does work.
3211 int x = pixcrd[0];
3212 int y = pixcrd[1];
3213 if (x > 0 && y > 0 && x < w && y < h)
3214 m_SkyObjects.append(new FITSSkyObject(object, x, y));
3215 }
3216 }
3217
3218 delete (num);
3219 return true;
3220}
3221
3222// See https://simbad.u-strasbg.fr/Pages/guide/sim-q.htx and
3223// https://simbad.u-strasbg.fr/Pages/guide/sim-url.htx
3224// for details of how to query Simbad and the various options
3225//
3226// Set the epoch from the datetime of the image and ensure that the reply is in this same epoch
3227bool FITSData::findSimbadObjectsInImage(SkyPoint searchCenter, double radius)
3228{
3229 m_CatObjects.clear();
3230 m_CatUpdateTable = false;
3231
3232 // Query Simbad
3233 QUrl simbadURL = QUrl("https://simbad.cds.unistra.fr/simbad/sim-coo");
3234 QUrlQuery simbadQuery(simbadURL.query());
3235
3236 QString coord = QString("%1 %2").arg(searchCenter.ra0().toHMSString(true, true))
3237 .arg(searchCenter.dec0().toDMSString(true, true, true));
3238
3239 QString radiusStr = QString("%1").arg(radius * 60.0, 0, 'f', 5);
3240 m_CatObjQuery = QString("%1 %2").arg(coord).arg(radiusStr);
3241
3242 coord.replace("+", "%2B"); // Need to replace the + otherwise it gets lost by URL processing
3243
3244 QString epoch = QString("J%1").arg(m_DateTime.epoch());
3245
3246 simbadQuery.addQueryItem("Coord", coord);
3247 simbadQuery.addQueryItem("Radius", radiusStr);
3248 simbadQuery.addQueryItem("Radius.unit", "arcsec");
3249 simbadQuery.addQueryItem("CooFrame", "ICRS");
3250 simbadQuery.addQueryItem("CooEpoch", epoch);
3251 simbadQuery.addQueryItem("output.format", "ASCII");
3252 simbadQuery.addQueryItem("output.max", QString("%1").arg(10000));
3253 simbadQuery.addQueryItem("list.otypesel", "on");
3254 simbadQuery.addQueryItem("otypedisp", "3"); // Use Simbad's 3 char object type
3255 simbadQuery.addQueryItem("list.coo1", "on"); // Display coord 1
3256 simbadQuery.addQueryItem("frame1", "ICRS");
3257 simbadQuery.addQueryItem("epoch1", "J2000");
3258 simbadQuery.addQueryItem("equi1", epoch);
3259 simbadQuery.addQueryItem("list.spsel", "off"); // Spectral Type
3260 simbadQuery.addQueryItem("list.sizesel", "on"); // Angular size
3261 simbadQuery.addQueryItem("list.bibsel", "off"); // Bibliography
3262
3263 simbadURL.setQuery(simbadQuery);
3264 qCDebug(KSTARS_FITS) << "Simbad query:" << simbadURL;
3265
3266 m_NetworkAccessManager.reset(new QNetworkAccessManager(this));
3267 connect(m_NetworkAccessManager.get(), &QNetworkAccessManager::finished, this, &FITSData::simbadResponseReady);
3268 QNetworkReply *response = m_NetworkAccessManager->get(QNetworkRequest(simbadURL));
3269 if (response->error() != QNetworkReply::NoError)
3270 {
3271 qCDebug(KSTARS_FITS) << "Error:" << response->errorString() << " occured querying SIMBAD with" << simbadURL;
3272 m_CatQueryInProgress = false;
3273 emit catalogQueryFailed(i18n("Error querying Simbad"));
3274 return false;
3275 }
3276
3277 m_CatQueryInProgress = true;
3278 m_CatQueryTimer.setSingleShot(true);
3279 connect(&m_CatQueryTimer, &QTimer::timeout, this, &FITSData::catTimeout);
3280 m_CatQueryTimer.start(60 * 1000);
3281 emit loadingCatalogData();
3282 return true;
3283}
3284#endif
3285
3286#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3287void FITSData::catTimeout()
3288{
3289 m_CatQueryInProgress = false;
3290 QString text = i18n("Simbad query timed out");
3291 qCDebug(KSTARS_FITS) << text;
3292 emit catalogQueryFailed(text);
3293}
3294#endif
3295
3296// There are 3 types of responses... zero objects, 1 object and >1 objects. The formats are as follows:
3297// (Note that Line xx isn't part of the reply, these are added for convenience)
3298//
3299// Zero Objects
3300// Line 0 "!! No astronomical object found : "
3301//
3302// 1 Object
3303// Line 0 "C.D.S. - SIMBAD4 rel 1.8 - 2024.09.30CEST14:07:24"
3304// Line 1 ""
3305// Line 2 "coord 21:36:49 +57:31:23 (ICRS, J2024.7, 2000.0), radius: 5.71083 arcsec"
3306// Line 3 "------------------------------------------------------------------------"
3307// Line 4 ""
3308// Line 5 "Object EM* LkHA 349C --- RS* --- OID=@157406 (@@122902,11) --- coobox=594"
3309// Line 6 ""
3310// 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"
3311// Line 8 "Coordinates(FK4,ep=B1950,eq=1950): 21 35 17.0484337158 +57 17 51.458438447"
3312// Line 9 "Coordinates(Gal,ep=J2000,eq=2000): 099.0986092951232 +03.9548388170815"
3313// Line 10 "hierarchy counts: #parents=2, #children=0, #siblings=0"
3314// Line 11 "Proper motions: -2.232 -5.540 [0.087 0.086 90] A 2020yCat.1350....0G"
3315// Line 12 "Parallax: 1.0825 [0.0718] A 2020yCat.1350....0G"
3316// Line 13 "Radial Velocity: ~ [~ ~ ] ~ ~ "
3317// Line 14 "Redshift: ~ [~ ~ ] ~ ~ "
3318// Line 15 "cz: ~ [~ ~ ] ~ ~ "
3319// Line 16 "Flux U : 19.219 [0.083] C 2010ApJ...710..597S"
3320// Line 17 "Flux V : 15.81 [0.01] B 2012MNRAS.426.2917G"
3321// Line 18 "Flux G : 14.948605 [0.007309] C 2020yCat.1350....0G"
3322// Line 19 "Flux R : 14.508 [0.030] C 2010ApJ...710..597S"
3323// Line 20 "Flux I : 13.094 [0.003] C 2010ApJ...710..597S"
3324// Line 21 "Flux J : 11.918 [0.022] C 2003yCat.2246....0C"
3325// Line 22 "Flux H : 10.914 [0.029] C 2003yCat.2246....0C"
3326// Line 23 "Flux K : 10.355 [0.021] C 2003yCat.2246....0C"
3327// Line 24 "Flux g : 16.702 [~] D 2020ApJS..249...18C"
3328// Line 25 "Flux r : 15.057 [~] D 2020ApJS..249...18C"
3329// Line 26 "Flux i : 14.84 [~] D 2012AJ....143...61N"
3330// Line 27 "Spectral type: K6 D 2006AJ....132.2135S"
3331// Line 28 "Morphological type: ~ ~ ~"
3332// Line 29 "Angular size: ~ ~ ~ (~) ~ ~"
3333// Line 30 ""
3334// Line 31 "Identifiers (12):"
3335// Line 32 " ZTF J213649.40+573121.9 Gaia DR3 2178442322028512512 TIC 469533750 "
3336// Line 33 " CoKu LkHA 349 c EM* LkHA 349C 2MASS J21364941+5731220 "
3337// Line 34 " [SHB2004] Trumpler 37 14-141 [MSR2009] IC1396A-29 [NSW2012] 125 "
3338// Line 35 " CXOIC1396A J213649.39+573122.1 [GFS2012] 173 Gaia DR2 2178442322028512512 "
3339// Line 36 ""
3340// Line 37 "Bibcodes 1850-2024 () (20):"
3341// Line 38 " 2024MNRAS.532.2108S 2020ApJS..249...18C 2019ApJ...878....7M 2018MNRAS.478.5091F"
3342// Line 39 " 2012AJ....143...61N 2012MNRAS.426.2917G 2011ApJ...742...39S 2010ApJ...710..597S"
3343// Line 40 " 2009ApJ...690..683R 2009ApJ...702.1507M 2006AJ....132.2135S 2006ApJ...638..897S"
3344// Line 41 " 2005A&A...443..535V 2005AJ....130..188S 2004AJ....128..805S 2004ApJS..154..385R"
3345// Line 42 " 1997A&A...325.1001S 1996ApJ...463L.105M 1995A&A...299..464H 1979ApJS...41..743C"
3346// Line 43 ""
3347// Line 44 "Measures (distance:1 PLX:2 PM:2 ROT:1 V*:3 velocities:1 ):"
3348// Line 45 "distance:1PLX:2PM:2ROT:1V*:3velocities:1"
3349// Line 46 ""
3350// Line 47 "Notes (0) :"
3351// Line 48 ""
3352// Line 49 "================================================================================"
3353// Line 50 ""
3354// Line 51 ""
3355//
3356// >1 Objects
3357// Line 0 "C.D.S. - SIMBAD4 rel 1.8 - 2024.09.30CEST14:13:23"
3358// Line 1 ""
3359// Line 2 "coord 21:36:51 +57:31:13 (ICRS, J2024.7, 2000.0), radius: 11.42166 arcsec"
3360// Line 3 "-------------------------------------------------------------------------"
3361// Line 4 ""
3362// Line 5 "Number of objects : 2"
3363// Line 6 ""
3364// Line 7 "#|dist(asec)| identifier |typ| coord1 (ICRS,J2000/J2024.7) |Mag U | Mag B | Mag V | Mag R | Mag I | ang. size |#not"
3365// Line 8 "-|----------|-----------------------------------|---|---------------------------------------|------|---------|---------|---------|---------|---------------|----"
3366// Line 9 "1| 2.62|DOBASHI 3187 |DNe|21 36 51.3 +57 31 12 | ~| ~ | ~ | ~ | ~ | ~ ~ | 1"
3367// 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"
3368// Line 11 "================================================================================"
3369// Line 12 ""
3370// Line 13 ""
3371
3372#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3373void FITSData::simbadResponseReady(QNetworkReply * simbadReply)
3374{
3375 m_CatQueryTimer.stop();
3376 m_CatQueryInProgress = false;
3377 if (simbadReply->error() != QNetworkReply::NoError)
3378 {
3379 qCDebug(KSTARS_FITS) << "Error:" << simbadReply->errorString() << " occured in SIMBAD query reply";
3380 m_CatQueryInProgress = false;
3381 emit catalogQueryFailed(i18n("Error querying Simbad"));
3382 return;
3383 }
3384
3385 auto data = simbadReply->readAll();
3386 QString dataStr = QString(data);
3387
3388 // Break the line into comma-separated components
3389 int numObjs = 0;
3390 bool ok;
3391 QString name, type, coord, sizeStr;
3392 double dist = 0, magnitude = 0;
3393 int num = 0;
3394 QStringList replyLines = dataStr.split(QLatin1Char('\n'));
3395 int dataRow = -1;
3396 for (int i = 0; i < replyLines.size(); i++)
3397 {
3398 // JEE qCDebug(KSTARS_FITS) << "Line" << i << replyLines[i];
3399 if (i == 0 && replyLines[i].contains("No astronomical object found : "))
3400 // No objects in the reply so we can just ignore the rest of the data
3401 break;
3402 if (i == 2)
3403 {
3404 // Check that the query coord & radius match the reply - otherwise ignore...
3405 QString replyStr;
3406 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3407 if (replyData.size() == 9 && replyData[0] == "coord" && replyData[6] == "radius:")
3408 replyStr = QString("%1 %2 %3").arg(replyData[1]).arg(replyData[2]).arg(replyData[7]);
3409 if (replyStr != m_CatObjQuery)
3410 {
3411 qCDebug(KSTARS_FITS) << "Simbad query:" << m_CatObjQuery << "Reply:" << replyStr << ". Ignoring...";
3412 break;
3413 }
3414 }
3415 else if (i == 5)
3416 {
3417 if (replyLines[i].contains("Number of objects : "))
3418 numObjs = QString(replyLines[i].mid(20)).toInt(&ok);
3419 else if (replyLines[i].startsWith("Object "))
3420 {
3421 // Dealing with a single object
3422 num = numObjs = 1;
3423 const int firstBreak = replyLines[i].indexOf("---");
3424 const int secondBreak = replyLines[i].indexOf("---", firstBreak + 3);
3425 name = replyLines[i].mid(7, firstBreak - 7).trimmed();
3426 type = replyLines[i].mid(firstBreak + 3, secondBreak - firstBreak - 3).trimmed();
3427 }
3428 else
3429 {
3430 qCDebug(KSTARS_FITS) << "Bad Simbad Reply, Ignoring...";
3431 break;
3432 }
3433 }
3434 else if (numObjs == 1 && i >= 7)
3435 {
3436 if (replyLines[i].startsWith("Coordinates(ICRS"))
3437 {
3438 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3439 if (replyData.size() >= 8)
3440 coord = replyData[1] + " " + replyData[2] + " " + replyData[3] + " " +
3441 replyData[5] + " " + replyData[6] + " " + replyData[7];
3442 }
3443 else if (replyLines[i].startsWith("Flux V :"))
3444 {
3445 QStringList replyData = replyLines[i].split(QLatin1Char(' '));
3446 if (replyData.size() >= 4)
3447 magnitude = QString(replyData[3]).toDouble(&ok);
3448 }
3449 else if (replyLines[i].startsWith("Angular size:"))
3450 sizeStr = replyLines[i].mid(13, replyLines[i].indexOf("~") - 13).trimmed();
3451 else if (replyLines[i].startsWith("=========================================="))
3452 {
3453 // End of data so add the object
3454 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3455 dataRow++;
3456 }
3457 }
3458 else if (numObjs > 1 && i >= 9)
3459 {
3460 // These are the data lines - 1 per object
3461 QStringList replyData = replyLines[i].split(QLatin1Char('|'));
3462 num = 0;
3463 name = type = coord = sizeStr = "";
3464 magnitude = dist = 0.0;
3465 for (int j = 0; j < replyData.size(); j++)
3466 {
3467 if (j == 0)
3468 num = QString(replyData[j]).toInt(&ok);
3469 else if (j == 1)
3470 dist = QString(replyData[j]).toDouble(&ok) / 60.0; // Convert from arcsecs to arcmins
3471 else if (j == 2)
3472 name = (QString(replyData[j]).trimmed());
3473 else if (j == 3)
3474 type = QString(replyData[j]).trimmed();
3475 else if (j == 4)
3476 coord = QString(replyData[j]).trimmed();
3477 else if (j == 7) // Mag V
3478 magnitude = QString(replyData[j]).toDouble(&ok);
3479 else if (j == 10) // Angular size
3480 sizeStr = QString(replyData[j]).trimmed();
3481 }
3482
3483 if (num == 0 || name.isEmpty() || type.isEmpty() || coord.isEmpty())
3484 continue;
3485
3486 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3487 dataRow++;
3488 }
3489 }
3490 if (numObjs != ++dataRow)
3491 qCDebug(KSTARS_FITS) << "Simbad Header rows:" << numObjs << ". Data rows:" << dataRow;
3492 // Signal fitsview to update based on new data. Fitsview will signal fitstab
3493 m_CatObjectsSearched = false;
3494 // Signal FITSView that there is new data to process
3495 emit dataChanged();
3496 // Signal FITSTab that there is new data to process
3497 emit loadedCatalogData();
3498 simbadReply->deleteLater();
3499}
3500#endif
3501
3502bool FITSData::addCatObject(const int num, const QString name, const QString type, const QString coord, const double dist,
3503 const double magnitude, const QString sizeStr)
3504{
3505 bool ok;
3506 dms r(0.0), d(0.0);
3507
3508 QStringList split = coord.split(QLatin1Char(' '));
3509 if (split.size() != 6)
3510 {
3511 qCDebug(KSTARS_FITS) << "Coordinates for " << name << "invalid: " << coord;
3512 return false;
3513 }
3514
3515 QString raStr = QString("%1 %2 %3").arg(split[0]).arg(split[1]).arg(split[2]);
3516 QString decStr = QString("%1 %2 %3").arg(split[3]).arg(split[4]).arg(split[5]);
3517 ok = r.setFromString(raStr, false); // angle in hours
3518 if (ok)
3519 ok = d.setFromString(decStr, true); // angle in degrees
3520 if (!ok)
3521 {
3522 qCDebug(KSTARS_FITS) << "Coordinates for " << name << "invalid: " << coord;
3523 return false;
3524 }
3525
3526 double size = 0.0;
3527 if (!sizeStr.contains("~"))
3528 {
3529 QStringList sizeSplit = sizeStr.split(QLatin1Char(' '));
3530 if (sizeSplit.size() >= 1)
3531 // Take the first (major axis) size. Note sometimes there's a double space between axes so
3532 // number of elements can be 3.
3533 size = QString(sizeSplit[0]).toDouble(&ok);
3534 else
3535 qCDebug(KSTARS_FITS) << "Angular size for " << name << "invalid: " << sizeStr;
3536 }
3537 CatObject catObject;
3538 catObject.num = num;
3539 catObject.catType = CAT_SIMBAD;
3540 catObject.name = name;
3541 catObject.typeCode = type;
3542 catObject.typeLabel = getCatObjectLabel(type);
3543 catObject.r = r;
3544 catObject.d = d;
3545 catObject.dist = dist;
3546 catObject.magnitude = magnitude;
3547 catObject.size = size;
3548 catObject.highlight = false;
3549 catObject.show = getCatObjectFilter(type);
3550
3551 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3552 int stat[2];
3553
3554 world[0] = r.Degrees();
3555 world[1] = d.Degrees();
3556
3557 int status = wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]);
3558 if (status != 0)
3559 {
3560 qCDebug(KSTARS_FITS) << QString("wcss2p error processing object %1 %2: %3.").arg(name).arg(status).arg(wcs_errmsg[status]);
3561 return false;
3562 }
3563
3564 //X and Y are set to the found position if in image
3565 double x = pixcrd[0];
3566 double y = pixcrd[1];
3567 if (x > 0 && y > 0 && x < width() && y < height())
3568 {
3569 catObject.x = x;
3570 catObject.y = y;
3571 m_CatObjects.append(catObject);
3572 }
3573 return true;
3574}
3575
3576void FITSData::setCatObjectsFilters(const QStringList filters)
3577{
3578 m_CatObjectsFilters = filters;
3579}
3580
3581#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3582void FITSData::setCatSearchROI(const QPoint searchCenter, const int radius)
3583{
3584 m_CatROIPt = searchCenter;
3585 m_CatROIRadius = radius;
3586 Q_UNUSED(searchCenter);
3587 Q_UNUSED(radius);
3588}
3589#endif
3590
3591#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3592QString FITSData::getCatObjectLabel(const QString code) const
3593{
3594 QString label;
3595 if (code != "")
3596 {
3597 for (int i = 0; i < MAX_CAT_OBJ_TYPES; i++)
3598 {
3599 if (code == catObjTypes[i].code)
3600 {
3601 label = catObjTypes[i].label;
3602 break;
3603 }
3604 else if (code == catObjTypes[i].candidateCode)
3605 {
3606 label = catObjTypes[i].label + "_Cand";
3607 break;
3608 }
3609 }
3610 }
3611 return label;
3612}
3613#endif
3614
3615#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3616bool FITSData::getCatObjectFilter(const QString type) const
3617{
3618 if (m_CatObjectsFilters.isEmpty() || m_CatObjectsFilters.contains(type))
3619 return true;
3620 return false;
3621}
3622#endif
3623
3624void FITSData::filterCatObjects()
3625{
3626#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3627 bool changed = false;
3628 for (int i = 0; i < m_CatObjects.size(); i++)
3629 {
3630 bool showState = getCatObjectFilter(m_CatObjects[i].typeCode);
3631 if (m_CatObjects[i].show != showState)
3632 {
3633 changed = true;
3634 m_CatObjects[i].show = showState;
3635 }
3636 }
3637 if (changed)
3638 emit dataChanged();
3639#endif
3640}
3641
3642bool FITSData::highlightCatObject(const int hilite, const int lolite)
3643{
3644#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3645 if (hilite < 0 || hilite >= m_CatObjects.size())
3646 return false;
3647
3648 // Lowlight previous highlighted catalog object if passed
3649 if (lolite >= 0 && lolite < m_CatObjects.size())
3650 m_CatObjects[lolite].highlight = false;
3651 else
3652 {
3653 // Loop through data to find current highlighted item
3654 for (int i = 0; i < m_CatObjects.size(); i++)
3655 {
3656 if (m_CatObjects[i].highlight)
3657 {
3658 m_CatObjects[i].highlight = false;
3659 break;
3660 }
3661 }
3662 }
3663
3664 // Highlight the catalog object
3665 m_CatObjects[hilite].highlight = true;
3666#endif
3667 return true;
3668}
3669
3670int FITSData::getFlipVCounter() const
3671{
3672 return flipVCounter;
3673}
3674
3675void FITSData::setFlipVCounter(int value)
3676{
3677 flipVCounter = value;
3678}
3679
3680int FITSData::getFlipHCounter() const
3681{
3682 return flipHCounter;
3683}
3684
3685void FITSData::setFlipHCounter(int value)
3686{
3687 flipHCounter = value;
3688}
3689
3690int FITSData::getRotCounter() const
3691{
3692 return rotCounter;
3693}
3694
3695void FITSData::setRotCounter(int value)
3696{
3697 rotCounter = value;
3698}
3699
3700/* Rotate an image by 90, 180, or 270 degrees, with an optional
3701 * reflection across the vertical or horizontal axis.
3702 * verbose generates extra info on stdout.
3703 * return nullptr if successful or rotated image.
3704 */
3705template <typename T>
3706bool FITSData::rotFITS(int rotate, int mirror)
3707{
3708 int ny, nx;
3709 int x1, y1, x2, y2;
3710 uint8_t * rotimage = nullptr;
3711 int offset = 0;
3712
3713 if (rotate == 1)
3714 rotate = 90;
3715 else if (rotate == 2)
3716 rotate = 180;
3717 else if (rotate == 3)
3718 rotate = 270;
3719 else if (rotate < 0)
3720 rotate = rotate + 360;
3721
3722 nx = m_Statistics.width;
3723 ny = m_Statistics.height;
3724
3725 int BBP = m_Statistics.bytesPerPixel;
3726
3727 /* Allocate buffer for rotated image */
3728 rotimage = new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
3729
3730 if (rotimage == nullptr)
3731 {
3732 qWarning() << "Unable to allocate memory for rotated image buffer!";
3733 return false;
3734 }
3735
3736 auto * rotBuffer = reinterpret_cast<T *>(rotimage);
3737 auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
3738
3739 /* Mirror image without rotation */
3740 if (rotate < 45 && rotate > -45)
3741 {
3742 if (mirror == 1)
3743 {
3744 for (int i = 0; i < m_Statistics.channels; i++)
3745 {
3746 offset = m_Statistics.samples_per_channel * i;
3747 for (x1 = 0; x1 < nx; x1++)
3748 {
3749 x2 = nx - x1 - 1;
3750 for (y1 = 0; y1 < ny; y1++)
3751 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3752 }
3753 }
3754 }
3755 else if (mirror == 2)
3756 {
3757 for (int i = 0; i < m_Statistics.channels; i++)
3758 {
3759 offset = m_Statistics.samples_per_channel * i;
3760 for (y1 = 0; y1 < ny; y1++)
3761 {
3762 y2 = ny - y1 - 1;
3763 for (x1 = 0; x1 < nx; x1++)
3764 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3765 }
3766 }
3767 }
3768 else
3769 {
3770 for (int i = 0; i < m_Statistics.channels; i++)
3771 {
3772 offset = m_Statistics.samples_per_channel * i;
3773 for (y1 = 0; y1 < ny; y1++)
3774 {
3775 for (x1 = 0; x1 < nx; x1++)
3776 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3777 }
3778 }
3779 }
3780 }
3781
3782 /* Rotate by 90 degrees */
3783 else if (rotate >= 45 && rotate < 135)
3784 {
3785 if (mirror == 1)
3786 {
3787 for (int i = 0; i < m_Statistics.channels; i++)
3788 {
3789 offset = m_Statistics.samples_per_channel * i;
3790 for (y1 = 0; y1 < ny; y1++)
3791 {
3792 x2 = ny - y1 - 1;
3793 for (x1 = 0; x1 < nx; x1++)
3794 {
3795 y2 = nx - x1 - 1;
3796 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3797 }
3798 }
3799 }
3800 }
3801 else if (mirror == 2)
3802 {
3803 for (int i = 0; i < m_Statistics.channels; i++)
3804 {
3805 offset = m_Statistics.samples_per_channel * i;
3806 for (y1 = 0; y1 < ny; y1++)
3807 {
3808 for (x1 = 0; x1 < nx; x1++)
3809 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3810 }
3811 }
3812 }
3813 else
3814 {
3815 for (int i = 0; i < m_Statistics.channels; i++)
3816 {
3817 offset = m_Statistics.samples_per_channel * i;
3818 for (y1 = 0; y1 < ny; y1++)
3819 {
3820 x2 = ny - y1 - 1;
3821 for (x1 = 0; x1 < nx; x1++)
3822 {
3823 y2 = x1;
3824 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3825 }
3826 }
3827 }
3828 }
3829
3830 m_Statistics.width = ny;
3831 m_Statistics.height = nx;
3832 }
3833
3834 /* Rotate by 180 degrees */
3835 else if (rotate >= 135 && rotate < 225)
3836 {
3837 if (mirror == 1)
3838 {
3839 for (int i = 0; i < m_Statistics.channels; i++)
3840 {
3841 offset = m_Statistics.samples_per_channel * i;
3842 for (y1 = 0; y1 < ny; y1++)
3843 {
3844 y2 = ny - y1 - 1;
3845 for (x1 = 0; x1 < nx; x1++)
3846 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3847 }
3848 }
3849 }
3850 else if (mirror == 2)
3851 {
3852 for (int i = 0; i < m_Statistics.channels; i++)
3853 {
3854 offset = m_Statistics.samples_per_channel * i;
3855 for (x1 = 0; x1 < nx; x1++)
3856 {
3857 x2 = nx - x1 - 1;
3858 for (y1 = 0; y1 < ny; y1++)
3859 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3860 }
3861 }
3862 }
3863 else
3864 {
3865 for (int i = 0; i < m_Statistics.channels; i++)
3866 {
3867 offset = m_Statistics.samples_per_channel * i;
3868 for (y1 = 0; y1 < ny; y1++)
3869 {
3870 y2 = ny - y1 - 1;
3871 for (x1 = 0; x1 < nx; x1++)
3872 {
3873 x2 = nx - x1 - 1;
3874 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3875 }
3876 }
3877 }
3878 }
3879 }
3880
3881 /* Rotate by 270 degrees */
3882 else if (rotate >= 225 && rotate < 315)
3883 {
3884 if (mirror == 1)
3885 {
3886 for (int i = 0; i < m_Statistics.channels; i++)
3887 {
3888 offset = m_Statistics.samples_per_channel * i;
3889 for (y1 = 0; y1 < ny; y1++)
3890 {
3891 for (x1 = 0; x1 < nx; x1++)
3892 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3893 }
3894 }
3895 }
3896 else if (mirror == 2)
3897 {
3898 for (int i = 0; i < m_Statistics.channels; i++)
3899 {
3900 offset = m_Statistics.samples_per_channel * i;
3901 for (y1 = 0; y1 < ny; y1++)
3902 {
3903 x2 = ny - y1 - 1;
3904 for (x1 = 0; x1 < nx; x1++)
3905 {
3906 y2 = nx - x1 - 1;
3907 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3908 }
3909 }
3910 }
3911 }
3912 else
3913 {
3914 for (int i = 0; i < m_Statistics.channels; i++)
3915 {
3916 offset = m_Statistics.samples_per_channel * i;
3917 for (y1 = 0; y1 < ny; y1++)
3918 {
3919 x2 = y1;
3920 for (x1 = 0; x1 < nx; x1++)
3921 {
3922 y2 = nx - x1 - 1;
3923 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3924 }
3925 }
3926 }
3927 }
3928
3929 m_Statistics.width = ny;
3930 m_Statistics.height = nx;
3931 }
3932
3933 /* If rotating by more than 315 degrees, assume top-bottom reflection */
3934 else if (rotate >= 315 && mirror)
3935 {
3936 for (int i = 0; i < m_Statistics.channels; i++)
3937 {
3938 offset = m_Statistics.samples_per_channel * i;
3939 for (y1 = 0; y1 < ny; y1++)
3940 {
3941 for (x1 = 0; x1 < nx; x1++)
3942 {
3943 x2 = y1;
3944 y2 = x1;
3945 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3946 }
3947 }
3948 }
3949 }
3950
3951 delete[] m_ImageBuffer;
3952 m_ImageBuffer = rotimage;
3953
3954 return true;
3955}
3956
3957void FITSData::rotWCSFITS(int angle, int mirror)
3958{
3959 int status = 0;
3960 char comment[100];
3961 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
3962 int WCS_DECIMALS = 6;
3963
3964 naxis1 = m_Statistics.width;
3965 naxis2 = m_Statistics.height;
3966
3967 if (fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
3968 {
3969 // No WCS keywords
3970 return;
3971 }
3972
3973 /* Reset CROTAn and CD matrix if axes have been exchanged */
3974 if (angle == 90)
3975 {
3976 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
3977 fits_update_key_dbl(fptr, "CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3978
3979 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
3980 fits_update_key_dbl(fptr, "CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3981 }
3982
3983 status = 0;
3984
3985 /* Negate rotation angle if mirrored */
3986 if (mirror != 0)
3987 {
3988 if (!fits_read_key_dbl(fptr, "CROTA1", &ctemp1, comment, &status))
3989 fits_update_key_dbl(fptr, "CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
3990
3991 if (!fits_read_key_dbl(fptr, "CROTA2", &ctemp1, comment, &status))
3992 fits_update_key_dbl(fptr, "CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
3993
3994 status = 0;
3995
3996 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
3997 fits_update_key_dbl(fptr, "LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
3998
3999 status = 0;
4000
4001 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
4002 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4003
4004 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp1, comment, &status))
4005 fits_update_key_dbl(fptr, "CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
4006
4007 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp1, comment, &status))
4008 fits_update_key_dbl(fptr, "CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
4009 }
4010
4011 status = 0;
4012
4013 /* Unbin CRPIX and CD matrix */
4014 if (!fits_read_key_dbl(fptr, "LTM1_1", &ctemp1, comment, &status))
4015 {
4016 if (ctemp1 != 1.0)
4017 {
4018 if (!fits_read_key_dbl(fptr, "LTM2_2", &ctemp2, comment, &status))
4019 if (ctemp1 == ctemp2)
4020 {
4021 double ltv1 = 0.0;
4022 double ltv2 = 0.0;
4023 status = 0;
4024 if (!fits_read_key_dbl(fptr, "LTV1", &ltv1, comment, &status))
4025 fits_delete_key(fptr, "LTV1", &status);
4026 if (!fits_read_key_dbl(fptr, "LTV2", &ltv2, comment, &status))
4027 fits_delete_key(fptr, "LTV2", &status);
4028
4029 status = 0;
4030
4031 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp3, comment, &status))
4032 fits_update_key_dbl(fptr, "CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
4033
4034 if (!fits_read_key_dbl(fptr, "CRPIX2", &ctemp3, comment, &status))
4035 fits_update_key_dbl(fptr, "CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
4036
4037 status = 0;
4038
4039 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp3, comment, &status))
4040 fits_update_key_dbl(fptr, "CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4041
4042 if (!fits_read_key_dbl(fptr, "CD1_2", &ctemp3, comment, &status))
4043 fits_update_key_dbl(fptr, "CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4044
4045 if (!fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status))
4046 fits_update_key_dbl(fptr, "CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4047
4048 if (!fits_read_key_dbl(fptr, "CD2_2", &ctemp3, comment, &status))
4049 fits_update_key_dbl(fptr, "CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4050
4051 status = 0;
4052
4053 fits_delete_key(fptr, "LTM1_1", &status);
4054 fits_delete_key(fptr, "LTM1_2", &status);
4055 }
4056 }
4057 }
4058
4059 status = 0;
4060
4061 /* Reset CRPIXn */
4062 if (!fits_read_key_dbl(fptr, "CRPIX1", &ctemp1, comment, &status) &&
4063 !fits_read_key_dbl(fptr, "CRPIX2", &ctemp2, comment, &status))
4064 {
4065 if (mirror != 0)
4066 {
4067 if (angle == 0)
4068 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4069 else if (angle == 90)
4070 {
4071 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4072 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4073 }
4074 else if (angle == 180)
4075 {
4076 fits_update_key_dbl(fptr, "CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
4077 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4078 }
4079 else if (angle == 270)
4080 {
4081 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4082 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4083 }
4084 }
4085 else
4086 {
4087 if (angle == 90)
4088 {
4089 fits_update_key_dbl(fptr, "CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4090 fits_update_key_dbl(fptr, "CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4091 }
4092 else if (angle == 180)
4093 {
4094 fits_update_key_dbl(fptr, "CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4095 fits_update_key_dbl(fptr, "CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4096 }
4097 else if (angle == 270)
4098 {
4099 fits_update_key_dbl(fptr, "CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4100 fits_update_key_dbl(fptr, "CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4101 }
4102 }
4103 }
4104
4105 status = 0;
4106
4107 /* Reset CDELTn (degrees per pixel) */
4108 if (!fits_read_key_dbl(fptr, "CDELT1", &ctemp1, comment, &status) &&
4109 !fits_read_key_dbl(fptr, "CDELT2", &ctemp2, comment, &status))
4110 {
4111 if (mirror != 0)
4112 {
4113 if (angle == 0)
4114 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4115 else if (angle == 90)
4116 {
4117 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4118 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4119 }
4120 else if (angle == 180)
4121 {
4122 fits_update_key_dbl(fptr, "CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
4123 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4124 }
4125 else if (angle == 270)
4126 {
4127 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4128 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4129 }
4130 }
4131 else
4132 {
4133 if (angle == 90)
4134 {
4135 fits_update_key_dbl(fptr, "CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4136 fits_update_key_dbl(fptr, "CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4137 }
4138 else if (angle == 180)
4139 {
4140 fits_update_key_dbl(fptr, "CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4141 fits_update_key_dbl(fptr, "CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4142 }
4143 else if (angle == 270)
4144 {
4145 fits_update_key_dbl(fptr, "CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4146 fits_update_key_dbl(fptr, "CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4147 }
4148 }
4149 }
4150
4151 /* Reset CD matrix, if present */
4152 ctemp1 = 0.0;
4153 ctemp2 = 0.0;
4154 ctemp3 = 0.0;
4155 ctemp4 = 0.0;
4156 status = 0;
4157 if (!fits_read_key_dbl(fptr, "CD1_1", &ctemp1, comment, &status))
4158 {
4159 fits_read_key_dbl(fptr, "CD1_2", &ctemp2, comment, &status);
4160 fits_read_key_dbl(fptr, "CD2_1", &ctemp3, comment, &status);
4161 fits_read_key_dbl(fptr, "CD2_2", &ctemp4, comment, &status);
4162 status = 0;
4163 if (mirror != 0)
4164 {
4165 if (angle == 0)
4166 {
4167 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4168 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4169 }
4170 else if (angle == 90)
4171 {
4172 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4173 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4174 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4175 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4176 }
4177 else if (angle == 180)
4178 {
4179 fits_update_key_dbl(fptr, "CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
4180 fits_update_key_dbl(fptr, "CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
4181 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4182 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4183 }
4184 else if (angle == 270)
4185 {
4186 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4187 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4188 fits_update_key_dbl(fptr, "CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
4189 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4190 }
4191 }
4192 else
4193 {
4194 if (angle == 90)
4195 {
4196 fits_update_key_dbl(fptr, "CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4197 fits_update_key_dbl(fptr, "CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4198 fits_update_key_dbl(fptr, "CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
4199 fits_update_key_dbl(fptr, "CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4200 }
4201 else if (angle == 180)
4202 {
4203 fits_update_key_dbl(fptr, "CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4204 fits_update_key_dbl(fptr, "CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4205 fits_update_key_dbl(fptr, "CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4206 fits_update_key_dbl(fptr, "CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4207 }
4208 else if (angle == 270)
4209 {
4210 fits_update_key_dbl(fptr, "CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4211 fits_update_key_dbl(fptr, "CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4212 fits_update_key_dbl(fptr, "CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4213 fits_update_key_dbl(fptr, "CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4214 }
4215 }
4216 }
4217
4218 /* Delete any polynomial solution */
4219 /* (These could maybe be switched, but I don't want to work them out yet */
4220 status = 0;
4221 if (!fits_read_key_dbl(fptr, "CO1_1", &ctemp1, comment, &status))
4222 {
4223 int i;
4224 char keyword[16];
4225
4226 for (i = 1; i < 13; i++)
4227 {
4228 sprintf(keyword, "CO1_%d", i);
4229 fits_delete_key(fptr, keyword, &status);
4230 }
4231 for (i = 1; i < 13; i++)
4232 {
4233 sprintf(keyword, "CO2_%d", i);
4234 fits_delete_key(fptr, keyword, &status);
4235 }
4236 }
4237
4238}
4239
4240uint8_t * FITSData::getWritableImageBuffer()
4241{
4242 return m_ImageBuffer;
4243}
4244
4245uint8_t const * FITSData::getImageBuffer() const
4246{
4247 return m_ImageBuffer;
4248}
4249
4250void FITSData::setImageBuffer(uint8_t * buffer)
4251{
4252 delete[] m_ImageBuffer;
4253 m_ImageBuffer = buffer;
4254}
4255
4256bool FITSData::checkDebayer()
4257{
4258 int status = 0;
4259 char bayerPattern[64], roworder[64];
4260
4261 // Let's search for BAYERPAT keyword, if it's not found we return as there is no bayer pattern in this image
4262 if (fits_read_keyword(fptr, "BAYERPAT", bayerPattern, nullptr, &status))
4263 return false;
4264
4265 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
4266 {
4267 m_LastError = i18n("Only 8 and 16 bits bayered images supported.");
4268 return false;
4269 }
4270 QString pattern(bayerPattern);
4271 pattern = pattern.remove('\'').trimmed();
4272
4273 QString order(roworder);
4274 order = order.remove('\'').trimmed();
4275
4276 if (order == "BOTTOM-UP" && !(m_Statistics.height % 2))
4277 {
4278 if (pattern == "RGGB")
4279 pattern = "GBRG";
4280 else if (pattern == "GBRG")
4281 pattern = "RGGB";
4282 else if (pattern == "GRBG")
4283 pattern = "BGGR";
4284 else if (pattern == "BGGR")
4285 pattern = "GRBG";
4286 else return false;
4287 }
4288
4289 if (pattern == "RGGB")
4290 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4291 else if (pattern == "GBRG")
4292 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4293 else if (pattern == "GRBG")
4294 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4295 else if (pattern == "BGGR")
4296 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4297 // We return unless we find a valid pattern
4298 else
4299 {
4300 m_LastError = i18n("Unsupported bayer pattern %1.", pattern);
4301 return false;
4302 }
4303
4304 fits_read_key(fptr, TINT, "XBAYROFF", &debayerParams.offsetX, nullptr, &status);
4305 fits_read_key(fptr, TINT, "YBAYROFF", &debayerParams.offsetY, nullptr, &status);
4306
4307 if (debayerParams.offsetX == 1)
4308 {
4309 // This may leave odd values in the 0th column if the color filter is not there
4310 // in the sensor, but otherwise should process the offset correctly.
4311 // Only offsets of 0 or 1 are implemented in debayer_8bit() and debayer_16bit().
4312 switch (debayerParams.filter)
4313 {
4314 case DC1394_COLOR_FILTER_RGGB:
4315 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4316 break;
4317 case DC1394_COLOR_FILTER_GBRG:
4318 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4319 break;
4320 case DC1394_COLOR_FILTER_GRBG:
4321 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4322 break;
4323 case DC1394_COLOR_FILTER_BGGR:
4324 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4325 break;
4326 }
4327 debayerParams.offsetX = 0;
4328 }
4329 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
4330 {
4331 m_LastError = i18n("Unsupported bayer offsets %1 %2.", debayerParams.offsetX, debayerParams.offsetY);
4332 return false;
4333 }
4334
4335 HasDebayer = true;
4336
4337 return true;
4338}
4339
4340void FITSData::getBayerParams(BayerParams * param)
4341{
4342 param->method = debayerParams.method;
4343 param->filter = debayerParams.filter;
4344 param->offsetX = debayerParams.offsetX;
4345 param->offsetY = debayerParams.offsetY;
4346}
4347
4348void FITSData::setBayerParams(BayerParams * param)
4349{
4350 debayerParams.method = param->method;
4351 debayerParams.filter = param->filter;
4352 debayerParams.offsetX = param->offsetX;
4353 debayerParams.offsetY = param->offsetY;
4354}
4355
4356bool FITSData::debayer(bool reload)
4357{
4358 if (reload)
4359 {
4360 int anynull = 0, status = 0;
4361
4362 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel, nullptr, m_ImageBuffer,
4363 &anynull, &status))
4364 {
4365 // char errmsg[512];
4366 // fits_get_errstatus(status, errmsg);
4367 // KSNotification::error(i18n("Error reading image: %1", QString(errmsg)), i18n("Debayer error"));
4368 return false;
4369 }
4370 }
4371
4372 switch (m_Statistics.dataType)
4373 {
4374 case TBYTE:
4375 return debayer_8bit();
4376
4377 case TUSHORT:
4378 return debayer_16bit();
4379
4380 default:
4381 return false;
4382 }
4383}
4384
4385bool FITSData::debayer_8bit()
4386{
4387 dc1394error_t error_code;
4388
4389 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4390 uint8_t * destinationBuffer = nullptr;
4391
4392 try
4393 {
4394 destinationBuffer = new uint8_t[rgb_size];
4395 }
4396 catch (const std::bad_alloc &e)
4397 {
4398 logOOMError(rgb_size);
4399 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4400 return false;
4401 }
4402
4403 auto * bayer_source_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
4404 auto * bayer_destination_buffer = reinterpret_cast<uint8_t *>(destinationBuffer);
4405
4406 if (bayer_destination_buffer == nullptr)
4407 {
4408 logOOMError(rgb_size);
4409 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer.");
4410 return false;
4411 }
4412
4413 int ds1394_height = m_Statistics.height;
4414 auto dc1394_source = bayer_source_buffer;
4415
4416 if (debayerParams.offsetY == 1)
4417 {
4418 dc1394_source += m_Statistics.width;
4419 ds1394_height--;
4420 }
4421 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
4422
4423 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4424 debayerParams.filter,
4425 debayerParams.method);
4426
4427 if (error_code != DC1394_SUCCESS)
4428 {
4429 m_LastError = i18n("Debayer failed (%1)", error_code);
4430 m_Statistics.channels = 1;
4431 delete[] destinationBuffer;
4432 return false;
4433 }
4434
4435 if (m_ImageBufferSize != rgb_size)
4436 {
4437 delete[] m_ImageBuffer;
4438 try
4439 {
4440 m_ImageBuffer = new uint8_t[rgb_size];
4441 }
4442 catch (const std::bad_alloc &e)
4443 {
4444 delete[] destinationBuffer;
4445 logOOMError(rgb_size);
4446 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4447 return false;
4448 }
4449
4450 m_ImageBufferSize = rgb_size;
4451 }
4452
4453 auto bayered_buffer = reinterpret_cast<uint8_t *>(m_ImageBuffer);
4454
4455 // Data in R1G1B1, we need to copy them into 3 layers for FITS
4456
4457 uint8_t * rBuff = bayered_buffer;
4458 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4459 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4460
4461 int imax = m_Statistics.samples_per_channel * 3 - 3;
4462 for (int i = 0; i <= imax; i += 3)
4463 {
4464 *rBuff++ = bayer_destination_buffer[i];
4465 *gBuff++ = bayer_destination_buffer[i + 1];
4466 *bBuff++ = bayer_destination_buffer[i + 2];
4467 }
4468
4469 // TODO Maybe all should be treated the same
4470 // Doing single channel saves lots of memory though for non-essential
4471 // frames
4472 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4473 m_Statistics.dataType = TBYTE;
4474 delete[] destinationBuffer;
4475 return true;
4476}
4477
4478bool FITSData::debayer_16bit()
4479{
4480 dc1394error_t error_code;
4481
4482 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4483 uint8_t *destinationBuffer = nullptr;
4484 try
4485 {
4486 destinationBuffer = new uint8_t[rgb_size];
4487 }
4488 catch (const std::bad_alloc &e)
4489 {
4490 logOOMError(rgb_size);
4491 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4492 return false;
4493 }
4494
4495 auto * bayer_source_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
4496 auto * bayer_destination_buffer = reinterpret_cast<uint16_t *>(destinationBuffer);
4497
4498 if (bayer_destination_buffer == nullptr)
4499 {
4500 logOOMError(rgb_size);
4501 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer.");
4502 return false;
4503 }
4504
4505 int ds1394_height = m_Statistics.height;
4506 auto dc1394_source = bayer_source_buffer;
4507
4508 if (debayerParams.offsetY == 1)
4509 {
4510 dc1394_source += m_Statistics.width;
4511 ds1394_height--;
4512 }
4513 // offsetX == 1 is handled in checkDebayer() and should be 0 here.
4514
4515 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4516 debayerParams.filter,
4517 debayerParams.method, 16);
4518
4519 if (error_code != DC1394_SUCCESS)
4520 {
4521 m_LastError = i18n("Debayer failed (%1)");
4522 m_Statistics.channels = 1;
4523 delete[] destinationBuffer;
4524 return false;
4525 }
4526
4527 if (m_ImageBufferSize != rgb_size)
4528 {
4529 delete[] m_ImageBuffer;
4530 try
4531 {
4532 m_ImageBuffer = new uint8_t[rgb_size];
4533 }
4534 catch (const std::bad_alloc &e)
4535 {
4536 logOOMError(rgb_size);
4537 delete[] destinationBuffer;
4538 m_LastError = i18n("Unable to allocate memory for temporary bayer buffer: %1", e.what());
4539 return false;
4540 }
4541
4542 m_ImageBufferSize = rgb_size;
4543 }
4544
4545 auto bayered_buffer = reinterpret_cast<uint16_t *>(m_ImageBuffer);
4546
4547 // Data in R1G1B1, we need to copy them into 3 layers for FITS
4548
4549 uint16_t * rBuff = bayered_buffer;
4550 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4551 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4552
4553 int imax = m_Statistics.samples_per_channel * 3 - 3;
4554 for (int i = 0; i <= imax; i += 3)
4555 {
4556 *rBuff++ = bayer_destination_buffer[i];
4557 *gBuff++ = bayer_destination_buffer[i + 1];
4558 *bBuff++ = bayer_destination_buffer[i + 2];
4559 }
4560
4561 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4562 m_Statistics.dataType = TUSHORT;
4563 delete[] destinationBuffer;
4564 return true;
4565}
4566
4567void FITSData::logOOMError(uint32_t requiredMemory)
4568{
4569 qCCritical(KSTARS_FITS) << "Debayed memory allocation failure. Required Memory:" << KFormat().formatByteSize(requiredMemory)
4570 << "Available system memory:" << KSUtils::getAvailableRAM();
4571}
4572
4573double FITSData::getADU() const
4574{
4575 double adu = 0;
4576 for (int i = 0; i < m_Statistics.channels; i++)
4577 adu += m_Statistics.mean[i];
4578
4579 return (adu / static_cast<double>(m_Statistics.channels));
4580}
4581
4582QString FITSData::getLastError() const
4583{
4584 return m_LastError;
4585}
4586
4587template <typename T>
4588void FITSData::convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image)
4589{
4590 const auto * buffer = reinterpret_cast<const T *>(getImageBuffer());
4591 const T limit = std::numeric_limits<T>::max();
4592 T bMin = dataMin < 0 ? 0 : dataMin;
4593 T bMax = dataMax > limit ? limit : dataMax;
4594 uint16_t w = width();
4595 uint16_t h = height();
4596 uint32_t size = w * h;
4597 double val;
4598
4599 if (channels() == 1)
4600 {
4601 /* Fill in pixel values using indexed map, linear scale */
4602 for (int j = 0; j < h; j++)
4603 {
4604 unsigned char * scanLine = image.scanLine(j);
4605
4606 for (int i = 0; i < w; i++)
4607 {
4608 val = qBound(bMin, buffer[j * w + i], bMax);
4609 val = val * scale + zero;
4610 scanLine[i] = qBound<unsigned char>(0, static_cast<uint8_t>(val), 255);
4611 }
4612 }
4613 }
4614 else
4615 {
4616 double rval = 0, gval = 0, bval = 0;
4617 QRgb value;
4618 /* Fill in pixel values using indexed map, linear scale */
4619 for (int j = 0; j < h; j++)
4620 {
4621 auto * scanLine = reinterpret_cast<QRgb *>((image.scanLine(j)));
4622
4623 for (int i = 0; i < w; i++)
4624 {
4625 rval = qBound(bMin, buffer[j * w + i], bMax);
4626 gval = qBound(bMin, buffer[j * w + i + size], bMax);
4627 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
4628
4629 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
4630
4631 scanLine[i] = value;
4632 }
4633 }
4634 }
4635}
4636
4637QImage FITSData::FITSToImage(const QString &filename)
4638{
4639 QImage fitsImage;
4640 double min, max;
4641
4642 FITSData data;
4643
4644 QFuture<bool> future = data.loadFromFile(filename);
4645
4646 // Wait synchronously
4647 future.waitForFinished();
4648 if (future.result() == false)
4649 return fitsImage;
4650
4651 data.getMinMax(&min, &max);
4652
4653 if (min == max)
4654 {
4655 fitsImage.fill(Qt::white);
4656 return fitsImage;
4657 }
4658
4659 if (data.channels() == 1)
4660 {
4661 fitsImage = QImage(data.width(), data.height(), QImage::Format_Indexed8);
4662
4663 fitsImage.setColorCount(256);
4664 for (int i = 0; i < 256; i++)
4665 fitsImage.setColor(i, qRgb(i, i, i));
4666 }
4667 else
4668 {
4669 fitsImage = QImage(data.width(), data.height(), QImage::Format_RGB32);
4670 }
4671
4672 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
4673 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
4674
4675 double bscale = 255. / (dataMax - dataMin);
4676 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
4677
4678 // Long way to do this since we do not want to use templated functions here
4679 switch (data.m_Statistics.dataType)
4680 {
4681 case TBYTE:
4682 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4683 break;
4684
4685 case TSHORT:
4686 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4687 break;
4688
4689 case TUSHORT:
4690 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4691 break;
4692
4693 case TLONG:
4694 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4695 break;
4696
4697 case TULONG:
4698 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4699 break;
4700
4701 case TFLOAT:
4702 data.convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
4703 break;
4704
4705 case TLONGLONG:
4706 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4707 break;
4708
4709 case TDOUBLE:
4710 data.convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
4711 break;
4712
4713 default:
4714 break;
4715 }
4716
4717 return fitsImage;
4718}
4719
4720bool FITSData::ImageToFITS(const QString &filename, const QString &format, QString &output)
4721{
4722 if (QImageReader::supportedImageFormats().contains(format.toLatin1()) == false)
4723 {
4724 qCCritical(KSTARS_FITS) << "Failed to convert" << filename << "to FITS since format" << format << "is not supported in Qt";
4725 return false;
4726 }
4727
4728 QImage input;
4729
4730 if (input.load(filename, format.toLatin1()) == false)
4731 {
4732 qCCritical(KSTARS_FITS) << "Failed to open image" << filename;
4733 return false;
4734 }
4735
4736 output = QDir(getTemporaryPath()).filePath(filename + ".fits");
4737
4738 //This section sets up the FITS File
4739 fitsfile *fptr = nullptr;
4740 int status = 0;
4741 long fpixel = 1, naxis = input.allGray() ? 2 : 3, nelements, exposure;
4742 long naxes[3] = { input.width(), input.height(), naxis == 3 ? 3 : 1 };
4743 char error_status[512] = {0};
4744
4745 if (fits_create_file(&fptr, QString('!' + output).toLocal8Bit().data(), &status))
4746 {
4747 fits_get_errstatus(status, error_status);
4748 qCCritical(KSTARS_FITS) << "Failed to create FITS file. Error:" << error_status;
4749 return false;
4750 }
4751
4752 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
4753 {
4754 qCWarning(KSTARS_FITS) << "fits_create_img failed:" << error_status;
4755 status = 0;
4756 fits_flush_file(fptr, &status);
4757 fits_close_file(fptr, &status);
4758 return false;
4759 }
4760
4761 exposure = 1;
4762 fits_update_key(fptr, TLONG, "EXPOSURE", &exposure, "Total Exposure Time", &status);
4763
4764 // Gray image
4765 if (naxis == 2)
4766 {
4767 nelements = naxes[0] * naxes[1];
4768 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.bits(), &status))
4769 {
4770 fits_get_errstatus(status, error_status);
4771 qCWarning(KSTARS_FITS) << "fits_write_img GRAY failed:" << error_status;
4772 status = 0;
4773 fits_flush_file(fptr, &status);
4774 fits_close_file(fptr, &status);
4775 return false;
4776 }
4777 }
4778 // RGB image, we have to convert from ARGB format to R G B for each plane
4779 else
4780 {
4781 nelements = naxes[0] * naxes[1] * 3;
4782
4783 uint8_t *srcBuffer = input.bits();
4784 // ARGB
4785 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
4786
4787 uint8_t *rgbBuffer = new uint8_t[nelements];
4788 if (rgbBuffer == nullptr)
4789 {
4790 qCWarning(KSTARS_FITS) << "Not enough memory for RGB buffer";
4791 fits_flush_file(fptr, &status);
4792 fits_close_file(fptr, &status);
4793 return false;
4794 }
4795
4796 uint8_t *subR = rgbBuffer;
4797 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
4798 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
4799 for (uint32_t i = 0; i < srcBytes; i += 4)
4800 {
4801 *subB++ = srcBuffer[i];
4802 *subG++ = srcBuffer[i + 1];
4803 *subR++ = srcBuffer[i + 2];
4804 }
4805
4806 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
4807 {
4808 fits_get_errstatus(status, error_status);
4809 qCWarning(KSTARS_FITS) << "fits_write_img RGB failed:" << error_status;
4810 status = 0;
4811 fits_flush_file(fptr, &status);
4812 fits_close_file(fptr, &status);
4813 delete [] rgbBuffer;
4814 return false;
4815 }
4816
4817 delete [] rgbBuffer;
4818 }
4819
4820 if (fits_flush_file(fptr, &status))
4821 {
4822 fits_get_errstatus(status, error_status);
4823 qCWarning(KSTARS_FITS) << "fits_flush_file failed:" << error_status;
4824 status = 0;
4825 fits_close_file(fptr, &status);
4826 return false;
4827 }
4828
4829 if (fits_close_file(fptr, &status))
4830 {
4831 fits_get_errstatus(status, error_status);
4832 qCWarning(KSTARS_FITS) << "fits_close_file failed:" << error_status;
4833 return false;
4834 }
4835
4836 return true;
4837}
4838
4839void FITSData::injectWCS(double orientation, double ra, double dec, double pixscale, bool eastToTheRight)
4840{
4841 int status = 0;
4842
4843 updateWCSHeaderData(orientation, ra, dec, pixscale, eastToTheRight);
4844
4845 if (fptr)
4846 {
4847 fits_update_key(fptr, TDOUBLE, "OBJCTRA", &ra, "Object RA", &status);
4848 fits_update_key(fptr, TDOUBLE, "OBJCTDEC", &dec, "Object DEC", &status);
4849
4850 int epoch = 2000;
4851
4852 fits_update_key(fptr, TINT, "EQUINOX", &epoch, "Equinox", &status);
4853
4854 fits_update_key(fptr, TDOUBLE, "CRVAL1", &ra, "CRVAL1", &status);
4855 fits_update_key(fptr, TDOUBLE, "CRVAL2", &dec, "CRVAL1", &status);
4856
4857 char radecsys[8] = "FK5";
4858 char ctype1[16] = "RA---TAN";
4859 char ctype2[16] = "DEC--TAN";
4860
4861 fits_update_key(fptr, TSTRING, "RADECSYS", radecsys, "RADECSYS", &status);
4862 fits_update_key(fptr, TSTRING, "CTYPE1", ctype1, "CTYPE1", &status);
4863 fits_update_key(fptr, TSTRING, "CTYPE2", ctype2, "CTYPE2", &status);
4864
4865 double crpix1 = width() / 2.0;
4866 double crpix2 = height() / 2.0;
4867
4868 fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix1, "CRPIX1", &status);
4869 fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix2, "CRPIX2", &status);
4870
4871 // Arcsecs per Pixel
4872 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4873 double secpix2 = pixscale;
4874
4875 fits_update_key(fptr, TDOUBLE, "SECPIX1", &secpix1, "SECPIX1", &status);
4876 fits_update_key(fptr, TDOUBLE, "SECPIX2", &secpix2, "SECPIX2", &status);
4877
4878 double degpix1 = secpix1 / 3600.0;
4879 double degpix2 = secpix2 / 3600.0;
4880
4881 fits_update_key(fptr, TDOUBLE, "CDELT1", &degpix1, "CDELT1", &status);
4882 fits_update_key(fptr, TDOUBLE, "CDELT2", &degpix2, "CDELT2", &status);
4883
4884 // Rotation is CW, we need to convert it to CCW per CROTA1 definition
4885 double rotation = 360 - orientation;
4886 if (rotation > 360)
4887 rotation -= 360;
4888
4889 fits_update_key(fptr, TDOUBLE, "CROTA1", &rotation, "CROTA1", &status);
4890 fits_update_key(fptr, TDOUBLE, "CROTA2", &rotation, "CROTA2", &status);
4891 }
4892
4893 emit headerChanged();
4894 m_WCSState = Idle;
4895}
4896
4897// Update header records based on plate solver solution
4898void FITSData::updateWCSHeaderData(const double orientation, const double ra, const double dec, const double pixscale,
4899 const bool eastToTheRight)
4900{
4901 updateRecordValue("OBJCTRA", ra, "Object RA");
4902 updateRecordValue("OBJCTDEC", dec, "Object DEC");
4903 updateRecordValue("EQUINOX", 2000, "Equinox");
4904 updateRecordValue("CRVAL1", ra, "CRVAL1");
4905 updateRecordValue("CRVAL2", dec, "CRVAL2");
4906
4907 updateRecordValue("RADECSYS", "'FK5'", "RADECSYS");
4908 updateRecordValue("CTYPE1", "'RA---TAN'", "CTYPE1");
4909 updateRecordValue("CTYPE2", "'DEC--TAN'", "CTYPE2");
4910
4911 updateRecordValue("CRPIX1", m_Statistics.width / 2.0, "CRPIX1");
4912 updateRecordValue("CRPIX2", m_Statistics.height / 2.0, "CRPIX2");
4913
4914 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4915 double secpix2 = pixscale;
4916
4917 double degpix1 = secpix1 / 3600.0;
4918 double degpix2 = secpix2 / 3600.0;
4919 updateRecordValue("CDELT1", degpix1, "CDELT1");
4920 updateRecordValue("CDELT2", degpix2, "CDELT2");
4921
4922 // Rotation is CW, we need to convert it to CCW per CROTA1 definition
4923 double rotation = 360 - orientation;
4924 if (rotation > 360)
4925 rotation -= 360;
4926 updateRecordValue("CROTA1", rotation, "CROTA1");
4927 updateRecordValue("CROTA2", rotation, "CROTA2");
4928}
4929
4930bool FITSData::contains(const QPointF &point) const
4931{
4932 return (point.x() >= 0 && point.y() >= 0 && point.x() <= m_Statistics.width && point.y() <= m_Statistics.height);
4933}
4934
4935void FITSData::saveStatistics(FITSImage::Statistic &other)
4936{
4937 other = m_Statistics;
4938}
4939
4940void FITSData::restoreStatistics(FITSImage::Statistic &other)
4941{
4942 m_Statistics = other;
4943
4944 emit dataChanged();
4945}
4946
4947void FITSData::constructHistogram()
4948{
4949 switch (m_Statistics.dataType)
4950 {
4951 case TBYTE:
4952 constructHistogramInternal<uint8_t>();
4953 break;
4954
4955 case TSHORT:
4956 constructHistogramInternal<int16_t>();
4957 break;
4958
4959 case TUSHORT:
4960 constructHistogramInternal<uint16_t>();
4961 break;
4962
4963 case TLONG:
4964 constructHistogramInternal<int32_t>();
4965 break;
4966
4967 case TULONG:
4968 constructHistogramInternal<uint32_t>();
4969 break;
4970
4971 case TFLOAT:
4972 constructHistogramInternal<float>();
4973 break;
4974
4975 case TLONGLONG:
4976 constructHistogramInternal<int64_t>();
4977 break;
4978
4979 case TDOUBLE:
4980 constructHistogramInternal<double>();
4981 break;
4982
4983 default:
4984 break;
4985 }
4986}
4987
4988template <typename T> int32_t FITSData::histogramBinInternal(T value, int channel) const
4989{
4990 return qMax(static_cast<T>(0), qMin(static_cast<T>(m_HistogramBinCount),
4991 static_cast<T>(rint((value - m_Statistics.min[channel]) / m_HistogramBinWidth[channel]))));
4992}
4993
4994template <typename T> int32_t FITSData::histogramBinInternal(int x, int y, int channel) const
4995{
4996 if (!m_ImageBuffer || !isHistogramConstructed())
4997 return 0;
4998 uint32_t samples = m_Statistics.width * m_Statistics.height;
4999 uint32_t offset = channel * samples;
5000 auto * const buffer = reinterpret_cast<T const *>(m_ImageBuffer);
5001 int index = y * m_Statistics.width + x;
5002 const T &sample = buffer[index + offset];
5003 return histogramBinInternal(sample, channel);
5004}
5005
5006int32_t FITSData::histogramBin(int x, int y, int channel) const
5007{
5008 switch (m_Statistics.dataType)
5009 {
5010 case TBYTE:
5011 return histogramBinInternal<uint8_t>(x, y, channel);
5012 break;
5013
5014 case TSHORT:
5015 return histogramBinInternal<int16_t>(x, y, channel);
5016 break;
5017
5018 case TUSHORT:
5019 return histogramBinInternal<uint16_t>(x, y, channel);
5020 break;
5021
5022 case TLONG:
5023 return histogramBinInternal<int32_t>(x, y, channel);
5024 break;
5025
5026 case TULONG:
5027 return histogramBinInternal<uint32_t>(x, y, channel);
5028 break;
5029
5030 case TFLOAT:
5031 return histogramBinInternal<float>(x, y, channel);
5032 break;
5033
5034 case TLONGLONG:
5035 return histogramBinInternal<int64_t>(x, y, channel);
5036 break;
5037
5038 case TDOUBLE:
5039 return histogramBinInternal<double>(x, y, channel);
5040 break;
5041
5042 default:
5043 return 0;
5044 break;
5045 }
5046}
5047
5048template <typename T> void FITSData::constructHistogramInternal()
5049{
5050 auto * const buffer = reinterpret_cast<T const *>(m_ImageBuffer);
5051 uint32_t samples = m_Statistics.width * m_Statistics.height;
5052 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
5053
5054 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
5055 if (m_HistogramBinCount <= 0)
5056 m_HistogramBinCount = 256;
5057
5058 for (int n = 0; n < m_Statistics.channels; n++)
5059 {
5060 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
5061 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
5062 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
5063 // Distinguish between 0-1.0 ranges and ranges with integer values.
5064 const double minBinSize = (m_Statistics.max[n] > 1.1) ? 1.0 : .0001;
5065 m_HistogramBinWidth[n] = qMax(minBinSize, (m_Statistics.max[n] - m_Statistics.min[n]) / m_HistogramBinCount);
5066 }
5067
5068 QVector<QFuture<void>> futures;
5069
5070 for (int n = 0; n < m_Statistics.channels; n++)
5071 {
5072 futures.append(QtConcurrent::run([ = ]()
5073 {
5074 for (int i = 0; i < m_HistogramBinCount; i++)
5075 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
5076 }));
5077 }
5078
5079 for (int n = 0; n < m_Statistics.channels; n++)
5080 {
5081 futures.append(QtConcurrent::run([ = ]()
5082 {
5083 uint32_t offset = n * samples;
5084
5085 for (uint32_t i = 0; i < samples; i += sampleBy)
5086 {
5087 int32_t id = histogramBinInternal<T>(buffer[i + offset], n);
5088 m_HistogramFrequency[n][id] += sampleBy;
5089 }
5090 }));
5091 }
5092
5093 for (QFuture<void> future : futures)
5094 future.waitForFinished();
5095
5096 futures.clear();
5097
5098 for (int n = 0; n < m_Statistics.channels; n++)
5099 {
5100 futures.append(QtConcurrent::run([ = ]()
5101 {
5102 uint32_t accumulator = 0;
5103 for (int i = 0; i < m_HistogramBinCount; i++)
5104 {
5105 accumulator += m_HistogramFrequency[n][i];
5106 m_CumulativeFrequency[n].replace(i, accumulator);
5107 }
5108 }));
5109 }
5110
5111 for (QFuture<void> future : futures)
5112 future.waitForFinished();
5113
5114 futures.clear();
5115
5116 // Custom index to indicate the overall contrast of the image
5117 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
5118 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] / static_cast<double>
5119 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
5120 4]);
5121 else
5122 m_JMIndex = 1;
5123
5124 qCDebug(KSTARS_FITS) << "FITHistogram: JMIndex " << m_JMIndex;
5125
5126 m_HistogramConstructed = true;
5127 emit histogramReady();
5128}
5129
5130void FITSData::recordLastError(int errorCode)
5131{
5132 char fitsErrorMessage[512] = {0};
5133 fits_get_errstatus(errorCode, fitsErrorMessage);
5134 m_LastError = fitsErrorMessage;
5135}
5136
5137double FITSData::getAverageMean(bool roi) const
5138{
5139 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5140 if (ptr->channels == 1)
5141 return ptr->mean[0];
5142 else
5143 return (ptr->mean[0] + ptr->mean[1] + ptr->mean[2]) / 3.0;
5144}
5145
5146double FITSData::getAverageMedian(bool roi) const
5147{
5148 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5149 if (ptr->channels == 1)
5150 return ptr->median[0];
5151 else
5152 return (ptr->median[0] + ptr->median[1] + ptr->median[2]) / 3.0;
5153}
5154
5155double FITSData::getAverageStdDev(bool roi) const
5156{
5157 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5158 if (ptr->channels == 1)
5159 return ptr->stddev[0];
5160 else
5161 return (ptr->stddev[0] + ptr->stddev[1] + ptr->stddev[2]) / 3.0;
5162}
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-2024 The KDE developers.
Generated on Sat Dec 21 2024 17:04:46 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.