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)
70 QString extension = info.completeSuffix().toLower();
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);
102 qRegisterMetaType<FITSMode>(
"FITSMode");
104 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
105 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
106 debayerParams.offsetX = debayerParams.offsetY = 0;
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);
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)
183 return privateLoad(buffer);
188 loadCommon(inFilename);
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;
219 return loadFITSImage(buffer);
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;
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())
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;
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)
640 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
644 m_Statistics.size = buffer.
size();
650 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
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.");
770 m_LastError =
i18n(
"Cannot open file %1: %2", m_Filename, libraw_strerror(ret));
771 RawProcessor.recycle();
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)
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;
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))
1252 tsString = tsString.remove(
'\'').trimmed();
1260 validObservationDate =
true;
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;
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;
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);
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();
1957 for (
int i = 0; i < nkeys; i++)
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;
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);
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)
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:
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:
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...";
2845 int nkeyrec = 1, nreject = 0;
2846 for(
auto &fitsKeyword : m_HeaderRecords)
2852 if (fitsKeyword.key.startsWith(
"PC1_") || fitsKeyword.key.startsWith(
"PC2_") ||
2853 fitsKeyword.key.startsWith(
"CD1_") || fitsKeyword.key.startsWith(
"CD2_") ||
2854 fitsKeyword.key.startsWith(
"PV1_") || fitsKeyword.key.startsWith(
"PV2_") ||
2855 fitsKeyword.key.startsWith(
"PS1_") || fitsKeyword.key.startsWith(
"PS2_"))
2858 rec.
append(fitsKeyword.key.leftJustified(8,
' ').toLatin1());
2860 rec.
append(fitsKeyword.value.toByteArray());
2862 rec.
append(fitsKeyword.comment.toLatin1());
2868 if ((status = wcspih(header_str.
data(), nkeyrec, WCSHDR_all, 0, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2870 wcsvfree(&m_nwcs, &m_WCSHandle);
2871 m_WCSHandle =
nullptr;
2873 m_LastError =
QString(
"wcspih ERROR %1: %2.").
arg(status).
arg(wcshdr_errmsg[status]);
2874 m_WCSState = Failure;
2878 if (m_WCSHandle ==
nullptr)
2880 m_LastError =
i18n(
"No world coordinate systems found.");
2881 m_WCSState = Failure;
2886 if (m_WCSHandle->crpix[0] == 0)
2888 wcsvfree(&m_nwcs, &m_WCSHandle);
2889 m_WCSHandle =
nullptr;
2891 m_LastError =
i18n(
"No world coordinate systems found.");
2892 m_WCSState = Failure;
2897 if ((status = wcsset(m_WCSHandle)) != 0)
2899 wcsvfree(&m_nwcs, &m_WCSHandle);
2900 m_WCSHandle =
nullptr;
2902 m_LastError =
QString(
"wcsset error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2903 m_WCSState = Failure;
2907 m_ObjectsSearched =
false;
2908 m_CatObjectsSearched =
false;
2909 m_WCSState = Success;
2912 qCDebug(KSTARS_FITS) <<
"Finished WCS Data processing...";
2922#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2924 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2926 if (m_WCSHandle ==
nullptr)
2928 m_LastError =
i18n(
"No world coordinate systems found.");
2932 worldcrd[0] = wcsCoord.
ra0().Degrees();
2933 worldcrd[1] = wcsCoord.
dec0().Degrees();
2935 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2937 m_LastError =
QString(
"wcss2p error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2941 wcsImagePoint.
setX(imgcrd[0]);
2942 wcsImagePoint.
setY(imgcrd[1]);
2944 wcsPixelPoint.
setX(pixcrd[0]);
2945 wcsPixelPoint.
setY(pixcrd[1]);
2950 Q_UNUSED(wcsPixelPoint);
2951 Q_UNUSED(wcsImagePoint);
2956bool FITSData::pixelToWCS(
const QPointF &wcsPixelPoint,
SkyPoint &wcsCoord)
2958#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2960 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2962 if (m_WCSHandle ==
nullptr)
2964 m_LastError =
i18n(
"No world coordinate systems found.");
2968 pixcrd[0] = wcsPixelPoint.
x();
2969 pixcrd[1] = wcsPixelPoint.
y();
2971 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2973 m_LastError =
QString(
"wcsp2s error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2978 wcsCoord.
setRA0(world[0] / 15.0);
2984 Q_UNUSED(wcsPixelPoint);
2990#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2991bool FITSData::searchObjects()
2993 return (Options::fitsCatalog() == CAT_SKYMAP) ? searchSkyMapObjects() : searchCatObjects();
2996bool FITSData::searchSkyMapObjects()
2998 if (m_ObjectsSearched)
3001 m_ObjectsSearched =
true;
3006 pixelToWCS(
QPointF(0, 0), startPoint);
3007 pixelToWCS(
QPointF(width() - 1, height() - 1), endPoint);
3009 return findObjectsInImage(startPoint, endPoint);
3012bool FITSData::searchCatObjects()
3014 if (m_CatObjectsSearched)
3017 m_CatObjectsSearched =
true;
3023 if (catROIRadius() > 0)
3029 ok = pixelToWCS(pt, searchCenter);
3031 ok = pixelToWCS(edgePt, searchEdge);
3037 if (getRecordValue(
"DATE-OBS", date))
3040 tsString = tsString.remove(
'\'').trimmed();
3052 num =
new KSNumbers(KStarsData::Instance()->ut().djd());
3065 pt =
QPoint((width() / 2) - 1, (height() / 2) - 1);
3066 ok = pixelToWCS(pt, searchCenter);
3070 double raEdge = searchCenter.
ra0().Degrees() + (radius / 60.0);
3071 double decEdge = searchCenter.
dec0().Degrees();
3072 SkyPoint searchEdge(raEdge / 15.0, decEdge);
3074 ok = wcsToPixel(searchEdge, pEdge, edgePoint);
3078 const double radiusPix = std::hypot((pt.
x() - pEdge.
x()), pt.
y() - pEdge.
y());
3079 setCatSearchROI(pt, radiusPix);
3085 qCDebug(KSTARS_FITS) <<
"Unable to process Catalog Object request...";
3088 if (Options::fitsCatalog() == CAT_SIMBAD)
3089 return findSimbadObjectsInImage(searchCenter, radius);
3094bool FITSData::findWCSBounds(
double &minRA,
double &maxRA,
double &minDec,
double &maxDec)
3096 if (m_WCSHandle ==
nullptr)
3098 m_LastError =
i18n(
"No world coordinate systems found.");
3107 auto updateMinMax = [&](
int x,
int y)
3110 double imgcrd[2], phi, pixcrd[2], theta, world[2];
3115 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
3118 minRA = std::min(minRA, world[0]);
3119 maxRA = std::max(maxRA, world[0]);
3120 minDec = std::min(minDec, world[1]);
3121 maxDec = std::max(maxDec, world[1]);
3125 for (
int y = 0; y < height(); y++)
3128 updateMinMax(width() - 1, y);
3131 for (
int x = 1; x < width() - 1; x++)
3134 updateMinMax(x, height() - 1);
3141 if (wcsToPixel(NCP, pPoint, imagePoint))
3143 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
3148 if (wcsToPixel(SCP, pPoint, imagePoint))
3150 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
3159#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3162 if (KStarsData::Instance() ==
nullptr)
3170 if (getRecordValue(
"DATE-OBS", date))
3173 tsString = tsString.remove(
'\'').trimmed();
3185 num =
new KSNumbers(KStarsData::Instance()->ut().djd());
3190 m_SkyObjects.
clear();
3195 int type = oneObject->type();
3196 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
3197 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
3198 type == SkyObject::SATELLITE);
3201 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3203 for (
auto &
object : list)
3205 world[0] =
object->ra0().Degrees();
3206 world[1] =
object->dec0().Degrees();
3208 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
3213 if (x > 0 && y > 0 && x < w && y < h)
3214 m_SkyObjects.
append(
new FITSSkyObject(
object, x, y));
3227bool FITSData::findSimbadObjectsInImage(
SkyPoint searchCenter,
double radius)
3229 m_CatObjects.
clear();
3230 m_CatUpdateTable =
false;
3233 QUrl simbadURL =
QUrl(
"https://simbad.cds.unistra.fr/simbad/sim-coo");
3237 .
arg(searchCenter.
dec0().toDMSString(
true,
true,
true));
3246 simbadQuery.addQueryItem(
"Coord", coord);
3247 simbadQuery.addQueryItem(
"Radius", radiusStr);
3248 simbadQuery.addQueryItem(
"Radius.unit",
"arcsec");
3249 simbadQuery.addQueryItem(
"CooFrame",
"ICRS");
3250 simbadQuery.addQueryItem(
"CooEpoch", epoch);
3251 simbadQuery.addQueryItem(
"output.format",
"ASCII");
3252 simbadQuery.addQueryItem(
"output.max",
QString(
"%1").arg(10000));
3253 simbadQuery.addQueryItem(
"list.otypesel",
"on");
3254 simbadQuery.addQueryItem(
"otypedisp",
"3");
3255 simbadQuery.addQueryItem(
"list.coo1",
"on");
3256 simbadQuery.addQueryItem(
"frame1",
"ICRS");
3257 simbadQuery.addQueryItem(
"epoch1",
"J2000");
3258 simbadQuery.addQueryItem(
"equi1", epoch);
3259 simbadQuery.addQueryItem(
"list.spsel",
"off");
3260 simbadQuery.addQueryItem(
"list.sizesel",
"on");
3261 simbadQuery.addQueryItem(
"list.bibsel",
"off");
3264 qCDebug(KSTARS_FITS) <<
"Simbad query:" << simbadURL;
3271 qCDebug(KSTARS_FITS) <<
"Error:" << response->
errorString() <<
" occured querying SIMBAD with" << simbadURL;
3272 m_CatQueryInProgress =
false;
3273 emit catalogQueryFailed(
i18n(
"Error querying Simbad"));
3277 m_CatQueryInProgress =
true;
3280 m_CatQueryTimer.
start(60 * 1000);
3281 emit loadingCatalogData();
3286#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3287void FITSData::catTimeout()
3289 m_CatQueryInProgress =
false;
3291 qCDebug(KSTARS_FITS) << text;
3292 emit catalogQueryFailed(text);
3372#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3373void FITSData::simbadResponseReady(
QNetworkReply * simbadReply)
3375 m_CatQueryTimer.
stop();
3376 m_CatQueryInProgress =
false;
3379 qCDebug(KSTARS_FITS) <<
"Error:" << simbadReply->
errorString() <<
" occured in SIMBAD query reply";
3380 m_CatQueryInProgress =
false;
3381 emit catalogQueryFailed(
i18n(
"Error querying Simbad"));
3385 auto data = simbadReply->
readAll();
3392 double dist = 0, magnitude = 0;
3396 for (
int i = 0; i < replyLines.
size(); i++)
3399 if (i == 0 && replyLines[i].contains(
"No astronomical object found : "))
3407 if (replyData.
size() == 9 && replyData[0] ==
"coord" && replyData[6] ==
"radius:")
3408 replyStr =
QString(
"%1 %2 %3").
arg(replyData[1]).
arg(replyData[2]).
arg(replyData[7]);
3409 if (replyStr != m_CatObjQuery)
3411 qCDebug(KSTARS_FITS) <<
"Simbad query:" << m_CatObjQuery <<
"Reply:" << replyStr <<
". Ignoring...";
3417 if (replyLines[i].contains(
"Number of objects : "))
3419 else if (replyLines[i].startsWith(
"Object "))
3423 const int firstBreak = replyLines[i].
indexOf(
"---");
3424 const int secondBreak = replyLines[i].
indexOf(
"---", firstBreak + 3);
3425 name = replyLines[i].
mid(7, firstBreak - 7).trimmed();
3426 type = replyLines[i].
mid(firstBreak + 3, secondBreak - firstBreak - 3).trimmed();
3430 qCDebug(KSTARS_FITS) <<
"Bad Simbad Reply, Ignoring...";
3434 else if (numObjs == 1 && i >= 7)
3436 if (replyLines[i].startsWith(
"Coordinates(ICRS"))
3439 if (replyData.
size() >= 8)
3440 coord = replyData[1] +
" " + replyData[2] +
" " + replyData[3] +
" " +
3441 replyData[5] +
" " + replyData[6] +
" " + replyData[7];
3443 else if (replyLines[i].startsWith(
"Flux V :"))
3446 if (replyData.
size() >= 4)
3449 else if (replyLines[i].startsWith(
"Angular size:"))
3450 sizeStr = replyLines[i].
mid(13, replyLines[i].indexOf(
"~") - 13).trimmed();
3451 else if (replyLines[i].startsWith(
"=========================================="))
3454 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3458 else if (numObjs > 1 && i >= 9)
3463 name =
type = coord = sizeStr =
"";
3464 magnitude = dist = 0.0;
3465 for (
int j = 0; j < replyData.
size(); j++)
3486 if (addCatObject(num, name, type, coord, dist, magnitude, sizeStr))
3490 if (numObjs != ++dataRow)
3491 qCDebug(KSTARS_FITS) <<
"Simbad Header rows:" << numObjs <<
". Data rows:" << dataRow;
3493 m_CatObjectsSearched =
false;
3497 emit loadedCatalogData();
3502bool FITSData::addCatObject(
const int num,
const QString name,
const QString type,
const QString coord,
const double dist,
3503 const double magnitude,
const QString sizeStr)
3509 if (split.
size() != 6)
3511 qCDebug(KSTARS_FITS) <<
"Coordinates for " <<
name <<
"invalid: " << coord;
3517 ok = r.setFromString(raStr,
false);
3519 ok = d.setFromString(decStr,
true);
3522 qCDebug(KSTARS_FITS) <<
"Coordinates for " <<
name <<
"invalid: " << coord;
3530 if (sizeSplit.
size() >= 1)
3535 qCDebug(KSTARS_FITS) <<
"Angular size for " <<
name <<
"invalid: " << sizeStr;
3537 CatObject catObject;
3538 catObject.num = num;
3539 catObject.catType = CAT_SIMBAD;
3540 catObject.name =
name;
3541 catObject.typeCode =
type;
3542 catObject.typeLabel = getCatObjectLabel(type);
3545 catObject.dist = dist;
3546 catObject.magnitude = magnitude;
3547 catObject.size = size;
3548 catObject.highlight =
false;
3549 catObject.show = getCatObjectFilter(type);
3551 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3554 world[0] = r.Degrees();
3555 world[1] = d.Degrees();
3557 int status = wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]);
3560 qCDebug(KSTARS_FITS) <<
QString(
"wcss2p error processing object %1 %2: %3.").
arg(name).
arg(status).
arg(wcs_errmsg[status]);
3565 double x = pixcrd[0];
3566 double y = pixcrd[1];
3567 if (x > 0 && y > 0 && x < width() && y < height())
3571 m_CatObjects.
append(catObject);
3576void FITSData::setCatObjectsFilters(
const QStringList filters)
3578 m_CatObjectsFilters = filters;
3581#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3582void FITSData::setCatSearchROI(
const QPoint searchCenter,
const int radius)
3584 m_CatROIPt = searchCenter;
3585 m_CatROIRadius = radius;
3586 Q_UNUSED(searchCenter);
3591#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3597 for (
int i = 0; i < MAX_CAT_OBJ_TYPES; i++)
3599 if (code == catObjTypes[i].code)
3601 label = catObjTypes[i].label;
3604 else if (code == catObjTypes[i].candidateCode)
3606 label = catObjTypes[i].label +
"_Cand";
3615#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3616bool FITSData::getCatObjectFilter(
const QString type)
const
3618 if (m_CatObjectsFilters.
isEmpty() || m_CatObjectsFilters.
contains(type))
3624void FITSData::filterCatObjects()
3626#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3627 bool changed =
false;
3628 for (
int i = 0; i < m_CatObjects.
size(); i++)
3630 bool showState = getCatObjectFilter(m_CatObjects[i].typeCode);
3631 if (m_CatObjects[i].show != showState)
3634 m_CatObjects[i].show = showState;
3642bool FITSData::highlightCatObject(
const int hilite,
const int lolite)
3644#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3645 if (hilite < 0 || hilite >= m_CatObjects.
size())
3649 if (lolite >= 0 && lolite < m_CatObjects.
size())
3650 m_CatObjects[lolite].highlight =
false;
3654 for (
int i = 0; i < m_CatObjects.
size(); i++)
3656 if (m_CatObjects[i].highlight)
3658 m_CatObjects[i].highlight =
false;
3665 m_CatObjects[hilite].highlight =
true;
3670int FITSData::getFlipVCounter()
const
3672 return flipVCounter;
3675void FITSData::setFlipVCounter(
int value)
3677 flipVCounter = value;
3680int FITSData::getFlipHCounter()
const
3682 return flipHCounter;
3685void FITSData::setFlipHCounter(
int value)
3687 flipHCounter = value;
3690int FITSData::getRotCounter()
const
3695void FITSData::setRotCounter(
int value)
3705template <
typename T>
3706bool FITSData::rotFITS(
int rotate,
int mirror)
3710 uint8_t * rotimage =
nullptr;
3715 else if (rotate == 2)
3717 else if (rotate == 3)
3719 else if (rotate < 0)
3720 rotate = rotate + 360;
3722 nx = m_Statistics.width;
3723 ny = m_Statistics.height;
3725 int BBP = m_Statistics.bytesPerPixel;
3728 rotimage =
new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
3730 if (rotimage ==
nullptr)
3732 qWarning() <<
"Unable to allocate memory for rotated image buffer!";
3736 auto * rotBuffer =
reinterpret_cast<T *
>(rotimage);
3737 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
3740 if (rotate < 45 && rotate > -45)
3744 for (
int i = 0; i < m_Statistics.channels; i++)
3746 offset = m_Statistics.samples_per_channel * i;
3747 for (x1 = 0; x1 < nx; x1++)
3750 for (y1 = 0; y1 < ny; y1++)
3751 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3755 else if (mirror == 2)
3757 for (
int i = 0; i < m_Statistics.channels; i++)
3759 offset = m_Statistics.samples_per_channel * i;
3760 for (y1 = 0; y1 < ny; y1++)
3763 for (x1 = 0; x1 < nx; x1++)
3764 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3770 for (
int i = 0; i < m_Statistics.channels; i++)
3772 offset = m_Statistics.samples_per_channel * i;
3773 for (y1 = 0; y1 < ny; y1++)
3775 for (x1 = 0; x1 < nx; x1++)
3776 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3783 else if (rotate >= 45 && rotate < 135)
3787 for (
int i = 0; i < m_Statistics.channels; i++)
3789 offset = m_Statistics.samples_per_channel * i;
3790 for (y1 = 0; y1 < ny; y1++)
3793 for (x1 = 0; x1 < nx; x1++)
3796 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3801 else if (mirror == 2)
3803 for (
int i = 0; i < m_Statistics.channels; i++)
3805 offset = m_Statistics.samples_per_channel * i;
3806 for (y1 = 0; y1 < ny; y1++)
3808 for (x1 = 0; x1 < nx; x1++)
3809 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3815 for (
int i = 0; i < m_Statistics.channels; i++)
3817 offset = m_Statistics.samples_per_channel * i;
3818 for (y1 = 0; y1 < ny; y1++)
3821 for (x1 = 0; x1 < nx; x1++)
3824 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3830 m_Statistics.width = ny;
3831 m_Statistics.height = nx;
3835 else if (rotate >= 135 && rotate < 225)
3839 for (
int i = 0; i < m_Statistics.channels; i++)
3841 offset = m_Statistics.samples_per_channel * i;
3842 for (y1 = 0; y1 < ny; y1++)
3845 for (x1 = 0; x1 < nx; x1++)
3846 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3850 else if (mirror == 2)
3852 for (
int i = 0; i < m_Statistics.channels; i++)
3854 offset = m_Statistics.samples_per_channel * i;
3855 for (x1 = 0; x1 < nx; x1++)
3858 for (y1 = 0; y1 < ny; y1++)
3859 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3865 for (
int i = 0; i < m_Statistics.channels; i++)
3867 offset = m_Statistics.samples_per_channel * i;
3868 for (y1 = 0; y1 < ny; y1++)
3871 for (x1 = 0; x1 < nx; x1++)
3874 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3882 else if (rotate >= 225 && rotate < 315)
3886 for (
int i = 0; i < m_Statistics.channels; i++)
3888 offset = m_Statistics.samples_per_channel * i;
3889 for (y1 = 0; y1 < ny; y1++)
3891 for (x1 = 0; x1 < nx; x1++)
3892 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3896 else if (mirror == 2)
3898 for (
int i = 0; i < m_Statistics.channels; i++)
3900 offset = m_Statistics.samples_per_channel * i;
3901 for (y1 = 0; y1 < ny; y1++)
3904 for (x1 = 0; x1 < nx; x1++)
3907 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3914 for (
int i = 0; i < m_Statistics.channels; i++)
3916 offset = m_Statistics.samples_per_channel * i;
3917 for (y1 = 0; y1 < ny; y1++)
3920 for (x1 = 0; x1 < nx; x1++)
3923 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3929 m_Statistics.width = ny;
3930 m_Statistics.height = nx;
3934 else if (rotate >= 315 && mirror)
3936 for (
int i = 0; i < m_Statistics.channels; i++)
3938 offset = m_Statistics.samples_per_channel * i;
3939 for (y1 = 0; y1 < ny; y1++)
3941 for (x1 = 0; x1 < nx; x1++)
3945 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3951 delete[] m_ImageBuffer;
3952 m_ImageBuffer = rotimage;
3957void FITSData::rotWCSFITS(
int angle,
int mirror)
3961 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
3962 int WCS_DECIMALS = 6;
3964 naxis1 = m_Statistics.width;
3965 naxis2 = m_Statistics.height;
3967 if (fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
3976 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3977 fits_update_key_dbl(fptr,
"CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3979 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
3980 fits_update_key_dbl(fptr,
"CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3988 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3989 fits_update_key_dbl(fptr,
"CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
3991 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
3992 fits_update_key_dbl(fptr,
"CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
3996 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
3997 fits_update_key_dbl(fptr,
"LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4001 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
4002 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4004 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp1, comment, &status))
4005 fits_update_key_dbl(fptr,
"CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
4007 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp1, comment, &status))
4008 fits_update_key_dbl(fptr,
"CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
4014 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
4018 if (!fits_read_key_dbl(fptr,
"LTM2_2", &ctemp2, comment, &status))
4019 if (ctemp1 == ctemp2)
4024 if (!fits_read_key_dbl(fptr,
"LTV1", <v1, comment, &status))
4025 fits_delete_key(fptr,
"LTV1", &status);
4026 if (!fits_read_key_dbl(fptr,
"LTV2", <v2, comment, &status))
4027 fits_delete_key(fptr,
"LTV2", &status);
4031 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp3, comment, &status))
4032 fits_update_key_dbl(fptr,
"CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
4034 if (!fits_read_key_dbl(fptr,
"CRPIX2", &ctemp3, comment, &status))
4035 fits_update_key_dbl(fptr,
"CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
4039 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp3, comment, &status))
4040 fits_update_key_dbl(fptr,
"CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4042 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp3, comment, &status))
4043 fits_update_key_dbl(fptr,
"CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4045 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status))
4046 fits_update_key_dbl(fptr,
"CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4048 if (!fits_read_key_dbl(fptr,
"CD2_2", &ctemp3, comment, &status))
4049 fits_update_key_dbl(fptr,
"CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
4053 fits_delete_key(fptr,
"LTM1_1", &status);
4054 fits_delete_key(fptr,
"LTM1_2", &status);
4062 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp1, comment, &status) &&
4063 !fits_read_key_dbl(fptr,
"CRPIX2", &ctemp2, comment, &status))
4068 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4069 else if (angle == 90)
4071 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4072 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4074 else if (angle == 180)
4076 fits_update_key_dbl(fptr,
"CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
4077 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4079 else if (angle == 270)
4081 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4082 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4089 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4090 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
4092 else if (angle == 180)
4094 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4095 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
4097 else if (angle == 270)
4099 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
4100 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
4108 if (!fits_read_key_dbl(fptr,
"CDELT1", &ctemp1, comment, &status) &&
4109 !fits_read_key_dbl(fptr,
"CDELT2", &ctemp2, comment, &status))
4114 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4115 else if (angle == 90)
4117 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4118 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4120 else if (angle == 180)
4122 fits_update_key_dbl(fptr,
"CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
4123 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4125 else if (angle == 270)
4127 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4128 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4135 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
4136 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
4138 else if (angle == 180)
4140 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
4141 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
4143 else if (angle == 270)
4145 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
4146 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
4157 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
4159 fits_read_key_dbl(fptr,
"CD1_2", &ctemp2, comment, &status);
4160 fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status);
4161 fits_read_key_dbl(fptr,
"CD2_2", &ctemp4, comment, &status);
4167 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4168 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4170 else if (angle == 90)
4172 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4173 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4174 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4175 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4177 else if (angle == 180)
4179 fits_update_key_dbl(fptr,
"CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
4180 fits_update_key_dbl(fptr,
"CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
4181 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4182 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4184 else if (angle == 270)
4186 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4187 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4188 fits_update_key_dbl(fptr,
"CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
4189 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4196 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
4197 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
4198 fits_update_key_dbl(fptr,
"CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
4199 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
4201 else if (angle == 180)
4203 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
4204 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
4205 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
4206 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
4208 else if (angle == 270)
4210 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
4211 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
4212 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
4213 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
4221 if (!fits_read_key_dbl(fptr,
"CO1_1", &ctemp1, comment, &status))
4226 for (i = 1; i < 13; i++)
4228 sprintf(keyword,
"CO1_%d", i);
4229 fits_delete_key(fptr, keyword, &status);
4231 for (i = 1; i < 13; i++)
4233 sprintf(keyword,
"CO2_%d", i);
4234 fits_delete_key(fptr, keyword, &status);
4240uint8_t * FITSData::getWritableImageBuffer()
4242 return m_ImageBuffer;
4245uint8_t
const * FITSData::getImageBuffer()
const
4247 return m_ImageBuffer;
4250void FITSData::setImageBuffer(uint8_t * buffer)
4252 delete[] m_ImageBuffer;
4253 m_ImageBuffer = buffer;
4256bool FITSData::checkDebayer()
4259 char bayerPattern[64], roworder[64];
4262 if (fits_read_keyword(fptr,
"BAYERPAT", bayerPattern,
nullptr, &status))
4265 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
4267 m_LastError =
i18n(
"Only 8 and 16 bits bayered images supported.");
4270 QString pattern(bayerPattern);
4271 pattern = pattern.remove(
'\'').trimmed();
4274 order = order.remove(
'\'').trimmed();
4276 if (order ==
"BOTTOM-UP" && !(m_Statistics.height % 2))
4278 if (pattern ==
"RGGB")
4280 else if (pattern ==
"GBRG")
4282 else if (pattern ==
"GRBG")
4284 else if (pattern ==
"BGGR")
4289 if (pattern ==
"RGGB")
4290 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4291 else if (pattern ==
"GBRG")
4292 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4293 else if (pattern ==
"GRBG")
4294 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4295 else if (pattern ==
"BGGR")
4296 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4300 m_LastError =
i18n(
"Unsupported bayer pattern %1.", pattern);
4304 fits_read_key(fptr, TINT,
"XBAYROFF", &debayerParams.offsetX,
nullptr, &status);
4305 fits_read_key(fptr, TINT,
"YBAYROFF", &debayerParams.offsetY,
nullptr, &status);
4307 if (debayerParams.offsetX == 1)
4312 switch (debayerParams.filter)
4314 case DC1394_COLOR_FILTER_RGGB:
4315 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
4317 case DC1394_COLOR_FILTER_GBRG:
4318 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
4320 case DC1394_COLOR_FILTER_GRBG:
4321 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
4323 case DC1394_COLOR_FILTER_BGGR:
4324 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
4327 debayerParams.offsetX = 0;
4329 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
4331 m_LastError =
i18n(
"Unsupported bayer offsets %1 %2.", debayerParams.offsetX, debayerParams.offsetY);
4340void FITSData::getBayerParams(BayerParams * param)
4342 param->method = debayerParams.method;
4343 param->filter = debayerParams.filter;
4344 param->offsetX = debayerParams.offsetX;
4345 param->offsetY = debayerParams.offsetY;
4348void FITSData::setBayerParams(BayerParams * param)
4350 debayerParams.method = param->method;
4351 debayerParams.filter = param->filter;
4352 debayerParams.offsetX = param->offsetX;
4353 debayerParams.offsetY = param->offsetY;
4356bool FITSData::debayer(
bool reload)
4360 int anynull = 0,
status = 0;
4362 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel,
nullptr, m_ImageBuffer,
4372 switch (m_Statistics.dataType)
4375 return debayer_8bit();
4378 return debayer_16bit();
4385bool FITSData::debayer_8bit()
4387 dc1394error_t error_code;
4389 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4390 uint8_t * destinationBuffer =
nullptr;
4394 destinationBuffer =
new uint8_t[rgb_size];
4396 catch (
const std::bad_alloc &e)
4398 logOOMError(rgb_size);
4399 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4403 auto * bayer_source_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
4404 auto * bayer_destination_buffer =
reinterpret_cast<uint8_t *
>(destinationBuffer);
4406 if (bayer_destination_buffer ==
nullptr)
4408 logOOMError(rgb_size);
4409 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
4413 int ds1394_height = m_Statistics.height;
4414 auto dc1394_source = bayer_source_buffer;
4416 if (debayerParams.offsetY == 1)
4418 dc1394_source += m_Statistics.width;
4423 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4424 debayerParams.filter,
4425 debayerParams.method);
4427 if (error_code != DC1394_SUCCESS)
4429 m_LastError =
i18n(
"Debayer failed (%1)", error_code);
4430 m_Statistics.channels = 1;
4431 delete[] destinationBuffer;
4435 if (m_ImageBufferSize != rgb_size)
4437 delete[] m_ImageBuffer;
4440 m_ImageBuffer =
new uint8_t[rgb_size];
4442 catch (
const std::bad_alloc &e)
4444 delete[] destinationBuffer;
4445 logOOMError(rgb_size);
4446 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4450 m_ImageBufferSize = rgb_size;
4453 auto bayered_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
4457 uint8_t * rBuff = bayered_buffer;
4458 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4459 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4461 int imax = m_Statistics.samples_per_channel * 3 - 3;
4462 for (
int i = 0; i <= imax; i += 3)
4464 *rBuff++ = bayer_destination_buffer[i];
4465 *gBuff++ = bayer_destination_buffer[i + 1];
4466 *bBuff++ = bayer_destination_buffer[i + 2];
4472 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4473 m_Statistics.dataType = TBYTE;
4474 delete[] destinationBuffer;
4478bool FITSData::debayer_16bit()
4480 dc1394error_t error_code;
4482 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
4483 uint8_t *destinationBuffer =
nullptr;
4486 destinationBuffer =
new uint8_t[rgb_size];
4488 catch (
const std::bad_alloc &e)
4490 logOOMError(rgb_size);
4491 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4495 auto * bayer_source_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
4496 auto * bayer_destination_buffer =
reinterpret_cast<uint16_t *
>(destinationBuffer);
4498 if (bayer_destination_buffer ==
nullptr)
4500 logOOMError(rgb_size);
4501 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
4505 int ds1394_height = m_Statistics.height;
4506 auto dc1394_source = bayer_source_buffer;
4508 if (debayerParams.offsetY == 1)
4510 dc1394_source += m_Statistics.width;
4515 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
4516 debayerParams.filter,
4517 debayerParams.method, 16);
4519 if (error_code != DC1394_SUCCESS)
4521 m_LastError =
i18n(
"Debayer failed (%1)");
4522 m_Statistics.channels = 1;
4523 delete[] destinationBuffer;
4527 if (m_ImageBufferSize != rgb_size)
4529 delete[] m_ImageBuffer;
4532 m_ImageBuffer =
new uint8_t[rgb_size];
4534 catch (
const std::bad_alloc &e)
4536 logOOMError(rgb_size);
4537 delete[] destinationBuffer;
4538 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
4542 m_ImageBufferSize = rgb_size;
4545 auto bayered_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
4549 uint16_t * rBuff = bayered_buffer;
4550 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
4551 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
4553 int imax = m_Statistics.samples_per_channel * 3 - 3;
4554 for (
int i = 0; i <= imax; i += 3)
4556 *rBuff++ = bayer_destination_buffer[i];
4557 *gBuff++ = bayer_destination_buffer[i + 1];
4558 *bBuff++ = bayer_destination_buffer[i + 2];
4561 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
4562 m_Statistics.dataType = TUSHORT;
4563 delete[] destinationBuffer;
4567void FITSData::logOOMError(uint32_t requiredMemory)
4569 qCCritical(KSTARS_FITS) <<
"Debayed memory allocation failure. Required Memory:" <<
KFormat().
formatByteSize(requiredMemory)
4570 <<
"Available system memory:" << KSUtils::getAvailableRAM();
4573double FITSData::getADU()
const
4576 for (
int i = 0; i < m_Statistics.channels; i++)
4577 adu += m_Statistics.mean[i];
4579 return (adu /
static_cast<double>(m_Statistics.channels));
4582QString FITSData::getLastError()
const
4587template <
typename T>
4588void FITSData::convertToQImage(
double dataMin,
double dataMax,
double scale,
double zero,
QImage &image)
4590 const auto * buffer =
reinterpret_cast<const T *
>(getImageBuffer());
4591 const T limit = std::numeric_limits<T>::max();
4592 T bMin = dataMin < 0 ? 0 : dataMin;
4593 T bMax = dataMax > limit ? limit : dataMax;
4594 uint16_t w = width();
4595 uint16_t h = height();
4596 uint32_t size = w * h;
4599 if (channels() == 1)
4602 for (
int j = 0; j < h; j++)
4604 unsigned char * scanLine = image.
scanLine(j);
4606 for (
int i = 0; i < w; i++)
4608 val = qBound(bMin, buffer[j * w + i], bMax);
4609 val = val * scale + zero;
4610 scanLine[i] = qBound<unsigned char>(0,
static_cast<uint8_t
>(val), 255);
4616 double rval = 0, gval = 0, bval = 0;
4619 for (
int j = 0; j < h; j++)
4621 auto * scanLine =
reinterpret_cast<QRgb *
>((image.
scanLine(j)));
4623 for (
int i = 0; i < w; i++)
4625 rval = qBound(bMin, buffer[j * w + i], bMax);
4626 gval = qBound(bMin, buffer[j * w + i + size], bMax);
4627 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
4629 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
4631 scanLine[i] = value;
4648 if (future.
result() ==
false)
4651 data.getMinMax(&min, &max);
4659 if (data.channels() == 1)
4664 for (
int i = 0; i < 256; i++)
4665 fitsImage.
setColor(i, qRgb(i, i, i));
4672 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
4673 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
4675 double bscale = 255. / (dataMax - dataMin);
4676 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
4679 switch (data.m_Statistics.dataType)
4682 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4686 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4690 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4694 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4698 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4702 data.convertToQImage<
float>(dataMin, dataMax, bscale, bzero, fitsImage);
4706 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4710 data.convertToQImage<
double>(dataMin, dataMax, bscale, bzero, fitsImage);
4724 qCCritical(KSTARS_FITS) <<
"Failed to convert" << filename <<
"to FITS since format" << format <<
"is not supported in Qt";
4732 qCCritical(KSTARS_FITS) <<
"Failed to open image" << filename;
4736 output =
QDir(getTemporaryPath()).
filePath(filename +
".fits");
4739 fitsfile *fptr =
nullptr;
4741 long fpixel = 1, naxis = input.
allGray() ? 2 : 3, nelements, exposure;
4742 long naxes[3] = { input.
width(), input.
height(), naxis == 3 ? 3 : 1 };
4743 char error_status[512] = {0};
4745 if (fits_create_file(&fptr,
QString(
'!' + output).toLocal8Bit().data(), &status))
4747 fits_get_errstatus(status, error_status);
4748 qCCritical(KSTARS_FITS) <<
"Failed to create FITS file. Error:" << error_status;
4752 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
4754 qCWarning(KSTARS_FITS) <<
"fits_create_img failed:" << error_status;
4756 fits_flush_file(fptr, &status);
4757 fits_close_file(fptr, &status);
4762 fits_update_key(fptr, TLONG,
"EXPOSURE", &exposure,
"Total Exposure Time", &status);
4767 nelements = naxes[0] * naxes[1];
4768 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.
bits(), &status))
4770 fits_get_errstatus(status, error_status);
4771 qCWarning(KSTARS_FITS) <<
"fits_write_img GRAY failed:" << error_status;
4773 fits_flush_file(fptr, &status);
4774 fits_close_file(fptr, &status);
4781 nelements = naxes[0] * naxes[1] * 3;
4783 uint8_t *srcBuffer = input.
bits();
4785 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
4787 uint8_t *rgbBuffer =
new uint8_t[nelements];
4788 if (rgbBuffer ==
nullptr)
4790 qCWarning(KSTARS_FITS) <<
"Not enough memory for RGB buffer";
4791 fits_flush_file(fptr, &status);
4792 fits_close_file(fptr, &status);
4796 uint8_t *subR = rgbBuffer;
4797 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
4798 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
4799 for (uint32_t i = 0; i < srcBytes; i += 4)
4801 *subB++ = srcBuffer[i];
4802 *subG++ = srcBuffer[i + 1];
4803 *subR++ = srcBuffer[i + 2];
4806 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
4808 fits_get_errstatus(status, error_status);
4809 qCWarning(KSTARS_FITS) <<
"fits_write_img RGB failed:" << error_status;
4811 fits_flush_file(fptr, &status);
4812 fits_close_file(fptr, &status);
4813 delete [] rgbBuffer;
4817 delete [] rgbBuffer;
4820 if (fits_flush_file(fptr, &status))
4822 fits_get_errstatus(status, error_status);
4823 qCWarning(KSTARS_FITS) <<
"fits_flush_file failed:" << error_status;
4825 fits_close_file(fptr, &status);
4829 if (fits_close_file(fptr, &status))
4831 fits_get_errstatus(status, error_status);
4832 qCWarning(KSTARS_FITS) <<
"fits_close_file failed:" << error_status;
4839void FITSData::injectWCS(
double orientation,
double ra,
double dec,
double pixscale,
bool eastToTheRight)
4843 updateWCSHeaderData(orientation, ra, dec, pixscale, eastToTheRight);
4847 fits_update_key(fptr, TDOUBLE,
"OBJCTRA", &ra,
"Object RA", &status);
4848 fits_update_key(fptr, TDOUBLE,
"OBJCTDEC", &dec,
"Object DEC", &status);
4852 fits_update_key(fptr, TINT,
"EQUINOX", &epoch,
"Equinox", &status);
4854 fits_update_key(fptr, TDOUBLE,
"CRVAL1", &ra,
"CRVAL1", &status);
4855 fits_update_key(fptr, TDOUBLE,
"CRVAL2", &dec,
"CRVAL1", &status);
4857 char radecsys[8] =
"FK5";
4858 char ctype1[16] =
"RA---TAN";
4859 char ctype2[16] =
"DEC--TAN";
4861 fits_update_key(fptr, TSTRING,
"RADECSYS", radecsys,
"RADECSYS", &status);
4862 fits_update_key(fptr, TSTRING,
"CTYPE1", ctype1,
"CTYPE1", &status);
4863 fits_update_key(fptr, TSTRING,
"CTYPE2", ctype2,
"CTYPE2", &status);
4865 double crpix1 = width() / 2.0;
4866 double crpix2 = height() / 2.0;
4868 fits_update_key(fptr, TDOUBLE,
"CRPIX1", &crpix1,
"CRPIX1", &status);
4869 fits_update_key(fptr, TDOUBLE,
"CRPIX2", &crpix2,
"CRPIX2", &status);
4872 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4873 double secpix2 = pixscale;
4875 fits_update_key(fptr, TDOUBLE,
"SECPIX1", &secpix1,
"SECPIX1", &status);
4876 fits_update_key(fptr, TDOUBLE,
"SECPIX2", &secpix2,
"SECPIX2", &status);
4878 double degpix1 = secpix1 / 3600.0;
4879 double degpix2 = secpix2 / 3600.0;
4881 fits_update_key(fptr, TDOUBLE,
"CDELT1", °pix1,
"CDELT1", &status);
4882 fits_update_key(fptr, TDOUBLE,
"CDELT2", °pix2,
"CDELT2", &status);
4885 double rotation = 360 - orientation;
4889 fits_update_key(fptr, TDOUBLE,
"CROTA1", &rotation,
"CROTA1", &status);
4890 fits_update_key(fptr, TDOUBLE,
"CROTA2", &rotation,
"CROTA2", &status);
4893 emit headerChanged();
4898void FITSData::updateWCSHeaderData(
const double orientation,
const double ra,
const double dec,
const double pixscale,
4899 const bool eastToTheRight)
4901 updateRecordValue(
"OBJCTRA", ra,
"Object RA");
4902 updateRecordValue(
"OBJCTDEC", dec,
"Object DEC");
4903 updateRecordValue(
"EQUINOX", 2000,
"Equinox");
4904 updateRecordValue(
"CRVAL1", ra,
"CRVAL1");
4905 updateRecordValue(
"CRVAL2", dec,
"CRVAL2");
4907 updateRecordValue(
"RADECSYS",
"'FK5'",
"RADECSYS");
4908 updateRecordValue(
"CTYPE1",
"'RA---TAN'",
"CTYPE1");
4909 updateRecordValue(
"CTYPE2",
"'DEC--TAN'",
"CTYPE2");
4911 updateRecordValue(
"CRPIX1", m_Statistics.width / 2.0,
"CRPIX1");
4912 updateRecordValue(
"CRPIX2", m_Statistics.height / 2.0,
"CRPIX2");
4914 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4915 double secpix2 = pixscale;
4917 double degpix1 = secpix1 / 3600.0;
4918 double degpix2 = secpix2 / 3600.0;
4919 updateRecordValue(
"CDELT1", degpix1,
"CDELT1");
4920 updateRecordValue(
"CDELT2", degpix2,
"CDELT2");
4923 double rotation = 360 - orientation;
4926 updateRecordValue(
"CROTA1", rotation,
"CROTA1");
4927 updateRecordValue(
"CROTA2", rotation,
"CROTA2");
4930bool FITSData::contains(
const QPointF &point)
const
4932 return (point.
x() >= 0 && point.
y() >= 0 && point.
x() <= m_Statistics.width && point.
y() <= m_Statistics.height);
4935void FITSData::saveStatistics(FITSImage::Statistic &other)
4937 other = m_Statistics;
4940void FITSData::restoreStatistics(FITSImage::Statistic &other)
4942 m_Statistics = other;
4947void FITSData::constructHistogram()
4949 switch (m_Statistics.dataType)
4952 constructHistogramInternal<uint8_t>();
4956 constructHistogramInternal<int16_t>();
4960 constructHistogramInternal<uint16_t>();
4964 constructHistogramInternal<int32_t>();
4968 constructHistogramInternal<uint32_t>();
4972 constructHistogramInternal<float>();
4976 constructHistogramInternal<int64_t>();
4980 constructHistogramInternal<double>();
4988template <
typename T> int32_t FITSData::histogramBinInternal(T value,
int channel)
const
4990 return qMax(
static_cast<T
>(0), qMin(
static_cast<T
>(m_HistogramBinCount),
4991 static_cast<T
>(rint((value - m_Statistics.min[channel]) / m_HistogramBinWidth[channel]))));
4994template <
typename T> int32_t FITSData::histogramBinInternal(
int x,
int y,
int channel)
const
4996 if (!m_ImageBuffer || !isHistogramConstructed())
4998 uint32_t samples = m_Statistics.width * m_Statistics.height;
4999 uint32_t offset = channel * samples;
5000 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
5001 int index = y * m_Statistics.width + x;
5002 const T &sample = buffer[index + offset];
5003 return histogramBinInternal(sample, channel);
5006int32_t FITSData::histogramBin(
int x,
int y,
int channel)
const
5008 switch (m_Statistics.dataType)
5011 return histogramBinInternal<uint8_t>(x, y, channel);
5015 return histogramBinInternal<int16_t>(x, y, channel);
5019 return histogramBinInternal<uint16_t>(x, y, channel);
5023 return histogramBinInternal<int32_t>(x, y, channel);
5027 return histogramBinInternal<uint32_t>(x, y, channel);
5031 return histogramBinInternal<float>(x, y, channel);
5035 return histogramBinInternal<int64_t>(x, y, channel);
5039 return histogramBinInternal<double>(x, y, channel);
5048template <
typename T>
void FITSData::constructHistogramInternal()
5050 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
5051 uint32_t samples = m_Statistics.width * m_Statistics.height;
5052 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
5054 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
5055 if (m_HistogramBinCount <= 0)
5056 m_HistogramBinCount = 256;
5058 for (
int n = 0; n < m_Statistics.channels; n++)
5060 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
5061 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
5062 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
5064 const double minBinSize = (m_Statistics.max[n] > 1.1) ? 1.0 : .0001;
5065 m_HistogramBinWidth[n] = qMax(minBinSize, (m_Statistics.max[n] - m_Statistics.min[n]) / m_HistogramBinCount);
5070 for (
int n = 0; n < m_Statistics.channels; n++)
5074 for (
int i = 0; i < m_HistogramBinCount; i++)
5075 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
5079 for (
int n = 0; n < m_Statistics.channels; n++)
5083 uint32_t offset = n * samples;
5085 for (uint32_t i = 0; i < samples; i += sampleBy)
5087 int32_t
id = histogramBinInternal<T>(buffer[i + offset], n);
5088 m_HistogramFrequency[n][id] += sampleBy;
5098 for (
int n = 0; n < m_Statistics.channels; n++)
5102 uint32_t accumulator = 0;
5103 for (
int i = 0; i < m_HistogramBinCount; i++)
5105 accumulator += m_HistogramFrequency[n][i];
5106 m_CumulativeFrequency[n].replace(i, accumulator);
5117 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
5118 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] /
static_cast<double>
5119 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
5124 qCDebug(KSTARS_FITS) <<
"FITHistogram: JMIndex " << m_JMIndex;
5126 m_HistogramConstructed =
true;
5127 emit histogramReady();
5130void FITSData::recordLastError(
int errorCode)
5132 char fitsErrorMessage[512] = {0};
5133 fits_get_errstatus(errorCode, fitsErrorMessage);
5134 m_LastError = fitsErrorMessage;
5137double FITSData::getAverageMean(
bool roi)
const
5139 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5140 if (ptr->channels == 1)
5141 return ptr->mean[0];
5143 return (ptr->mean[0] + ptr->mean[1] + ptr->mean[2]) / 3.0;
5146double FITSData::getAverageMedian(
bool roi)
const
5148 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5149 if (ptr->channels == 1)
5150 return ptr->median[0];
5152 return (ptr->median[0] + ptr->median[1] + ptr->median[2]) / 3.0;
5155double FITSData::getAverageStdDev(
bool roi)
const
5157 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
5158 if (ptr->channels == 1)
5159 return ptr->stddev[0];
5161 return (ptr->stddev[0] + ptr->stddev[1] + ptr->stddev[2]) / 3.0;
There are several time-dependent values used in position calculations, that are not specific to an ob...
SkyMapComposite * skyComposite()
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day,...
static KStarsDateTime currentDateTimeUtc()
QList< SkyObject * > findObjectsInArea(const SkyPoint &p1, const SkyPoint &p2)
Provides all necessary information about an object in the sky: its coordinates, name(s),...
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.
An angle, stored as degrees, but expressible in many ways.
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
const double & Degrees() const
QString i18n(const char *text, const TYPE &arg...)
char * toString(const EngineQuery &query)
QString name(GameStandardAction id)
QAction * end(const QObject *recvr, const char *slot, QObject *parent)
KIOCORE_EXPORT StatJob * stat(const QUrl &url, JobFlags flags=DefaultFlags)
KIOCORE_EXPORT QString number(KIO::filesize_t size)
void error(QWidget *parent, const QString &text, const QString &title, const KGuiItem &buttonOk, Options options=Notify)
VehicleSection::Type type(QStringView coachNumber, QStringView coachClassification)
KIOCORE_EXPORT QStringList list(const QString &fileClass)
QString label(StandardShortcut id)
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
QByteArray & append(QByteArrayView data)
const_iterator constBegin() const const
const char * constData() const const
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 filePath(const QString &fileName) const const
QString path() const const
virtual qint64 size() const const override
virtual void close() override
qint64 size() const const
QString suffix() const const
bool isRunning() 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
qint64 write(const QByteArray &data)
void append(QList< T > &&value)
const_reference at(qsizetype i) const const
bool contains(const AT &value) const const
qsizetype count() const const
iterator erase(const_iterator begin, const_iterator end)
bool isEmpty() const const
QList< T > mid(qsizetype pos, qsizetype length) const const
qsizetype size() const const
void finished(QNetworkReply *reply)
NetworkError error() const const
QMetaObject::Connection connect(const QObject *sender, PointerToMemberFunction signal, Functor functor)
void 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
bool startsWith(QChar c, Qt::CaseSensitivity cs) const const
double toDouble(bool *ok) const const
int toInt(bool *ok, int base) const const
QByteArray toLatin1() const const
QByteArray toLocal8Bit() const const
QString trimmed() const const
bool contains(QLatin1StringView str, Qt::CaseSensitivity cs) const const
qsizetype indexOf(const QRegularExpression &re, qsizetype from) const const
QTextStream & center(QTextStream &stream)
QFuture< void > map(Iterator begin, Iterator end, MapFunctor &&function)
QFuture< T > run(Function function,...)
virtual QString fileName() const const override
void setFileTemplate(const QString &name)
void setSingleShot(bool singleShot)
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.