• Skip to content
  • Skip to link menu
KDE API Reference
  • KDE API Reference
  • edu API Reference
  • KDE Home
  • Contact Us
 

kstars

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

KDE's Doxygen guidelines are available online.

kstars

Skip menu "kstars"
  • Main Page
  • Namespace List
  • Namespace Members
  • Alphabetical List
  • Class List
  • Class Hierarchy
  • Class Members
  • File List
  • File Members
  • Modules
  • Related Pages

edu API Reference

Skip menu "edu API Reference"
  •     core
  • kstars

Search



Report problems with this website to our bug tracking system.
Contact the specific authors with questions and comments about the page contents.

KDE® and the K Desktop Environment® logo are registered trademarks of KDE e.V. | Legal