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

KDE's Doxygen guidelines are available online.