Kstars

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

KDE's Doxygen guidelines are available online.