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#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
38#if !defined(KSTARS_LITE) && defined(HAVE_LIBRAW)
39#include <libraw/libraw.h>
49#include <fits_debug.h>
51#define ZOOM_DEFAULT 100.0
54#define ZOOM_LOW_INCR 10
55#define ZOOM_HIGH_INCR 50
63const QStringList RAWFormats = {
"cr2",
"cr3",
"crw",
"nef",
"raf",
"dng",
"arw",
"orf" };
65bool FITSData::readableFilename(
const QString &filename)
68 QString extension = info.completeSuffix().toLower();
75FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)
79 qRegisterMetaType<FITSMode>(
"FITSMode");
81 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
82 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
83 debayerParams.offsetX = debayerParams.offsetY = 0;
86 m_CumulativeFrequency.resize(3);
87 m_HistogramBinWidth.resize(3);
88 m_HistogramFrequency.resize(3);
89 m_HistogramIntensity.resize(3);
100 qRegisterMetaType<FITSMode>(
"FITSMode");
102 debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
103 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
104 debayerParams.offsetX = debayerParams.offsetY = 0;
108 this->m_Mode = other->m_Mode;
109 this->m_Statistics.channels = other->m_Statistics.channels;
110 memcpy(&m_Statistics, &(other->m_Statistics),
sizeof(m_Statistics));
111 m_ImageBuffer =
new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel];
112 memcpy(m_ImageBuffer, other->m_ImageBuffer,
113 m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel);
131 if (m_WCSHandle !=
nullptr)
133 wcsvfree(&m_nwcs, &m_WCSHandle);
134 m_WCSHandle =
nullptr;
139 if (starCenters.
count() > 0)
140 qDeleteAll(starCenters);
143 if (m_SkyObjects.
count() > 0)
144 qDeleteAll(m_SkyObjects);
145 m_SkyObjects.
clear();
149 fits_flush_file(fptr, &status);
150 fits_close_file(fptr, &status);
152 m_PackBuffer =
nullptr;
157void FITSData::loadCommon(
const QString &inFilename)
160 qDeleteAll(starCenters);
165 fits_flush_file(fptr, &status);
166 fits_close_file(fptr, &status);
168 m_PackBuffer =
nullptr;
172 m_Filename = inFilename;
175bool FITSData::loadFromBuffer(
const QByteArray &buffer)
179 return privateLoad(buffer);
184 loadCommon(inFilename);
186 m_Extension = info.completeSuffix().toLower();
187 qCDebug(KSTARS_FITS) <<
"Loading file " << m_Filename;
194QString fitsErrorToString(
int status)
196 char error_status[512] = {0};
197 fits_report_error(stderr, status);
198 fits_get_errstatus(status, error_status);
199 QString message = error_status;
204bool FITSData::privateLoad(
const QByteArray &buffer)
208 cacheEccentricity = -1;
211 return loadFITSImage(buffer);
213 return loadXISFImage(buffer);
215 return loadCanonicalImage(buffer);
216 else if (RAWFormats.
contains(m_Extension))
217 return loadRAWImage(buffer);
222bool FITSData::loadFITSImage(
const QByteArray &buffer,
const bool isCompressed)
224 int status = 0, anynull = 0;
227 m_HistogramConstructed =
false;
229 if (m_Extension.
contains(
".fz") || isCompressed)
238 m_compressedFilename = m_Filename;
243 rc = fp_unpack_file_to_fits(m_Filename.
toLocal8Bit().
data(), &fptr, fpvar) == 0;
246 m_Filename = uncompressedFile;
251 size_t m_PackBufferSize = 100000;
253 m_PackBuffer = (uint8_t *)malloc(m_PackBufferSize);
254 rc = fp_unpack_data_to_data(buffer.
data(), buffer.
size(), &m_PackBuffer, &m_PackBufferSize, fpvar) == 0;
258 void *data =
reinterpret_cast<void *
>(m_PackBuffer);
259 if (fits_open_memfile(&fptr, m_Filename.
toLocal8Bit().
data(), READONLY, &data, &m_PackBufferSize, 0,
262 m_LastError =
i18n(
"Error reading fits buffer: %1.", fitsErrorToString(status));
266 m_Statistics.size = m_PackBufferSize;
274 m_PackBuffer =
nullptr;
275 m_LastError =
i18n(
"Failed to unpack compressed fits");
276 qCCritical(KSTARS_FITS) << m_LastError;
280 m_isTemporary =
true;
281 m_isCompressed =
true;
282 m_Statistics.size = fptr->Fptr->logfilesize;
289 if (fits_open_diskfile(&fptr, m_Filename.
toLocal8Bit(), READONLY, &status))
291 m_LastError =
i18n(
"Error opening fits file %1 : %2", m_Filename, fitsErrorToString(status));
294 m_Statistics.size =
QFile(m_Filename).
size();
299 void *temp_buffer =
const_cast<void *
>(
reinterpret_cast<const void *
>(buffer.
data()));
300 size_t temp_size = buffer.
size();
301 if (fits_open_memfile(&fptr, m_Filename.
toLocal8Bit().
data(), READONLY,
302 &temp_buffer, &temp_size, 0,
nullptr, &
status))
304 m_LastError =
i18n(
"Error reading fits buffer: %1", fitsErrorToString(status));
308 m_Statistics.size = temp_size;
311 if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status))
315 m_PackBuffer =
nullptr;
316 m_LastError =
i18n(
"Could not locate image HDU: %1", fitsErrorToString(status));
319 if (fits_get_img_param(fptr, 3, &m_FITSBITPIX, &(m_Statistics.ndim), naxes, &status))
322 m_PackBuffer =
nullptr;
323 m_LastError =
i18n(
"FITS file open error (fits_get_img_param): %1", fitsErrorToString(status));
328 if ((fits_is_compressed_image(fptr, &status) || m_Statistics.ndim <= 0) && !isCompressed)
330 loadCommon(m_Filename);
331 qCDebug(KSTARS_FITS) <<
"Image is compressed. Reloading...";
332 return loadFITSImage(buffer,
true);
335 if (m_Statistics.ndim < 2)
337 m_LastError =
i18n(
"1D FITS images are not supported in KStars.");
338 qCCritical(KSTARS_FITS) << m_LastError;
340 m_PackBuffer =
nullptr;
344 switch (m_FITSBITPIX)
347 m_Statistics.dataType = TBYTE;
348 m_Statistics.bytesPerPixel =
sizeof(uint8_t);
352 m_FITSBITPIX = USHORT_IMG;
353 m_Statistics.dataType = TUSHORT;
354 m_Statistics.bytesPerPixel =
sizeof(int16_t);
357 m_Statistics.dataType = TUSHORT;
358 m_Statistics.bytesPerPixel =
sizeof(uint16_t);
362 m_FITSBITPIX = ULONG_IMG;
363 m_Statistics.dataType = TULONG;
364 m_Statistics.bytesPerPixel =
sizeof(int32_t);
367 m_Statistics.dataType = TULONG;
368 m_Statistics.bytesPerPixel =
sizeof(uint32_t);
371 m_Statistics.dataType = TFLOAT;
372 m_Statistics.bytesPerPixel =
sizeof(float);
375 m_Statistics.dataType = TLONGLONG;
376 m_Statistics.bytesPerPixel =
sizeof(int64_t);
379 m_Statistics.dataType = TDOUBLE;
380 m_Statistics.bytesPerPixel =
sizeof(double);
383 m_LastError =
i18n(
"Bit depth %1 is not supported.", m_FITSBITPIX);
384 qCCritical(KSTARS_FITS) << m_LastError;
388 if (m_Statistics.ndim < 3)
391 if (naxes[0] == 0 || naxes[1] == 0)
393 m_LastError =
i18n(
"Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
394 qCCritical(KSTARS_FITS) << m_LastError;
396 m_PackBuffer =
nullptr;
400 m_Statistics.width = naxes[0];
401 m_Statistics.height = naxes[1];
402 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
403 roiCenter.
setX(m_Statistics.width / 2);
404 roiCenter.
setY(m_Statistics.height / 2);
405 if(m_Statistics.width % 2)
406 roiCenter.
setX(roiCenter.
x() + 1);
407 if(m_Statistics.height % 2)
408 roiCenter.
setY(roiCenter.
y() + 1);
412 m_Statistics.channels = naxes[2];
416 if ( (m_Mode != FITS_NORMAL && m_Mode != FITS_CALIBRATE) || !Options::auto3DCube())
417 m_Statistics.channels = 1;
419 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
420 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
421 if (m_ImageBuffer ==
nullptr)
423 qCWarning(KSTARS_FITS) <<
"FITSData: Not enough memory for image_buffer channel. Requested: "
424 << m_ImageBufferSize <<
" bytes.";
427 m_PackBuffer =
nullptr;
434 long nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
436 if (fits_read_img(fptr, m_Statistics.dataType, 1, nelements,
nullptr, m_ImageBuffer, &anynull, &status))
438 m_LastError =
i18n(
"Error reading image: %1", fitsErrorToString(status));
446 if (getRecordValue(
"DATE-OBS", value) && value.
isValid())
454 if (naxes[2] == 1 && m_Statistics.channels == 1 && Options::autoDebayer() && checkDebayer())
457 if (m_isTemporary && m_TemporaryDataFile.
open())
459 m_TemporaryDataFile.
write(buffer);
460 m_TemporaryDataFile.
close();
461 m_Filename = m_TemporaryDataFile.
fileName();
465 calculateStats(
false,
false);
468 calculateStats(
false,
false);
470 if (m_Mode == FITS_NORMAL || m_Mode == FITS_ALIGN)
473 starsSearched =
false;
478bool FITSData::loadXISFImage(
const QByteArray &buffer)
480 m_HistogramConstructed =
false;
486 LibXISF::XISFReader xisfReader;
493 LibXISF::ByteArray byteArray(buffer.
constData(), buffer.
size());
494 xisfReader.open(byteArray);
497 if (xisfReader.imagesCount() == 0)
499 m_LastError =
i18n(
"File contain no images");
503 const LibXISF::Image &image = xisfReader.getImage(0);
505 switch (image.sampleFormat())
507 case LibXISF::Image::UInt8:
508 m_Statistics.dataType = TBYTE;
509 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt8);
510 m_FITSBITPIX = TBYTE;
512 case LibXISF::Image::UInt16:
513 m_Statistics.dataType = TUSHORT;
514 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt16);
515 m_FITSBITPIX = TUSHORT;
517 case LibXISF::Image::UInt32:
518 m_Statistics.dataType = TULONG;
519 m_Statistics.bytesPerPixel =
sizeof(LibXISF::UInt32);
520 m_FITSBITPIX = TULONG;
522 case LibXISF::Image::Float32:
523 m_Statistics.dataType = TFLOAT;
524 m_Statistics.bytesPerPixel =
sizeof(LibXISF::Float32);
525 m_FITSBITPIX = TFLOAT;
528 m_LastError =
i18n(
"Sample format %1 is not supported.", LibXISF::Image::sampleFormatString(image.sampleFormat()).c_str());
529 qCCritical(KSTARS_FITS) << m_LastError;
533 m_Statistics.width = image.width();
534 m_Statistics.height = image.height();
535 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
536 m_Statistics.channels = image.channelCount();
537 m_Statistics.size = buffer.
size();
538 roiCenter.
setX(m_Statistics.width / 2);
539 roiCenter.
setY(m_Statistics.height / 2);
540 if(m_Statistics.width % 2)
541 roiCenter.
setX(roiCenter.
x() + 1);
542 if(m_Statistics.height % 2)
543 roiCenter.
setY(roiCenter.
y() + 1);
545 m_HeaderRecords.clear();
546 auto &fitsKeywords = image.fitsKeywords();
547 for(
auto &fitsKeyword : fitsKeywords)
551 if (getRecordValue(
"DATE-OBS", value) && value.
isValid())
557 m_ImageBufferSize = image.imageDataSize();
558 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
559 std::memcpy(m_ImageBuffer, image.imageData(), m_ImageBufferSize);
561 calculateStats(
false,
false);
564 catch (LibXISF::Error &error)
566 m_LastError =
i18n(
"XISF file open error: ") +
error.what();
577bool FITSData::saveXISFImage(
const QString &newFilename)
582 LibXISF::XISFWriter xisfWriter;
583 LibXISF::Image image;
584 image.setGeometry(m_Statistics.width, m_Statistics.height, m_Statistics.channels);
585 if (m_Statistics.channels > 1)
586 image.setColorSpace(LibXISF::Image::RGB);
588 switch (m_FITSBITPIX)
591 image.setSampleFormat(LibXISF::Image::UInt8);
594 image.setSampleFormat(LibXISF::Image::UInt16);
597 image.setSampleFormat(LibXISF::Image::UInt32);
600 image.setSampleFormat(LibXISF::Image::Float32);
603 m_LastError =
i18n(
"Bit depth %1 is not supported.", m_FITSBITPIX);
604 qCCritical(KSTARS_FITS) << m_LastError;
608 std::memcpy(image.imageData(), m_ImageBuffer, m_ImageBufferSize);
609 for (
auto &fitsKeyword : m_HeaderRecords)
610 image.addFITSKeyword({fitsKeyword.key.toUtf8().data(), fitsKeyword.value.toString().toUtf8().data(), fitsKeyword.comment.toUtf8().data()});
612 xisfWriter.writeImage(image);
614 m_Filename = newFilename;
616 catch (LibXISF::Error &err)
618 m_LastError =
i18n(
"Error saving XISF image") + err.what();
623 Q_UNUSED(newFilename)
628bool FITSData::loadCanonicalImage(
const QByteArray &buffer)
635 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
639 m_Statistics.size = buffer.
size();
645 qCCritical(KSTARS_FITS) <<
"Failed to open image.";
657 m_FITSBITPIX = BYTE_IMG;
658 switch (m_FITSBITPIX)
661 m_Statistics.dataType = TBYTE;
662 m_Statistics.bytesPerPixel =
sizeof(uint8_t);
666 m_Statistics.dataType = TUSHORT;
667 m_Statistics.bytesPerPixel =
sizeof(int16_t);
670 m_Statistics.dataType = TUSHORT;
671 m_Statistics.bytesPerPixel =
sizeof(uint16_t);
675 m_Statistics.dataType = TULONG;
676 m_Statistics.bytesPerPixel =
sizeof(int32_t);
679 m_Statistics.dataType = TULONG;
680 m_Statistics.bytesPerPixel =
sizeof(uint32_t);
683 m_Statistics.dataType = TFLOAT;
684 m_Statistics.bytesPerPixel =
sizeof(float);
687 m_Statistics.dataType = TLONGLONG;
688 m_Statistics.bytesPerPixel =
sizeof(int64_t);
691 m_Statistics.dataType = TDOUBLE;
692 m_Statistics.bytesPerPixel =
sizeof(double);
695 m_LastError =
QString(
"Bit depth %1 is not supported.").
arg(m_FITSBITPIX);
696 qCCritical(KSTARS_FITS) << m_LastError;
700 m_Statistics.width =
static_cast<uint16_t
>(imageFromFile.
width());
701 m_Statistics.height =
static_cast<uint16_t
>(imageFromFile.
height());
703 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
705 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels *
static_cast<uint16_t
>
706 (m_Statistics.bytesPerPixel);
707 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
708 if (m_ImageBuffer ==
nullptr)
710 m_LastError =
i18n(
"FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
711 qCCritical(KSTARS_FITS) << m_LastError;
716 if (m_Statistics.channels == 1)
718 memcpy(m_ImageBuffer, imageFromFile.
bits(), m_ImageBufferSize);
723 auto debayered_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
724 auto * original_bayered_buffer =
reinterpret_cast<uint8_t *
>(imageFromFile.
bits());
727 uint8_t * rBuff = debayered_buffer;
728 uint8_t * gBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height);
729 uint8_t * bBuff = debayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
731 int imax = m_Statistics.samples_per_channel * 4 - 4;
732 for (
int i = 0; i <= imax; i += 4)
734 *rBuff++ = original_bayered_buffer[i + 2];
735 *gBuff++ = original_bayered_buffer[i + 1];
736 *bBuff++ = original_bayered_buffer[i + 0];
740 calculateStats(
false,
false);
744bool FITSData::loadRAWImage(
const QByteArray &buffer)
746#if !defined(KSTARS_LITE) && !defined(HAVE_LIBRAW)
747 m_LastError =
i18n(
"Unable to find dcraw and cjpeg. Please install the required tools to convert CR2/NEF to JPEG.");
761 m_LastError =
i18n(
"Cannot open file %1: %2", m_Filename, libraw_strerror(ret));
762 RawProcessor.recycle();
771 if ((ret = RawProcessor.open_buffer(
const_cast<void *
>(
reinterpret_cast<const void *
>(buffer.
data())),
772 buffer.
size())) != LIBRAW_SUCCESS)
774 m_LastError =
i18n(
"Cannot open buffer: %1", libraw_strerror(ret));
775 RawProcessor.recycle();
779 m_Statistics.size = buffer.
size();
783 if ((ret = RawProcessor.unpack()) != LIBRAW_SUCCESS)
785 m_LastError =
i18n(
"Cannot unpack_thumb: %1", libraw_strerror(ret));
786 RawProcessor.recycle();
790 if ((ret = RawProcessor.dcraw_process()) != LIBRAW_SUCCESS)
792 m_LastError =
i18n(
"Cannot dcraw_process: %1", libraw_strerror(ret));
793 RawProcessor.recycle();
797 libraw_processed_image_t *image = RawProcessor.dcraw_make_mem_image(&ret);
798 if (ret != LIBRAW_SUCCESS)
800 m_LastError =
i18n(
"Cannot load to memory: %1", libraw_strerror(ret));
801 RawProcessor.recycle();
805 RawProcessor.recycle();
807 m_Statistics.bytesPerPixel = image->bits / 8;
809 if (m_Statistics.bytesPerPixel == 1)
810 m_Statistics.dataType = TBYTE;
812 m_Statistics.dataType = TUSHORT;
813 m_Statistics.width = image->width;
814 m_Statistics.height = image->height;
815 m_Statistics.channels = image->colors;
816 m_Statistics.samples_per_channel = m_Statistics.width * m_Statistics.height;
818 m_ImageBufferSize = m_Statistics.samples_per_channel * m_Statistics.channels * m_Statistics.bytesPerPixel;
819 m_ImageBuffer =
new uint8_t[m_ImageBufferSize];
820 if (m_ImageBuffer ==
nullptr)
822 m_LastError =
i18n(
"FITSData: Not enough memory for image_buffer channel. Requested: %1 bytes ", m_ImageBufferSize);
823 qCCritical(KSTARS_FITS) << m_LastError;
824 libraw_dcraw_clear_mem(image);
829 auto destination_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
830 auto source_buffer =
reinterpret_cast<uint8_t *
>(image->data);
833 if (image->colors == 1)
835 memcpy(destination_buffer, source_buffer, m_ImageBufferSize);
840 uint8_t * rBuff = destination_buffer;
841 uint8_t * gBuff = destination_buffer + (m_Statistics.width * m_Statistics.height);
842 uint8_t * bBuff = destination_buffer + (m_Statistics.width * m_Statistics.height * 2);
844 int imax = m_Statistics.samples_per_channel * 3 - 3;
845 for (
int i = 0; i <= imax; i += 3)
847 *rBuff++ = source_buffer[i + 0];
848 *gBuff++ = source_buffer[i + 1];
849 *bBuff++ = source_buffer[i + 2];
852 libraw_dcraw_clear_mem(image);
854 calculateStats(
false,
false);
859bool FITSData::saveImage(
const QString &newFilename)
861 if (newFilename == m_Filename)
866 if (ext ==
"jpg" || ext ==
"png")
869 QImage fitsImage = FITSToImage(newFilename);
870 getMinMax(&min, &max);
876 else if (channels() == 1)
881 for (
int i = 0; i < 256; i++)
882 fitsImage.
setColor(i, qRgb(i, i, i));
889 double dataMin = m_Statistics.mean[0] - m_Statistics.stddev[0];
890 double dataMax = m_Statistics.mean[0] + m_Statistics.stddev[0] * 3;
892 double bscale = 255. / (dataMax - dataMin);
893 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
896 switch (m_Statistics.dataType)
899 convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
903 convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
907 convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
911 convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
915 convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
919 convertToQImage<float>(dataMin, dataMax, bscale, bzero, fitsImage);
923 convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
927 convertToQImage<double>(dataMin, dataMax, bscale, bzero, fitsImage);
936 m_Filename = newFilename;
937 qCInfo(KSTARS_FITS) <<
"Saved image file:" << m_Filename;
943 m_LastError =
i18n(
"Saving compressed files is not supported.");
955 fits_close_file(fptr, &status);
958 rotCounter = flipHCounter = flipVCounter = 0;
959 return saveXISFImage(newFilename);
998 nelements = m_Statistics.samples_per_channel * m_Statistics.channels;
1001 if (fptr && fits_close_file(fptr, &status))
1003 m_LastError =
i18n(
"Failed to close file: %1", fitsErrorToString(status));
1008 if (fits_create_file(&new_fptr,
QString(
"!%1").arg(newFilename).toLocal8Bit(), &status))
1010 m_LastError =
i18n(
"Failed to create file: %1", fitsErrorToString(status));
1019 long naxis = m_Statistics.channels == 1 ? 2 : 3;
1020 long naxes[3] = {m_Statistics.width, m_Statistics.height, naxis};
1023 if (fits_create_img(fptr, m_FITSBITPIX, naxis, naxes, &status))
1025 m_LastError =
i18n(
"Failed to create image: %1", fitsErrorToString(status));
1032 if (fits_update_key(fptr, TDOUBLE,
"DATAMIN", &(m_Statistics.min),
"Minimum value", &status))
1034 m_LastError =
i18n(
"Failed to update key: %1", fitsErrorToString(status));
1039 if (fits_update_key(fptr, TDOUBLE,
"DATAMAX", &(m_Statistics.max),
"Maximum value", &status))
1041 m_LastError =
i18n(
"Failed to update key: %1", fitsErrorToString(status));
1046 fits_write_key(fptr, TDOUBLE,
"MIN1", &m_Statistics.min[0],
"Min Channel 1", &status);
1049 fits_write_key(fptr, TDOUBLE,
"MIN2", &m_Statistics.min[1],
"Min Channel 2", &status);
1050 fits_write_key(fptr, TDOUBLE,
"MIN3", &m_Statistics.min[2],
"Min Channel 3", &status);
1054 fits_write_key(fptr, TDOUBLE,
"MAX1", &m_Statistics.max[0],
"Max Channel 1", &status);
1057 fits_write_key(fptr, TDOUBLE,
"MAX2", &m_Statistics.max[1],
"Max Channel 2", &status);
1058 fits_write_key(fptr, TDOUBLE,
"MAX3", &m_Statistics.max[2],
"Max Channel 3", &status);
1062 if (m_Statistics.mean[0] > 0)
1064 fits_write_key(fptr, TDOUBLE,
"MEAN1", &m_Statistics.mean[0],
"Mean Channel 1", &status);
1067 fits_write_key(fptr, TDOUBLE,
"MEAN2", &m_Statistics.mean[1],
"Mean Channel 2", &status);
1068 fits_write_key(fptr, TDOUBLE,
"MEAN3", &m_Statistics.mean[2],
"Mean Channel 3", &status);
1073 if (m_Statistics.median[0] > 0)
1075 fits_write_key(fptr, TDOUBLE,
"MEDIAN1", &m_Statistics.median[0],
"Median Channel 1", &status);
1078 fits_write_key(fptr, TDOUBLE,
"MEDIAN2", &m_Statistics.median[1],
"Median Channel 2", &status);
1079 fits_write_key(fptr, TDOUBLE,
"MEDIAN3", &m_Statistics.median[2],
"Median Channel 3", &status);
1084 if (m_Statistics.stddev[0] > 0)
1086 fits_write_key(fptr, TDOUBLE,
"STDDEV1", &m_Statistics.stddev[0],
"Standard Deviation Channel 1", &status);
1089 fits_write_key(fptr, TDOUBLE,
"STDDEV2", &m_Statistics.stddev[1],
"Standard Deviation Channel 2", &status);
1090 fits_write_key(fptr, TDOUBLE,
"STDDEV3", &m_Statistics.stddev[2],
"Standard Deviation Channel 3", &status);
1095 for (
int i = 10; i < m_HeaderRecords.count(); i++)
1097 QString key = m_HeaderRecords[i].key;
1098 const char *comment = m_HeaderRecords[i].comment.toLatin1().constBegin();
1099 QVariant value = m_HeaderRecords[i].value;
1101 switch (value.
type())
1106 fits_write_key(fptr, TINT, key.
toLatin1().
constData(), &number, comment, &status);
1113 fits_write_key(fptr, TDOUBLE, key.
toLatin1().
constData(), &number, comment, &status);
1120 char valueBuffer[256] = {0};
1122 fits_write_key(fptr, TSTRING, key.
toLatin1().
constData(), valueBuffer, comment, &status);
1128 if (fits_write_date(fptr, &status))
1130 m_LastError =
i18n(
"Failed to update date: %1", fitsErrorToString(status));
1137 if (fits_write_history(fptr, history.
toLatin1(), &status))
1139 m_LastError =
i18n(
"Failed to update history: %1", fitsErrorToString(status));
1143 int rot = 0, mirror = 0;
1146 rot = (90 * rotCounter) % 360;
1150 if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0)
1153 if ((rot != 0) || (mirror != 0))
1154 rotWCSFITS(rot, mirror);
1156 rotCounter = flipHCounter = flipVCounter = 0;
1159 if (fits_write_img(fptr, m_Statistics.dataType, 1, nelements, m_ImageBuffer, &status))
1161 m_LastError =
i18n(
"Failed to write image: %1", fitsErrorToString(status));
1165 m_Filename = newFilename;
1167 fits_flush_file(fptr, &status);
1169 qCInfo(KSTARS_FITS) <<
"Saved FITS file:" << m_Filename;
1174void FITSData::clearImageBuffers()
1176 delete[] m_ImageBuffer;
1177 m_ImageBuffer =
nullptr;
1178 if(m_ImageRoiBuffer !=
nullptr )
1180 delete[] m_ImageRoiBuffer;
1181 m_ImageRoiBuffer =
nullptr;
1187void FITSData::makeRoiBuffer(
QRect roi)
1192 uint32_t channelSize = roi.
height() * roi.
width();
1193 if(channelSize > m_Statistics.samples_per_channel || channelSize == 1)
1197 if(m_ImageRoiBuffer !=
nullptr )
1199 delete[] m_ImageRoiBuffer;
1200 m_ImageRoiBuffer =
nullptr;
1202 int xoffset = roi.
topLeft().
x() - 1;
1203 int yoffset = roi.
topLeft().
y() - 1;
1204 uint32_t bpp = m_Statistics.bytesPerPixel;
1205 m_ImageRoiBuffer =
new uint8_t[roi.
height()*roi.
width()*m_Statistics.channels * m_Statistics.bytesPerPixel]();
1206 for(
int n = 0 ; n < m_Statistics.channels ; n++)
1208 for(
int i = 0; i < roi.
height(); i++)
1210 size_t i1 = n * channelSize * bpp + i * roi.
width() * bpp;
1211 size_t i2 = n * m_Statistics.samples_per_channel * bpp + (yoffset + i) * width() * bpp + xoffset * bpp;
1212 memcpy(&m_ImageRoiBuffer[i1],
1218 memcpy(&m_ROIStatistics, &m_Statistics,
sizeof(FITSImage::Statistic));
1219 m_ROIStatistics.samples_per_channel = roi.
height() * roi.
width();
1220 m_ROIStatistics.width = roi.
width();
1221 m_ROIStatistics.height = roi.
height();
1222 calculateStats(
false,
true);
1224void FITSData::calculateStats(
bool refresh,
bool roi)
1229 calculateMinMax(refresh);
1230 calculateMedian(refresh);
1233 if (refresh ==
false && fptr)
1237 if (fits_read_key_dbl(fptr,
"MEAN1", &m_Statistics.mean[0],
nullptr, &status) == 0)
1240 fits_read_key_dbl(fptr,
"MEAN2", & m_Statistics.mean[1],
nullptr, &status);
1241 fits_read_key_dbl(fptr,
"MEAN3", &m_Statistics.mean[2],
nullptr, &status);
1244 if (fits_read_key_dbl(fptr,
"STDDEV1", &m_Statistics.stddev[0],
nullptr, &status) == 0)
1247 fits_read_key_dbl(fptr,
"STDDEV2", &m_Statistics.stddev[1],
nullptr, &status);
1248 fits_read_key_dbl(fptr,
"STDDEV3", &m_Statistics.stddev[2],
nullptr, &status);
1256 switch (m_Statistics.dataType)
1259 calculateStdDev<uint8_t>();
1263 calculateStdDev<int16_t>();
1267 calculateStdDev<uint16_t>();
1271 calculateStdDev<int32_t>();
1275 calculateStdDev<uint32_t>();
1279 calculateStdDev<float>();
1283 calculateStdDev<int64_t>();
1287 calculateStdDev<double>();
1295 m_Statistics.SNR = m_Statistics.mean[0] / m_Statistics.stddev[0];
1299 calculateMinMax(refresh, roi);
1300 calculateMedian(refresh, roi);
1302 switch (m_ROIStatistics.dataType)
1305 calculateStdDev<uint8_t>(roi);
1309 calculateStdDev<int16_t>(roi);
1313 calculateStdDev<uint16_t>(roi);
1317 calculateStdDev<int32_t>(roi);
1321 calculateStdDev<uint32_t>(roi);
1325 calculateStdDev<float>(roi);
1329 calculateStdDev<int64_t>(roi);
1333 calculateStdDev<double>(roi);
1342void FITSData::calculateMinMax(
bool refresh,
bool roi)
1352 if (fptr !=
nullptr && !refresh)
1357 if (fits_read_key_dbl(fptr,
"DATAMIN", &(m_Statistics.min[0]),
nullptr, &status) == 0)
1359 else if (fits_read_key_dbl(fptr,
"MIN1", &(m_Statistics.min[0]),
nullptr, &status) == 0)
1363 fits_read_key_dbl(fptr,
"MIN2", &m_Statistics.min[1],
nullptr, &status);
1364 fits_read_key_dbl(fptr,
"MIN3", &m_Statistics.min[2],
nullptr, &status);
1368 if (fits_read_key_dbl(fptr,
"DATAMAX", &(m_Statistics.max[0]),
nullptr, &status) == 0)
1370 else if (fits_read_key_dbl(fptr,
"MAX1", &(m_Statistics.max[0]),
nullptr, &status) == 0)
1374 fits_read_key_dbl(fptr,
"MAX2", &m_Statistics.max[1],
nullptr, &status);
1375 fits_read_key_dbl(fptr,
"MAX3", &m_Statistics.max[2],
nullptr, &status);
1378 if (nfound == 2 && !(m_Statistics.min[0] == 0 && m_Statistics.max[0] == 0))
1382 m_Statistics.min[0] = 1.0E30;
1383 m_Statistics.max[0] = -1.0E30;
1385 m_Statistics.min[1] = 1.0E30;
1386 m_Statistics.max[1] = -1.0E30;
1388 m_Statistics.min[2] = 1.0E30;
1389 m_Statistics.max[2] = -1.0E30;
1391 switch (m_Statistics.dataType)
1394 calculateMinMax<uint8_t>();
1398 calculateMinMax<int16_t>();
1402 calculateMinMax<uint16_t>();
1406 calculateMinMax<int32_t>();
1410 calculateMinMax<uint32_t>();
1414 calculateMinMax<float>();
1418 calculateMinMax<int64_t>();
1422 calculateMinMax<double>();
1431 m_ROIStatistics.min[0] = 1.0E30;
1432 m_ROIStatistics.max[0] = -1.0E30;
1434 m_ROIStatistics.min[1] = 1.0E30;
1435 m_ROIStatistics.max[1] = -1.0E30;
1437 m_ROIStatistics.min[2] = 1.0E30;
1438 m_ROIStatistics.max[2] = -1.0E30;
1440 switch (m_Statistics.dataType)
1443 calculateMinMax<uint8_t>(roi);
1447 calculateMinMax<int16_t>(roi);
1451 calculateMinMax<uint16_t>(roi);
1455 calculateMinMax<int32_t>(roi);
1459 calculateMinMax<uint32_t>(roi);
1463 calculateMinMax<float>(roi);
1467 calculateMinMax<int64_t>(roi);
1471 calculateMinMax<double>(roi);
1481void FITSData::calculateMedian(
bool refresh,
bool roi)
1491 if (fptr !=
nullptr && !refresh)
1494 if (fits_read_key_dbl(fptr,
"MEDIAN1", &m_Statistics.median[0],
nullptr, &status) == 0)
1498 fits_read_key_dbl(fptr,
"MEDIAN2", &m_Statistics.median[1],
nullptr, &status);
1499 fits_read_key_dbl(fptr,
"MEDIAN3", &m_Statistics.median[2],
nullptr, &status);
1505 m_Statistics.median[RED_CHANNEL] = 0;
1506 m_Statistics.median[GREEN_CHANNEL] = 0;
1507 m_Statistics.median[BLUE_CHANNEL] = 0;
1509 switch (m_Statistics.dataType)
1512 calculateMedian<uint8_t>();
1516 calculateMedian<int16_t>();
1520 calculateMedian<uint16_t>();
1524 calculateMedian<int32_t>();
1528 calculateMedian<uint32_t>();
1532 calculateMedian<float>();
1536 calculateMedian<int64_t>();
1540 calculateMedian<double>();
1549 m_ROIStatistics.median[RED_CHANNEL] = 0;
1550 m_ROIStatistics.median[GREEN_CHANNEL] = 0;
1551 m_ROIStatistics.median[BLUE_CHANNEL] = 0;
1553 switch (m_ROIStatistics.dataType)
1556 calculateMedian<uint8_t>(roi);
1560 calculateMedian<int16_t>(roi);
1564 calculateMedian<uint16_t>(roi);
1568 calculateMedian<int32_t>(roi);
1572 calculateMedian<uint32_t>(roi);
1576 calculateMedian<float>(roi);
1580 calculateMedian<int64_t>(roi);
1584 calculateMedian<double>(roi);
1594template <
typename T>
1595void FITSData::calculateMedian(
bool roi)
1597 auto * buffer =
reinterpret_cast<T *
>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1598 const uint32_t maxMedianSize = 500000;
1599 uint32_t medianSize = roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel;
1600 uint8_t downsample = 1;
1601 if (medianSize > maxMedianSize)
1603 downsample = (
static_cast<double>(medianSize) / maxMedianSize) + 0.999;
1604 medianSize /= downsample;
1610 std::vector<int32_t> samples;
1611 samples.reserve(medianSize);
1613 for (uint8_t n = 0; n < m_Statistics.channels; n++)
1615 auto *oneChannel = buffer + n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1616 for (uint32_t upto = 0; upto < (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1619 auto median = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_MEDIAN, samples);
1620 roi ? m_ROIStatistics.median[n] = median : m_Statistics.median[n] = median;
1624template <
typename T>
1625QPair<T, T> FITSData::getParitionMinMax(uint32_t start, uint32_t stride,
bool roi)
1627 auto * buffer =
reinterpret_cast<T *
>(roi ? m_ImageRoiBuffer : m_ImageBuffer);
1628 T min = std::numeric_limits<T>::max();
1629 T max = std::numeric_limits<T>::min();
1631 uint32_t
end = start + stride;
1633 for (uint32_t i = start; i <
end; i++)
1635 min = qMin(buffer[i], min);
1636 max = qMax(buffer[i], max);
1643 return qMakePair(min, max);
1646template <
typename T>
1647void FITSData::calculateMinMax(
bool roi)
1649 T min = std::numeric_limits<T>::max();
1650 T max = std::numeric_limits<T>::min();
1654 const uint8_t nThreads = 16;
1656 for (
int n = 0; n < m_Statistics.channels; n++)
1658 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1661 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1664 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1665 (tStride * nThreads));
1668 uint32_t tStart = cStart;
1673 for (
int i = 0; i < nThreads; i++)
1676 futures.
append(
QtConcurrent::run(
this, &FITSData::getParitionMinMax<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride,
1682 for (
int i = 0; i < nThreads; i++)
1684 QPair<T, T> result = futures[i].result();
1685 min = qMin(result.first, min);
1686 max = qMax(result.second, max);
1691 m_Statistics.min[n] = min;
1692 m_Statistics.max[n] = max;
1696 m_ROIStatistics.min[n] = min;
1697 m_ROIStatistics.max[n] = max;
1710 SumData(
double s,
double sq,
int n) : sum(s), squaredSum(sq), numSamples(n) {}
1711 SumData() : sum(0), squaredSum(0), numSamples(0) {}
1714template <
typename T>
1715SumData getSumAndSquaredSum(uint32_t start, uint32_t stride, uint8_t *buff)
1717 auto * buffer =
reinterpret_cast<T *
>(buff);
1718 const uint32_t
end = start + stride;
1720 double squaredSum = 0;
1721 for (uint32_t i = start; i <
end; i++)
1723 double sample = buffer[i];
1725 squaredSum += sample * sample;
1727 const double numSamples =
end - start;
1728 return SumData(sum, squaredSum, numSamples);
1731template <
typename T>
1732void FITSData::calculateStdDev(
bool roi )
1735 const uint8_t nThreads = 16;
1737 for (
int n = 0; n < m_Statistics.channels; n++)
1739 uint32_t cStart = n * (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel);
1742 uint32_t tStride = (roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) / nThreads;
1745 uint32_t fStride = tStride + ((roi ? m_ROIStatistics.samples_per_channel : m_Statistics.samples_per_channel) -
1746 (tStride * nThreads));
1749 uint32_t tStart = cStart;
1754 for (
int i = 0; i < nThreads; i++)
1757 uint8_t *buff = roi ? m_ImageRoiBuffer : m_ImageBuffer;
1759 (i == (nThreads - 1)) ? fStride : tStride, buff));
1764 double sum = 0, squared_sum = 0;
1765 double numSamples = 0;
1766 for (
int i = 0; i < nThreads; i++)
1768 const SumData result = futures[i].result();
1770 squared_sum += result.squaredSum;
1771 numSamples += result.numSamples;
1774 if (numSamples <= 0)
continue;
1775 const double mean = sum / numSamples;
1776 const double variance = squared_sum / numSamples - mean * mean;
1779 m_Statistics.mean[n] = mean;
1780 m_Statistics.stddev[n] = sqrt(variance);
1784 m_ROIStatistics.mean[n] = mean;
1785 m_ROIStatistics.stddev[n] = sqrt(variance);
1793 kernel.fill(0.0, size * size);
1795 double kernelSum = 0.0;
1796 int fOff = (size - 1) / 2;
1797 double normal = 1.0 / (2.0 * M_PI * sigma * sigma);
1798 for (
int y = -fOff; y <= fOff; y++)
1800 for (
int x = -fOff; x <= fOff; x++)
1802 double distance = ((y * y) + (x * x)) / (2.0 * sigma * sigma);
1803 int index = (y + fOff) * size + (x + fOff);
1804 kernel[index] = normal * qExp(-distance);
1805 kernelSum += kernel.at(index);
1808 for (
int y = 0; y < size; y++)
1810 for (
int x = 0; x < size; x++)
1812 int index = y * size + x;
1813 kernel[index] = kernel.at(index) * 1.0 / kernelSum;
1820template <
typename T>
1821void FITSData::convolutionFilter(
const QVector<double> &kernel,
int kernelSize)
1823 T * imagePtr =
reinterpret_cast<T *
>(m_ImageBuffer);
1829 int fOff = (kernelSize - 1) / 2;
1833 for (
int offsetY = 0; offsetY < m_Statistics.height; offsetY++)
1835 for (
int offsetX = 0; offsetX < m_Statistics.width; offsetX++)
1840 int byteOffset = offsetY * m_Statistics.width + offsetX;
1843 for (
int filterY = -fOff; filterY <= fOff; filterY++)
1845 for (
int filterX = -fOff; filterX <= fOff; filterX++)
1847 if ((offsetY + filterY) >= 0 && (offsetY + filterY) < m_Statistics.height
1848 && ((offsetX + filterX) >= 0 && (offsetX + filterX) < m_Statistics.width ))
1851 int calcOffset = byteOffset + filterX + filterY * m_Statistics.width;
1852 int index = (filterY + fOff) * kernelSize + (filterX + fOff);
1853 double kernelValue = kernel.
at(index);
1854 gt += (imagePtr[calcOffset]) * kernelValue;
1860 imagePtr[byteOffset] = gt;
1865template <
typename T>
1866void FITSData::gaussianBlur(
int kernelSize,
double sigma)
1869 if (kernelSize % 2 == 0)
1872 qCInfo(KSTARS_FITS) <<
"Warning, size must be an odd number, correcting size to " << kernelSize;
1880 QVector<double> gaussianKernel = createGaussianKernel(kernelSize, sigma);
1881 convolutionFilter<T>(gaussianKernel, kernelSize);
1884void FITSData::setMinMax(
double newMin,
double newMax, uint8_t channel)
1886 m_Statistics.min[channel] = newMin;
1887 m_Statistics.max[channel] = newMax;
1890bool FITSData::parseHeader()
1892 char * header =
nullptr;
1893 int status = 0, nkeys = 0;
1895 if (fits_hdr2str(fptr, 0,
nullptr, 0, &header, &nkeys, &status))
1897 fits_report_error(stderr, status);
1902 m_HeaderRecords.clear();
1905 for (
int i = 0; i < nkeys; i++)
1915 oneRecord.comment =
properties[0].mid(8).simplified();
1920 oneRecord.value =
properties[1].simplified();
1922 oneRecord.comment =
properties[2].simplified();
1929 oneRecord.value.toInt(&ok);
1935 oneRecord.value.toDouble(&ok);
1941 m_HeaderRecords.append(oneRecord);
1949bool FITSData::getRecordValue(
const QString &key,
QVariant &value)
const
1951 auto result = std::find_if(m_HeaderRecords.begin(), m_HeaderRecords.end(), [&key](
const Record & oneRecord)
1953 return (oneRecord.key == key && oneRecord.value.isValid());
1956 if (result != m_HeaderRecords.end())
1958 value = (*result).
value;
1964bool FITSData::parseSolution(FITSImage::Solution &solution)
const
1967 bool raOK =
false, deOK =
false, coordOK =
false, scaleOK =
false;
1971 solution.fieldWidth = solution.fieldHeight = solution.pixscale = solution.ra = solution.dec = -1;
1974 if (getRecordValue(
"OBJCTRA", value))
1977 solution.ra = angleValue.
Degrees();
1980 else if (getRecordValue(
"RA", value))
1982 solution.ra = value.
toDouble(&raOK);
1986 if (getRecordValue(
"OBJCTDEC", value))
1989 solution.dec = angleValue.
Degrees();
1992 else if (getRecordValue(
"DEC", value))
1994 solution.dec = value.
toDouble(&deOK);
1997 coordOK = raOK && deOK;
2001 if (getRecordValue(
"SCALE", value))
2006 double focal_length = -1;
2007 if (getRecordValue(
"FOCALLEN", value))
2012 double pixsize1 = -1, pixsize2 = -1;
2014 if (getRecordValue(
"PIXSIZE1", value))
2019 if (getRecordValue(
"PIXSIZE2", value))
2024 int binx = 1, biny = 1;
2026 if (getRecordValue(
"XBINNING", value))
2031 if (getRecordValue(
"YBINNING", value))
2036 if (pixsize1 > 0 && pixsize2 > 0)
2042 solution.pixscale = scale;
2044 solution.fieldWidth = (m_Statistics.width * scale) / 60.0;
2046 solution.fieldHeight = (m_Statistics.height * scale * (pixsize2 / pixsize1)) / 60.0;
2048 else if (focal_length > 0)
2051 solution.fieldWidth = ((206264.8062470963552 * m_Statistics.width * (pixsize1 / 1000.0)) / (focal_length * binx)) / 60.0;
2053 solution.fieldHeight = ((206264.8062470963552 * m_Statistics.height * (pixsize2 / 1000.0)) / (focal_length * biny)) / 60.0;
2055 solution.pixscale = (206264.8062470963552 * (pixsize1 / 1000.0)) / (focal_length * binx);
2061 return (coordOK || scaleOK);
2069 starAlgorithm = algorithm;
2070 qDeleteAll(starCenters);
2071 starCenters.
clear();
2072 starsSearched =
true;
2078 m_StarDetector.
reset(
new FITSSEPDetector(
this));
2079 m_StarDetector->setSettings(m_SourceExtractorSettings);
2080 if (m_Mode == FITS_NORMAL && trackingBox.
isNull())
2082 if (Options::quickHFR())
2085 const int w = getStatistics().width;
2086 const int h = getStatistics().height;
2087 QRect middle(
static_cast<int>(w * 0.25),
static_cast<int>(h * 0.25), w / 2, h / 2);
2088 m_StarFindFuture = m_StarDetector->findSources(middle);
2089 return m_StarFindFuture;
2092 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2093 return m_StarFindFuture;
2096 case ALGORITHM_GRADIENT:
2099 m_StarDetector.
reset(
new FITSGradientDetector(
this));
2100 m_StarDetector->setSettings(m_SourceExtractorSettings);
2101 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2102 return m_StarFindFuture;
2105 case ALGORITHM_CENTROID:
2108 m_StarDetector.
reset(
new FITSCentroidDetector(
this));
2109 m_StarDetector->setSettings(m_SourceExtractorSettings);
2111 if (!isHistogramConstructed())
2112 constructHistogram();
2113 m_StarDetector->configure(
"JMINDEX", m_JMIndex);
2114 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2115 return m_StarFindFuture;
2119 m_StarDetector.
reset(
new FITSCentroidDetector(
this));
2120 m_StarFindFuture = starDetector->findSources(trackingBox);
2121 return m_StarFindFuture;
2125 case ALGORITHM_THRESHOLD:
2127 m_StarDetector.
reset(
new FITSThresholdDetector(
this));
2128 m_StarDetector->setSettings(m_SourceExtractorSettings);
2129 m_StarDetector->configure(
"THRESHOLD_PERCENTAGE", Options::focusThreshold());
2130 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2131 return m_StarFindFuture;
2134 case ALGORITHM_BAHTINOV:
2136 m_StarDetector.
reset(
new FITSBahtinovDetector(
this));
2137 m_StarDetector->setSettings(m_SourceExtractorSettings);
2138 m_StarDetector->configure(
"NUMBER_OF_AVERAGE_ROWS", Options::focusMultiRowAverage());
2139 m_StarFindFuture = m_StarDetector->findSources(trackingBox);
2140 return m_StarFindFuture;
2147 if (mask.
isNull() ==
false)
2149 starCenters.
erase(std::remove_if(starCenters.
begin(), starCenters.
end(),
2152 return (mask->isVisible(edge->x, edge->y) == false);
2153 }), starCenters.
end());
2156 return starCenters.
count();
2160double FITSData::getHFR(HFRType type)
2162 if (starCenters.
empty())
2165 if (cacheHFR >= 0 && cacheHFRType == type)
2168 m_SelectedHFRStar.invalidate();
2170 if (type == HFR_MAX)
2175 for (
int i = 0; i < starCenters.
count(); i++)
2177 if (starCenters[i]->HFR > maxVal)
2180 maxVal = starCenters[i]->HFR;
2184 m_SelectedHFRStar = *starCenters[maxIndex];
2185 cacheHFR = starCenters[maxIndex]->HFR;
2186 cacheHFRType =
type;
2189 else if (type == HFR_HIGH)
2192 int minX = width() / 10;
2193 int minY = height() / 10;
2194 int maxX = width() - minX;
2195 int maxY = height() - minY;
2196 starCenters.
erase(std::remove_if(starCenters.
begin(), starCenters.
end(), [minX, minY, maxX, maxY](Edge * oneStar)
2198 return (oneStar->x < minX || oneStar->x > maxX || oneStar->y < minY || oneStar->y > maxY);
2199 }), starCenters.
end());
2201 if (starCenters.
empty())
2204 m_SelectedHFRStar = *starCenters[
static_cast<int>(starCenters.
size() * 0.05)];
2205 cacheHFR = m_SelectedHFRStar.HFR;
2206 cacheHFRType =
type;
2209 else if (type == HFR_MEDIAN)
2211 std::nth_element(starCenters.
begin(), starCenters.
begin() + starCenters.
size() / 2, starCenters.
end());
2212 m_SelectedHFRStar = *starCenters[starCenters.
size() / 2];
2214 cacheHFR = m_SelectedHFRStar.HFR;
2215 cacheHFRType =
type;
2223 const int saturationValue = m_Statistics.dataType == TBYTE ? 250 : 50000;
2224 int numSaturated = 0;
2225 for (
auto center : starCenters)
2226 if (
center->val > saturationValue)
2228 const bool removeSaturatedStars = starCenters.
size() - numSaturated > 20;
2229 if (removeSaturatedStars && numSaturated > 0)
2230 qCDebug(KSTARS_FITS) <<
"Removing " << numSaturated <<
" stars from HFR calculation";
2232 std::vector<double> HFRs;
2234 for (
auto center : starCenters)
2236 if (removeSaturatedStars &&
center->val > saturationValue)
continue;
2238 if (type == HFR_AVERAGE)
2239 HFRs.push_back(
center->HFR);
2245 if (m_SkyBackground.mean <= 0.0 ||
center->val < m_SkyBackground.mean)
2247 HFRs.push_back(
center->HFR);
2248 qCDebug(KSTARS_FITS) <<
"HFR Adj, sky background " << m_SkyBackground.mean <<
" star peak " <<
center->val <<
2253 const double a_div_b =
center->val / m_SkyBackground.mean;
2254 const double factor = erf(sqrt(log(a_div_b))) / (1 - (1 / a_div_b));
2255 HFRs.push_back(
center->HFR * factor);
2256 qCDebug(KSTARS_FITS) <<
"HFR Adj, brightness adjusted from " <<
center->HFR <<
" to " <<
center->HFR * factor;
2261 auto m = Mathematics::RobustStatistics::ComputeLocation(Mathematics::RobustStatistics::LOCATION_SIGMACLIPPING,
2265 cacheHFRType = HFR_AVERAGE;
2269double FITSData::getHFR(
int x,
int y,
double scale)
2271 if (starCenters.
empty())
2274 for (
int i = 0; i < starCenters.
count(); i++)
2276 const int maxDist = std::max(2,
static_cast<int>(0.5 + 2 * starCenters[i]->width / scale));
2277 const int dx = std::fabs(starCenters[i]->x - x);
2278 const int dy = std::fabs(starCenters[i]->y - y);
2279 if (dx <= maxDist && dy <= maxDist)
2281 m_SelectedHFRStar = *starCenters[i];
2282 return starCenters[i]->HFR;
2289double FITSData::getEccentricity()
2291 if (starCenters.
empty())
2293 if (cacheEccentricity >= 0)
2294 return cacheEccentricity;
2295 std::vector<float> eccs;
2296 for (
const auto &s : starCenters)
2297 eccs.push_back(s->ellipticity);
2298 int middle = eccs.size() / 2;
2299 std::nth_element(eccs.begin(), eccs.begin() + middle, eccs.end());
2300 float medianEllipticity = eccs[middle];
2306 const float eccentricity = sqrt(medianEllipticity * (2 - medianEllipticity));
2307 cacheEccentricity = eccentricity;
2308 return eccentricity;
2313 if (type == FITS_NONE)
2326 case FITS_AUTO_STRETCH:
2328 for (
int i = 0; i < 3; i++)
2330 dataMin[i] = m_Statistics.mean[i] - m_Statistics.stddev[i];
2331 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2336 case FITS_HIGH_CONTRAST:
2338 for (
int i = 0; i < 3; i++)
2340 dataMin[i] = m_Statistics.mean[i] + m_Statistics.stddev[i];
2341 dataMax[i] = m_Statistics.mean[i] + m_Statistics.stddev[i] * 3;
2346 case FITS_HIGH_PASS:
2348 for (
int i = 0; i < 3; i++)
2350 dataMin[i] = m_Statistics.mean[i];
2359 switch (m_Statistics.dataType)
2363 for (
int i = 0; i < 3; i++)
2365 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2366 dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i];
2368 applyFilter<uint8_t>(type, image, &dataMin, &dataMax);
2374 for (
int i = 0; i < 3; i++)
2376 dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i];
2377 dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i];
2379 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2386 for (
int i = 0; i < 3; i++)
2388 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2389 dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i];
2391 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2397 for (
int i = 0; i < 3; i++)
2399 dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i];
2400 dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i];
2402 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2408 for (
int i = 0; i < 3; i++)
2410 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
2411 dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i];
2413 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
2419 for (
int i = 0; i < 3; i++)
2421 dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i];
2422 dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i];
2424 applyFilter<float>(type, image, &dataMin, &dataMax);
2430 for (
int i = 0; i < 3; i++)
2432 dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i];
2433 dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i];
2436 applyFilter<long>(type, image, &dataMin, &dataMax);
2442 for (
int i = 0; i < 3; i++)
2444 dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i];
2445 dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i];
2447 applyFilter<double>(type, image, &dataMin, &dataMax);
2464template <
typename T>
2467 bool calcStats =
false;
2468 T * image =
nullptr;
2471 image =
reinterpret_cast<T *
>(targetImage);
2474 image =
reinterpret_cast<T *
>(m_ImageBuffer);
2479 for (
int i = 0; i < 3; i++)
2481 min[i] = (*targetMin)[i] < std::numeric_limits<T>::min() ? std::numeric_limits<T>::min() : (*targetMin)[i];
2482 max[i] = (*targetMax)[i] > std::numeric_limits<T>::max() ? std::numeric_limits<T>::max() : (*targetMax)[i];
2487 const uint8_t nThreads = 16;
2489 uint32_t width = m_Statistics.width;
2490 uint32_t height = m_Statistics.height;
2498 case FITS_AUTO_STRETCH:
2499 case FITS_HIGH_CONTRAST:
2502 case FITS_HIGH_PASS:
2508 if (type == FITS_LOG)
2510 for (
int i = 0; i < 3; i++)
2511 coeff[i] = max[i] / std::log(1 + max[i]);
2513 else if (type == FITS_SQRT)
2515 for (
int i = 0; i < 3; i++)
2516 coeff[i] = max[i] / sqrt(max[i]);
2519 for (
int n = 0; n < m_Statistics.channels; n++)
2521 if (type == FITS_HIGH_PASS)
2522 min[n] = m_Statistics.mean[n];
2524 uint32_t cStart = n * m_Statistics.samples_per_channel;
2527 uint32_t tStride = m_Statistics.samples_per_channel / nThreads;
2530 uint32_t fStride = tStride + (m_Statistics.samples_per_channel - (tStride * nThreads));
2532 T * runningBuffer = image + cStart;
2534 if (type == FITS_LOG)
2536 for (
int i = 0; i < nThreads; i++)
2539 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2542 a = qBound(min[n],
static_cast<T
>(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]);
2545 runningBuffer += tStride;
2548 else if (type == FITS_SQRT)
2550 for (
int i = 0; i < nThreads; i++)
2553 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2556 a = qBound(min[n],
static_cast<T
>(round(coeff[n] * a)), max[n]);
2560 runningBuffer += tStride;
2564 for (
int i = 0; i < nThreads; i++)
2567 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max,
2570 a = qBound(min[n], a, max[n]);
2572 runningBuffer += tStride;
2577 for (
int i = 0; i < nThreads * m_Statistics.channels; i++)
2578 futures[i].waitForFinished();
2582 for (
int i = 0; i < 3; i++)
2584 m_Statistics.min[i] = min[i];
2585 m_Statistics.max[i] = max[i];
2587 calculateStdDev<T>();
2595 if (!isHistogramConstructed())
2596 constructHistogram();
2600 double coeff = 255.0 / (height * width);
2604 for (
int i = 0; i < m_Statistics.channels; i++)
2606 uint32_t offset = i * m_Statistics.samples_per_channel;
2607 for (uint32_t j = 0; j < height; j++)
2609 row = offset + j * width;
2610 for (uint32_t k = 0; k < width; k++)
2613 bufferVal = (image[index] - min[i]) / m_HistogramBinWidth[i];
2615 if (bufferVal >= m_CumulativeFrequency[i].size())
2616 bufferVal = m_CumulativeFrequency[i].size() - 1;
2618 image[index] = qBound(min[i],
static_cast<T
>(round(coeff * m_CumulativeFrequency[i][bufferVal])), max[i]);
2625 calculateStats(
true,
false);
2631 uint8_t BBP = m_Statistics.bytesPerPixel;
2632 auto * extension =
new T[(width + 2) * (height + 2)];
2637 for (uint32_t ch = 0; ch < m_Statistics.channels; ch++)
2639 uint32_t offset = ch * m_Statistics.samples_per_channel;
2640 uint32_t N = width, M = height;
2642 for (uint32_t i = 0; i < M; ++i)
2644 memcpy(extension + (N + 2) * (i + 1) + 1, image + (N * i) + offset, N * BBP);
2645 extension[(N + 2) * (i + 1)] = image[N * i + offset];
2646 extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1 + offset];
2649 memcpy(extension, extension + N + 2, (N + 2) * BBP);
2651 memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP);
2658 for (uint32_t m = 1; m < M - 1; ++m)
2659 for (uint32_t n = 1; n < N - 1; ++n)
2665 memset(&window[0], 0, 9 *
sizeof(
float));
2666 for (uint32_t j = m - 1; j < m + 2; ++j)
2667 for (uint32_t i = n - 1; i < n + 2; ++i)
2668 window[k++] = extension[j * N + i];
2670 for (uint32_t j = 0; j < 5; ++j)
2674 for (uint32_t l = j + 1; l < 9; ++l)
2675 if (window[l] < window[mine])
2678 const float temp =
window[j];
2683 image[(m - 1) * (N - 2) + n - 1 + offset] =
window[4];
2691 calculateStdDev<T>();
2696 gaussianBlur<T>(Options::focusGaussianKernelSize(), Options::focusGaussianSigma());
2698 calculateStats(
true,
false);
2701 case FITS_ROTATE_CW:
2706 case FITS_ROTATE_CCW:
2711 case FITS_MOUNT_FLIP_H:
2716 case FITS_MOUNT_FLIP_V:
2729 for (
int i = 0; i < starCenters.
count(); i++)
2731 int x =
static_cast<int>(starCenters[i]->x);
2732 int y =
static_cast<int>(starCenters[i]->y);
2735 starCentersInSubFrame.
append(starCenters[i]);
2738 return starCentersInSubFrame;
2741bool FITSData::loadWCS()
2743#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2745 if (m_WCSState == Success)
2747 qCWarning(KSTARS_FITS) <<
"WCS data already loaded";
2751 if (m_WCSHandle !=
nullptr)
2753 wcsvfree(&m_nwcs, &m_WCSHandle);
2755 m_WCSHandle =
nullptr;
2758 qCDebug(KSTARS_FITS) <<
"Started WCS Data Processing...";
2762 int nkeyrec = 0, nreject = 0;
2765 char *header =
nullptr;
2766 if (fits_hdr2str(fptr, 1,
nullptr, 0, &header, &nkeyrec, &status))
2769 fits_get_errstatus(status, errmsg);
2770 m_LastError = errmsg;
2771 m_WCSState = Failure;
2775 fits_free_memory(header, &status);
2780 for(
auto &fitsKeyword : m_HeaderRecords)
2783 rec.
append(fitsKeyword.key.leftJustified(8,
' ').toLatin1());
2785 rec.
append(fitsKeyword.value.toByteArray());
2787 rec.
append(fitsKeyword.comment.toLatin1());
2794 if ((status = wcspih(header_str.
data(), nkeyrec, WCSHDR_all, 0, &nreject, &m_nwcs, &m_WCSHandle)) != 0)
2796 wcsvfree(&m_nwcs, &m_WCSHandle);
2797 m_WCSHandle =
nullptr;
2799 m_LastError =
QString(
"wcspih ERROR %1: %2.").
arg(status).
arg(wcshdr_errmsg[status]);
2800 m_WCSState = Failure;
2804 if (m_WCSHandle ==
nullptr)
2806 m_LastError =
i18n(
"No world coordinate systems found.");
2807 m_WCSState = Failure;
2812 if (m_WCSHandle->crpix[0] == 0)
2814 wcsvfree(&m_nwcs, &m_WCSHandle);
2815 m_WCSHandle =
nullptr;
2817 m_LastError =
i18n(
"No world coordinate systems found.");
2818 m_WCSState = Failure;
2823 if ((status = wcsset(m_WCSHandle)) != 0)
2825 wcsvfree(&m_nwcs, &m_WCSHandle);
2826 m_WCSHandle =
nullptr;
2828 m_LastError =
QString(
"wcsset error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2829 m_WCSState = Failure;
2833 m_ObjectsSearched =
false;
2834 m_WCSState = Success;
2837 qCDebug(KSTARS_FITS) <<
"Finished WCS Data processing...";
2847#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2849 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, worldcrd[NWCSFIX];
2851 if (m_WCSHandle ==
nullptr)
2853 m_LastError =
i18n(
"No world coordinate systems found.");
2860 if ((status = wcss2p(m_WCSHandle, 1, 2, worldcrd, &phi, &theta, imgcrd, pixcrd, stat)) != 0)
2862 m_LastError =
QString(
"wcss2p error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2866 wcsImagePoint.
setX(imgcrd[0]);
2867 wcsImagePoint.
setY(imgcrd[1]);
2869 wcsPixelPoint.
setX(pixcrd[0]);
2870 wcsPixelPoint.
setY(pixcrd[1]);
2875 Q_UNUSED(wcsPixelPoint);
2876 Q_UNUSED(wcsImagePoint);
2881bool FITSData::pixelToWCS(
const QPointF &wcsPixelPoint,
SkyPoint &wcsCoord)
2883#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2885 double imgcrd[NWCSFIX], phi, pixcrd[NWCSFIX], theta, world[NWCSFIX];
2887 if (m_WCSHandle ==
nullptr)
2889 m_LastError =
i18n(
"No world coordinate systems found.");
2893 pixcrd[0] = wcsPixelPoint.
x();
2894 pixcrd[1] = wcsPixelPoint.
y();
2896 if ((status = wcsp2s(m_WCSHandle, 1, 2, pixcrd, imgcrd, &phi, &theta, world, stat)) != 0)
2898 m_LastError =
QString(
"wcsp2s error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2903 wcsCoord.
setRA0(world[0] / 15.0);
2909 Q_UNUSED(wcsPixelPoint);
2915#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2916bool FITSData::searchObjects()
2918 if (m_ObjectsSearched)
2921 m_ObjectsSearched =
true;
2926 pixelToWCS(
QPointF(0, 0), startPoint);
2927 pixelToWCS(
QPointF(width() - 1, height() - 1), endPoint);
2929 return findObjectsInImage(startPoint, endPoint);
2932bool FITSData::findWCSBounds(
double &minRA,
double &maxRA,
double &minDec,
double &maxDec)
2934 if (m_WCSHandle ==
nullptr)
2936 m_LastError =
i18n(
"No world coordinate systems found.");
2945 auto updateMinMax = [&](
int x,
int y)
2948 double imgcrd[2], phi, pixcrd[2], theta, world[2];
2953 if (wcsp2s(m_WCSHandle, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0]))
2956 minRA = std::min(minRA, world[0]);
2957 maxRA = std::max(maxRA, world[0]);
2958 minDec = std::min(minDec, world[1]);
2959 maxDec = std::max(maxDec, world[1]);
2963 for (
int y = 0; y < height(); y++)
2966 updateMinMax(width() - 1, y);
2969 for (
int x = 1; x < width() - 1; x++)
2972 updateMinMax(x, height() - 1);
2979 if (wcsToPixel(NCP, pPoint, imagePoint))
2981 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
2986 if (wcsToPixel(SCP, pPoint, imagePoint))
2988 if (pPoint.
x() > 0 && pPoint.
x() < width() && pPoint.
y() > 0 && pPoint.
y() < height())
2997#if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
3000 if (KStarsData::Instance() ==
nullptr)
3008 if (getRecordValue(
"DATE-OBS", date))
3011 tsString = tsString.remove(
'\'').trimmed();
3023 num =
new KSNumbers(KStarsData::Instance()->ut().djd());
3028 m_SkyObjects.
clear();
3033 int type = oneObject->type();
3034 return (type == SkyObject::STAR || type == SkyObject::PLANET || type == SkyObject::ASTEROID ||
3035 type == SkyObject::COMET || type == SkyObject::SUPERNOVA || type == SkyObject::MOON ||
3036 type == SkyObject::SATELLITE);
3039 double world[2], phi, theta, imgcrd[2], pixcrd[2];
3041 for (
auto &
object :
list)
3043 world[0] =
object->ra0().Degrees();
3044 world[1] =
object->dec0().Degrees();
3046 if (wcss2p(m_WCSHandle, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0]) == 0)
3051 if (x > 0 && y > 0 && x < w && y < h)
3052 m_SkyObjects.
append(
new FITSSkyObject(
object, x, y));
3061int FITSData::getFlipVCounter()
const
3063 return flipVCounter;
3066void FITSData::setFlipVCounter(
int value)
3068 flipVCounter = value;
3071int FITSData::getFlipHCounter()
const
3073 return flipHCounter;
3076void FITSData::setFlipHCounter(
int value)
3078 flipHCounter = value;
3081int FITSData::getRotCounter()
const
3086void FITSData::setRotCounter(
int value)
3096template <
typename T>
3097bool FITSData::rotFITS(
int rotate,
int mirror)
3101 uint8_t * rotimage =
nullptr;
3106 else if (rotate == 2)
3108 else if (rotate == 3)
3110 else if (rotate < 0)
3111 rotate = rotate + 360;
3113 nx = m_Statistics.width;
3114 ny = m_Statistics.height;
3116 int BBP = m_Statistics.bytesPerPixel;
3119 rotimage =
new uint8_t[m_Statistics.samples_per_channel * m_Statistics.channels * BBP];
3121 if (rotimage ==
nullptr)
3123 qWarning() <<
"Unable to allocate memory for rotated image buffer!";
3127 auto * rotBuffer =
reinterpret_cast<T *
>(rotimage);
3128 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
3131 if (rotate < 45 && rotate > -45)
3135 for (
int i = 0; i < m_Statistics.channels; i++)
3137 offset = m_Statistics.samples_per_channel * i;
3138 for (x1 = 0; x1 < nx; x1++)
3141 for (y1 = 0; y1 < ny; y1++)
3142 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3146 else if (mirror == 2)
3148 for (
int i = 0; i < m_Statistics.channels; i++)
3150 offset = m_Statistics.samples_per_channel * i;
3151 for (y1 = 0; y1 < ny; y1++)
3154 for (x1 = 0; x1 < nx; x1++)
3155 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3161 for (
int i = 0; i < m_Statistics.channels; i++)
3163 offset = m_Statistics.samples_per_channel * i;
3164 for (y1 = 0; y1 < ny; y1++)
3166 for (x1 = 0; x1 < nx; x1++)
3167 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3174 else if (rotate >= 45 && rotate < 135)
3178 for (
int i = 0; i < m_Statistics.channels; i++)
3180 offset = m_Statistics.samples_per_channel * i;
3181 for (y1 = 0; y1 < ny; y1++)
3184 for (x1 = 0; x1 < nx; x1++)
3187 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3192 else if (mirror == 2)
3194 for (
int i = 0; i < m_Statistics.channels; i++)
3196 offset = m_Statistics.samples_per_channel * i;
3197 for (y1 = 0; y1 < ny; y1++)
3199 for (x1 = 0; x1 < nx; x1++)
3200 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3206 for (
int i = 0; i < m_Statistics.channels; i++)
3208 offset = m_Statistics.samples_per_channel * i;
3209 for (y1 = 0; y1 < ny; y1++)
3212 for (x1 = 0; x1 < nx; x1++)
3215 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3221 m_Statistics.width = ny;
3222 m_Statistics.height = nx;
3226 else if (rotate >= 135 && rotate < 225)
3230 for (
int i = 0; i < m_Statistics.channels; i++)
3232 offset = m_Statistics.samples_per_channel * i;
3233 for (y1 = 0; y1 < ny; y1++)
3236 for (x1 = 0; x1 < nx; x1++)
3237 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
3241 else if (mirror == 2)
3243 for (
int i = 0; i < m_Statistics.channels; i++)
3245 offset = m_Statistics.samples_per_channel * i;
3246 for (x1 = 0; x1 < nx; x1++)
3249 for (y1 = 0; y1 < ny; y1++)
3250 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3256 for (
int i = 0; i < m_Statistics.channels; i++)
3258 offset = m_Statistics.samples_per_channel * i;
3259 for (y1 = 0; y1 < ny; y1++)
3262 for (x1 = 0; x1 < nx; x1++)
3265 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3273 else if (rotate >= 225 && rotate < 315)
3277 for (
int i = 0; i < m_Statistics.channels; i++)
3279 offset = m_Statistics.samples_per_channel * i;
3280 for (y1 = 0; y1 < ny; y1++)
3282 for (x1 = 0; x1 < nx; x1++)
3283 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
3287 else if (mirror == 2)
3289 for (
int i = 0; i < m_Statistics.channels; i++)
3291 offset = m_Statistics.samples_per_channel * i;
3292 for (y1 = 0; y1 < ny; y1++)
3295 for (x1 = 0; x1 < nx; x1++)
3298 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3305 for (
int i = 0; i < m_Statistics.channels; i++)
3307 offset = m_Statistics.samples_per_channel * i;
3308 for (y1 = 0; y1 < ny; y1++)
3311 for (x1 = 0; x1 < nx; x1++)
3314 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3320 m_Statistics.width = ny;
3321 m_Statistics.height = nx;
3325 else if (rotate >= 315 && mirror)
3327 for (
int i = 0; i < m_Statistics.channels; i++)
3329 offset = m_Statistics.samples_per_channel * i;
3330 for (y1 = 0; y1 < ny; y1++)
3332 for (x1 = 0; x1 < nx; x1++)
3336 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
3342 delete[] m_ImageBuffer;
3343 m_ImageBuffer = rotimage;
3348void FITSData::rotWCSFITS(
int angle,
int mirror)
3352 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
3353 int WCS_DECIMALS = 6;
3355 naxis1 = m_Statistics.width;
3356 naxis2 = m_Statistics.height;
3358 if (fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
3367 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3368 fits_update_key_dbl(fptr,
"CROTA1", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3370 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
3371 fits_update_key_dbl(fptr,
"CROTA2", ctemp1 + 90.0, WCS_DECIMALS, comment, &status);
3379 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
3380 fits_update_key_dbl(fptr,
"CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
3382 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
3383 fits_update_key_dbl(fptr,
"CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
3387 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
3388 fits_update_key_dbl(fptr,
"LTM1_1", -ctemp1, WCS_DECIMALS, comment, &status);
3392 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
3393 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
3395 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp1, comment, &status))
3396 fits_update_key_dbl(fptr,
"CD1_2", -ctemp1, WCS_DECIMALS, comment, &status);
3398 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp1, comment, &status))
3399 fits_update_key_dbl(fptr,
"CD2_1", -ctemp1, WCS_DECIMALS, comment, &status);
3405 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
3409 if (!fits_read_key_dbl(fptr,
"LTM2_2", &ctemp2, comment, &status))
3410 if (ctemp1 == ctemp2)
3415 if (!fits_read_key_dbl(fptr,
"LTV1", <v1, comment, &status))
3416 fits_delete_key(fptr,
"LTV1", &status);
3417 if (!fits_read_key_dbl(fptr,
"LTV2", <v2, comment, &status))
3418 fits_delete_key(fptr,
"LTV2", &status);
3422 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp3, comment, &status))
3423 fits_update_key_dbl(fptr,
"CRPIX1", (ctemp3 - ltv1) / ctemp1, WCS_DECIMALS, comment, &status);
3425 if (!fits_read_key_dbl(fptr,
"CRPIX2", &ctemp3, comment, &status))
3426 fits_update_key_dbl(fptr,
"CRPIX2", (ctemp3 - ltv2) / ctemp1, WCS_DECIMALS, comment, &status);
3430 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp3, comment, &status))
3431 fits_update_key_dbl(fptr,
"CD1_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3433 if (!fits_read_key_dbl(fptr,
"CD1_2", &ctemp3, comment, &status))
3434 fits_update_key_dbl(fptr,
"CD1_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3436 if (!fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status))
3437 fits_update_key_dbl(fptr,
"CD2_1", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3439 if (!fits_read_key_dbl(fptr,
"CD2_2", &ctemp3, comment, &status))
3440 fits_update_key_dbl(fptr,
"CD2_2", ctemp3 / ctemp1, WCS_DECIMALS, comment, &status);
3444 fits_delete_key(fptr,
"LTM1_1", &status);
3445 fits_delete_key(fptr,
"LTM1_2", &status);
3453 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp1, comment, &status) &&
3454 !fits_read_key_dbl(fptr,
"CRPIX2", &ctemp2, comment, &status))
3459 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3460 else if (angle == 90)
3462 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3463 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3465 else if (angle == 180)
3467 fits_update_key_dbl(fptr,
"CRPIX1", ctemp1, WCS_DECIMALS, comment, &status);
3468 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3470 else if (angle == 270)
3472 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
3473 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
3480 fits_update_key_dbl(fptr,
"CRPIX1", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3481 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
3483 else if (angle == 180)
3485 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3486 fits_update_key_dbl(fptr,
"CRPIX2", naxis2 - ctemp2, WCS_DECIMALS, comment, &status);
3488 else if (angle == 270)
3490 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
3491 fits_update_key_dbl(fptr,
"CRPIX2", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3499 if (!fits_read_key_dbl(fptr,
"CDELT1", &ctemp1, comment, &status) &&
3500 !fits_read_key_dbl(fptr,
"CDELT2", &ctemp2, comment, &status))
3505 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3506 else if (angle == 90)
3508 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3509 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3511 else if (angle == 180)
3513 fits_update_key_dbl(fptr,
"CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
3514 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3516 else if (angle == 270)
3518 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3519 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3526 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3527 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3529 else if (angle == 180)
3531 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3532 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3534 else if (angle == 270)
3536 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3537 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3548 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
3550 fits_read_key_dbl(fptr,
"CD1_2", &ctemp2, comment, &status);
3551 fits_read_key_dbl(fptr,
"CD2_1", &ctemp3, comment, &status);
3552 fits_read_key_dbl(fptr,
"CD2_2", &ctemp4, comment, &status);
3558 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
3559 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3561 else if (angle == 90)
3563 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
3564 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
3565 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
3566 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
3568 else if (angle == 180)
3570 fits_update_key_dbl(fptr,
"CD1_1", ctemp1, WCS_DECIMALS, comment, &status);
3571 fits_update_key_dbl(fptr,
"CD1_2", ctemp2, WCS_DECIMALS, comment, &status);
3572 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3573 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
3575 else if (angle == 270)
3577 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
3578 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
3579 fits_update_key_dbl(fptr,
"CD2_1", ctemp3, WCS_DECIMALS, comment, &status);
3580 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
3587 fits_update_key_dbl(fptr,
"CD1_1", -ctemp4, WCS_DECIMALS, comment, &status);
3588 fits_update_key_dbl(fptr,
"CD1_2", -ctemp3, WCS_DECIMALS, comment, &status);
3589 fits_update_key_dbl(fptr,
"CD2_1", ctemp1, WCS_DECIMALS, comment, &status);
3590 fits_update_key_dbl(fptr,
"CD2_2", ctemp1, WCS_DECIMALS, comment, &status);
3592 else if (angle == 180)
3594 fits_update_key_dbl(fptr,
"CD1_1", -ctemp1, WCS_DECIMALS, comment, &status);
3595 fits_update_key_dbl(fptr,
"CD1_2", -ctemp2, WCS_DECIMALS, comment, &status);
3596 fits_update_key_dbl(fptr,
"CD2_1", -ctemp3, WCS_DECIMALS, comment, &status);
3597 fits_update_key_dbl(fptr,
"CD2_2", -ctemp4, WCS_DECIMALS, comment, &status);
3599 else if (angle == 270)
3601 fits_update_key_dbl(fptr,
"CD1_1", ctemp4, WCS_DECIMALS, comment, &status);
3602 fits_update_key_dbl(fptr,
"CD1_2", ctemp3, WCS_DECIMALS, comment, &status);
3603 fits_update_key_dbl(fptr,
"CD2_1", -ctemp2, WCS_DECIMALS, comment, &status);
3604 fits_update_key_dbl(fptr,
"CD2_2", -ctemp1, WCS_DECIMALS, comment, &status);
3612 if (!fits_read_key_dbl(fptr,
"CO1_1", &ctemp1, comment, &status))
3617 for (i = 1; i < 13; i++)
3619 sprintf(keyword,
"CO1_%d", i);
3620 fits_delete_key(fptr, keyword, &status);
3622 for (i = 1; i < 13; i++)
3624 sprintf(keyword,
"CO2_%d", i);
3625 fits_delete_key(fptr, keyword, &status);
3631uint8_t * FITSData::getWritableImageBuffer()
3633 return m_ImageBuffer;
3636uint8_t
const * FITSData::getImageBuffer()
const
3638 return m_ImageBuffer;
3641void FITSData::setImageBuffer(uint8_t * buffer)
3643 delete[] m_ImageBuffer;
3644 m_ImageBuffer = buffer;
3647bool FITSData::checkDebayer()
3650 char bayerPattern[64], roworder[64];
3653 if (fits_read_keyword(fptr,
"BAYERPAT", bayerPattern,
nullptr, &status))
3656 if (m_Statistics.dataType != TUSHORT && m_Statistics.dataType != TBYTE)
3658 m_LastError =
i18n(
"Only 8 and 16 bits bayered images supported.");
3661 QString pattern(bayerPattern);
3662 pattern = pattern.remove(
'\'').trimmed();
3665 order = order.remove(
'\'').trimmed();
3667 if (order ==
"BOTTOM-UP" && !(m_Statistics.height % 2))
3669 if (pattern ==
"RGGB")
3671 else if (pattern ==
"GBRG")
3673 else if (pattern ==
"GRBG")
3675 else if (pattern ==
"BGGR")
3680 if (pattern ==
"RGGB")
3681 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
3682 else if (pattern ==
"GBRG")
3683 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
3684 else if (pattern ==
"GRBG")
3685 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
3686 else if (pattern ==
"BGGR")
3687 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
3691 m_LastError =
i18n(
"Unsupported bayer pattern %1.", pattern);
3695 fits_read_key(fptr, TINT,
"XBAYROFF", &debayerParams.offsetX,
nullptr, &status);
3696 fits_read_key(fptr, TINT,
"YBAYROFF", &debayerParams.offsetY,
nullptr, &status);
3698 if (debayerParams.offsetX == 1)
3703 switch (debayerParams.filter)
3705 case DC1394_COLOR_FILTER_RGGB:
3706 debayerParams.filter = DC1394_COLOR_FILTER_GRBG;
3708 case DC1394_COLOR_FILTER_GBRG:
3709 debayerParams.filter = DC1394_COLOR_FILTER_BGGR;
3711 case DC1394_COLOR_FILTER_GRBG:
3712 debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
3714 case DC1394_COLOR_FILTER_BGGR:
3715 debayerParams.filter = DC1394_COLOR_FILTER_GBRG;
3718 debayerParams.offsetX = 0;
3720 if (debayerParams.offsetX != 0 || debayerParams.offsetY > 1 || debayerParams.offsetY < 0)
3722 m_LastError =
i18n(
"Unsupported bayer offsets %1 %2.", debayerParams.offsetX, debayerParams.offsetY);
3731void FITSData::getBayerParams(BayerParams * param)
3733 param->method = debayerParams.method;
3734 param->filter = debayerParams.filter;
3735 param->offsetX = debayerParams.offsetX;
3736 param->offsetY = debayerParams.offsetY;
3739void FITSData::setBayerParams(BayerParams * param)
3741 debayerParams.method = param->method;
3742 debayerParams.filter = param->filter;
3743 debayerParams.offsetX = param->offsetX;
3744 debayerParams.offsetY = param->offsetY;
3747bool FITSData::debayer(
bool reload)
3751 int anynull = 0,
status = 0;
3753 if (fits_read_img(fptr, m_Statistics.dataType, 1, m_Statistics.samples_per_channel,
nullptr, m_ImageBuffer,
3763 switch (m_Statistics.dataType)
3766 return debayer_8bit();
3769 return debayer_16bit();
3776bool FITSData::debayer_8bit()
3778 dc1394error_t error_code;
3780 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
3781 uint8_t * destinationBuffer =
nullptr;
3785 destinationBuffer =
new uint8_t[rgb_size];
3787 catch (
const std::bad_alloc &e)
3789 logOOMError(rgb_size);
3790 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
3794 auto * bayer_source_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
3795 auto * bayer_destination_buffer =
reinterpret_cast<uint8_t *
>(destinationBuffer);
3797 if (bayer_destination_buffer ==
nullptr)
3799 logOOMError(rgb_size);
3800 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
3804 int ds1394_height = m_Statistics.height;
3805 auto dc1394_source = bayer_source_buffer;
3807 if (debayerParams.offsetY == 1)
3809 dc1394_source += m_Statistics.width;
3814 error_code = dc1394_bayer_decoding_8bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
3815 debayerParams.filter,
3816 debayerParams.method);
3818 if (error_code != DC1394_SUCCESS)
3820 m_LastError =
i18n(
"Debayer failed (%1)", error_code);
3821 m_Statistics.channels = 1;
3822 delete[] destinationBuffer;
3826 if (m_ImageBufferSize != rgb_size)
3828 delete[] m_ImageBuffer;
3831 m_ImageBuffer =
new uint8_t[rgb_size];
3833 catch (
const std::bad_alloc &e)
3835 delete[] destinationBuffer;
3836 logOOMError(rgb_size);
3837 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
3841 m_ImageBufferSize = rgb_size;
3844 auto bayered_buffer =
reinterpret_cast<uint8_t *
>(m_ImageBuffer);
3848 uint8_t * rBuff = bayered_buffer;
3849 uint8_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
3850 uint8_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
3852 int imax = m_Statistics.samples_per_channel * 3 - 3;
3853 for (
int i = 0; i <= imax; i += 3)
3855 *rBuff++ = bayer_destination_buffer[i];
3856 *gBuff++ = bayer_destination_buffer[i + 1];
3857 *bBuff++ = bayer_destination_buffer[i + 2];
3863 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
3864 m_Statistics.dataType = TBYTE;
3865 delete[] destinationBuffer;
3869bool FITSData::debayer_16bit()
3871 dc1394error_t error_code;
3873 uint32_t rgb_size = m_Statistics.samples_per_channel * 3 * m_Statistics.bytesPerPixel;
3874 uint8_t *destinationBuffer =
nullptr;
3877 destinationBuffer =
new uint8_t[rgb_size];
3879 catch (
const std::bad_alloc &e)
3881 logOOMError(rgb_size);
3882 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
3886 auto * bayer_source_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
3887 auto * bayer_destination_buffer =
reinterpret_cast<uint16_t *
>(destinationBuffer);
3889 if (bayer_destination_buffer ==
nullptr)
3891 logOOMError(rgb_size);
3892 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer.");
3896 int ds1394_height = m_Statistics.height;
3897 auto dc1394_source = bayer_source_buffer;
3899 if (debayerParams.offsetY == 1)
3901 dc1394_source += m_Statistics.width;
3906 error_code = dc1394_bayer_decoding_16bit(dc1394_source, bayer_destination_buffer, m_Statistics.width, ds1394_height,
3907 debayerParams.filter,
3908 debayerParams.method, 16);
3910 if (error_code != DC1394_SUCCESS)
3912 m_LastError =
i18n(
"Debayer failed (%1)");
3913 m_Statistics.channels = 1;
3914 delete[] destinationBuffer;
3918 if (m_ImageBufferSize != rgb_size)
3920 delete[] m_ImageBuffer;
3923 m_ImageBuffer =
new uint8_t[rgb_size];
3925 catch (
const std::bad_alloc &e)
3927 logOOMError(rgb_size);
3928 delete[] destinationBuffer;
3929 m_LastError =
i18n(
"Unable to allocate memory for temporary bayer buffer: %1", e.what());
3933 m_ImageBufferSize = rgb_size;
3936 auto bayered_buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
3940 uint16_t * rBuff = bayered_buffer;
3941 uint16_t * gBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height);
3942 uint16_t * bBuff = bayered_buffer + (m_Statistics.width * m_Statistics.height * 2);
3944 int imax = m_Statistics.samples_per_channel * 3 - 3;
3945 for (
int i = 0; i <= imax; i += 3)
3947 *rBuff++ = bayer_destination_buffer[i];
3948 *gBuff++ = bayer_destination_buffer[i + 1];
3949 *bBuff++ = bayer_destination_buffer[i + 2];
3952 m_Statistics.channels = (m_Mode == FITS_NORMAL || m_Mode == FITS_CALIBRATE) ? 3 : 1;
3953 m_Statistics.dataType = TUSHORT;
3954 delete[] destinationBuffer;
3958void FITSData::logOOMError(uint32_t requiredMemory)
3960 qCCritical(KSTARS_FITS) <<
"Debayed memory allocation failure. Required Memory:" <<
KFormat().
formatByteSize(requiredMemory)
3961 <<
"Available system memory:" << KSUtils::getAvailableRAM();
3964double FITSData::getADU()
const
3967 for (
int i = 0; i < m_Statistics.channels; i++)
3968 adu += m_Statistics.mean[i];
3970 return (adu /
static_cast<double>(m_Statistics.channels));
3973QString FITSData::getLastError()
const
3978template <
typename T>
3979void FITSData::convertToQImage(
double dataMin,
double dataMax,
double scale,
double zero,
QImage &image)
3981 const auto * buffer =
reinterpret_cast<const T *
>(getImageBuffer());
3982 const T limit = std::numeric_limits<T>::max();
3983 T bMin = dataMin < 0 ? 0 : dataMin;
3984 T bMax = dataMax > limit ? limit : dataMax;
3985 uint16_t w = width();
3986 uint16_t h = height();
3987 uint32_t size = w * h;
3990 if (channels() == 1)
3993 for (
int j = 0; j < h; j++)
3995 unsigned char * scanLine = image.
scanLine(j);
3997 for (
int i = 0; i < w; i++)
3999 val = qBound(bMin, buffer[j * w + i], bMax);
4000 val = val * scale + zero;
4001 scanLine[i] = qBound<unsigned char>(0,
static_cast<uint8_t
>(val), 255);
4007 double rval = 0, gval = 0, bval = 0;
4010 for (
int j = 0; j < h; j++)
4012 auto * scanLine =
reinterpret_cast<QRgb *
>((image.
scanLine(j)));
4014 for (
int i = 0; i < w; i++)
4016 rval = qBound(bMin, buffer[j * w + i], bMax);
4017 gval = qBound(bMin, buffer[j * w + i + size], bMax);
4018 bval = qBound(bMin, buffer[j * w + i + size * 2], bMax);
4020 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
4022 scanLine[i] = value;
4039 if (future.
result() ==
false)
4042 data.getMinMax(&min, &max);
4050 if (data.channels() == 1)
4055 for (
int i = 0; i < 256; i++)
4056 fitsImage.
setColor(i, qRgb(i, i, i));
4063 double dataMin = data.m_Statistics.mean[0] - data.m_Statistics.stddev[0];
4064 double dataMax = data.m_Statistics.mean[0] + data.m_Statistics.stddev[0] * 3;
4066 double bscale = 255. / (dataMax - dataMin);
4067 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
4070 switch (data.m_Statistics.dataType)
4073 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4077 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4081 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4085 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4089 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4093 data.convertToQImage<
float>(dataMin, dataMax, bscale, bzero, fitsImage);
4097 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
4101 data.convertToQImage<
double>(dataMin, dataMax, bscale, bzero, fitsImage);
4115 qCCritical(KSTARS_FITS) <<
"Failed to convert" << filename <<
"to FITS since format" << format <<
"is not supported in Qt";
4123 qCCritical(KSTARS_FITS) <<
"Failed to open image" << filename;
4127 output =
QDir(getTemporaryPath()).
filePath(filename +
".fits");
4130 fitsfile *fptr =
nullptr;
4132 long fpixel = 1, naxis = input.
allGray() ? 2 : 3, nelements, exposure;
4133 long naxes[3] = { input.
width(), input.
height(), naxis == 3 ? 3 : 1 };
4134 char error_status[512] = {0};
4136 if (fits_create_file(&fptr,
QString(
'!' + output).toLocal8Bit().data(), &status))
4138 fits_get_errstatus(status, error_status);
4139 qCCritical(KSTARS_FITS) <<
"Failed to create FITS file. Error:" << error_status;
4143 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
4145 qCWarning(KSTARS_FITS) <<
"fits_create_img failed:" << error_status;
4147 fits_flush_file(fptr, &status);
4148 fits_close_file(fptr, &status);
4153 fits_update_key(fptr, TLONG,
"EXPOSURE", &exposure,
"Total Exposure Time", &status);
4158 nelements = naxes[0] * naxes[1];
4159 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.
bits(), &status))
4161 fits_get_errstatus(status, error_status);
4162 qCWarning(KSTARS_FITS) <<
"fits_write_img GRAY failed:" << error_status;
4164 fits_flush_file(fptr, &status);
4165 fits_close_file(fptr, &status);
4172 nelements = naxes[0] * naxes[1] * 3;
4174 uint8_t *srcBuffer = input.
bits();
4176 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
4178 uint8_t *rgbBuffer =
new uint8_t[nelements];
4179 if (rgbBuffer ==
nullptr)
4181 qCWarning(KSTARS_FITS) <<
"Not enough memory for RGB buffer";
4182 fits_flush_file(fptr, &status);
4183 fits_close_file(fptr, &status);
4187 uint8_t *subR = rgbBuffer;
4188 uint8_t *subG = rgbBuffer + naxes[0] * naxes[1];
4189 uint8_t *subB = rgbBuffer + naxes[0] * naxes[1] * 2;
4190 for (uint32_t i = 0; i < srcBytes; i += 4)
4192 *subB++ = srcBuffer[i];
4193 *subG++ = srcBuffer[i + 1];
4194 *subR++ = srcBuffer[i + 2];
4197 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
4199 fits_get_errstatus(status, error_status);
4200 qCWarning(KSTARS_FITS) <<
"fits_write_img RGB failed:" << error_status;
4202 fits_flush_file(fptr, &status);
4203 fits_close_file(fptr, &status);
4204 delete [] rgbBuffer;
4208 delete [] rgbBuffer;
4211 if (fits_flush_file(fptr, &status))
4213 fits_get_errstatus(status, error_status);
4214 qCWarning(KSTARS_FITS) <<
"fits_flush_file failed:" << error_status;
4216 fits_close_file(fptr, &status);
4220 if (fits_close_file(fptr, &status))
4222 fits_get_errstatus(status, error_status);
4223 qCWarning(KSTARS_FITS) <<
"fits_close_file failed:" << error_status;
4230void FITSData::injectWCS(
double orientation,
double ra,
double dec,
double pixscale,
bool eastToTheRight)
4234 if (fptr ==
nullptr)
4237 fits_update_key(fptr, TDOUBLE,
"OBJCTRA", &ra,
"Object RA", &status);
4238 fits_update_key(fptr, TDOUBLE,
"OBJCTDEC", &dec,
"Object DEC", &status);
4242 fits_update_key(fptr, TINT,
"EQUINOX", &epoch,
"Equinox", &status);
4244 fits_update_key(fptr, TDOUBLE,
"CRVAL1", &ra,
"CRVAL1", &status);
4245 fits_update_key(fptr, TDOUBLE,
"CRVAL2", &dec,
"CRVAL1", &status);
4247 char radecsys[8] =
"FK5";
4248 char ctype1[16] =
"RA---TAN";
4249 char ctype2[16] =
"DEC--TAN";
4251 fits_update_key(fptr, TSTRING,
"RADECSYS", radecsys,
"RADECSYS", &status);
4252 fits_update_key(fptr, TSTRING,
"CTYPE1", ctype1,
"CTYPE1", &status);
4253 fits_update_key(fptr, TSTRING,
"CTYPE2", ctype2,
"CTYPE2", &status);
4255 double crpix1 = width() / 2.0;
4256 double crpix2 = height() / 2.0;
4258 fits_update_key(fptr, TDOUBLE,
"CRPIX1", &crpix1,
"CRPIX1", &status);
4259 fits_update_key(fptr, TDOUBLE,
"CRPIX2", &crpix2,
"CRPIX2", &status);
4262 double secpix1 = eastToTheRight ? pixscale : -pixscale;
4263 double secpix2 = pixscale;
4265 fits_update_key(fptr, TDOUBLE,
"SECPIX1", &secpix1,
"SECPIX1", &status);
4266 fits_update_key(fptr, TDOUBLE,
"SECPIX2", &secpix2,
"SECPIX2", &status);
4268 double degpix1 = secpix1 / 3600.0;
4269 double degpix2 = secpix2 / 3600.0;
4271 fits_update_key(fptr, TDOUBLE,
"CDELT1", °pix1,
"CDELT1", &status);
4272 fits_update_key(fptr, TDOUBLE,
"CDELT2", °pix2,
"CDELT2", &status);
4275 double rotation = 360 - orientation;
4279 fits_update_key(fptr, TDOUBLE,
"CROTA1", &rotation,
"CROTA1", &status);
4280 fits_update_key(fptr, TDOUBLE,
"CROTA2", &rotation,
"CROTA2", &status);
4285bool FITSData::contains(
const QPointF &point)
const
4287 return (point.
x() >= 0 && point.
y() >= 0 && point.
x() <= m_Statistics.width && point.
y() <= m_Statistics.height);
4290void FITSData::saveStatistics(FITSImage::Statistic &other)
4292 other = m_Statistics;
4295void FITSData::restoreStatistics(FITSImage::Statistic &other)
4297 m_Statistics = other;
4302void FITSData::constructHistogram()
4304 switch (m_Statistics.dataType)
4307 constructHistogramInternal<uint8_t>();
4311 constructHistogramInternal<int16_t>();
4315 constructHistogramInternal<uint16_t>();
4319 constructHistogramInternal<int32_t>();
4323 constructHistogramInternal<uint32_t>();
4327 constructHistogramInternal<float>();
4331 constructHistogramInternal<int64_t>();
4335 constructHistogramInternal<double>();
4343template <
typename T> int32_t FITSData::histogramBinInternal(T value,
int channel)
const
4345 return qMax(
static_cast<T
>(0), qMin(
static_cast<T
>(m_HistogramBinCount),
4346 static_cast<T
>(rint((value - m_Statistics.min[channel]) / m_HistogramBinWidth[channel]))));
4349template <
typename T> int32_t FITSData::histogramBinInternal(
int x,
int y,
int channel)
const
4351 if (!m_ImageBuffer || !isHistogramConstructed())
4353 uint32_t samples = m_Statistics.width * m_Statistics.height;
4354 uint32_t offset = channel * samples;
4355 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
4356 int index = y * m_Statistics.width + x;
4357 const T &sample = buffer[index + offset];
4358 return histogramBinInternal(sample, channel);
4361int32_t FITSData::histogramBin(
int x,
int y,
int channel)
const
4363 switch (m_Statistics.dataType)
4366 return histogramBinInternal<uint8_t>(x, y, channel);
4370 return histogramBinInternal<int16_t>(x, y, channel);
4374 return histogramBinInternal<uint16_t>(x, y, channel);
4378 return histogramBinInternal<int32_t>(x, y, channel);
4382 return histogramBinInternal<uint32_t>(x, y, channel);
4386 return histogramBinInternal<float>(x, y, channel);
4390 return histogramBinInternal<int64_t>(x, y, channel);
4394 return histogramBinInternal<double>(x, y, channel);
4403template <
typename T>
void FITSData::constructHistogramInternal()
4405 auto *
const buffer =
reinterpret_cast<T
const *
>(m_ImageBuffer);
4406 uint32_t samples = m_Statistics.width * m_Statistics.height;
4407 const uint32_t sampleBy = samples > 500000 ? samples / 500000 : 1;
4409 m_HistogramBinCount = qMax(0., qMin(m_Statistics.max[0] - m_Statistics.min[0], 256.0));
4410 if (m_HistogramBinCount <= 0)
4411 m_HistogramBinCount = 256;
4413 for (
int n = 0; n < m_Statistics.channels; n++)
4415 m_HistogramIntensity[n].fill(0, m_HistogramBinCount + 1);
4416 m_HistogramFrequency[n].fill(0, m_HistogramBinCount + 1);
4417 m_CumulativeFrequency[n].fill(0, m_HistogramBinCount + 1);
4419 const double minBinSize = (m_Statistics.max[n] > 1.1) ? 1.0 : .0001;
4420 m_HistogramBinWidth[n] = qMax(minBinSize, (m_Statistics.max[n] - m_Statistics.min[n]) / m_HistogramBinCount);
4425 for (
int n = 0; n < m_Statistics.channels; n++)
4429 for (
int i = 0; i < m_HistogramBinCount; i++)
4430 m_HistogramIntensity[n][i] = m_Statistics.min[n] + (m_HistogramBinWidth[n] * i);
4434 for (
int n = 0; n < m_Statistics.channels; n++)
4438 uint32_t offset = n * samples;
4440 for (uint32_t i = 0; i < samples; i += sampleBy)
4442 int32_t
id = histogramBinInternal<T>(buffer[i + offset], n);
4443 m_HistogramFrequency[n][id] += sampleBy;
4449 future.waitForFinished();
4453 for (
int n = 0; n < m_Statistics.channels; n++)
4457 uint32_t accumulator = 0;
4458 for (
int i = 0; i < m_HistogramBinCount; i++)
4460 accumulator += m_HistogramFrequency[n][i];
4461 m_CumulativeFrequency[n].replace(i, accumulator);
4467 future.waitForFinished();
4472 if (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 4] > 0)
4473 m_JMIndex = m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount / 8] /
static_cast<double>
4474 (m_CumulativeFrequency[RED_CHANNEL][m_HistogramBinCount /
4479 qCDebug(KSTARS_FITS) <<
"FITHistogram: JMIndex " << m_JMIndex;
4481 m_HistogramConstructed =
true;
4482 emit histogramReady();
4485void FITSData::recordLastError(
int errorCode)
4487 char fitsErrorMessage[512] = {0};
4488 fits_get_errstatus(errorCode, fitsErrorMessage);
4489 m_LastError = fitsErrorMessage;
4492double FITSData::getAverageMean(
bool roi)
const
4494 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
4495 if (ptr->channels == 1)
4496 return ptr->mean[0];
4498 return (ptr->mean[0] + ptr->mean[1] + ptr->mean[2]) / 3.0;
4501double FITSData::getAverageMedian(
bool roi)
const
4503 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
4504 if (ptr->channels == 1)
4505 return ptr->median[0];
4507 return (ptr->median[0] + ptr->median[1] + ptr->median[2]) / 3.0;
4510double FITSData::getAverageStdDev(
bool roi)
const
4512 const FITSImage::Statistic* ptr = (roi ? &m_ROIStatistics : &m_Statistics);
4513 if (ptr->channels == 1)
4514 return ptr->stddev[0];
4516 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,...
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,...
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)
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)
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
KOSM_EXPORT double distance(const std::vector< const OSM::Node * > &path, Coordinate coord)
QByteArray & append(QByteArrayView data)
const char * constData() const const
bool isEmpty() const const
QByteArray leftJustified(qsizetype width, char fill, bool truncate) const const
void push_back(QByteArrayView str)
qsizetype size() const const
QDateTime currentDateTime()
QDateTime fromString(QStringView string, QStringView format, QCalendar cal)
bool isValid() const const
QString 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)
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)
qsizetype size() const const
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 mid(qsizetype position, qsizetype n) const const
QString & remove(QChar ch, Qt::CaseSensitivity cs)
QStringList split(QChar sep, Qt::SplitBehavior behavior, Qt::CaseSensitivity cs) const const
bool startsWith(QChar c, Qt::CaseSensitivity cs) const const
QByteArray toLatin1() const const
QByteArray toLocal8Bit() const const
bool contains(QLatin1StringView str, Qt::CaseSensitivity cs) 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)
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