11#include "fitsbahtinovdetector.h"
12#include "fitsthresholddetector.h"
13#include "fitsgradientdetector.h"
14#include "fitscentroiddetector.h"
15#include "fitssepdetector.h"
19#include "kstarsdata.h"
23#include "skymapcomposite.h"
24#include "auxiliary/ksnotification.h"
25#include "auxiliary/robuststatistics.h"
28#include <QApplication>
30#include <QtConcurrent>
31#include <QImageReader>
33#include <QNetworkAccessManager>
35#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
40#if !defined(KSTARS_LITE) && defined(HAVE_LIBRAW)
41#include <libraw/libraw.h>
51#include <fits_debug.h>
53#define ZOOM_DEFAULT 100.0
56#define ZOOM_LOW_INCR 10
57#define ZOOM_HIGH_INCR 50
65const QStringList RAWFormats = {
"cr2",
"cr3",
"crw",
"nef",
"raf",
"dng",
"arw",
"orf" };
67bool FITSData::readableFilename(
const QString &filename)
69 QFileInfo info(filename);
70 QString extension = info.completeSuffix().toLower();
71 return extension.contains(
"fit") || extension.contains(
"fz") ||
72 extension.contains(
"xisf") ||
77FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)
81 qRegisterMetaType<FITSMode>(
"FITSMode");
83 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
84 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
85 debayerParams.offsetX = debayerParams.offsetY = 0;
88 m_CumulativeFrequency.resize(3);
89 m_HistogramBinWidth.resize(3);
90 m_HistogramFrequency.resize(3);
91 m_HistogramIntensity.resize(3);
101 static const QRegularExpression re(
"[-{}]");
102 qRegisterMetaType<FITSMode>(
"FITSMode");
104 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
105 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
106 debayerParams.offsetX = debayerParams.offsetY = 0;
108 m_TemporaryDataFile.setFileTemplate(
"fits_memory_XXXXXX");
110 this->m_Mode = other->m_Mode;
111 this->m_Statistics.channels = other->m_Statistics.channels;
112 memcpy(&m_Statistics, &(other->m_Statistics),
sizeof(m_Statistics));
113 m_ImageBuffer =
new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel];
114 memcpy(m_ImageBuffer, other->m_ImageBuffer,
115 m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel);
127 if (m_StarFindFuture.isRunning())
128 m_StarFindFuture.waitForFinished();
133 if (m_WCSHandle !=
nullptr)
135 wcsvfree(&m_nwcs, &m_WCSHandle);
136 m_WCSHandle =
nullptr;
141 if (starCenters.count() > 0)
142 qDeleteAll(starCenters);
145 if (m_SkyObjects.count() > 0)
146 qDeleteAll(m_SkyObjects);
147 m_SkyObjects.clear();
149 m_CatObjects.clear();
153 fits_flush_file(fptr, &status);
154 fits_close_file(fptr, &status);
156 m_PackBuffer =
nullptr;
161void FITSData::loadCommon(
const QString &inFilename)
164 qDeleteAll(starCenters);
169 fits_flush_file(fptr, &status);
170 fits_close_file(fptr, &status);
172 m_PackBuffer =
nullptr;
176 m_Filename = inFilename;
179bool FITSData::loadFromBuffer(
const QByteArray &buffer)
182 qCDebug(KSTARS_FITS) <<
"Reading file buffer (" << KFormat().formatByteSize(buffer.
size()) <<
")";
183 return privateLoad(buffer);
188 loadCommon(inFilename);
189 QFileInfo info(m_Filename);
190 m_Extension = info.completeSuffix().toLower();
191 qCDebug(KSTARS_FITS) <<
"Loading file " << m_Filename;
192#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
202QString fitsErrorToString(
int status)
204 char error_status[512] = {0};
205 fits_report_error(stderr, status);
206 fits_get_errstatus(status, error_status);
207 QString message = error_status;
212bool FITSData::privateLoad(
const QByteArray &buffer)
216 cacheEccentricity = -1;
218 if (m_Extension.contains(
"fit") || m_Extension.contains(
"fz"))
219 return loadFITSImage(buffer);
220 if (m_Extension.contains(
"xisf"))
221 return loadXISFImage(buffer);
223 return loadCanonicalImage(buffer);
224 else if (RAWFormats.
contains(m_Extension))
225 return loadRAWImage(buffer);
230bool FITSData::loadFITSImage(
const QByteArray &buffer,
const bool isCompressed)
232 int status = 0, anynull = 0;
235 m_HistogramConstructed =
false;
237 if (m_Extension.contains(
".fz") || isCompressed)
246 m_compressedFilename = m_Filename;
249 QRegularExpression(
"[-{}]")));
251 rc = fp_unpack_file_to_fits(m_Filename.toLocal8Bit().data(), &fptr, fpvar) == 0;
254 m_Filename = uncompressedFile;
259 size_t m_PackBufferSize = 100000;
261 m_PackBuffer = (uint8_t *)malloc(m_PackBufferSize);
262 rc = fp_unpack_data_to_data(buffer.
data(), buffer.
size(), &m_PackBuffer, &m_PackBufferSize, fpvar) == 0;
266 void *data =
reinterpret_cast<void *
>(m_PackBuffer);
267 if (fits_open_memfile(&fptr, m_Filename.toLocal8Bit().data(), READONLY, &data, &m_PackBufferSize, 0,
270 m_LastError =
i18n(
"Error reading fits buffer: %1.", fitsErrorToString(status));
274 m_Statistics.size = m_PackBufferSize;
282 m_PackBuffer =
nullptr;
283 m_LastError =
i18n(
"Failed to unpack compressed fits");
284 qCCritical(KSTARS_FITS) << m_LastError;
288 m_isTemporary =
true;
289 m_isCompressed =
true;
290 m_Statistics.size = fptr->Fptr->logfilesize;
297 if (fits_open_diskfile(&fptr, m_Filename.toLocal8Bit(), READONLY, &status))
299 m_LastError =
i18n(
"Error opening fits file %1 : %2", m_Filename, fitsErrorToString(status));
302 m_Statistics.size = QFile(m_Filename).size();
307 void *temp_buffer =
const_cast<void *
>(
reinterpret_cast<const void *
>(buffer.
data()));
308 size_t temp_size = buffer.
size();
309 if (fits_open_memfile(&fptr, m_Filename.toLocal8Bit().data(), READONLY,
310 &temp_buffer, &temp_size, 0,
nullptr, &status))
312 m_LastError =
i18n(
"Error reading fits buffer: %1", fitsErrorToString(status));
316 m_Statistics.size = temp_size;
319 if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status))
323 m_PackBuffer =
nullptr;
324 m_LastError =
i18n(
"Could not locate image HDU: %1", fitsErrorToString(status));
327 if (fits_get_img_param(fptr, 3, &m_FITSBITPIX, &(m_Statistics.ndim), naxes, &status))
330 m_PackBuffer =
nullptr;
331 m_LastError =
i18n(
"FITS file open error (fits_get_img_param): %1", fitsErrorToString(status));
336 if ((fits_is_compressed_image(fptr, &status) || m_Statistics.ndim <= 0) && !isCompressed)
338 loadCommon(m_Filename);
339 qCDebug(KSTARS_FITS) <<
"Image is compressed. Reloading...";
340 return loadFITSImage(buffer,
true);
343 if (m_Statistics.ndim < 2)
345 m_LastError =
i18n(
"1D FITS images are not supported in KStars.");
346 qCCritical(KSTARS_FITS) << m_LastError;
348 m_PackBuffer =
nullptr;
352 switch (m_FITSBITPIX)
355 m_Statistics.dataType = TBYTE;
356 m_Statistics.bytesPerPixel =
sizeof(uint8_t);
360 m_FITSBITPIX = USHORT_IMG;
361 m_Statistics.dataType = TUSHORT;
362 m_Statistics.bytesPerPixel =
sizeof(int16_t);
365 m_Statistics.dataType = TUSHORT;
366 m_Statistics.bytesPerPixel =
sizeof(uint16_t);
370 m_FITSBITPIX = ULONG_IMG;
371 m_Statistics.dataType = TULONG;
372 m_Statistics.bytesPerPixel =
sizeof(int32_t);
375 m_Statistics.dataType = TULONG;
376 m_Statistics.bytesPerPixel =
sizeof(uint32_t);
379 m_Statistics.dataType = TFLOAT;
380 m_Statistics.bytesPerPixel =
sizeof(float);
383 m_Statistics.dataType = TLONGLONG;
384 m_Statistics.bytesPerPixel =
sizeof(int64_t);
387 m_Statistics.dataType = TDOUBLE;
388 m_Statistics.bytesPerPixel =
sizeof(double);
391 m_LastError =
i18n(
"Bit depth %1 is not supported.", m_FITSBITPIX);
392 qCCritical(KSTARS_FITS) << m_LastError;
396 if (m_Statistics.ndim < 3)
399 if (naxes[0] == 0 || naxes[1] == 0)
401 m_LastError =
i18n(
"Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
402 qCCritical(KSTARS_FITS) << m_LastError;
404 m_PackBuffer =
nullptr;
408 m_Statistics.width = naxes[0];
409 m_Statistics.height = naxes[1];
410 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
411 roiCenter.setX(m_Statistics.width / 2);
412 roiCenter.setY(m_Statistics.height / 2);
413 if(m_Statistics.width % 2)
414 roiCenter.setX(roiCenter.x() + 1);
415 if(m_Statistics.height % 2)
416 roiCenter.setY(roiCenter.y() + 1);
420 m_Statistics.channels = naxes[2];
424 if ( (m_Mode != FITS_NORMAL && m_Mode != FITS_CALIBRATE) || !Options::auto3DCube())
425 m_Statistics.channels = 1;
427 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
428 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
429 if (m_ImageBuffer ==
nullptr)
431 qCWarning(KSTARS_FITS) <<
"FITSData: Not enough memory for image_buffer channel. Requested: "
432 << m_ImageBufferSize <<
" bytes.";
435 m_PackBuffer =
nullptr;
442 long nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
444 if (fits_read_img(fptr, m_Statistics.dataType, 1, nelements,
nullptr, m_ImageBuffer, &anynull, &status))
446 m_LastError =
i18n(
"Error reading image: %1", fitsErrorToString(status));
454 if (getRecordValue(
"DATE-OBS", value) && value.
isValid())
457 m_DateTime = KStarsDateTime(ts.
date(), ts.
time());
462 if (naxes[2] == 1 && m_Statistics.channels == 1 && Options::autoDebayer() && checkDebayer())
465 if (m_isTemporary && m_TemporaryDataFile.open())
467 m_TemporaryDataFile.write(buffer);
468 m_TemporaryDataFile.close();
469 m_Filename = m_TemporaryDataFile.fileName();
473 calculateStats(
false,
false);
476 calculateStats(
false,
false);
478 if (m_Mode == FITS_NORMAL || m_Mode == FITS_ALIGN)
481 starsSearched =
false;
486bool FITSData::loadXISFImage(
const QByteArray &buffer)
488 m_HistogramConstructed =
false;
494 LibXISF::XISFReader xisfReader;
497 xisfReader.open(m_Filename.toLocal8Bit().data());
501 LibXISF::ByteArray byteArray(buffer.
constData(), buffer.
size());
502 xisfReader.open(byteArray);
505 if (xisfReader.imagesCount() == 0)
507 m_LastError =
i18n(
"File contain no images");
511 const LibXISF::Image &image = xisfReader.getImage(0);
513 switch (image.sampleFormat())
515 case LibXISF::Image::UInt8:
516 m_Statistics.dataType = TBYTE;
517 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt8);
518 m_FITSBITPIX = TBYTE;
520 case LibXISF::Image::UInt16:
521 m_Statistics.dataType = TUSHORT;
522 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt16);
523 m_FITSBITPIX = TUSHORT;
525 case LibXISF::Image::UInt32:
526 m_Statistics.dataType = TULONG;
527 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt32);
528 m_FITSBITPIX = TULONG;
530 case LibXISF::Image::Float32:
531 m_Statistics.dataType = TFLOAT;
532 m_Statistics.bytesPerPixel =
sizeof(LibXISF::Float32);
533 m_FITSBITPIX = TFLOAT;
536 m_LastError =
i18n(
"Sample format %1 is not supported.", LibXISF::Image::sampleFormatString(image.sampleFormat()).c_str());
537 qCCritical(KSTARS_FITS) << m_LastError;
541 m_Statistics.width = image.width();
542 m_Statistics.height = image.height();
543 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
544 m_Statistics.channels = image.channelCount();
545 m_Statistics.size = buffer.
size();
546 roiCenter.setX(m_Statistics.width / 2);
547 roiCenter.setY(m_Statistics.height / 2);
548 if(m_Statistics.width % 2)
549 roiCenter.setX(roiCenter.x() + 1);
550 if(m_Statistics.height % 2)
551 roiCenter.setY(roiCenter.y() + 1);
553 m_HeaderRecords.clear();
554 auto &fitsKeywords = image.fitsKeywords();
555 for(
auto &fitsKeyword : fitsKeywords)
562 m_ImageBufferSize = image.imageDataSize();
563 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
564 std::memcpy(m_ImageBuffer, image.imageData(), m_ImageBufferSize);
566 calculateStats(
false,
false);
569 catch (LibXISF::Error &error)
571 m_LastError =
i18n(
"XISF file open error: ") +
error.what();
582bool FITSData::saveXISFImage(
const QString &newFilename)
587 LibXISF::XISFWriter xisfWriter;
588 LibXISF::Image image;
589 image.setGeometry(m_Statistics.width, m_Statistics.height, m_Statistics.channels);
590 if (m_Statistics.channels > 1)
591 image.setColorSpace(LibXISF::Image::RGB);
593 switch (m_FITSBITPIX)
596 image.setSampleFormat(LibXISF::Image::UInt8);
599 image.setSampleFormat(LibXISF::Image::UInt16);
602 image.setSampleFormat(LibXISF::Image::UInt32);
605 image.setSampleFormat(LibXISF::Image::Float32);
608 m_LastError =
i18n(
"Bit depth %1 is not supported.", m_FITSBITPIX);
609 qCCritical(KSTARS_FITS) << m_LastError;
613 std::memcpy(image.imageData(), m_ImageBuffer, m_ImageBufferSize);
614 for (
auto &fitsKeyword : m_HeaderRecords)
615 image.addFITSKeyword({fitsKeyword.key.toUtf8().data(), fitsKeyword.value.toString().toUtf8().data(), fitsKeyword.comment.toUtf8().data()});
617 xisfWriter.writeImage(image);
619 m_Filename = newFilename;
621 catch (LibXISF::Error &err)
623 m_LastError =
i18n(
"Error saving XISF image") + err.what();
628 Q_UNUSED(newFilename)
633bool FITSData::loadCanonicalImage(
const QByteArray &buffer)
635 QImage imageFromFile;
640 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
644 m_Statistics.size = buffer.
size();
648 if(!imageFromFile.
load(m_Filename.toLocal8Bit()))
650 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
654 m_Statistics.size = QFileInfo(m_Filename).size();
662 m_FITSBITPIX = BYTE_IMG;
663 switch (m_FITSBITPIX)
666 m_Statistics.dataType = TBYTE;
667 m_Statistics.bytesPerPixel =
sizeof(uint8_t);
671 m_Statistics.dataType = TUSHORT;
672 m_Statistics.bytesPerPixel =
sizeof(int16_t);
675 m_Statistics.dataType = TUSHORT;
676 m_Statistics.bytesPerPixel =
sizeof(uint16_t);
680 m_Statistics.dataType = TULONG;
681 m_Statistics.bytesPerPixel =
sizeof(int32_t);
684 m_Statistics.dataType = TULONG;
685 m_Statistics.bytesPerPixel =
sizeof(uint32_t);
688 m_Statistics.dataType = TFLOAT;
689 m_Statistics.bytesPerPixel =
sizeof(float);
692 m_Statistics.dataType = TLONGLONG;
693 m_Statistics.bytesPerPixel =
sizeof(int64_t);
696 m_Statistics.dataType = TDOUBLE;
697 m_Statistics.bytesPerPixel =
sizeof(double);
700 m_LastError = QString(
"Bit depth %1 is not supported.").arg(m_FITSBITPIX);
701 qCCritical(KSTARS_FITS) << m_LastError;
705 m_Statistics.width =
static_cast<uint16_t
>(imageFromFile.
width());
706 m_Statistics.height =
static_cast<uint16_t
>(imageFromFile.
height());
708 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
710 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels *
static_cast<uint16_t
>
711 (m_Statistics.bytesPerPixel);
712 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
713 if (m_ImageBuffer ==
nullptr)
715 m_LastError =
i18n(
"FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
716 qCCritical(KSTARS_FITS) << m_LastError;
721 if (m_Statistics.channels == 1)
723 memcpy(m_ImageBuffer, imageFromFile.
bits(), m_ImageBufferSize);
728 auto debayered_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
729 auto * original_bayered_buffer =
reinterpret_cast<uint8_t *
>(imageFromFile.
bits());
732 uint8_t * rBuff = debayered_buffer;
733 uint8_t * gBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height);
734 uint8_t * bBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
736 int imax = m_Statistics.samples_per_channel * 4 - 4;
737 for (
int i = 0; i <= imax; i += 4)
739 *rBuff++ = original_bayered_buffer[i + 2];
740 *gBuff++ = original_bayered_buffer[i + 1];
741 *bBuff++ = original_bayered_buffer[i + 0];
745 m_HeaderRecords.clear();
748 calculateStats(
false,
false);
753bool FITSData::loadRAWImage(
const QByteArray &buffer)
755#if !defined(KSTARS_LITE) && !defined(HAVE_LIBRAW)
756 m_LastError =
i18n(
"Unable to find dcraw and cjpeg. Please install the required tools to convert CR2/NEF to JPEG.");
768 if ((ret = RawProcessor.open_file(m_Filename.toLocal8Bit().constData())) != LIBRAW_SUCCESS)
770 m_LastError =
i18n(
"Cannot open file %1: %2", m_Filename, libraw_strerror(ret));
771 RawProcessor.recycle();
775 m_Statistics.size = QFileInfo(m_Filename).size();
780 if ((ret = RawProcessor.open_buffer(
const_cast<void *
>(
reinterpret_cast<const void *
>(buffer.
data())),
781 buffer.
size())) != LIBRAW_SUCCESS)
783 m_LastError =
i18n(
"Cannot open buffer: %1", libraw_strerror(ret));
784 RawProcessor.recycle();
788 m_Statistics.size = buffer.
size();
792 if ((ret = RawProcessor.unpack()) != LIBRAW_SUCCESS)
794 m_LastError =
i18n(
"Cannot unpack_thumb: %1", libraw_strerror(ret));
795 RawProcessor.recycle();
799 if ((ret = RawProcessor.dcraw_process()) != LIBRAW_SUCCESS)
801 m_LastError =
i18n(
"Cannot dcraw_process: %1", libraw_strerror(ret));
802 RawProcessor.recycle();
806 libraw_processed_image_t *image = RawProcessor.dcraw_make_mem_image(&ret);
807 if (ret != LIBRAW_SUCCESS)
809 m_LastError =
i18n(
"Cannot load to memory: %1", libraw_strerror(ret));
810 RawProcessor.recycle();
814 RawProcessor.recycle();
816 m_Statistics.bytesPerPixel = image->bits / 8;
818 if (m_Statistics.bytesPerPixel == 1)
819 m_Statistics.dataType = TBYTE;
821 m_Statistics.dataType = TUSHORT;
822 m_Statistics.width = image->width;
823 m_Statistics.height = image->height;
824 m_Statistics.channels = image->colors;
825 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
827 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
828 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
829 if (m_ImageBuffer ==
nullptr)
831 m_LastError =
i18n(
"FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
832 qCCritical(KSTARS_FITS) << m_LastError;
833 libraw_dcraw_clear_mem(image);
838 auto destination_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
839 auto source_buffer =
reinterpret_cast<uint8_t *
>(image->data);
842 if (image->colors == 1)
844 memcpy(destination_buffer, source_buffer, m_ImageBufferSize);
849 uint8_t * rBuff = destination_buffer;
850 uint8_t * gBuff = destination_buffer + (m_Statistics.width * m_Statistics.height);
851 uint8_t * bBuff = destination_buffer + (m_Statistics.width * m_Statistics.height * 2);
853 int imax = m_Statistics.samples_per_channel * 3 - 3;
854 for (
int i = 0; i <= imax; i += 3)
856 *rBuff++ = source_buffer[i + 0];
857 *gBuff++ = source_buffer[i + 1];
858 *bBuff++ = source_buffer[i + 2];
861 libraw_dcraw_clear_mem(image);
863 m_HeaderRecords.clear();
866 calculateStats(
false,
false);
872bool FITSData::saveImage(
const QString &newFilename)
874 if (newFilename == m_Filename)
877 const QString ext = QFileInfo(newFilename).suffix();
879 if (ext ==
"jpg" || ext ==
"png")
882 QImage fitsImage = FITSToImage(newFilename);
883 getMinMax(&min, &max);
889 else if (channels() == 1)
894 for (
int i = 0; i < 256; i++)
895 fitsImage.
setColor(i, qRgb(i, i, i));
902 double dataMin = m_Statistics.mean[0] - m_Statistics.stddev[0];
903 double dataMax = m_Statistics.mean[0] + m_Statistics.stddev[0] * 3;
905 double bscale = 255. / (dataMax - dataMin);
906 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
909 switch (m_Statistics.dataType)
912 convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
916 convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
920 convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
924 convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
928 convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
932 convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
936 convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
940 convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
949 m_Filename = newFilename;
950 qCInfo(KSTARS_FITS) <<
"Saved image file:" << m_Filename;
956 m_LastError =
i18n(
"Saving compressed files is not supported.");
968 fits_close_file(fptr, &status);
971 rotCounter = flipHCounter = flipVCounter = 0;
972 return saveXISFImage(newFilename);
1011 nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
1014 if (fptr && fits_close_file(fptr, &status))
1016 m_LastError =
i18n(
"Failed to close file: %1", fitsErrorToString(status));
1021 if (fits_create_file(&new_fptr, QString(
"!%1").arg(newFilename).toLocal8Bit(), &status))
1023 m_LastError =
i18n(
"Failed to create file: %1", fitsErrorToString(status));
1032 long naxis = m_Statistics.channels == 1 ? 2 : 3;
1033 long naxes[3] = {m_Statistics.width, m_Statistics.height, naxis};
1036 if (fits_create_img(fptr, m_FITSBITPIX, naxis, naxes, &status))
1038 m_LastError =
i18n(
"Failed to create image: %1", fitsErrorToString(status));
1045 if (fits_update_key(fptr, TDOUBLE,
"DATAMIN", &(m_Statistics.min),
"Minimum value", &status))
1047 m_LastError =
i18n(
"Failed to update key: %1", fitsErrorToString(status));
1052 if (fits_update_key(fptr, TDOUBLE,
"DATAMAX", &(m_Statistics.max),
"Maximum value", &status))
1054 m_LastError =
i18n(
"Failed to update key: %1", fitsErrorToString(status));
1059 fits_write_key(fptr, TDOUBLE,
"MIN1", &m_Statistics.min[0],
"Min Channel 1", &status);
1062 fits_write_key(fptr, TDOUBLE,
"MIN2", &m_Statistics.min[1],
"Min Channel 2", &status);
1063 fits_write_key(fptr, TDOUBLE,
"MIN3", &m_Statistics.min[2],
"Min Channel 3", &status);
1067 fits_write_key(fptr, TDOUBLE,
"MAX1", &m_Statistics.max[0],
"Max Channel 1", &status);
1070 fits_write_key(fptr, TDOUBLE,
"MAX2", &m_Statistics.max[1],
"Max Channel 2", &status);
1071 fits_write_key(fptr, TDOUBLE,
"MAX3", &m_Statistics.max[2],
"Max Channel 3", &status);
1075 if (m_Statistics.mean[0] > 0)
1077 fits_write_key(fptr, TDOUBLE,
"MEAN1", &m_Statistics.mean[0],
"Mean Channel 1", &status);
1080 fits_write_key(fptr, TDOUBLE,
"MEAN2", &m_Statistics.mean[1],
"Mean Channel 2", &status);
1081 fits_write_key(fptr, TDOUBLE,
"MEAN3", &m_Statistics.mean[2],
"Mean Channel 3", &status);
1086 if (m_Statistics.median[0] > 0)
1088 fits_write_key(fptr, TDOUBLE,
"MEDIAN1", &m_Statistics.median[0],
"Median Channel 1", &status);
1091 fits_write_key(fptr, TDOUBLE,
"MEDIAN2", &m_Statistics.median[1],
"Median Channel 2", &status);
1092 fits_write_key(fptr, TDOUBLE,
"MEDIAN3", &m_Statistics.median[2],
"Median Channel 3", &status);
1097 if (m_Statistics.stddev[0] > 0)
1099 fits_write_key(fptr, TDOUBLE,
"STDDEV1", &m_Statistics.stddev[0],
"Standard Deviation Channel 1", &status);
1102 fits_write_key(fptr, TDOUBLE,
"STDDEV2", &m_Statistics.stddev[1],
"Standard Deviation Channel 2", &status);
1103 fits_write_key(fptr, TDOUBLE,
"STDDEV3", &m_Statistics.stddev[2],
"Standard Deviation Channel 3", &status);
1108 for (
int i = 10; i < m_HeaderRecords.count(); i++)
1110 QString key = m_HeaderRecords[i].key;
1111 const char *comment = m_HeaderRecords[i].comment.toLatin1().constBegin();
1112 QVariant value = m_HeaderRecords[i].
value;
1114 switch (value.
type())
1119 fits_write_key(fptr, TINT, key.
toLatin1().
constData(), &number, comment, &status);
1126 fits_write_key(fptr, TDOUBLE, key.
toLatin1().
constData(), &number, comment, &status);
1133 char valueBuffer[256] = {0};
1135 fits_write_key(fptr, TSTRING, key.
toLatin1().
constData(), valueBuffer, comment, &status);
1141 if (fits_write_date(fptr, &status))
1143 m_LastError =
i18n(
"Failed to update date: %1", fitsErrorToString(status));
1150 if (fits_write_history(fptr, history.
toLatin1(), &status))
1152 m_LastError =
i18n(
"Failed to update history: %1", fitsErrorToString(status));
1156 int rot = 0, mirror = 0;
1159 rot = (90 * rotCounter) % 360;
1163 if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0)
1166 if ((rot != 0) || (mirror != 0))
1167 rotWCSFITS(rot, mirror);
1169 rotCounter = flipHCounter = flipVCounter = 0;
1172 if (fits_write_img(fptr, m_Statistics.dataType, 1, nelements, m_ImageBuffer, &status))
1174 m_LastError =
i18n(
"Failed to write image: %1", fitsErrorToString(status));
1178 m_Filename = newFilename;
1180 fits_flush_file(fptr, &status);
1182 qCInfo(KSTARS_FITS) <<
"Saved FITS file:" << m_Filename;
1187void FITSData::clearImageBuffers()
1189 delete[] m_ImageBuffer;
1190 m_ImageBuffer =
nullptr;
1191 if(m_ImageRoiBuffer !=
nullptr )
1193 delete[] m_ImageRoiBuffer;
1194 m_ImageRoiBuffer =
nullptr;
1200void FITSData::makeRoiBuffer(
QRect roi)
1205 uint32_t channelSize = roi.
height() * roi.
width();
1206 if(channelSize > m_Statistics.samples_per_channel || channelSize == 1)
1210 if(m_ImageRoiBuffer !=
nullptr )
1212 delete[] m_ImageRoiBuffer;
1213 m_ImageRoiBuffer =
nullptr;
1215 int xoffset = roi.
topLeft().
x() - 1;
1216 int yoffset = roi.
topLeft().
y() - 1;
1217 uint32_t bpp = m_Statistics.bytesPerPixel;
1218 m_ImageRoiBuffer =
new uint8_t[roi.
height()*roi.
width()*m_Statistics.channels * m_Statistics.bytesPerPixel]();
1219 for(
int n = 0 ; n < m_Statistics.channels ; n++)
1221 for(
int i = 0; i < roi.
height(); i++)
1223 size_t i1 = n * channelSize * bpp + i * roi.
width() * bpp;
1224 size_t i2 = n * m_Statistics.samples_per_channel * bpp + (yoffset + i) * width() * bpp + xoffset * bpp;
1225 memcpy(&m_ImageRoiBuffer[i1],
1231 memcpy(&m_ROIStatistics, &m_Statistics,
sizeof(FITSImage::Statistic));
1232 m_ROIStatistics.samples_per_channel = roi.
height() * roi.
width();
1233 m_ROIStatistics.width = roi.
width();
1234 m_ROIStatistics.height = roi.
height();
1235 calculateStats(
false,
true);
1238void FITSData::setupWCSParams()
1240 FITSImage::Solution solution;
1241 if (parseSolution(solution))
1243 const bool eastToTheRight = solution.parity == FITSImage::POSITIVE ? false :
true;
1244 updateWCSHeaderData(solution.orientation, solution.ra, solution.dec, solution.pixscale, eastToTheRight);
1247 bool validObservationDate =
false;
1249 if (getRecordValue(
"DATE-OBS", value))
1251 QString tsString(value.
toString());
1252 tsString = tsString.remove(
'\'').trimmed();
1260 validObservationDate =
true;
1261 m_DateTime = KStarsDateTime(ts.
date(), ts.
time());
1266 if (!validObservationDate)
1271void FITSData::calculateStats(
bool refresh,
bool roi)
1276 calculateMinMax(refresh);
1277 calculateMedian(refresh);
1280 if (refresh ==
false && fptr)
1284 if (fits_read_key_dbl(fptr,
"MEAN1", &m_Statistics.mean[0],
nullptr, &status) == 0)
1287 fits_read_key_dbl(fptr,
"MEAN2", & m_Statistics.mean[1],
nullptr, &status);
1288 fits_read_key_dbl(fptr,
"MEAN3", &m_Statistics.mean[2],
nullptr, &status);
1291 if (fits_read_key_dbl(fptr,
"STDDEV1", &m_Statistics.stddev[0],
nullptr, &status) == 0)
1294 fits_read_key_dbl(fptr,
"STDDEV2", &m_Statistics.stddev[1],
nullptr, &status);
1295 fits_read_key_dbl(fptr,
"STDDEV3", &m_Statistics.stddev[2],
nullptr, &status);
1303 switch (m_Statistics.dataType)
1306 calculateStdDev<uint8_t>();
1310 calculateStdDev<int16_t>();
1314 calculateStdDev<uint16_t>();
1318 calculateStdDev<int32_t>();
1322 calculateStdDev<uint32_t>();
1326 calculateStdDev<float>();
1330 calculateStdDev<int64_t>();
1334 calculateStdDev<double>();
1342 m_Statistics.SNR = m_Statistics.mean[0] / m_Statistics.stddev[0];
1346 calculateMinMax(refresh, roi);
1347 calculateMedian(refresh, roi);
1349 switch (m_ROIStatistics.dataType)
1352 calculateStdDev<uint8_t>(roi);
1356 calculateStdDev<int16_t>(roi);
1360 calculateStdDev<uint16_t>(roi);
1364 calculateStdDev<int32_t>(roi);
1368 calculateStdDev<uint32_t>(roi);
1372 calculateStdDev<float>(roi);
1376 calculateStdDev<int64_t>(roi);
1380 calculateStdDev<double>(roi);
1389void FITSData::calculateMinMax(
bool refresh,
bool roi)
1399 if (fptr !=
nullptr && !refresh)
1404 if (fits_read_key_dbl(fptr,
"DATAMIN", &(m_Statistics.min[0]),
nullptr, &status) == 0)
1406 else if (fits_read_key_dbl(fptr,
"MIN1", &(m_Statistics.min[0]),
nullptr, &status) == 0)
1410 fits_read_key_dbl(fptr,
"MIN2", &m_Statistics.min[1],
nullptr, &status);
1411 fits_read_key_dbl(fptr,
"MIN3", &m_Statistics.min[2],
nullptr, &status);
1415 if (fits_read_key_dbl(fptr,
"DATAMAX", &(m_Statistics.max[0]),
nullptr, &status) == 0)
1417 else if (fits_read_key_dbl(fptr,
"MAX1", &(m_Statistics.max[0]),
nullptr, &status) == 0)
1421 fits_read_key_dbl(fptr,
"MAX2", &m_Statistics.max[1],
nullptr, &status);
1422 fits_read_key_dbl(fptr,
"MAX3", &m_Statistics.max[2],
nullptr, &status);
1425 if (nfound == 2 && !(m_Statistics.min[0] == 0 && m_Statistics.max[0] == 0))
1429 m_Statistics.min[0] = 1.0E30;
1430 m_Statistics.max[0] = -1.0E30;
1432 m_Statistics.min[1] = 1.0E30;
1433 m_Statistics.max[1] = -1.0E30;
1435 m_Statistics.min[2] = 1.0E30;
1436 m_Statistics.max[2] = -1.0E30;
1438 switch (m_Statistics.dataType)
1441 calculateMinMax<uint8_t>();
1445 calculateMinMax<int16_t>();
1449 calculateMinMax<uint16_t>();
1453 calculateMinMax<int32_t>();
1457 calculateMinMax<uint32_t>();
1461 calculateMinMax<float>();
1465 calculateMinMax<int64_t>();
1469 calculateMinMax<double>();
1478 m_ROIStatistics.min[0] = 1.0E30;
1479 m_ROIStatistics.max[0] = -1.0E30;
1481 m_ROIStatistics.min[1] = 1.0E30;
1482 m_ROIStatistics.max[1] = -1.0E30;
1484 m_ROIStatistics.min[2] = 1.0E30;
1485 m_ROIStatistics.max[2] = -1.0E30;
1487 switch (m_Statistics.dataType)
1490 calculateMinMax<uint8_t>(roi);
1494 calculateMinMax<int16_t>(roi);
1498 calculateMinMax<uint16_t>(roi);
1502 calculateMinMax<int32_t>(roi);
1506 calculateMinMax<uint32_t>(roi);
1510 calculateMinMax<float>(roi);
1514 calculateMinMax<int64_t>(roi);
1518 calculateMinMax<double>(roi);
1528void FITSData::calculateMedian(
bool refresh,
bool roi)
1538 if (fptr !=
nullptr && !refresh)
1541 if (fits_read_key_dbl(fptr,
"MEDIAN1", &m_Statistics.median[0],
nullptr, &status) == 0)
1545 fits_read_key_dbl(fptr,
"MEDIAN2", &m_Statistics.median[1],
nullptr, &status);
1546 fits_read_key_dbl(fptr,
"MEDIAN3", &m_Statistics.median[2],
nullptr, &status);
1552 m_Statistics.median[RED_CHANNEL] = 0;
1553 m_Statistics.median[GREEN_CHANNEL] = 0;
1554 m_Statistics.median[BLUE_CHANNEL] = 0;
1556 switch (m_Statistics.dataType)
1559 calculateMedian<uint8_t>();
1563 calculateMedian<int16_t>();
1567 calculateMedian<uint16_t>();
1571 calculateMedian<int32_t>();
1575 calculateMedian<uint32_t>();
1579 calculateMedian<float>();
1583 calculateMedian<int64_t>();
1587 calculateMedian<double>();
1596 m_ROIStatistics.median[RED_CHANNEL] = 0;
1597 m_ROIStatistics.median[GREEN_CHANNEL] = 0;
1598 m_ROIStatistics.median[BLUE_CHANNEL] = 0;
1600 switch (m_ROIStatistics.dataType)
1603 calculateMedian<uint8_t>(roi);
1607 calculateMedian<int16_t>(roi);
1611 calculateMedian<uint16_t>(roi);
1615 calculateMedian<int32_t>(roi);
1619 calculateMedian<uint32_t>(roi);
1623 calculateMedian<float>(roi);
1627 calculateMedian<int64_t>(roi);
1631 calculateMedian<double>(roi);
1641template <
typename T>
1642void FITSData::calculateMedian(
bool roi)
1644 auto * buffer =
reinterpret_cast<T *
>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1645 const uint32_t maxMedianSize = 500000;
1646 uint32_t medianSize = roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel;
1647 uint8_t downsample = 1;
1648 if (medianSize > maxMedianSize)
1650 downsample = (
static_cast<double>(medianSize) / maxMedianSize) + 0.999;
1651 medianSize /= downsample;
1657 std::vector<int32_t> samples;
1658 samples.reserve(medianSize);
1660 for (uint8_t n = 0; n < m_Statistics.channels; n++)
1662 auto *oneChannel = buffer + n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1663 for (uint32_t upto = 0; upto < (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1666 auto median = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_MEDIAN, samples);
1667 roi ? m_ROIStatistics.median[n] = median : m_Statistics.median[n] = median;
1671template <
typename T>
1672QPair<T, T> FITSData::getParitionMinMax(uint32_t start, uint32_t stride,
bool roi)
1674 auto * buffer =
reinterpret_cast<T *
>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1675 T min = std::numeric_limits<T>::max();
1676 T max = std::numeric_limits<T>::min();
1678 uint32_t
end = start + stride;
1680 for (uint32_t i = start; i <
end; i++)
1682 min = qMin(buffer[i], min);
1683 max = qMax(buffer[i], max);
1690 return qMakePair(min, max);
1693template <
typename T>
1694void FITSData::calculateMinMax(
bool roi)
1696 T min = std::numeric_limits<T>::max();
1697 T max = std::numeric_limits<T>::min();
1701 const uint8_t nThreads = 16;
1703 for (
int n = 0; n < m_Statistics.channels; n++)
1705 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1708 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1711 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1712 (tStride * nThreads));
1715 uint32_t tStart = cStart;
1718 QList<QFuture<QPair<T, T>>> futures;
1720 for (
int i = 0; i < nThreads; i++)
1723#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
1724 futures.
append(
QtConcurrent::run(&FITSData::getParitionMinMax<T>,
this, tStart, (i == (nThreads - 1)) ? fStride : tStride,
1727 futures.
append(
QtConcurrent::run(
this, &FITSData::getParitionMinMax<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride,
1734 for (
int i = 0; i < nThreads; i++)
1736 QPair<T, T> result = futures[i].result();
1737 min = qMin(result.first, min);
1738 max = qMax(result.second, max);
1743 m_Statistics.min[n] = min;
1744 m_Statistics.max[n] = max;
1748 m_ROIStatistics.min[n] = min;
1749 m_ROIStatistics.max[n] = max;
1762 SumData(
double s,
double sq,
int n) : sum(s), squaredSum(sq), numSamples(n) {}
1763 SumData() : sum(0), squaredSum(0), numSamples(0) {}
1766template <
typename T>
1767SumData getSumAndSquaredSum(uint32_t start, uint32_t stride, uint8_t *buff)
1769 auto * buffer =
reinterpret_cast<T *
>(buff);
1770 const uint32_t
end = start + stride;
1772 double squaredSum = 0;
1773 for (uint32_t i = start; i <
end; i++)
1775 double sample = buffer[i];
1777 squaredSum += sample * sample;
1779 const double numSamples =
end - start;
1780 return SumData(sum, squaredSum, numSamples);
1783template <
typename T>
1784void FITSData::calculateStdDev(
bool roi )
1787 const uint8_t nThreads = 16;
1789 for (
int n = 0; n < m_Statistics.channels; n++)
1791 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1794 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1797 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1798 (tStride * nThreads));
1801 uint32_t tStart = cStart;
1804 QList<QFuture<SumData>> futures;
1806 for (
int i = 0; i < nThreads; i++)
1809 uint8_t *buff = roi ? m_ImageRoiBuffer : m_ImageBuffer;
1811 (i == (nThreads - 1)) ? fStride : tStride, buff));
1816 double sum = 0, squared_sum = 0;
1817 double numSamples = 0;
1818 for (
int i = 0; i < nThreads; i++)
1820 const SumData result = futures[i].result();
1822 squared_sum += result.squaredSum;
1823 numSamples += result.numSamples;
1826 if (numSamples <= 0)
continue;
1827 const double mean = sum / numSamples;
1828 const double variance = squared_sum / numSamples - mean * mean;
1831 m_Statistics.mean[n] = mean;
1832 m_Statistics.stddev[n] = sqrt(variance);
1836 m_ROIStatistics.mean[n] = mean;
1837 m_ROIStatistics.stddev[n] = sqrt(variance);
1844 QVector<double> kernel(size * size);
1845 kernel.fill(0.0, size * size);
1847 double kernelSum = 0.0;
1848 int fOff = (size - 1) / 2;
1849 double normal = 1.0 / (2.0 * M_PI * sigma * sigma);
1850 for (
int y = -fOff; y <= fOff; y++)
1852 for (
int x = -fOff; x <= fOff; x++)
1854 double distance = ((y * y) + (x * x)) / (2.0 * sigma * sigma);
1855 int index = (y + fOff) * size + (x + fOff);
1856 kernel[index] = normal * qExp(-distance);
1857 kernelSum += kernel.at(index);
1860 for (
int y = 0; y < size; y++)
1862 for (
int x = 0; x < size; x++)
1864 int index = y * size + x;
1865 kernel[index] = kernel.at(index) * 1.0 / kernelSum;
1872template <
typename T>
1873void FITSData::convolutionFilter(
const QVector<double> &kernel,
int kernelSize)
1875 T * imagePtr =
reinterpret_cast<T *
>(m_ImageBuffer);
1881 int fOff = (kernelSize - 1) / 2;
1885 for (
int offsetY = 0; offsetY < m_Statistics.height; offsetY++)
1887 for (
int offsetX = 0; offsetX < m_Statistics.width; offsetX++)
1892 int byteOffset = offsetY * m_Statistics.width + offsetX;
1895 for (
int filterY = -fOff; filterY <= fOff; filterY++)
1897 for (
int filterX = -fOff; filterX <= fOff; filterX++)
1899 if ((offsetY + filterY) >= 0 && (offsetY + filterY) < m_Statistics.height
1900 && ((offsetX + filterX) >= 0 && (offsetX + filterX) < m_Statistics.width ))
1903 int calcOffset = byteOffset + filterX + filterY * m_Statistics.width;
1904 int index = (filterY + fOff) * kernelSize + (filterX + fOff);
1905 double kernelValue = kernel.
at(index);
1906 gt += (imagePtr[calcOffset]) * kernelValue;
1912 imagePtr[byteOffset] = gt;
1917template <
typename T>
1918void FITSData::gaussianBlur(
int kernelSize,
double sigma)
1921 if (kernelSize % 2 == 0)
1924 qCInfo(KSTARS_FITS) <<
"Warning, size must be an odd number, correcting size to " << kernelSize;
1932 QVector<double> gaussianKernel = createGaussianKernel(kernelSize, sigma);
1933 convolutionFilter<T>(gaussianKernel, kernelSize);
1936void FITSData::setMinMax(
double newMin,
double newMax, uint8_t channel)
1938 m_Statistics.min[channel] = newMin;
1939 m_Statistics.max[channel] = newMax;
1942bool FITSData::parseHeader()
1944 char * header =
nullptr;
1945 int status = 0, nkeys = 0;
1947 if (fits_hdr2str(fptr, 0,
nullptr, 0, &header, &nkeys, &status))
1949 fits_report_error(stderr, status);
1954 m_HeaderRecords.clear();
1955 QString recordList = QString(header);
1957 for (
int i = 0; i < nkeys; i++)
1961 QString record = recordList.
mid(i * 80, 80).
remove(
"'");
1967 oneRecord.comment =
properties[0].mid(8).simplified();
1972 oneRecord.value =
properties[1].simplified();
1974 oneRecord.comment =
properties[2].simplified();
1981 oneRecord.value.toInt(&ok);
1987 oneRecord.value.toDouble(&ok);
1993 m_HeaderRecords.
append(oneRecord);
2001bool FITSData::getRecordValue(
const QString &key,
QVariant &value)
const
2003 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](
const Record & oneRecord)
2005 return (oneRecord.key == key && oneRecord.value.isValid());
2008 if (result != m_HeaderRecords.end())
2010 value = (*result).
value;
2018 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](
const Record & oneRecord)
2020 return (oneRecord.key == key && oneRecord.value.isValid());
2023 if (result != m_HeaderRecords.end())
2025 (*result).value = value;
2030 FITSData::Record record = {key, value.
toString(), comment};
2031#if QT_VERSION >= QT_VERSION_CHECK(6, 0, 0)
2032 m_HeaderRecords.
insert(std::max(0LL, m_HeaderRecords.size() - 1), record);
2034 m_HeaderRecords.
insert(std::max(0, m_HeaderRecords.size() - 1), record);
2039bool FITSData::parseSolution(FITSImage::Solution &solution)
const
2042 bool raOK =
false, deOK =
false, coordOK =
false, scaleOK =
false;
2046 solution.fieldWidth = solution.fieldHeight = solution.pixscale = solution.ra = solution.dec = -1;
2049 if (getRecordValue(
"OBJCTRA", value))
2052 solution.ra = angleValue.
Degrees();
2055 else if (getRecordValue(
"RA", value))
2057 solution.ra = value.
toDouble(&raOK);
2061 if (getRecordValue(
"OBJCTDEC", value))
2064 solution.dec = angleValue.
Degrees();
2067 else if (getRecordValue(
"DEC", value))
2069 solution.dec = value.
toDouble(&deOK);
2072 coordOK = raOK && deOK;
2076 if (getRecordValue(
"SCALE", value))
2081 double focal_length = -1;
2082 if (getRecordValue(
"FOCALLEN", value))
2087 double pixsize1 = -1, pixsize2 = -1;
2089 if (getRecordValue(
"PIXSIZE1", value))
2093 else if (getRecordValue(
"XPIXSZ", value))
2098 if (getRecordValue(
"PIXSIZE2", value))
2102 else if (getRecordValue(
"YPIXSZ", value))
2107 int binx = 1, biny = 1;
2109 if (getRecordValue(
"XBINNING", value))
2114 if (getRecordValue(
"YBINNING", value))
2119 if (pixsize1 > 0 && pixsize2 > 0)
2125 solution.pixscale = scale;
2127 solution.fieldWidth = (m_Statistics.width * scale) / 60.0;
2129 solution.fieldHeight = (m_Statistics.height * scale * (pixsize2 / pixsize1)) / 60.0;
2131 else if (focal_length > 0)
2134 solution.fieldWidth = ((206264.8062470963552 * m_Statistics.width * (pixsize1 / 1000.0)) / (focal_length * binx)) / 60.0;
2136 solution.fieldHeight = ((206264.8062470963552 * m_Statistics.height * (pixsize2 / 1000.0)) / (focal_length * biny)) / 60.0;
2138 solution.pixscale = (206264.8062470963552 * (pixsize1 / 1000.0)) / (focal_length * binx);
2144 return (coordOK || scaleOK);
2149 if (m_StarFindFuture.isRunning())
2152 starAlgorithm = algorithm;
2153 qDeleteAll(starCenters);
2154 starCenters.clear();
2155 starsSearched =
true;
2161 m_StarDetector.reset(
new FITSSEPDetector(
this));
2162 m_StarDetector->setSettings(m_SourceExtractorSettings);
2163 if (m_Mode == FITS_NORMAL && trackingBox.
isNull())
2165 if (Options::quickHFR())
2168 const int w = getStatistics().width;
2169 const int h = getStatistics().height;
2170 QRect middle(
static_cast<int>(w * 0.25),
static_cast<int>(h * 0.25), w / 2, h / 2);
2171 m_StarFindFuture = m_StarDetector->findSources(middle);
2172 return m_StarFindFuture;
2175 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2176 return m_StarFindFuture;
2179 case ALGORITHM_GRADIENT:
2182 m_StarDetector.reset(
new FITSGradientDetector(
this));
2183 m_StarDetector->setSettings(m_SourceExtractorSettings);
2184 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2185 return m_StarFindFuture;
2188 case ALGORITHM_CENTROID:
2191 m_StarDetector.reset(
new FITSCentroidDetector(
this));
2192 m_StarDetector->setSettings(m_SourceExtractorSettings);
2194 if (!isHistogramConstructed())
2195 constructHistogram();
2196 m_StarDetector->configure(
"JMINDEX", m_JMIndex);
2197 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2198 return m_StarFindFuture;
2202 m_StarDetector.reset(
new FITSCentroidDetector(
this));
2203 m_StarFindFuture = starDetector->findSources(trackingBox);
2204 return m_StarFindFuture;
2208 case ALGORITHM_THRESHOLD:
2210 m_StarDetector.reset(
new FITSThresholdDetector(
this));
2211 m_StarDetector->setSettings(m_SourceExtractorSettings);
2212 m_StarDetector->configure(
"THRESHOLD_PERCENTAGE", Options::focusThreshold());
2213 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2214 return m_StarFindFuture;
2217 case ALGORITHM_BAHTINOV:
2219 m_StarDetector.reset(
new FITSBahtinovDetector(
this));
2220 m_StarDetector->setSettings(m_SourceExtractorSettings);
2221 m_StarDetector->configure(
"NUMBER_OF_AVERAGE_ROWS", Options::focusMultiRowAverage());
2222 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2223 return m_StarFindFuture;
2230 if (mask.
isNull() ==
false)
2232 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(),
2235 return (mask->isVisible(edge->x, edge->y) == false);
2236 }), starCenters.end());
2239 return starCenters.count();
2243double FITSData::getHFR(HFRType type)
2245 if (starCenters.empty())
2248 if (cacheHFR >= 0 && cacheHFRType == type)
2251 m_SelectedHFRStar.invalidate();
2253 if (type == HFR_MAX)
2258 for (
int i = 0; i < starCenters.count(); i++)
2260 if (starCenters[i]->HFR > maxVal)
2263 maxVal = starCenters[i]->HFR;
2267 m_SelectedHFRStar = *starCenters[maxIndex];
2268 cacheHFR = starCenters[maxIndex]->HFR;
2269 cacheHFRType =
type;
2272 else if (type == HFR_HIGH)
2275 int minX = width() / 10;
2276 int minY = height() / 10;
2277 int maxX = width() - minX;
2278 int maxY = height() - minY;
2279 starCenters.erase(std::remove_if(starCenters.begin(), starCenters.end(), [minX, minY, maxX, maxY](Edge * oneStar)
2281 return (oneStar->x < minX || oneStar->x > maxX || oneStar->y < minY || oneStar->y > maxY);
2282 }), starCenters.end());
2284 if (starCenters.empty())
2287 m_SelectedHFRStar = *starCenters[
static_cast<int>(starCenters.size() * 0.05)];
2288 cacheHFR = m_SelectedHFRStar.HFR;
2289 cacheHFRType =
type;
2292 else if (type == HFR_MEDIAN)
2294 std::nth_element(starCenters.begin(), starCenters.begin() + starCenters.size() / 2, starCenters.end());
2295 m_SelectedHFRStar = *starCenters[starCenters.size() / 2];
2297 cacheHFR = m_SelectedHFRStar.HFR;
2298 cacheHFRType =
type;
2306 const int saturationValue = m_Statistics.dataType == TBYTE ? 250 : 50000;
2307 int numSaturated = 0;
2308 for (
auto center : starCenters)
2309 if (
center->val > saturationValue)
2311 const bool removeSaturatedStars = starCenters.size() - numSaturated > 20;
2312 if (removeSaturatedStars && numSaturated > 0)
2313 qCDebug(KSTARS_FITS) <<
"Removing " << numSaturated <<
" stars from HFR calculation";
2315 std::vector<double> HFRs;
2317 for (
auto center : starCenters)
2319 if (removeSaturatedStars &&
center->val > saturationValue)
continue;
2321 if (type == HFR_AVERAGE)
2322 HFRs.push_back(
center->HFR);
2328 if (m_SkyBackground.mean <= 0.0 ||
center->val < m_SkyBackground.mean)
2330 HFRs.push_back(
center->HFR);
2331 qCDebug(KSTARS_FITS) <<
"HFR Adj, sky background " << m_SkyBackground.mean <<
" star peak " <<
center->val <<
2336 const double a_div_b =
center->val / m_SkyBackground.mean;
2337 const double factor = erf(sqrt(log(a_div_b))) / (1 - (1 / a_div_b));
2338 HFRs.push_back(
center->HFR * factor);
2339 qCDebug(KSTARS_FITS) <<
"HFR Adj, brightness adjusted from " <<
center->HFR <<
" to " <<
center->HFR * factor;
2344 auto m = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_SIGMACLIPPING,
2348 cacheHFRType = HFR_AVERAGE;
2352double FITSData::getHFR(
int x,
int y,
double scale)
2354 if (starCenters.empty())
2357 for (
int i = 0; i < starCenters.count(); i++)
2359 const int maxDist = std::max(2,
static_cast<int>(0.5 + 2 * starCenters[i]->width / scale));
2360 const int dx = std::fabs(starCenters[i]->x - x);
2361 const int dy = std::fabs(starCenters[i]->y - y);
2362 if (dx <= maxDist && dy <= maxDist)
2364 m_SelectedHFRStar = *starCenters[i];
2365 return starCenters[i]->HFR;
2372double FITSData::getEccentricity()
2374 if (starCenters.empty())
2376 if (cacheEccentricity >= 0)
2377 return cacheEccentricity;
2378 std::vector<float> eccs;
2379 for (
const auto &s : starCenters)
2380 eccs.push_back(s->ellipticity);
2381 int middle = eccs.size() / 2;
2382 std::nth_element(eccs.begin(), eccs.begin() + middle, eccs.end());
2383 float medianEllipticity = eccs[middle];
2389 const float eccentricity = sqrt(medianEllipticity * (2 - medianEllipticity));
2390 cacheEccentricity = eccentricity;
2391 return eccentricity;
2396 if (type == FITS_NONE)
2399 QVector<double> dataMin(3);
2400 QVector<double> dataMax(3);
2409 case FITS_AUTO_STRETCH:
2411 for (
int i = 0; i < 3; i++)
2413 dataMin[i] = m_Statistics.mean[i] - m_Statistics.stddev[i];
2414 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2419 case FITS_HIGH_CONTRAST:
2421 for (
int i = 0; i < 3; i++)
2423 dataMin[i] = m_Statistics.mean[i] + m_Statistics.stddev[i];
2424 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2429 case FITS_HIGH_PASS:
2431 for (
int i = 0; i < 3; i++)
2433 dataMin[i] = m_Statistics.mean[i];
2442 switch (m_Statistics.dataType)
2446 for (
int i = 0; i < 3; i++)
2448 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2449 dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i];
2451 applyFilter<uint8_t>(type, image, &dataMin, &dataMax);
2457 for (
int i = 0; i < 3; i++)
2459 dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i];
2460 dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i];
2462 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2469 for (
int i = 0; i < 3; i++)
2471 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2472 dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i];
2474 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2480 for (
int i = 0; i < 3; i++)
2482 dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i];
2483 dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i];
2485 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2491 for (
int i = 0; i < 3; i++)
2493 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2494 dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i];
2496 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2502 for (
int i = 0; i < 3; i++)
2504 dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i];
2505 dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i];
2507 applyFilter<float>(type, image, &dataMin, &dataMax);
2513 for (
int i = 0; i < 3; i++)
2515 dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i];
2516 dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i];
2519 applyFilter<long>(type, image, &dataMin, &dataMax);
2525 for (
int i = 0; i < 3; i++)
2527 dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i];
2528 dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i];
2530 applyFilter<double>(type, image, &dataMin, &dataMax);
2547template <
typename T>
2550 bool calcStats =
false;
2551 T * image =
nullptr;
2554 image =
reinterpret_cast<T *
>(targetImage);
2557 image =
reinterpret_cast<T *
>(m_ImageBuffer);
2562 for (
int i = 0; i < 3; i++)
2564 min[i] = (*targetMin)[i] < std::numeric_limits<T>::min() ? std::numeric_limits<T>::min() : (*targetMin)[i];
2565 max[i] = (*targetMax)[i] > std::numeric_limits<T>::max() ? std::numeric_limits<T>::max() : (*targetMax)[i];
2570 const uint8_t nThreads = 16;
2572 uint32_t width = m_Statistics.width;
2573 uint32_t height = m_Statistics.height;
2581 case FITS_AUTO_STRETCH:
2582 case FITS_HIGH_CONTRAST:
2585 case FITS_HIGH_PASS:
2588 QList<QFuture<void>> futures;
2589 QVector<double> coeff(3);
2591 if (type == FITS_LOG)
2593 for (
int i = 0; i < 3; i++)
2594 coeff[i] = max[i] / std::log(1 + max[i]);
2596 else if (type == FITS_SQRT)
2598 for (
int i = 0; i < 3; i++)
2599 coeff[i] = max[i] / sqrt(max[i]);
2602 for (
int n = 0; n < m_Statistics.channels; n++)
2604 if (type == FITS_HIGH_PASS)
2605 min[n] = m_Statistics.mean[n];
2607 uint32_t cStart = n * m_Statistics.samples_per_channel;
2610 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
2613 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
2615 T * runningBuffer = image + cStart;
2617 if (type == FITS_LOG)
2619 for (
int i = 0; i < nThreads; i++)
2622 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2625 a = qBound(min[n],
static_cast<T
>(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]);
2628 runningBuffer += tStride;
2631 else if (type == FITS_SQRT)
2633 for (
int i = 0; i < nThreads; i++)
2636 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2639 a = qBound(min[n],
static_cast<T
>(round(coeff[n] * a)), max[n]);
2643 runningBuffer += tStride;
2647 for (
int i = 0; i < nThreads; i++)
2650 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2653 a = qBound(min[n], a, max[n]);
2655 runningBuffer += tStride;
2660 for (
int i = 0; i < nThreads * m_Statistics.channels; i++)
2661 futures[i].waitForFinished();
2665 for (
int i = 0; i < 3; i++)
2667 m_Statistics.min[i] = min[i];
2668 m_Statistics.max[i] = max[i];
2670 calculateStdDev<T>();
2678 if (!isHistogramConstructed())
2679 constructHistogram();
2683 double coeff = 255.0 / (height * width);
2687 for (
int i = 0; i < m_Statistics.channels; i++)
2689 uint32_t offset = i * m_Statistics.samples_per_channel;
2690 for (uint32_t j = 0; j < height; j++)
2692 row = offset + j * width;
2693 for (uint32_t k = 0; k < width; k++)
2696 bufferVal = (image[index] - min[i]) / m_HistogramBinWidth[i];
2698 if (bufferVal >= m_CumulativeFrequency[i].size())
2699 bufferVal = m_CumulativeFrequency[i].size() - 1;
2701 image[index] = qBound(min[i],
static_cast<T
>(round(coeff * m_CumulativeFrequency[i][bufferVal])), max[i]);
2708 calculateStats(
true,
false);
2714 uint8_t BBP = m_Statistics.bytesPerPixel;
2715 auto * extension =
new T[(width + 2) * (height + 2)];
2720 for (uint32_t ch = 0; ch < m_Statistics.channels; ch++)
2722 uint32_t offset = ch * m_Statistics.samples_per_channel;
2723 uint32_t N = width, M = height;
2725 for (uint32_t i = 0; i < M; ++i)
2727 memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP);
2728 extension[(N + 2) * (i + 1)] = image[N * i + offset];
2729 extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset];
2732 memcpy(extension, extension + N + 2, (N + 2) * BBP);
2734 memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP);
2741 for (uint32_t m = 1; m < M - 1; ++m)
2742 for (uint32_t n = 1; n < N - 1; ++n)
2748 memset(&window[0], 0, 9 *
sizeof(
float));
2749 for (uint32_t j = m - 1; j < m + 2; ++j)
2750 for (uint32_t i = n - 1; i < n + 2; ++i)
2751 window[k++] = extension[j * N + i];
2753 for (uint32_t j = 0; j < 5; ++j)
2757 for (uint32_t l = j + 1; l < 9; ++l)
2758 if (window[l] < window[mine])
2761 const float temp =
window[j];
2766 image[(m - 1) * (N - 2) + n - 1 + offset] =
window[4];
2774 calculateStdDev<T>();
2779 gaussianBlur<T>(Options::focusGaussianKernelSize(), Options::focusGaussianSigma());
2781 calculateStats(
true,
false);
2784 case FITS_ROTATE_CW:
2789 case FITS_ROTATE_CCW:
2794 case FITS_MOUNT_FLIP_H:
2799 case FITS_MOUNT_FLIP_V:
2811 QList<Edge *> starCentersInSubFrame;
2812 for (
int i = 0; i < starCenters.count(); i++)
2814 int x =
static_cast<int>(starCenters[i]->x);
2815 int y =
static_cast<int>(starCenters[i]->y);
2818 starCentersInSubFrame.
append(starCenters[i]);
2821 return starCentersInSubFrame;
2824bool FITSData::loadWCS()
2826#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2828 if (m_WCSState == Success)
2830 qCWarning(KSTARS_FITS) <<
"WCS data already loaded";
2834 if (m_WCSHandle !=
nullptr)
2836 wcsvfree(&m_nwcs, &m_WCSHandle);
2838 m_WCSHandle =
nullptr;
2841 qCDebug(KSTARS_FITS) <<
"Started WCS Data Processing...";
2843 QByteArray header_str;
2845 int nkeyrec = 0, nreject = 0;
2848 char *header =
nullptr;
2849 if (fits_hdr2str(fptr, 1,
nullptr, 0, &header, &nkeyrec, &status))
2852 fits_get_errstatus(status, errmsg);
2853 m_LastError = errmsg;
2854 m_WCSState = Failure;
2857 header_str = QByteArray(header);
2858 fits_free_memory(header, &status);
2863 for(
auto &fitsKeyword : m_HeaderRecords)
2866 rec.
append(fitsKeyword.key.leftJustified(8,
' ').toLatin1());
2868 rec.
append(fitsKeyword.value.toByteArray());
2870 rec.
append(fitsKeyword.comment.toLatin1());
2874 header_str.
append(QByteArray(
"END").leftJustified(80));
2877 if ((status = wcspih(header_str.
data(), nkeyrec, WCSHDR_all, 0, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2879 wcsvfree(&m_nwcs, &m_WCSHandle);
2880 m_WCSHandle =
nullptr;
2882 m_LastError = QString(
"wcspih ERROR %1: %2.").arg(status).arg(wcshdr_errmsg[status]);
2883 m_WCSState = Failure;
2887 if (m_WCSHandle ==
nullptr)
2889 m_LastError =
i18n(
"No world coordinate systems found.");
2890 m_WCSState = Failure;
2895 if (m_WCSHandle->crpix[0] == 0)
2897 wcsvfree(&m_nwcs, &m_WCSHandle);
2898 m_WCSHandle =
nullptr;
2900 m_LastError =
i18n(
"No world coordinate systems found.");
2901 m_WCSState = Failure;
2906 if ((status = wcsset(m_WCSHandle)) != 0)
2908 wcsvfree(&m_nwcs, &m_WCSHandle);
2909 m_WCSHandle =
nullptr;
2911 m_LastError = QString(
"wcsset error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2912 m_WCSState = Failure;
2916 m_ObjectsSearched =
false;
2917 m_CatObjectsSearched =
false;
2918 m_WCSState = Success;
2921 qCDebug(KSTARS_FITS) <<
"Finished WCS Data processing...";
2931#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2933 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2935 if (m_WCSHandle ==
nullptr)
2937 m_LastError =
i18n(
"No world coordinate systems found.");
2944 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2946 m_LastError = QString(
"wcss2p error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2950 wcsImagePoint.
setX(imgcrd[0]);
2951 wcsImagePoint.
setY(imgcrd[1]);
2953 wcsPixelPoint.
setX(pixcrd[0]);
2954 wcsPixelPoint.
setY(pixcrd[1]);
2959 Q_UNUSED(wcsPixelPoint);
2960 Q_UNUSED(wcsImagePoint);
2965bool FITSData::pixelToWCS(
const QPointF &wcsPixelPoint,
SkyPoint &wcsCoord)
2967#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2969 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2971 if (m_WCSHandle ==
nullptr)
2973 m_LastError =
i18n(
"No world coordinate systems found.");
2977 pixcrd[0] = wcsPixelPoint.
x();
2978 pixcrd[1] = wcsPixelPoint.
y();
2980 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2982 m_LastError = QString(
"wcsp2s error %1: %2.").arg(status).arg(wcs_errmsg[status]);
2987 wcsCoord.
setRA0(world[0] / 15.0);
2993 Q_UNUSED(wcsPixelPoint);
2999#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3000bool FITSData::searchObjects()
3002 return (Options::fitsCatalog() == CAT_SKYMAP) ? searchSkyMapObjects() : searchCatObjects();
3005bool FITSData::searchSkyMapObjects()
3007 if (m_ObjectsSearched)
3010 m_ObjectsSearched =
true;
3012 SkyPoint startPoint;
3015 pixelToWCS(QPointF(0, 0), startPoint);
3016 pixelToWCS(QPointF(width() - 1, height() - 1), endPoint);
3018 return findObjectsInImage(startPoint, endPoint);
3021bool FITSData::searchCatObjects()
3023 if (m_CatObjectsSearched)
3026 m_CatObjectsSearched =
true;
3028 SkyPoint searchCenter;
3032 if (catROIRadius() > 0)
3036 QPoint edgePt = QPoint(pt.
x() + catROIRadius(), pt.
y());
3037 SkyPoint searchEdge;
3038 ok = pixelToWCS(pt, searchCenter);
3040 ok = pixelToWCS(edgePt, searchEdge);
3044 KSNumbers * num =
nullptr;
3046 if (getRecordValue(
"DATE-OBS", date))
3049 tsString = tsString.remove(
'\'').trimmed();
3056 num =
new KSNumbers(KStarsDateTime(ts).djd());
3061 num =
new KSNumbers(KStarsData::Instance()->ut().djd());
3074 pt = QPoint((width() / 2) - 1, (height() / 2) - 1);
3075 ok = pixelToWCS(pt, searchCenter);
3079 double raEdge = searchCenter.
ra0().
Degrees() + (radius / 60.0);
3081 SkyPoint searchEdge(raEdge / 15.0, decEdge);
3082 QPointF edgePoint, pEdge;
3083 ok = wcsToPixel(searchEdge, pEdge, edgePoint);
3087 const double radiusPix = std::hypot((pt.
x() - pEdge.
x()), pt.
y() - pEdge.
y());
3088 setCatSearchROI(pt, radiusPix);
3094 qCDebug(KSTARS_FITS) <<
"Unable to process Catalog Object request...";
3097 if (Options::fitsCatalog() == CAT_SIMBAD)
3098 return findSimbadObjectsInImage(searchCenter, radius);
3103bool FITSData::findWCSBounds(
double &minRA,
double &maxRA,
double &minDec,
double &maxDec)
3105 if (m_WCSHandle ==
nullptr)
3107 m_LastError =
i18n(
"No world coordinate systems found.");
3116 auto updateMinMax = [&](
int x,
int y)
3119 double imgcrd[2], phi, pixcrd[2], theta, world[2];
3124 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
3127 minRA = std::min(minRA, world[0]);
3128 maxRA = std::max(maxRA, world[0]);
3129 minDec = std::min(minDec, world[1]);
3130 maxDec = std::max(maxDec, world[1]);
3134 for (
int y = 0; y < height(); y++)
3137 updateMinMax(width() - 1, y);
3140 for (
int x = 1; x < width() - 1; x++)
3143 updateMinMax(x, height() - 1);
3147 SkyPoint NCP(0, 90);
3148 SkyPoint SCP(0, -90);
3149 QPointF imagePoint, pPoint;
3150 if (wcsToPixel(NCP, pPoint, imagePoint))
3152 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
3157 if (wcsToPixel(SCP, pPoint, imagePoint))
3159 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
3168#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3171 if (KStarsData::Instance() ==
nullptr)
3177 KSNumbers * num =
nullptr;
3179 if (getRecordValue(
"DATE-OBS", date))
3182 tsString = tsString.remove(
'\'').trimmed();
3189 num =
new KSNumbers(KStarsDateTime(ts).djd());
3194 num =
new KSNumbers(KStarsData::Instance()->ut().djd());
3199 m_SkyObjects.clear();
3204 int type = oneObject->type();
3205 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
3206 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
3207 type == SkyObject::SATELLITE);
3210 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3212 for (
auto &
object : list)
3214 world[0] =
object->ra0().Degrees();
3215 world[1] =
object->dec0().Degrees();
3217 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
3222 if (x > 0 && y > 0 && x < w && y < h)
3223 m_SkyObjects.append(
new FITSSkyObject(
object, x, y));
3236bool FITSData::findSimbadObjectsInImage(
SkyPoint searchCenter,
double radius)
3238 m_CatObjects.clear();
3239 m_CatUpdateTable =
false;
3242 QUrl simbadURL = QUrl(
"https://simbad.cds.unistra.fr/simbad/sim-coo");
3243 QUrlQuery simbadQuery(simbadURL.
query());
3245 QString coord = QString(
"%1 %2").arg(searchCenter.
ra0().
toHMSString(
true,
true))
3248 QString radiusStr = QString(
"%1").arg(radius * 60.0, 0,
'f', 5);
3249 m_CatObjQuery = QString(
"%1 %2").arg(coord).arg(radiusStr);
3253 QString epoch = QString(
"J%1").
arg(m_DateTime.epoch());
3255 simbadQuery.addQueryItem(
"Coord", coord);
3256 simbadQuery.addQueryItem(
"Radius", radiusStr);
3257 simbadQuery.addQueryItem(
"Radius.unit",
"arcsec");
3258 simbadQuery.addQueryItem(
"CooFrame",
"ICRS");
3259 simbadQuery.addQueryItem(
"CooEpoch", epoch);
3260 simbadQuery.addQueryItem(
"output.format",
"ASCII");
3261 simbadQuery.addQueryItem(
"output.max", QString(
"%1").arg(10000));
3262 simbadQuery.addQueryItem(
"list.otypesel",
"on");
3263 simbadQuery.addQueryItem(
"otypedisp",
"3");
3264 simbadQuery.addQueryItem(
"list.coo1",
"on");
3265 simbadQuery.addQueryItem(
"frame1",
"ICRS");
3266 simbadQuery.addQueryItem(
"epoch1",
"J2000");
3267 simbadQuery.addQueryItem(
"equi1", epoch);
3268 simbadQuery.addQueryItem(
"list.spsel",
"off");
3269 simbadQuery.addQueryItem(
"list.sizesel",
"on");
3270 simbadQuery.addQueryItem(
"list.bibsel",
"off");
3273 qCDebug(KSTARS_FITS) <<
"Simbad query:" << simbadURL;
3275 m_NetworkAccessManager.reset(
new QNetworkAccessManager(
this));
3277 QNetworkReply *response = m_NetworkAccessManager->get(QNetworkRequest(simbadURL));
3280 qCDebug(KSTARS_FITS) <<
"Error:" << response->
errorString() <<
" occured querying SIMBAD with" << simbadURL;
3281 m_CatQueryInProgress =
false;
3282 emit catalogQueryFailed(
i18n(
"Error querying Simbad"));
3286 m_CatQueryInProgress =
true;
3287 m_CatQueryTimer.setSingleShot(
true);
3289 m_CatQueryTimer.start(60 * 1000);
3290 emit loadingCatalogData();
3295#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3296void FITSData::catTimeout()
3298 m_CatQueryInProgress =
false;
3299 QString text =
i18n(
"Simbad query timed out");
3300 qCDebug(KSTARS_FITS) << text;
3301 emit catalogQueryFailed(text);
3381#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3382void FITSData::simbadResponseReady(
QNetworkReply * simbadReply)
3384 m_CatQueryTimer.stop();
3385 m_CatQueryInProgress =
false;
3388 qCDebug(KSTARS_FITS) <<
"Error:" << simbadReply->
errorString() <<
" occured in SIMBAD query reply";
3389 m_CatQueryInProgress =
false;
3390 emit catalogQueryFailed(
i18n(
"Error querying Simbad"));
3394 auto data = simbadReply->
readAll();
3395 QString dataStr = QString(data);
3400 QString
name,
type, coord, sizeStr;
3401 double dist = 0, magnitude = 0;
3403 QStringList replyLines = dataStr.
split(QLatin1Char(
'\n'));
3405 for (
int i = 0; i < replyLines.
size(); i++)
3408 if (i == 0 && replyLines[i].contains(
"No astronomical object found : "))
3415 QStringList replyData = replyLines[i].split(QLatin1Char(
' '));
3416 if (replyData.
size() == 9 && replyData[0] ==
"coord" && replyData[6] ==
"radius:")
3417 replyStr = QString(
"%1 %2 %3").
arg(replyData[1]).
arg(replyData[2]).
arg(replyData[7]);
3418 if (replyStr != m_CatObjQuery)
3420 qCDebug(KSTARS_FITS) <<
"Simbad query:" << m_CatObjQuery <<
"Reply:" << replyStr <<
". Ignoring...";
3426 if (replyLines[i].contains(
"Number of objects : "))
3427 numObjs = QString(replyLines[i].mid(20)).toInt(&ok);
3428 else if (replyLines[i].startsWith(
"Object "))
3432 const int firstBreak = replyLines[i].
indexOf(
"---");
3433 const int secondBreak = replyLines[i].
indexOf(
"---", firstBreak + 3);
3434 name = replyLines[i].
mid(7, firstBreak - 7).trimmed();
3435 type = replyLines[i].
mid(firstBreak + 3, secondBreak - firstBreak - 3).trimmed();
3439 qCDebug(KSTARS_FITS) <<
"Bad Simbad Reply, Ignoring...";
3443 else if (numObjs == 1 && i >= 7)
3445 if (replyLines[i].startsWith(
"Coordinates(ICRS"))
3447 QStringList replyData = replyLines[i].split(QLatin1Char(
' '));
3448 if (replyData.
size() >= 8)
3449 coord = replyData[1] +
" " + replyData[2] +
" " + replyData[3] +
" " +
3450 replyData[5] +
" " + replyData[6] +
" " + replyData[7];
3452 else if (replyLines[i].startsWith(
"Flux V :"))
3454 QStringList replyData = replyLines[i].split(QLatin1Char(
' '));
3455 if (replyData.
size() >= 4)
3456 magnitude = QString(replyData[3]).toDouble(&ok);
3458 else if (replyLines[i].startsWith(
"Angular size:"))
3459 sizeStr = replyLines[i].
mid(13, replyLines[i].indexOf(
"~") - 13).trimmed();
3460 else if (replyLines[i].startsWith(
"=========================================="))
3463 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3467 else if (numObjs > 1 && i >= 9)
3470 QStringList replyData = replyLines[i].split(QLatin1Char(
'|'));
3472 name =
type = coord = sizeStr =
"";
3473 magnitude = dist = 0.0;
3474 for (
int j = 0; j < replyData.
size(); j++)
3477 num = QString(replyData[j]).toInt(&ok);
3479 dist = QString(replyData[j]).toDouble(&ok) / 60.0;
3481 name = (QString(replyData[j]).trimmed());
3483 type = QString(replyData[j]).trimmed();
3485 coord = QString(replyData[j]).
trimmed();
3487 magnitude = QString(replyData[j]).
toDouble(&ok);
3489 sizeStr = QString(replyData[j]).
trimmed();
3495 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3499 if (numObjs != ++dataRow)
3500 qCDebug(KSTARS_FITS) <<
"Simbad Header rows:" << numObjs <<
". Data rows:" << dataRow;
3502 m_CatObjectsSearched =
false;
3506 emit loadedCatalogData();
3511bool FITSData::addCatObject(
const int num,
const QString name,
const QString type,
const QString coord,
const double dist,
3512 const double magnitude,
const QString sizeStr)
3517 QStringList split = coord.
split(QLatin1Char(
' '));
3518 if (split.
size() != 6)
3520 qCDebug(KSTARS_FITS) <<
"Coordinates for " <<
name <<
"invalid: " << coord;
3524 QString raStr = QString(
"%1 %2 %3").arg(split[0]).arg(split[1]).arg(split[2]);
3525 QString decStr = QString(
"%1 %2 %3").arg(split[3]).arg(split[4]).arg(split[5]);
3526 ok = r.setFromString(raStr,
false);
3528 ok = d.setFromString(decStr,
true);
3531 qCDebug(KSTARS_FITS) <<
"Coordinates for " <<
name <<
"invalid: " << coord;
3538 QStringList sizeSplit = sizeStr.
split(QLatin1Char(
' '));
3539 if (sizeSplit.
size() >= 1)
3542 size = QString(sizeSplit[0]).toDouble(&ok);
3544 qCDebug(KSTARS_FITS) <<
"Angular size for " <<
name <<
"invalid: " << sizeStr;
3546 CatObject catObject;
3547 catObject.num = num;
3548 catObject.catType = CAT_SIMBAD;
3549 catObject.name =
name;
3550 catObject.typeCode =
type;
3551 catObject.typeLabel = getCatObjectLabel(type);
3554 catObject.dist = dist;
3555 catObject.magnitude = magnitude;
3556 catObject.size = size;
3557 catObject.highlight =
false;
3558 catObject.show = getCatObjectFilter(type);
3560 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3563 world[0] = r.Degrees();
3564 world[1] = d.Degrees();
3566 int status = wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]);
3569 qCDebug(KSTARS_FITS) << QString(
"wcss2p error processing object %1 %2: %3.").arg(name).arg(status).arg(wcs_errmsg[status]);
3574 double x = pixcrd[0];
3575 double y = pixcrd[1];
3576 if (x > 0 && y > 0 && x < width() && y < height())
3580 m_CatObjects.append(catObject);
3585void FITSData::setCatObjectsFilters(
const QStringList filters)
3587 m_CatObjectsFilters = filters;
3590#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3591void FITSData::setCatSearchROI(
const QPoint searchCenter,
const int radius)
3593 m_CatROIPt = searchCenter;
3594 m_CatROIRadius = radius;
3595 Q_UNUSED(searchCenter);
3600#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3606 for (
int i = 0; i < MAX_CAT_OBJ_TYPES; i++)
3608 if (code == catObjTypes[i].code)
3610 label = catObjTypes[i].label;
3613 else if (code == catObjTypes[i].candidateCode)
3615 label = catObjTypes[i].label +
"_Cand";
3624#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3625bool FITSData::getCatObjectFilter(
const QString type)
const
3627 if (m_CatObjectsFilters.isEmpty() || m_CatObjectsFilters.contains(type))
3633void FITSData::filterCatObjects()
3635#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3636 bool changed =
false;
3637 for (
int i = 0; i < m_CatObjects.size(); i++)
3639 bool showState = getCatObjectFilter(m_CatObjects[i].typeCode);
3640 if (m_CatObjects[i].show != showState)
3643 m_CatObjects[i].show = showState;
3651bool FITSData::highlightCatObject(
const int hilite,
const int lolite)
3653#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3654 if (hilite < 0 || hilite >= m_CatObjects.size())
3658 if (lolite >= 0 && lolite < m_CatObjects.size())
3659 m_CatObjects[lolite].highlight =
false;
3663 for (
int i = 0; i < m_CatObjects.size(); i++)
3665 if (m_CatObjects[i].highlight)
3667 m_CatObjects[i].highlight =
false;
3674 m_CatObjects[hilite].highlight =
true;
3679int FITSData::getFlipVCounter()
const
3681 return flipVCounter;
3684void FITSData::setFlipVCounter(
int value)
3686 flipVCounter = value;
3689int FITSData::getFlipHCounter()
const
3691 return flipHCounter;
3694void FITSData::setFlipHCounter(
int value)
3696 flipHCounter = value;
3699int FITSData::getRotCounter()
const
3704void FITSData::setRotCounter(
int value)
3714template <
typename T>
3715bool FITSData::rotFITS(
int rotate,
int mirror)
3719 uint8_t * rotimage =
nullptr;
3724 else if (rotate == 2)
3726 else if (rotate == 3)
3728 else if (rotate < 0)
3729 rotate = rotate + 360;
3731 nx = m_Statistics.width;
3732 ny = m_Statistics.height;
3734 int BBP = m_Statistics.bytesPerPixel;
3737 rotimage =
new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
3739 if (rotimage ==
nullptr)
3741 qWarning() <<
"Unable to allocate memory for rotated image buffer!";
3745 auto * rotBuffer =
reinterpret_cast<T *
>(rotimage);
3746 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
3749 if (rotate < 45 && rotate > -45)
3753 for (
int i = 0; i < m_Statistics.channels; i++)
3755 offset = m_Statistics.samples_per_channel * i;
3756 for (x1 = 0; x1 < nx; x1++)
3759 for (y1 = 0; y1 < ny; y1++)
3760 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3764 else if (mirror == 2)
3766 for (
int i = 0; i < m_Statistics.channels; i++)
3768 offset = m_Statistics.samples_per_channel * i;
3769 for (y1 = 0; y1 < ny; y1++)
3772 for (x1 = 0; x1 < nx; x1++)
3773 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3779 for (
int i = 0; i < m_Statistics.channels; i++)
3781 offset = m_Statistics.samples_per_channel * i;
3782 for (y1 = 0; y1 < ny; y1++)
3784 for (x1 = 0; x1 < nx; x1++)
3785 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3792 else if (rotate >= 45 && rotate < 135)
3796 for (
int i = 0; i < m_Statistics.channels; i++)
3798 offset = m_Statistics.samples_per_channel * i;
3799 for (y1 = 0; y1 < ny; y1++)
3802 for (x1 = 0; x1 < nx; x1++)
3805 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3810 else if (mirror == 2)
3812 for (
int i = 0; i < m_Statistics.channels; i++)
3814 offset = m_Statistics.samples_per_channel * i;
3815 for (y1 = 0; y1 < ny; y1++)
3817 for (x1 = 0; x1 < nx; x1++)
3818 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3824 for (
int i = 0; i < m_Statistics.channels; i++)
3826 offset = m_Statistics.samples_per_channel * i;
3827 for (y1 = 0; y1 < ny; y1++)
3830 for (x1 = 0; x1 < nx; x1++)
3833 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3839 m_Statistics.width = ny;
3840 m_Statistics.height = nx;
3844 else if (rotate >= 135 && rotate < 225)
3848 for (
int i = 0; i < m_Statistics.channels; i++)
3850 offset = m_Statistics.samples_per_channel * i;
3851 for (y1 = 0; y1 < ny; y1++)
3854 for (x1 = 0; x1 < nx; x1++)
3855 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3859 else if (mirror == 2)
3861 for (
int i = 0; i < m_Statistics.channels; i++)
3863 offset = m_Statistics.samples_per_channel * i;
3864 for (x1 = 0; x1 < nx; x1++)
3867 for (y1 = 0; y1 < ny; y1++)
3868 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3874 for (
int i = 0; i < m_Statistics.channels; i++)
3876 offset = m_Statistics.samples_per_channel * i;
3877 for (y1 = 0; y1 < ny; y1++)
3880 for (x1 = 0; x1 < nx; x1++)
3883 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3891 else if (rotate >= 225 && rotate < 315)
3895 for (
int i = 0; i < m_Statistics.channels; i++)
3897 offset = m_Statistics.samples_per_channel * i;
3898 for (y1 = 0; y1 < ny; y1++)
3900 for (x1 = 0; x1 < nx; x1++)
3901 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3905 else if (mirror == 2)
3907 for (
int i = 0; i < m_Statistics.channels; i++)
3909 offset = m_Statistics.samples_per_channel * i;
3910 for (y1 = 0; y1 < ny; y1++)
3913 for (x1 = 0; x1 < nx; x1++)
3916 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3923 for (
int i = 0; i < m_Statistics.channels; i++)
3925 offset = m_Statistics.samples_per_channel * i;
3926 for (y1 = 0; y1 < ny; y1++)
3929 for (x1 = 0; x1 < nx; x1++)
3932 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3938 m_Statistics.width = ny;
3939 m_Statistics.height = nx;
3943 else if (rotate >= 315 && mirror)
3945 for (
int i = 0; i < m_Statistics.channels; i++)
3947 offset = m_Statistics.samples_per_channel * i;
3948 for (y1 = 0; y1 < ny; y1++)
3950 for (x1 = 0; x1 < nx; x1++)
3954 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3960 delete[] m_ImageBuffer;
3961 m_ImageBuffer = rotimage;
3966void FITSData::rotWCSFITS(
int angle,
int mirror)
3970 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
3971 int WCS_DECIMALS = 6;
3973 naxis1 = m_Statistics.width;
3974 naxis2 = m_Statistics.height;
3976 if (fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
3985 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3986 fits_update_key_dbl(fptr,
"CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3988 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
3989 fits_update_key_dbl(fptr,
"CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3997 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3998 fits_update_key_dbl(fptr,
"CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
4000 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
4001 fits_update_key_dbl(fptr,
"CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
4005 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
4006 fits_update_key_dbl(fptr,
"LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4010 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
4011 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4013 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp1, comment, &status))
4014 fits_update_key_dbl(fptr,
"CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
4016 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp1, comment, &status))
4017 fits_update_key_dbl(fptr,
"CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
4023 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
4027 if (!fits_read_key_dbl(fptr,
"LTM2_2", &ctemp2, comment, &status))
4028 if (ctemp1 == ctemp2)
4033 if (!fits_read_key_dbl(fptr,
"LTV1", <v1, comment, &status))
4034 fits_delete_key(fptr,
"LTV1", &status);
4035 if (!fits_read_key_dbl(fptr,
"LTV2", <v2, comment, &status))
4036 fits_delete_key(fptr,
"LTV2", &status);
4040 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp3, comment, &status))
4041 fits_update_key_dbl(fptr,
"CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
4043 if (!fits_read_key_dbl(fptr,
"CRPIX2", &ctemp3, comment, &status))
4044 fits_update_key_dbl(fptr,
"CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
4048 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp3, comment, &status))
4049 fits_update_key_dbl(fptr,
"CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4051 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp3, comment, &status))
4052 fits_update_key_dbl(fptr,
"CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4054 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status))
4055 fits_update_key_dbl(fptr,
"CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4057 if (!fits_read_key_dbl(fptr,
"CD2_2", &ctemp3, comment, &status))
4058 fits_update_key_dbl(fptr,
"CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4062 fits_delete_key(fptr,
"LTM1_1", &status);
4063 fits_delete_key(fptr,
"LTM1_2", &status);
4071 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp1, comment, &status) &&
4072 !fits_read_key_dbl(fptr,
"CRPIX2", &ctemp2, comment, &status))
4077 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4078 else if (angle == 90)
4080 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4081 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4083 else if (angle == 180)
4085 fits_update_key_dbl(fptr,
"CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
4086 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4088 else if (angle == 270)
4090 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4091 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4098 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4099 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4101 else if (angle == 180)
4103 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4104 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4106 else if (angle == 270)
4108 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4109 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4117 if (!fits_read_key_dbl(fptr,
"CDELT1", &ctemp1, comment, &status) &&
4118 !fits_read_key_dbl(fptr,
"CDELT2", &ctemp2, comment, &status))
4123 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4124 else if (angle == 90)
4126 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4127 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4129 else if (angle == 180)
4131 fits_update_key_dbl(fptr,
"CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
4132 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4134 else if (angle == 270)
4136 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4137 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4144 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4145 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4147 else if (angle == 180)
4149 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4150 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4152 else if (angle == 270)
4154 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4155 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4166 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
4168 fits_read_key_dbl(fptr,
"CD1_2", &ctemp2, comment, &status);
4169 fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status);
4170 fits_read_key_dbl(fptr,
"CD2_2", &ctemp4, comment, &status);
4176 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4177 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4179 else if (angle == 90)
4181 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4182 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4183 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4184 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4186 else if (angle == 180)
4188 fits_update_key_dbl(fptr,
"CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
4189 fits_update_key_dbl(fptr,
"CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
4190 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4191 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4193 else if (angle == 270)
4195 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4196 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4197 fits_update_key_dbl(fptr,
"CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
4198 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4205 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4206 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4207 fits_update_key_dbl(fptr,
"CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
4208 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4210 else if (angle == 180)
4212 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4213 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4214 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4215 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4217 else if (angle == 270)
4219 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4220 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4221 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4222 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4230 if (!fits_read_key_dbl(fptr,
"CO1_1", &ctemp1, comment, &status))
4235 for (i = 1; i < 13; i++)
4237 sprintf(keyword,
"CO1_%d", i);
4238 fits_delete_key(fptr, keyword, &status);
4240 for (i = 1; i < 13; i++)
4242 sprintf(keyword,
"CO2_%d", i);
4243 fits_delete_key(fptr, keyword, &status);
4249uint8_t * FITSData::getWritableImageBuffer()
4251 return m_ImageBuffer;
4254uint8_t
const * FITSData::getImageBuffer()
const
4256 return m_ImageBuffer;
4259void FITSData::setImageBuffer(uint8_t * buffer)
4261 delete[] m_ImageBuffer;
4262 m_ImageBuffer = buffer;
4265bool FITSData::checkDebayer()
4268 char bayerPattern[64], roworder[64];
4271 if (fits_read_keyword(fptr,
"BAYERPAT", bayerPattern,
nullptr, &status))
4274 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
4276 m_LastError =
i18n(
"Only 8 and 16 bits bayered images supported.");
4279 QString pattern(bayerPattern);
4280 pattern = pattern.remove(
'\'').trimmed();
4282 QString order(roworder);
4283 order = order.remove(
'\'').trimmed();
4285 if (order ==
"BOTTOM-UP" && !(m_Statistics.height % 2))
4287 if (pattern ==
"RGGB")
4289 else if (pattern ==
"GBRG")
4291 else if (pattern ==
"GRBG")
4293 else if (pattern ==
"BGGR")
4298 if (pattern ==
"RGGB")
4299 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4300 else if (pattern ==
"GBRG")
4301 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4302 else if (pattern ==
"GRBG")
4303 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4304 else if (pattern ==
"BGGR")
4305 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4309 m_LastError =
i18n(
"Unsupported bayer pattern %1.", pattern);
4313 fits_read_key(fptr, TINT,
"XBAYROFF", &debayerParams.offsetX,
nullptr, &status);
4314 fits_read_key(fptr, TINT,
"YBAYROFF", &debayerParams.offsetY,
nullptr, &status);
4316 if (debayerParams.offsetX == 1)
4321 switch (debayerParams.filter)
4323 case DC1394_COLOR_FILTER_RGGB:
4324 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4326 case DC1394_COLOR_FILTER_GBRG:
4327 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4329 case DC1394_COLOR_FILTER_GRBG:
4330 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4332 case DC1394_COLOR_FILTER_BGGR:
4333 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4336 debayerParams.offsetX = 0;
4338 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
4340 m_LastError =
i18n(
"Unsupported bayer offsets %1 %2.", debayerParams.offsetX, debayerParams.offsetY);
4349void FITSData::getBayerParams(BayerParams * param)
4351 param->method = debayerParams.method;
4352 param->filter = debayerParams.filter;
4353 param->offsetX = debayerParams.offsetX;
4354 param->offsetY = debayerParams.offsetY;
4357void FITSData::setBayerParams(BayerParams * param)
4359 debayerParams.method = param->method;
4360 debayerParams.filter = param->filter;
4361 debayerParams.offsetX = param->offsetX;
4362 debayerParams.offsetY = param->offsetY;
4365bool FITSData::debayer(
bool reload)
4369 int anynull = 0,
status = 0;
4371 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel,
nullptr, m_ImageBuffer,
4381 switch (m_Statistics.dataType)
4384 return debayer_8bit();
4387 return debayer_16bit();
4394bool FITSData::debayer_8bit()
4396 dc1394error_t error_code;
4398 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4399 uint8_t * destinationBuffer =
nullptr;
4403 destinationBuffer =
new uint8_t[rgb_size];
4405 catch (
const std::bad_alloc &e)
4407 logOOMError(rgb_size);
4408 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4412 auto * bayer_source_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
4413 auto * bayer_destination_buffer =
reinterpret_cast<uint8_t *
>(destinationBuffer);
4415 if (bayer_destination_buffer ==
nullptr)
4417 logOOMError(rgb_size);
4418 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
4422 int ds1394_height = m_Statistics.height;
4423 auto dc1394_source = bayer_source_buffer;
4425 if (debayerParams.offsetY == 1)
4427 dc1394_source += m_Statistics.width;
4432 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4433 debayerParams.filter,
4434 debayerParams.method);
4436 if (error_code != DC1394_SUCCESS)
4438 m_LastError =
i18n(
"Debayer failed (%1)", error_code);
4439 m_Statistics.channels = 1;
4440 delete[] destinationBuffer;
4444 if (m_ImageBufferSize != rgb_size)
4446 delete[] m_ImageBuffer;
4449 m_ImageBuffer =
new uint8_t[rgb_size];
4451 catch (
const std::bad_alloc &e)
4453 delete[] destinationBuffer;
4454 logOOMError(rgb_size);
4455 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4459 m_ImageBufferSize = rgb_size;
4462 auto bayered_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
4466 uint8_t * rBuff = bayered_buffer;
4467 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4468 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4470 int imax = m_Statistics.samples_per_channel * 3 - 3;
4471 for (
int i = 0; i <= imax; i += 3)
4473 *rBuff++ = bayer_destination_buffer[i];
4474 *gBuff++ = bayer_destination_buffer[i + 1];
4475 *bBuff++ = bayer_destination_buffer[i + 2];
4481 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4482 m_Statistics.dataType = TBYTE;
4483 delete[] destinationBuffer;
4487bool FITSData::debayer_16bit()
4489 dc1394error_t error_code;
4491 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4492 uint8_t *destinationBuffer =
nullptr;
4495 destinationBuffer =
new uint8_t[rgb_size];
4497 catch (
const std::bad_alloc &e)
4499 logOOMError(rgb_size);
4500 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4504 auto * bayer_source_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
4505 auto * bayer_destination_buffer =
reinterpret_cast<uint16_t *
>(destinationBuffer);
4507 if (bayer_destination_buffer ==
nullptr)
4509 logOOMError(rgb_size);
4510 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
4514 int ds1394_height = m_Statistics.height;
4515 auto dc1394_source = bayer_source_buffer;
4517 if (debayerParams.offsetY == 1)
4519 dc1394_source += m_Statistics.width;
4524 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4525 debayerParams.filter,
4526 debayerParams.method, 16);
4528 if (error_code != DC1394_SUCCESS)
4530 m_LastError =
i18n(
"Debayer failed (%1)");
4531 m_Statistics.channels = 1;
4532 delete[] destinationBuffer;
4536 if (m_ImageBufferSize != rgb_size)
4538 delete[] m_ImageBuffer;
4541 m_ImageBuffer =
new uint8_t[rgb_size];
4543 catch (
const std::bad_alloc &e)
4545 logOOMError(rgb_size);
4546 delete[] destinationBuffer;
4547 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4551 m_ImageBufferSize = rgb_size;
4554 auto bayered_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
4558 uint16_t * rBuff = bayered_buffer;
4559 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4560 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4562 int imax = m_Statistics.samples_per_channel * 3 - 3;
4563 for (
int i = 0; i <= imax; i += 3)
4565 *rBuff++ = bayer_destination_buffer[i];
4566 *gBuff++ = bayer_destination_buffer[i + 1];
4567 *bBuff++ = bayer_destination_buffer[i + 2];
4570 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4571 m_Statistics.dataType = TUSHORT;
4572 delete[] destinationBuffer;
4576void FITSData::logOOMError(uint32_t requiredMemory)
4578 qCCritical(KSTARS_FITS) <<
"Debayed memory allocation failure. Required Memory:" << KFormat().formatByteSize(requiredMemory)
4579 <<
"Available system memory:" << KSUtils::getAvailableRAM();
4582double FITSData::getADU()
const
4585 for (
int i = 0; i < m_Statistics.channels; i++)
4586 adu += m_Statistics.mean[i];
4588 return (adu /
static_cast<double>(m_Statistics.channels));
4591QString FITSData::getLastError()
const
4596template <
typename T>
4597void FITSData::convertToQImage(
double dataMin,
double dataMax,
double scale,
double zero,
QImage &image)
4599 const auto * buffer =
reinterpret_cast<const T *
>(getImageBuffer());
4600 const T limit = std::numeric_limits<T>::max();
4601 T bMin = dataMin < 0 ? 0 : dataMin;
4602 T bMax = dataMax > limit ? limit : dataMax;
4603 uint16_t w = width();
4604 uint16_t h = height();
4605 uint32_t size = w * h;
4608 if (channels() == 1)
4611 for (
int j = 0; j < h; j++)
4613 unsigned char * scanLine = image.
scanLine(j);
4615 for (
int i = 0; i < w; i++)
4617 val = qBound(bMin, buffer[j * w + i], bMax);
4618 val = val * scale + zero;
4619 scanLine[i] = qBound<unsigned char>(0,
static_cast<uint8_t
>(val), 255);
4625 double rval = 0, gval = 0, bval = 0;
4628 for (
int j = 0; j < h; j++)
4630 auto * scanLine =
reinterpret_cast<QRgb *
>((image.
scanLine(j)));
4632 for (
int i = 0; i < w; i++)
4634 rval = qBound(bMin, buffer[j * w + i], bMax);
4635 gval = qBound(bMin, buffer[j * w + i + size], bMax);
4636 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
4638 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
4640 scanLine[i] = value;
4653 QFuture<bool> future = data.loadFromFile(filename);
4657 if (future.
result() ==
false)
4660 data.getMinMax(&min, &max);
4668 if (data.channels() == 1)
4673 for (
int i = 0; i < 256; i++)
4674 fitsImage.
setColor(i, qRgb(i, i, i));
4681 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
4682 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
4684 double bscale = 255. / (dataMax - dataMin);
4685 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
4688 switch (data.m_Statistics.dataType)
4691 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4695 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4699 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4703 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4707 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4711 data.convertToQImage<
float>(dataMin, dataMax, bscale, bzero, fitsImage);
4715 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4719 data.convertToQImage<
double>(dataMin, dataMax, bscale, bzero, fitsImage);
4733 qCCritical(KSTARS_FITS) <<
"Failed to convert" << filename <<
"to FITS since format" << format <<
"is not supported in Qt";
4741 qCCritical(KSTARS_FITS) <<
"Failed to open image" << filename;
4745 output = QDir(getTemporaryPath()).filePath(filename +
".fits");
4748 fitsfile *fptr =
nullptr;
4750 long fpixel = 1, naxis = input.
allGray() ? 2 : 3, nelements, exposure;
4751 long naxes[3] = { input.
width(), input.
height(), naxis == 3 ? 3 : 1 };
4752 char error_status[512] = {0};
4754 if (fits_create_file(&fptr, QString(
'!' + output).toLocal8Bit().data(), &status))
4756 fits_get_errstatus(status, error_status);
4757 qCCritical(KSTARS_FITS) <<
"Failed to create FITS file. Error:" << error_status;
4761 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
4763 qCWarning(KSTARS_FITS) <<
"fits_create_img failed:" << error_status;
4765 fits_flush_file(fptr, &status);
4766 fits_close_file(fptr, &status);
4771 fits_update_key(fptr, TLONG,
"EXPOSURE", &exposure,
"Total Exposure Time", &status);
4776 nelements = naxes[0] * naxes[1];
4777 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.
bits(), &status))
4779 fits_get_errstatus(status, error_status);
4780 qCWarning(KSTARS_FITS) <<
"fits_write_img GRAY failed:" << error_status;
4782 fits_flush_file(fptr, &status);
4783 fits_close_file(fptr, &status);
4790 nelements = naxes[0] * naxes[1] * 3;
4792 uint8_t *srcBuffer = input.
bits();
4794 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
4796 uint8_t *rgbBuffer =
new uint8_t[nelements];
4797 if (rgbBuffer ==
nullptr)
4799 qCWarning(KSTARS_FITS) <<
"Not enough memory for RGB buffer";
4800 fits_flush_file(fptr, &status);
4801 fits_close_file(fptr, &status);
4805 uint8_t *subR = rgbBuffer;
4806 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
4807 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
4808 for (uint32_t i = 0; i < srcBytes; i += 4)
4810 *subB++ = srcBuffer[i];
4811 *subG++ = srcBuffer[i + 1];
4812 *subR++ = srcBuffer[i + 2];
4815 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
4817 fits_get_errstatus(status, error_status);
4818 qCWarning(KSTARS_FITS) <<
"fits_write_img RGB failed:" << error_status;
4820 fits_flush_file(fptr, &status);
4821 fits_close_file(fptr, &status);
4822 delete [] rgbBuffer;
4826 delete [] rgbBuffer;
4829 if (fits_flush_file(fptr, &status))
4831 fits_get_errstatus(status, error_status);
4832 qCWarning(KSTARS_FITS) <<
"fits_flush_file failed:" << error_status;
4834 fits_close_file(fptr, &status);
4838 if (fits_close_file(fptr, &status))
4840 fits_get_errstatus(status, error_status);
4841 qCWarning(KSTARS_FITS) <<
"fits_close_file failed:" << error_status;
4848void FITSData::injectWCS(
double orientation,
double ra,
double dec,
double pixscale,
bool eastToTheRight)
4852 updateWCSHeaderData(orientation, ra, dec, pixscale, eastToTheRight);
4856 fits_update_key(fptr, TDOUBLE,
"OBJCTRA", &ra,
"Object RA", &status);
4857 fits_update_key(fptr, TDOUBLE,
"OBJCTDEC", &dec,
"Object DEC", &status);
4861 fits_update_key(fptr, TINT,
"EQUINOX", &epoch,
"Equinox", &status);
4863 fits_update_key(fptr, TDOUBLE,
"CRVAL1", &ra,
"CRVAL1", &status);
4864 fits_update_key(fptr, TDOUBLE,
"CRVAL2", &dec,
"CRVAL1", &status);
4866 char radecsys[8] =
"FK5";
4867 char ctype1[16] =
"RA---TAN";
4868 char ctype2[16] =
"DEC--TAN";
4870 fits_update_key(fptr, TSTRING,
"RADECSYS", radecsys,
"RADECSYS", &status);
4871 fits_update_key(fptr, TSTRING,
"CTYPE1", ctype1,
"CTYPE1", &status);
4872 fits_update_key(fptr, TSTRING,
"CTYPE2", ctype2,
"CTYPE2", &status);
4874 double crpix1 = width() / 2.0;
4875 double crpix2 = height() / 2.0;
4877 fits_update_key(fptr, TDOUBLE,
"CRPIX1", &crpix1,
"CRPIX1", &status);
4878 fits_update_key(fptr, TDOUBLE,
"CRPIX2", &crpix2,
"CRPIX2", &status);
4881 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4882 double secpix2 = pixscale;
4884 fits_update_key(fptr, TDOUBLE,
"SECPIX1", &secpix1,
"SECPIX1", &status);
4885 fits_update_key(fptr, TDOUBLE,
"SECPIX2", &secpix2,
"SECPIX2", &status);
4887 double degpix1 = secpix1 / 3600.0;
4888 double degpix2 = secpix2 / 3600.0;
4890 fits_update_key(fptr, TDOUBLE,
"CDELT1", °pix1,
"CDELT1", &status);
4891 fits_update_key(fptr, TDOUBLE,
"CDELT2", °pix2,
"CDELT2", &status);
4894 double rotation = 360 - orientation;
4898 fits_update_key(fptr, TDOUBLE,
"CROTA1", &rotation,
"CROTA1", &status);
4899 fits_update_key(fptr, TDOUBLE,
"CROTA2", &rotation,
"CROTA2", &status);
4902 emit headerChanged();
4907void FITSData::updateWCSHeaderData(
const double orientation,
const double ra,
const double dec,
const double pixscale,
4908 const bool eastToTheRight)
4910 updateRecordValue(
"OBJCTRA", ra,
"Object RA");
4911 updateRecordValue(
"OBJCTDEC", dec,
"Object DEC");
4912 updateRecordValue(
"EQUINOX", 2000,
"Equinox");
4913 updateRecordValue(
"CRVAL1", ra,
"CRVAL1");
4914 updateRecordValue(
"CRVAL2", dec,
"CRVAL2");
4916 updateRecordValue(
"RADECSYS",
"'FK5'",
"RADECSYS");
4917 updateRecordValue(
"CTYPE1",
"'RA---TAN'",
"CTYPE1");
4918 updateRecordValue(
"CTYPE2",
"'DEC--TAN'",
"CTYPE2");
4920 updateRecordValue(
"CRPIX1", m_Statistics.width / 2.0,
"CRPIX1");
4921 updateRecordValue(
"CRPIX2", m_Statistics.height / 2.0,
"CRPIX2");
4923 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4924 double secpix2 = pixscale;
4926 double degpix1 = secpix1 / 3600.0;
4927 double degpix2 = secpix2 / 3600.0;
4928 updateRecordValue(
"CDELT1", degpix1,
"CDELT1");
4929 updateRecordValue(
"CDELT2", degpix2,
"CDELT2");
4932 double rotation = 360 - orientation;
4935 updateRecordValue(
"CROTA1", rotation,
"CROTA1");
4936 updateRecordValue(
"CROTA2", rotation,
"CROTA2");
4939bool FITSData::contains(
const QPointF &point)
const
4941 return (point.
x() >= 0 && point.
y() >= 0 && point.
x() <= m_Statistics.width && point.
y() <= m_Statistics.height);
4944void FITSData::saveStatistics(FITSImage::Statistic &other)
4946 other = m_Statistics;
4949void FITSData::restoreStatistics(FITSImage::Statistic &other)
4951 m_Statistics = other;
4956void FITSData::constructHistogram()
4958 switch (m_Statistics.dataType)
4961 constructHistogramInternal<uint8_t>();
4965 constructHistogramInternal<int16_t>();
4969 constructHistogramInternal<uint16_t>();
4973 constructHistogramInternal<int32_t>();
4977 constructHistogramInternal<uint32_t>();
4981 constructHistogramInternal<float>();
4985 constructHistogramInternal<int64_t>();
4989 constructHistogramInternal<double>();
4997template <
typename T> int32_t FITSData::histogramBinInternal(T value,
int channel)
const
4999 return qMax(
static_cast<T
>(0), qMin(
static_cast<T
>(m_HistogramBinCount),
5000 static_cast<T
>(rint((value - m_Statistics.min[channel]) / m_HistogramBinWidth[channel]))));
5003template <
typename T> int32_t FITSData::histogramBinInternal(
int x,
int y,
int channel)
const
5005 if (!m_ImageBuffer || !isHistogramConstructed())
5007 uint32_t samples = m_Statistics.width * m_Statistics.height;
5008 uint32_t offset = channel * samples;
5009 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
5010 int index = y * m_Statistics.width + x;
5011 const T &sample = buffer[index + offset];
5012 return histogramBinInternal(sample, channel);
5015int32_t FITSData::histogramBin(
int x,
int y,
int channel)
const
5017 switch (m_Statistics.dataType)
5020 return histogramBinInternal<uint8_t>(x, y, channel);
5024 return histogramBinInternal<int16_t>(x, y, channel);
5028 return histogramBinInternal<uint16_t>(x, y, channel);
5032 return histogramBinInternal<int32_t>(x, y, channel);
5036 return histogramBinInternal<uint32_t>(x, y, channel);
5040 return histogramBinInternal<float>(x, y, channel);
5044 return histogramBinInternal<int64_t>(x, y, channel);
5048 return histogramBinInternal<double>(x, y, channel);
5057template <
typename T>
void FITSData::constructHistogramInternal()
5059 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
5060 uint32_t samples = m_Statistics.width * m_Statistics.height;
5061 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
5063 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
5064 if (m_HistogramBinCount <= 0)
5065 m_HistogramBinCount = 256;
5067 for (
int n = 0; n < m_Statistics.channels; n++)
5069 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
5070 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
5071 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
5073 const double minBinSize = (m_Statistics.max[n] > 1.1) ? 1.0 : .0001;
5074 m_HistogramBinWidth[n] = qMax(minBinSize, (m_Statistics.max[n] - m_Statistics.min[n]) / m_HistogramBinCount);
5077 QVector<QFuture<void>> futures;
5079 for (
int n = 0; n < m_Statistics.channels; n++)
5083 for (
int i = 0; i < m_HistogramBinCount; i++)
5084 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
5088 for (
int n = 0; n < m_Statistics.channels; n++)
5092 uint32_t offset = n * samples;
5094 for (uint32_t i = 0; i < samples; i += sampleBy)
5096 int32_t
id = histogramBinInternal<T>(buffer[i + offset], n);
5097 m_HistogramFrequency[n][id] += sampleBy;
5102 for (QFuture<void> future : futures)
5107 for (
int n = 0; n < m_Statistics.channels; n++)
5111 uint32_t accumulator = 0;
5112 for (
int i = 0; i < m_HistogramBinCount; i++)
5114 accumulator += m_HistogramFrequency[n][i];
5115 m_CumulativeFrequency[n].replace(i, accumulator);
5120 for (QFuture<void> future : futures)
5126 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
5127 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] /
static_cast<double>
5128 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
5133 qCDebug(KSTARS_FITS) <<
"FITHistogram: JMIndex " << m_JMIndex;
5135 m_HistogramConstructed =
true;
5136 emit histogramReady();
5139void FITSData::recordLastError(
int errorCode)
5141 char fitsErrorMessage[512] = {0};
5142 fits_get_errstatus(errorCode, fitsErrorMessage);
5143 m_LastError = fitsErrorMessage;
5146double FITSData::getAverageMean(
bool roi)
const
5148 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5149 if (ptr->channels == 1)
5150 return ptr->mean[0];
5152 return (ptr->mean[0] + ptr->mean[1] + ptr->mean[2]) / 3.0;
5155double FITSData::getAverageMedian(
bool roi)
const
5157 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5158 if (ptr->channels == 1)
5159 return ptr->median[0];
5161 return (ptr->median[0] + ptr->median[1] + ptr->median[2]) / 3.0;
5164double FITSData::getAverageStdDev(
bool roi)
const
5166 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5167 if (ptr->channels == 1)
5168 return ptr->stddev[0];
5170 return (ptr->stddev[0] + ptr->stddev[1] + ptr->stddev[2]) / 3.0;
SkyMapComposite * skyComposite()
static KStarsDateTime currentDateTimeUtc()
QList< SkyObject * > findObjectsInArea(const SkyPoint &p1, const SkyPoint &p2)
The sky coordinates of a point in the sky.
const CachingDms & ra0() const
virtual void updateCoordsNow(const KSNumbers *num)
updateCoordsNow Shortcut for updateCoords( const KSNumbers *num, false, nullptr, nullptr,...
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
const CachingDms & dec0() const
void setDec0(dms d)
Sets Dec0, the catalog Declination.
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
const double & Degrees() const
QString i18n(const char *text, const TYPE &arg...)
Type type(const QSqlDatabase &db)
bool remove(const QString &column, const QVariant &value)
char * toString(const EngineQuery &query)
KIOCORE_EXPORT StatJob * stat(const QUrl &url, JobFlags flags=DefaultFlags)
KIOCORE_EXPORT QString number(KIO::filesize_t size)
KGUIADDONS_EXPORT QWindow * window(QObject *job)
void error(QWidget *parent, const QString &text, const QString &title, const KGuiItem &buttonOk, Options options=Notify)
KIOCORE_EXPORT QStringList list(const QString &fileClass)
QString name(StandardAction id)
QString label(StandardShortcut id)
const QList< QKeySequence > & end()
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
QByteArray & append(QByteArrayView data)
const char * constData() const const
bool isEmpty() const const
QByteArray leftJustified(qsizetype width, char fill, bool truncate) const const
void push_back(QByteArrayView str)
qsizetype size() const const
QDateTime currentDateTime()
QDateTime fromString(QStringView string, QStringView format, QCalendar cal)
bool isValid() const const
QString path() const const
bool allGray() const const
void fill(Qt::GlobalColor color)
bool load(QIODevice *device, const char *format)
bool loadFromData(QByteArrayView data, const char *format)
bool save(QIODevice *device, const char *format, int quality) const const
void setColor(int index, QRgb colorValue)
void setColorCount(int colorCount)
QString errorString() const const
void append(QList< T > &&value)
const_reference at(qsizetype i) const const
bool contains(const AT &value) const const
iterator erase(const_iterator begin, const_iterator end)
QList< T > mid(qsizetype pos, qsizetype length) const const
qsizetype size() const const
void finished(QNetworkReply *reply)
NetworkError error() const const
QMetaObject::Connection connect(const QObject *sender, PointerToMemberFunction signal, Functor functor)
void setObjectName(QAnyStringView name)
bool contains(const QPoint &point, bool proper) const const
bool isNull() const const
bool isValid() const const
QPoint topLeft() const const
bool isNull() const const
QString arg(Args &&... args) const const
bool contains(QChar ch, Qt::CaseSensitivity cs) const const
QString fromStdString(const std::string &str)
QString & insert(qsizetype position, QChar ch)
bool isEmpty() const const
QString mid(qsizetype position, qsizetype n) const const
QString & remove(QChar ch, Qt::CaseSensitivity cs)
QString & replace(QChar before, QChar after, Qt::CaseSensitivity cs)
QStringList split(QChar sep, Qt::SplitBehavior behavior, Qt::CaseSensitivity cs) const const
double toDouble(bool *ok) const const
QByteArray toLatin1() const const
QByteArray toLocal8Bit() const const
QString trimmed() const const
bool contains(QLatin1StringView str, Qt::CaseSensitivity cs) const const
qsizetype indexOf(const QRegularExpression &re, qsizetype from) const const
QTextStream & center(QTextStream &stream)
QFuture< void > map(Iterator begin, Iterator end, MapFunctor &&function)
QFuture< T > run(Function function,...)
QString query(ComponentFormattingOptions options) const const
void setQuery(const QString &query, ParsingMode mode)
QString toString(StringFormat mode) 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
Object to hold FITS Header records.