32 #include <QApplication>
34 #include <QtConcurrent>
35 #include <QImageReader>
37 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
49 #include <fits_debug.h>
51 #define ZOOM_DEFAULT 100.0
54 #define ZOOM_LOW_INCR 10
55 #define ZOOM_HIGH_INCR 50
58 const QString FITSData::m_TemporaryPath = QStandardPaths::writableLocation(QStandardPaths::TempLocation);
60 #define DIFFUSE_THRESHOLD 0.15
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
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));
100 if (starCenters.
count() > 0)
101 qDeleteAll(starCenters);
110 fits_close_file(fptr, &status);
113 if (m_isTemporary && autoRemoveTemporaryFITS)
120 void FITSData::loadCommon(
const QString &inFilename)
123 qDeleteAll(starCenters);
128 fits_close_file(fptr, &status);
136 if (m_isTemporary && autoRemoveTemporaryFITS && inFilename != m_Filename)
140 m_Filename = inFilename;
144 size_t fits_buffer_size,
bool silent)
146 loadCommon(inFilename);
147 qCInfo(KSTARS_FITS) <<
"Reading FITS file buffer ";
148 return privateLoad(fits_buffer, fits_buffer_size, silent);
153 loadCommon(inFilename);
154 qCInfo(KSTARS_FITS) <<
"Loading FITS file " << m_Filename;
156 this, &FITSData::privateLoad,
nullptr, 0, silent);
164 bool fitsOpenError(
int status,
const QString &message,
bool silent)
166 char error_status[512];
167 fits_report_error(stderr, status);
168 fits_get_errstatus(status, error_status);
173 qCCritical(KSTARS_FITS) << errMessage;
178 bool FITSData::privateLoad(
void *fits_buffer,
size_t fits_buffer_size,
bool silent)
180 int status = 0, anynull = 0;
184 m_isTemporary = m_Filename.
startsWith(m_TemporaryPath);
186 if (fits_buffer ==
nullptr && m_Filename.
endsWith(
".fz"))
189 m_compressedFilename = m_Filename;
196 errMessage = i18n(
"Failed to unpack compressed fits");
199 qCCritical(KSTARS_FITS) << errMessage;
204 if (m_isTemporary && autoRemoveTemporaryFITS)
207 m_Filename = uncompressedFile;
208 m_isTemporary =
true;
209 m_isCompressed =
true;
212 if (fits_buffer ==
nullptr)
216 if (fits_open_diskfile(&fptr, m_Filename.
toLatin1(), READONLY, &status))
217 return fitsOpenError(status, i18n(
"Error opening fits file %1", m_Filename), silent);
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);
230 stats.
size = fits_buffer_size;
233 if (fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status))
234 return fitsOpenError(status, i18n(
"Could not locate image HDU."), silent);
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);
241 errMessage = i18n(
"1D FITS images are not supported in KStars.");
244 qCCritical(KSTARS_FITS) << errMessage;
256 m_DataType = TUSHORT;
260 m_DataType = TUSHORT;
277 m_DataType = TLONGLONG;
281 m_DataType = TDOUBLE;
285 errMessage = i18n(
"Bit depth %1 is not supported.", stats.
bitpix);
288 qCCritical(KSTARS_FITS) << errMessage;
295 if (naxes[0] == 0 || naxes[1] == 0)
297 errMessage = i18n(
"Image has invalid dimensions %1x%2", naxes[0], naxes[1]);
300 qCCritical(KSTARS_FITS) << errMessage;
304 stats.
width = naxes[0];
310 m_Channels = naxes[2];
318 if (m_ImageBuffer ==
nullptr)
320 qCWarning(KSTARS_FITS) <<
"FITSData: Not enough memory for image_buffer channel. Requested: "
331 if (fits_read_img(fptr, m_DataType, 1, nelements,
nullptr, m_ImageBuffer, &anynull, &status))
332 return fitsOpenError(status, i18n(
"Error reading image."), silent);
338 bayerBuffer = m_ImageBuffer;
350 starsSearched =
false;
357 if (newFilename == m_Filename)
366 int status = 0, exttype = 0;
373 if (fits_close_file(fptr, &status))
375 fits_report_error(stderr, status);
380 QString finalFileName(newFilename);
382 finalFileName.
remove(
'!');
389 qCCritical(KSTARS_FITS()) <<
"FITS: Failed to copy " << m_Filename <<
" to " << finalFileName;
394 if (m_isTemporary && autoRemoveTemporaryFITS)
397 m_isTemporary =
false;
400 m_Filename = finalFileName;
404 fits_open_diskfile(&fptr, m_Filename.
toLatin1(), READONLY, &status);
405 fits_movabs_hdu(fptr, 1, IMAGE_HDU, &status);
413 if (fits_create_file(&new_fptr, newFilename.
toLatin1(), &status))
415 fits_report_error(stderr, status);
425 if (fits_copy_header(fptr, new_fptr, &status))
427 fits_report_error(stderr, status);
432 if (fits_close_file(fptr, &status))
434 fits_report_error(stderr, status);
442 if (fits_movabs_hdu(fptr, 1, &exttype, &status))
444 fits_report_error(stderr, status);
449 if (fits_write_img(fptr, m_DataType, 1, nelements, m_ImageBuffer, &status))
451 fits_report_error(stderr, status);
458 if (fits_update_key(fptr, TDOUBLE,
"DATAMIN", &(stats.
min),
"Minimum value", &status))
460 fits_report_error(stderr, status);
465 if (fits_update_key(fptr, TDOUBLE,
"DATAMAX", &(stats.
max),
"Maximum value", &status))
467 fits_report_error(stderr, status);
472 if (fits_update_key(fptr, TUSHORT,
"NAXIS1", &(stats.
width),
"length of data axis 1", &status))
474 fits_report_error(stderr, status);
479 if (fits_update_key(fptr, TUSHORT,
"NAXIS2", &(stats.
height),
"length of data axis 2", &status))
481 fits_report_error(stderr, status);
486 if (fits_write_date(fptr, &status))
488 fits_report_error(stderr, status);
495 if (fits_write_history(fptr, history.
toLatin1(), &status))
497 fits_report_error(stderr, status);
501 int rot = 0, mirror = 0;
504 rot = (90 * rotCounter) % 360;
508 if (flipHCounter % 2 != 0 || flipVCounter % 2 != 0)
511 if ((rot != 0) || (mirror != 0))
512 rotWCSFITS(rot, mirror);
514 rotCounter = flipHCounter = flipVCounter = 0;
516 if (m_isTemporary && autoRemoveTemporaryFITS)
519 m_isTemporary =
false;
522 m_Filename = newFilename;
524 fits_flush_file(fptr, &status);
526 qCInfo(KSTARS_FITS) <<
"Saved FITS file:" << m_Filename;
533 delete[] m_ImageBuffer;
534 m_ImageBuffer =
nullptr;
535 bayerBuffer =
nullptr;
541 calculateMinMax(refresh);
547 runningAverageStdDev<uint8_t>();
551 runningAverageStdDev<int16_t>();
555 runningAverageStdDev<uint16_t>();
559 runningAverageStdDev<int32_t>();
563 runningAverageStdDev<uint32_t>();
567 runningAverageStdDev<float>();
571 runningAverageStdDev<int64_t>();
575 runningAverageStdDev<double>();
585 if (refresh && markStars)
587 starsSearched =
false;
590 int FITSData::calculateMinMax(
bool refresh)
592 int status, nfound = 0;
596 if ((fptr !=
nullptr) && !refresh)
598 if (fits_read_key_dbl(fptr,
"DATAMIN", &(stats.
min[0]),
nullptr, &status) == 0)
601 if (fits_read_key_dbl(fptr,
"DATAMAX", &(stats.
max[0]),
nullptr, &status) == 0)
605 if (nfound == 2 && !(stats.
min[0] == 0 && stats.
max[0] == 0))
609 stats.
min[0] = 1.0E30;
610 stats.
max[0] = -1.0E30;
612 stats.
min[1] = 1.0E30;
613 stats.
max[1] = -1.0E30;
615 stats.
min[2] = 1.0E30;
616 stats.
max[2] = -1.0E30;
621 calculateMinMax<uint8_t>();
625 calculateMinMax<int16_t>();
629 calculateMinMax<uint16_t>();
633 calculateMinMax<int32_t>();
637 calculateMinMax<uint32_t>();
641 calculateMinMax<float>();
645 calculateMinMax<int64_t>();
649 calculateMinMax<double>();
660 template <
typename T>
661 QPair<T, T> FITSData::getParitionMinMax(uint32_t start, uint32_t stride)
663 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
664 T min = std::numeric_limits<T>::max();
665 T max = std::numeric_limits<T>::min();
667 uint32_t end = start + stride;
669 for (uint32_t i = start; i < end; i++)
673 else if (buffer[i] > max)
677 return qMakePair(min, max);
680 template <
typename T>
681 void FITSData::calculateMinMax()
683 T min = std::numeric_limits<T>::max();
684 T max = std::numeric_limits<T>::min();
687 const uint8_t nThreads = 16;
689 for (
int n = 0; n < m_Channels; n++)
700 uint32_t tStart = cStart;
705 for (
int i = 0; i < nThreads; i++)
708 futures.
append(
QtConcurrent::run(
this, &FITSData::getParitionMinMax<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride));
713 for (
int i = 0; i < nThreads; i++)
716 if (result.first < min)
718 if (result.second > max)
727 template <
typename T>
731 double m_oldM = 0, m_newM = 0, m_oldS = 0, m_newS = 0;
733 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
734 uint32_t end = start + stride;
736 for (uint32_t i = start; i < end; i++)
738 m_newM = m_oldM + (buffer[i] - m_oldM) / m_n;
739 m_newS = m_oldS + (buffer[i] - m_oldM) * (buffer[i] - m_newM);
746 return qMakePair<double, double>(m_newM, m_newS);
749 template <
typename T>
750 void FITSData::runningAverageStdDev()
753 const uint8_t nThreads = 16;
755 for (
int n = 0; n < m_Channels; n++)
766 uint32_t tStart = cStart;
771 for (
int i = 0; i < nThreads; i++)
774 futures.
append(
QtConcurrent::run(
this, &FITSData::getSquaredSumAndMean<T>, tStart, (i == (nThreads - 1)) ? fStride : tStride));
778 double mean = 0, squared_sum = 0;
781 for (
int i = 0; i < nThreads; i++)
784 mean += result.first;
785 squared_sum += result.second;
790 stats.
mean[n] = mean / nThreads;
791 stats.
stddev[n] = sqrt(variance);
797 stats.
min[channel] = newMin;
798 stats.
max[channel] = newMax;
801 bool FITSData::parseHeader()
803 char * header =
nullptr;
804 int status = 0, nkeys = 0;
806 if (fits_hdr2str(fptr, 0,
nullptr, 0, &header, &nkeys, &status))
808 fits_report_error(stderr, status);
815 for (
int i = 0; i < nkeys; i++)
817 Record * oneRecord =
new Record;
822 if (properties.
size() == 1)
824 oneRecord->key = properties[0].
mid(0, 7);
825 oneRecord->comment = properties[0].
mid(8).simplified();
829 oneRecord->key = properties[0].simplified();
830 oneRecord->
value = properties[1].simplified();
831 if (properties.
size() > 2)
832 oneRecord->comment = properties[2].simplified();
839 oneRecord->
value.toInt(&ok);
841 oneRecord->value.convert(QMetaType::Int);
845 oneRecord->value.toDouble(&ok);
847 oneRecord->value.convert(QMetaType::Double);
851 records.append(oneRecord);
861 for (
Record * oneRecord : records)
863 if (oneRecord->key == key)
865 value = oneRecord->
value;
873 bool FITSData::checkCollision(
Edge * s1,
Edge * s2)
877 int diff_x = s1->
x - s2->
x;
878 int diff_y = s1->
y - s2->
y;
880 dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y));
881 dis -= s1->
width / 2;
882 dis -= s2->
width / 2;
896 return FITSData::findCannyStar<uint8_t>(data, boundary);
899 return FITSData::findCannyStar<int16_t>(data, boundary);
902 return FITSData::findCannyStar<uint16_t>(data, boundary);
905 return FITSData::findCannyStar<int32_t>(data, boundary);
908 return FITSData::findCannyStar<uint16_t>(data, boundary);
911 return FITSData::findCannyStar<float>(data, boundary);
914 return FITSData::findCannyStar<int64_t>(data, boundary);
917 return FITSData::findCannyStar<double>(data, boundary);
929 starAlgorithm = algorithm;
931 qDeleteAll(starCenters);
945 count = findCentroid(trackingBox);
953 starsSearched =
true;
961 long const sqInnerRadius = std::lround(sqDiagonal * innerRadius * innerRadius);
962 long const sqOuterRadius = std::lround(sqDiagonal * outerRadius * outerRadius);
964 starCenters.
erase(std::remove_if(starCenters.
begin(), starCenters.
end(),
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());
973 return starCenters.
count();
976 template <
typename T>
979 int subX = qMax(0, boundary.
isNull() ? 0 : boundary.
x());
980 int subY = qMax(0, boundary.
isNull() ? 0 : boundary.
y());
986 uint16_t dataWidth = data->
width();
989 uint32_t
size = subW * subH;
990 uint32_t offset = subX + subY * dataWidth;
993 auto * buffer =
new uint8_t[size * BBP];
999 uint8_t * dataPtr = buffer;
1001 uint32_t lineOffset = 0;
1005 lineOffset = (subX +
height * dataWidth) * BBP;
1006 memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
1007 dataPtr += (subW * BBP);
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;
1020 boundedImage->setProperty(
"dataType", data->
property(
"dataType"));
1023 boundedImage->setImageBuffer(buffer);
1025 boundedImage->calculateStats(
true);
1037 boundedImage->sobel<T>(gradients, directions);
1041 int maxID = boundedImage->partition(subW, subH, gradients, ids);
1044 delete boundedImage;
1053 float totalMass = 0;
1059 for (
int y = 0; y < subH; y++)
1061 for (
int x = 0; x < subW; x++)
1063 int index = x + y * subW;
1065 int regionID = ids[index];
1068 float pixel = gradients[index];
1070 masses[regionID].totalMass += pixel;
1071 masses[regionID].massX += x * pixel;
1072 masses[regionID].massY += y * pixel;
1078 int maxRegionID = 1;
1079 int maxTotalMass = masses[1].totalMass;
1080 double totalMassRatio = 1e6;
1081 for (
auto key : masses.
keys())
1083 massInfo oneMass = masses.
value(key);
1084 if (oneMass.totalMass > maxTotalMass)
1086 totalMassRatio = oneMass.totalMass / maxTotalMass;
1087 maxTotalMass = oneMass.totalMass;
1094 if (maxID > 10 && totalMassRatio < 1.5)
1097 auto * center =
new Edge;
1099 center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
1100 center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
1104 int maxR = qMin(subW - 1, subH - 1) / 2;
1106 for (
int r = maxR; r > 1; r--)
1110 for (
float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
1112 int testX = center->x + std::cos(theta) * r;
1113 int testY = center->y + std::sin(theta) * r;
1116 if (testX < 0 || testX >= subW || testY < 0 || testY >= subH)
1119 if (gradients[testX + testY * subW] > 0)
1124 center->width = r * 2;
1133 qCDebug(KSTARS_FITS) <<
"FITS: Weighted Center is X: " << center->x <<
" Y: " << center->y <<
" Width: " << center->width;
1136 if (center->width == -1)
1145 double FSum = 0, HF = 0, TF = 0;
1146 const double resolution = 1.0 / 20.0;
1148 int cen_y = qRound(center->y);
1150 double rightEdge = center->x + center->width / 2.0;
1151 double leftEdge = center->x - center->width / 2.0;
1154 subPixels.
reserve(center->width / resolution);
1156 const T * origBuffer =
reinterpret_cast<T *
>(data->
getImageBuffer()) + offset;
1158 for (
double x = leftEdge; x <= rightEdge; x += resolution)
1160 double slice = resolution * (origBuffer[
static_cast<int>(floor(x)) + cen_y * dataWidth]);
1168 int subPixelCenter = (center->width / resolution) / 2;
1171 TF = subPixels[subPixelCenter];
1175 for (
int k = 1; k < subPixelCenter; k++)
1177 TF += subPixels[subPixelCenter + k];
1178 TF += subPixels[subPixelCenter - k];
1185 center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
1201 qCDebug(KSTARS_FITS) <<
"Flux: " << FSum <<
" Half-Flux: " << HF <<
" HFR: " << center->HFR;
1211 return findOneStar<uint8_t>(boundary);
1214 return findOneStar<int16_t>(boundary);
1217 return findOneStar<uint16_t>(boundary);
1220 return findOneStar<int32_t>(boundary);
1223 return findOneStar<uint32_t>(boundary);
1226 return findOneStar<float>(boundary);
1229 return findOneStar<int64_t>(boundary);
1232 return findOneStar<double>(boundary);
1241 template <
typename T>
1247 int subX = boundary.
x();
1248 int subY = boundary.
y();
1249 int subW = subX + boundary.
width();
1250 int subH = subY + boundary.
height();
1252 float massX = 0, massY = 0, totalMass = 0;
1254 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
1259 for (
int y = subY; y < subH; y++)
1261 for (
int x = subX; x < subW; x++)
1263 T pixel = buffer[x + y * stats.
width];
1264 if (pixel > threshold)
1273 qCDebug(KSTARS_FITS) <<
"FITS: Weighted Center is X: " << massX / totalMass <<
" Y: " << massY / totalMass;
1275 auto * center =
new Edge;
1277 center->x = massX / totalMass + 0.5;
1278 center->y = massY / totalMass + 0.5;
1282 int maxR = qMin(subW - 1, subH - 1) / 2;
1285 double critical_threshold = threshold * 0.7;
1286 double running_threshold = threshold;
1288 while (running_threshold >= critical_threshold)
1290 for (
int r = maxR; r > 1; r--)
1294 for (
float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0)
1296 int testX = center->x + std::cos(theta) * r;
1297 int testY = center->y + std::sin(theta) * r;
1300 if (testX < subX || testX > subW || testY < subY || testY > subH)
1303 if (buffer[testX + testY * stats.
width] > running_threshold)
1311 center->width = r * 2;
1316 if (center->width > 0)
1320 running_threshold -= running_threshold * 0.1;
1324 if (center->width == -1)
1333 starCenters.
append(center);
1335 double FSum = 0, HF = 0, TF = 0, min = stats.
min[0];
1336 const double resolution = 1.0 / 20.0;
1338 int cen_y = qRound(center->y);
1340 double rightEdge = center->x + center->width / 2.0;
1341 double leftEdge = center->x - center->width / 2.0;
1344 subPixels.
reserve(center->width / resolution);
1346 for (
double x = leftEdge; x <= rightEdge; x += resolution)
1349 double slice = resolution * (buffer[
static_cast<int>(floor(x)) + cen_y * stats.
width] - min);
1358 int subPixelCenter = (center->width / resolution) / 2;
1361 TF = subPixels[subPixelCenter];
1365 for (
int k = 1; k < subPixelCenter; k++)
1367 TF += subPixels[subPixelCenter + k];
1368 TF += subPixels[subPixelCenter - k];
1376 center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
1390 int FITSData::findCentroid(
const QRect &boundary,
int initStdDev,
int minEdgeWidth)
1395 return findCentroid<uint8_t>(boundary, initStdDev, minEdgeWidth);
1398 return findCentroid<int16_t>(boundary, initStdDev, minEdgeWidth);
1401 return findCentroid<uint16_t>(boundary, initStdDev, minEdgeWidth);
1404 return findCentroid<int32_t>(boundary, initStdDev, minEdgeWidth);
1407 return findCentroid<uint32_t>(boundary, initStdDev, minEdgeWidth);
1410 return findCentroid<float>(boundary, initStdDev, minEdgeWidth);
1413 return findCentroid<int64_t>(boundary, initStdDev, minEdgeWidth);
1416 return findCentroid<double>(boundary, initStdDev, minEdgeWidth);
1423 template <
typename T>
1424 int FITSData::findCentroid(
const QRect &boundary,
int initStdDev,
int minEdgeWidth)
1426 double threshold = 0, sum = 0, avg = 0, min = 0;
1427 int starDiameter = 0;
1431 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
1433 double JMIndex = 100;
1443 float dispersion_ratio = 1.5;
1449 minEdgeWidth = JMIndex * 35 + 1;
1450 minimumEdgeCount = minEdgeWidth - 1;
1455 minimumEdgeCount = 4;
1458 while (initStdDev >= 1)
1463 minEdgeWidth = qMax(3, minEdgeWidth);
1464 minimumEdgeCount = qMax(3, minimumEdgeCount);
1472 if (threshold - min < 0)
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;
1494 int subX, subY, subW, subH;
1501 subX = round(stats.
width * 0.15);
1502 subY = round(stats.
height * 0.15);
1503 subW = stats.
width - subX;
1504 subH = stats.
height - subY;
1517 subX = boundary.
x();
1518 subY = boundary.
y();
1519 subW = subX + boundary.
width();
1520 subH = subY + boundary.
height();
1524 for (
int i = subY; i < subH; i++)
1528 for (
int j = subX; j < subW; j++)
1530 pixVal = buffer[j + (i * stats.
width)] - min;
1533 if (pixVal >= threshold)
1543 if (starDiameter >= minEdgeWidth)
1545 float center = avg / sum + 0.5;
1548 int i_center = std::floor(center);
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) >=
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
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;
1566 auto * newEdge =
new Edge();
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;
1582 avg = sum = starDiameter = 0;
1587 qCDebug(KSTARS_FITS) <<
"Total number of edges found is: " << edges.
count();
1590 if (edges.
count() == 1 && initStdDev > 1)
1598 qCWarning(KSTARS_FITS) <<
"Too many edges, aborting... " << edges.
count();
1603 if (edges.
count() >= minimumEdgeCount)
1622 for (
int i = 0; i < edges.
count(); i++)
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;
1628 if (edges[i]->scanned == 1)
1630 qCDebug(KSTARS_FITS) <<
"Skipping check for center " << i <<
" because it was already counted";
1634 qCDebug(KSTARS_FITS) <<
"Investigating edge # " << i <<
" now ...";
1637 cen_x = edges[i]->x;
1638 cen_y = edges[i]->y;
1639 cen_v = edges[i]->sum;
1640 cen_w = edges[i]->width;
1649 for (
int j = 0; j < edges.
count(); j++)
1651 if (edges[j]->scanned)
1654 if (checkCollision(edges[j], edges[i]))
1656 if (edges[j]->sum >= cen_v)
1658 cen_v = edges[j]->sum;
1659 cen_w = edges[j]->width;
1662 edges[j]->scanned = 1;
1665 avg_x += edges[j]->x * edges[j]->val;
1666 avg_y += edges[j]->y * edges[j]->val;
1667 sum += edges[j]->val;
1683 qCDebug(KSTARS_FITS) <<
"center_count: " << cen_count <<
" and initstdDev= " << initStdDev <<
" and limit is "
1691 if (cen_count >= cen_limit)
1694 auto * rCenter =
new Edge();
1696 rCenter->x = avg_x / sum;
1697 rCenter->y = avg_y / sum;
1698 width_sum += rCenter->width;
1699 rCenter->width = cen_w;
1701 qCDebug(KSTARS_FITS) <<
"Found a real center with number with (" << rCenter->x <<
"," << rCenter->y <<
")";
1708 cen_x = (int)std::floor(rCenter->x);
1709 cen_y = (int)std::floor(rCenter->y);
1711 if (cen_x < 0 || cen_x > stats.
width || cen_y < 0 || cen_y > stats.
height)
1719 for (
int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--)
1721 FSum += buffer[cen_x - k + (cen_y * stats.
width)] - min;
1729 TF = buffer[cen_y * stats.
width + cen_x] - min;
1731 int pixelCounter = 1;
1734 for (
int k = 1; k < rCenter->width / 2; k++)
1738 qCDebug(KSTARS_FITS) <<
"Stopping at TF " << TF <<
" after #" << k <<
" pixels.";
1742 TF += buffer[cen_y * stats.
width + cen_x + k] - min;
1743 TF += buffer[cen_y * stats.
width + cen_x - k] - min;
1749 rCenter->HFR = pixelCounter * (HF / TF);
1751 rCenter->val = FSum;
1753 qCDebug(KSTARS_FITS) <<
"HFR for this center is " << rCenter->HFR <<
" pixels and the total flux is " << FSum;
1755 starCenters.
append(rCenter);
1761 float width_avg = (float)width_sum / starCenters.
count();
1762 float lsum = 0, sdev = 0;
1764 for (
auto ¢er : starCenters)
1765 lsum += (center->width - width_avg) * (center->width - width_avg);
1767 sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4;
1770 foreach (
Edge * center, starCenters)
1771 if (center->
width > sdev)
1772 starCenters.removeOne(center);
1781 return starCenters.count();
1791 if (starCenters.empty())
1796 maxHFRStar =
nullptr;
1799 for (
int i = 0; i < starCenters.count(); i++)
1801 if (starCenters[i]->HFR > maxVal)
1804 maxVal = starCenters[i]->
HFR;
1808 maxHFRStar = starCenters[maxIndex];
1809 return static_cast<double>(starCenters[maxIndex]->HFR);
1813 for (
auto center : starCenters)
1814 HFRs << center->
HFR;
1815 std::sort(HFRs.
begin(), HFRs.
end());
1817 double sum = std::accumulate(HFRs.
begin(), HFRs.
end(), 0.0);
1818 double m = sum / HFRs.
size();
1820 if (HFRs.
size() > 3)
1823 std::for_each (HFRs.
begin(), HFRs.
end(), [&](
const double d)
1825 accum += (
d - m) * (
d - m);
1827 double stddev = sqrt(accum / (HFRs.
size() - 1));
1830 auto end1 = std::remove_if(HFRs.
begin(), HFRs.
end(), [m, stddev](
const double hfr)
1832 return hfr > (m + stddev * 2);
1834 auto end2 = std::remove_if(HFRs.
begin(), end1, [m, stddev](
const double hfr)
1836 return hfr < (m - stddev * 2);
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;
1850 if (starCenters.empty())
1853 for (
int i = 0; i < starCenters.count(); i++)
1855 if (std::fabs(starCenters[i]->x - x) <= starCenters[i]->width / 2 &&
1856 std::fabs(starCenters[i]->y - y) <= starCenters[i]->width / 2)
1858 return starCenters[i]->HFR;
1882 for (
int i = 0; i < 3; i++)
1884 dataMin[i] = stats.
mean[i] - stats.
stddev[i];
1885 dataMax[i] = stats.
mean[i] + stats.
stddev[i] * 3;
1892 for (
int i = 0; i < 3; i++)
1894 dataMin[i] = stats.
mean[i] + stats.
stddev[i];
1895 dataMax[i] = stats.
mean[i] + stats.
stddev[i] * 3;
1902 for (
int i = 0; i < 3; i++)
1904 dataMin[i] = stats.
mean[i];
1917 for (
int i = 0; i < 3; i++)
1919 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1920 dataMax[i] = dataMax[i] > UINT8_MAX ? UINT8_MAX : dataMax[i];
1922 applyFilter<uint8_t>(type, image, &dataMin, &dataMax);
1928 for (
int i = 0; i < 3; i++)
1930 dataMin[i] = dataMin[i] < INT16_MIN ? INT16_MIN : dataMin[i];
1931 dataMax[i] = dataMax[i] > INT16_MAX ? INT16_MAX : dataMax[i];
1933 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1940 for (
int i = 0; i < 3; i++)
1942 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1943 dataMax[i] = dataMax[i] > UINT16_MAX ? UINT16_MAX : dataMax[i];
1945 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1951 for (
int i = 0; i < 3; i++)
1953 dataMin[i] = dataMin[i] < INT_MIN ? INT_MIN : dataMin[i];
1954 dataMax[i] = dataMax[i] > INT_MAX ? INT_MAX : dataMax[i];
1956 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1962 for (
int i = 0; i < 3; i++)
1964 dataMin[i] = dataMin[i] < 0 ? 0 : dataMin[i];
1965 dataMax[i] = dataMax[i] > UINT_MAX ? UINT_MAX : dataMax[i];
1967 applyFilter<uint16_t>(type, image, &dataMin, &dataMax);
1973 for (
int i = 0; i < 3; i++)
1975 dataMin[i] = dataMin[i] < FLT_MIN ? FLT_MIN : dataMin[i];
1976 dataMax[i] = dataMax[i] > FLT_MAX ? FLT_MAX : dataMax[i];
1978 applyFilter<float>(type, image, &dataMin, &dataMax);
1984 for (
int i = 0; i < 3; i++)
1986 dataMin[i] = dataMin[i] < LLONG_MIN ? LLONG_MIN : dataMin[i];
1987 dataMax[i] = dataMax[i] > LLONG_MAX ? LLONG_MAX : dataMax[i];
1990 applyFilter<long>(type, image, &dataMin, &dataMax);
1996 for (
int i = 0; i < 3; i++)
1998 dataMin[i] = dataMin[i] < DBL_MIN ? DBL_MIN : dataMin[i];
1999 dataMax[i] = dataMax[i] > DBL_MAX ? DBL_MAX : dataMax[i];
2001 applyFilter<double>(type, image, &dataMin, &dataMax);
2016 template <
typename T>
2019 bool calcStats =
false;
2020 T * image =
nullptr;
2023 image =
reinterpret_cast<T *
>(targetImage);
2026 image =
reinterpret_cast<T *
>(m_ImageBuffer);
2031 for (
int i = 0; i < 3; i++)
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];
2039 const uint8_t nThreads = 16;
2062 for (
int i = 0; i < 3; i++)
2063 coeff[i] = max[i] / std::log(1 + max[i]);
2067 for (
int i = 0; i < 3; i++)
2068 coeff[i] = max[i] / sqrt(max[i]);
2071 for (
int n = 0; n < m_Channels; n++)
2074 min[n] = stats.
mean[n];
2084 T * runningBuffer = image + cStart;
2088 for (
int i = 0; i < nThreads; i++)
2091 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, coeff, n](T & a)
2093 a = qBound(min[n], static_cast<T>(round(coeff[n] * std::log(1 + qBound(min[n], a, max[n])))), max[n]);
2096 runningBuffer += tStride;
2101 for (
int i = 0; i < nThreads; i++)
2104 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, coeff, n](T & a)
2106 a = qBound(min[n], static_cast<T>(round(coeff[n] * a)), max[n]);
2110 runningBuffer += tStride;
2114 for (
int i = 0; i < nThreads; i++)
2117 futures.
append(
QtConcurrent::map(runningBuffer, (runningBuffer + ((i == (nThreads - 1)) ? fStride : tStride)), [min, max, n](T & a)
2119 a = qBound(min[n], a, max[n]);
2121 runningBuffer += tStride;
2126 for (
int i = 0; i < nThreads * m_Channels; i++)
2127 futures[i].waitForFinished();
2131 for (
int i = 0; i < 3; i++)
2133 stats.
min[i] = min[i];
2134 stats.
max[i] = max[i];
2137 runningAverageStdDev<T>();
2146 if (histogram ==
nullptr)
2155 double coeff = 255.0 / (height *
width);
2159 for (
int i = 0; i < m_Channels; i++)
2162 for (uint32_t j = 0; j <
height; j++)
2164 row = offset + j *
width;
2165 for (uint32_t k = 0; k <
width; k++)
2168 bufferVal = (image[index] - min[i]) / histogram->
getBinWidth(i);
2170 if (bufferVal >= cumulativeFreq.
size())
2171 bufferVal = cumulativeFreq.
size() - 1;
2173 image[index] = qBound(min[i], static_cast<T>(round(coeff * cumulativeFreq[bufferVal])), max[i]);
2187 auto * extension =
new T[(width + 2) * (height + 2)];
2192 for (uint32_t ch = 0; ch < m_Channels; ch++)
2197 for (uint32_t i = 0; i < M; ++i)
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];
2204 memcpy(extension, extension + N + 2, (N + 2) * BBP);
2206 memcpy(extension + (N + 2) * (M + 1), extension + (N + 2) * M, (N + 2) * BBP);
2213 for (uint32_t m = 1; m < M - 1; ++m)
2214 for (uint32_t n = 1; n < N - 1; ++n)
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];
2225 for (uint32_t j = 0; j < 5; ++j)
2229 for (uint32_t l = j + 1; l < 9; ++l)
2230 if (window[l] < window[mine])
2233 const float temp = window[j];
2234 window[j] = window[mine];
2235 window[mine] = temp;
2238 image[(m - 1) * (N - 2) + n - 1 + offset] = window[4];
2246 runningAverageStdDev<T>();
2278 for (
int i = 0; i < starCenters.count(); i++)
2280 int x =
static_cast<int>(starCenters[i]->x);
2281 int y =
static_cast<int>(starCenters[i]->y);
2284 starCentersInSubFrame.
append(starCenters[i]);
2287 return starCentersInSubFrame;
2292 if (starCenters.count() == 0)
2295 auto * pEdge =
new Edge();
2300 foreach (
Edge * center, starCenters)
2301 if (checkCollision(pEdge, center))
2303 *x =
static_cast<int>(center->
x);
2304 *y =
static_cast<int>(center->
y);
2318 int nkeyrec, nreject, nwcs;
2320 if (fits_hdr2str(fptr, 1,
nullptr, 0, &header, &nkeyrec, &status))
2323 fits_get_errstatus(status, errmsg);
2328 if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &wcs)) != 0)
2331 lastError =
QString(
"wcspih ERROR %1: %2.").
arg(status).
arg(wcshdr_errmsg[status]);
2340 lastError = i18n(
"No world coordinate systems found.");
2345 if (wcs->crpix[0] == 0)
2347 lastError = i18n(
"No world coordinate systems found.");
2351 if ((status = wcsset(wcs)) != 0)
2353 lastError =
QString(
"wcsset error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2365 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2369 qWarning() <<
"WCS data already loaded";
2373 qCDebug(KSTARS_FITS) <<
"Started WCS Data Processing...";
2377 int nkeyrec, nreject, nwcs, stat[2];
2378 double imgcrd[2], phi = 0, pixcrd[2], theta = 0, world[2];
2382 if (fits_hdr2str(fptr, 1,
nullptr, 0, &header, &nkeyrec, &status))
2385 fits_get_errstatus(status, errmsg);
2390 if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs, &wcs)) != 0)
2393 lastError =
QString(
"wcspih ERROR %1: %2.").
arg(status).
arg(wcshdr_errmsg[status]);
2402 lastError = i18n(
"No world coordinate systems found.");
2407 if (wcs->crpix[0] == 0)
2409 lastError = i18n(
"No world coordinate systems found.");
2413 if ((status = wcsset(wcs)) != 0)
2415 lastError =
QString(
"wcsset error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2423 if (wcs_coord ==
nullptr)
2425 lastError =
"Not enough memory for WCS data!";
2431 for (
int i = 0; i < h; i++)
2433 for (
int j = 0; j < w; j++)
2438 if ((status = wcsp2s(wcs, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0])) != 0)
2440 lastError =
QString(
"wcsp2s error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2452 findObjectsInImage(&world[0], phi, theta, &imgcrd[0], &pixcrd[0], &stat[0]);
2457 qCDebug(KSTARS_FITS) <<
"Finished WCS Data processing...";
2467 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2470 double imgcrd[2], worldcrd[2], pixcrd[2], phi[2], theta[2];
2474 lastError = i18n(
"No world coordinate systems found.");
2481 if ((status = wcss2p(wcs, 1, 2, &worldcrd[0], &phi[0], &theta[0], &imgcrd[0], &pixcrd[0], &stat[0])) != 0)
2483 lastError =
QString(
"wcss2p error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2487 wcsImagePoint.
setX(imgcrd[0]);
2488 wcsImagePoint.
setY(imgcrd[1]);
2490 wcsPixelPoint.
setX(pixcrd[0]);
2491 wcsPixelPoint.
setY(pixcrd[1]);
2496 Q_UNUSED(wcsPixelPoint);
2497 Q_UNUSED(wcsImagePoint);
2504 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2507 double imgcrd[2], phi, pixcrd[2], theta, world[2];
2511 lastError = i18n(
"No world coordinate systems found.");
2515 pixcrd[0] = wcsPixelPoint.
x();
2516 pixcrd[1] = wcsPixelPoint.
y();
2518 if ((status = wcsp2s(wcs, 1, 2, &pixcrd[0], &imgcrd[0], &phi, &theta, &world[0], &stat[0])) != 0)
2520 lastError =
QString(
"wcsp2s error %1: %2.").
arg(status).
arg(wcs_errmsg[status]);
2525 wcsCoord.
setRA0(world[0] / 15.0);
2531 Q_UNUSED(wcsPixelPoint);
2537 #if !defined(KSTARS_LITE) && defined(HAVE_WCSLIB)
2538 void FITSData::findObjectsInImage(
double world[],
double phi,
double theta,
double imgcrd[],
double pixcrd[],
2547 if (fits_read_keyword(fptr,
"DATE-OBS", date,
nullptr, &status) == 0)
2550 tsString = tsString.remove(
'\'').trimmed();
2565 if (wcs_coord !=
nullptr)
2583 int type =
object->type();
2594 world[0] =
object->ra0().Degrees();
2595 world[1] =
object->dec0().Degrees();
2597 if ((status = wcss2p(wcs, 1, 2, &world[0], &phi, &theta, &imgcrd[0], &pixcrd[0], &stat[0])) != 0)
2599 fprintf(stderr,
"wcsp2s ERROR %d: %s.\n", status, wcs_errmsg[status]);
2607 if (x > 0 && y > 0 && x < w && y < h)
2623 skyObjectStored = object;
2630 return skyObjectStored;
2655 return flipVCounter;
2660 flipVCounter = value;
2665 return flipHCounter;
2670 flipHCounter = value;
2688 template <
typename T>
2689 bool FITSData::rotFITS(
int rotate,
int mirror)
2693 uint8_t * rotimage =
nullptr;
2698 else if (rotate == 2)
2700 else if (rotate == 3)
2702 else if (rotate < 0)
2703 rotate = rotate + 360;
2713 if (rotimage ==
nullptr)
2715 qWarning() <<
"Unable to allocate memory for rotated image buffer!";
2719 auto * rotBuffer =
reinterpret_cast<T *
>(rotimage);
2720 auto * buffer =
reinterpret_cast<T *
>(m_ImageBuffer);
2723 if (rotate < 45 && rotate > -45)
2727 for (
int i = 0; i < m_Channels; i++)
2730 for (x1 = 0; x1 < nx; x1++)
2733 for (y1 = 0; y1 < ny; y1++)
2734 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2738 else if (mirror == 2)
2740 for (
int i = 0; i < m_Channels; i++)
2743 for (y1 = 0; y1 < ny; y1++)
2746 for (x1 = 0; x1 < nx; x1++)
2747 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2753 for (
int i = 0; i < m_Channels; i++)
2756 for (y1 = 0; y1 < ny; y1++)
2758 for (x1 = 0; x1 < nx; x1++)
2759 rotBuffer[(y1 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2766 else if (rotate >= 45 && rotate < 135)
2770 for (
int i = 0; i < m_Channels; i++)
2773 for (y1 = 0; y1 < ny; y1++)
2776 for (x1 = 0; x1 < nx; x1++)
2779 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2784 else if (mirror == 2)
2786 for (
int i = 0; i < m_Channels; i++)
2789 for (y1 = 0; y1 < ny; y1++)
2791 for (x1 = 0; x1 < nx; x1++)
2792 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
2798 for (
int i = 0; i < m_Channels; i++)
2801 for (y1 = 0; y1 < ny; y1++)
2804 for (x1 = 0; x1 < nx; x1++)
2807 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2818 else if (rotate >= 135 && rotate < 225)
2822 for (
int i = 0; i < m_Channels; i++)
2825 for (y1 = 0; y1 < ny; y1++)
2828 for (x1 = 0; x1 < nx; x1++)
2829 rotBuffer[(y2 * nx) + x1 + offset] = buffer[(y1 * nx) + x1 + offset];
2833 else if (mirror == 2)
2835 for (
int i = 0; i < m_Channels; i++)
2838 for (x1 = 0; x1 < nx; x1++)
2841 for (y1 = 0; y1 < ny; y1++)
2842 rotBuffer[(y1 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2848 for (
int i = 0; i < m_Channels; i++)
2851 for (y1 = 0; y1 < ny; y1++)
2854 for (x1 = 0; x1 < nx; x1++)
2857 rotBuffer[(y2 * nx) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2865 else if (rotate >= 225 && rotate < 315)
2869 for (
int i = 0; i < m_Channels; i++)
2872 for (y1 = 0; y1 < ny; y1++)
2874 for (x1 = 0; x1 < nx; x1++)
2875 rotBuffer[(x1 * ny) + y1 + offset] = buffer[(y1 * nx) + x1 + offset];
2879 else if (mirror == 2)
2881 for (
int i = 0; i < m_Channels; i++)
2884 for (y1 = 0; y1 < ny; y1++)
2887 for (x1 = 0; x1 < nx; x1++)
2890 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2897 for (
int i = 0; i < m_Channels; i++)
2900 for (y1 = 0; y1 < ny; y1++)
2903 for (x1 = 0; x1 < nx; x1++)
2906 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2917 else if (rotate >= 315 && mirror)
2919 for (
int i = 0; i < m_Channels; i++)
2922 for (y1 = 0; y1 < ny; y1++)
2924 for (x1 = 0; x1 < nx; x1++)
2928 rotBuffer[(y2 * ny) + x2 + offset] = buffer[(y1 * nx) + x1 + offset];
2934 delete[] m_ImageBuffer;
2935 m_ImageBuffer = rotimage;
2940 void FITSData::rotWCSFITS(
int angle,
int mirror)
2944 double ctemp1, ctemp2, ctemp3, ctemp4, naxis1, naxis2;
2945 int WCS_DECIMALS = 6;
2947 naxis1 = stats.
width;
2950 if (fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
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);
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);
2971 if (!fits_read_key_dbl(fptr,
"CROTA1", &ctemp1, comment, &status))
2972 fits_update_key_dbl(fptr,
"CROTA1", -ctemp1, WCS_DECIMALS, comment, &status);
2974 if (!fits_read_key_dbl(fptr,
"CROTA2", &ctemp1, comment, &status))
2975 fits_update_key_dbl(fptr,
"CROTA2", -ctemp1, WCS_DECIMALS, comment, &status);
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);
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);
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);
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);
2997 if (!fits_read_key_dbl(fptr,
"LTM1_1", &ctemp1, comment, &status))
3001 if (!fits_read_key_dbl(fptr,
"LTM2_2", &ctemp2, comment, &status))
3002 if (ctemp1 == ctemp2)
3007 if (!fits_read_key_dbl(fptr,
"LTV1", <v1, comment, &status))
3008 fits_delete_key(fptr,
"LTV1", &status);
3009 if (!fits_read_key_dbl(fptr,
"LTV2", <v2, comment, &status))
3010 fits_delete_key(fptr,
"LTV2", &status);
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);
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);
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);
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);
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);
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);
3036 fits_delete_key(fptr,
"LTM1_1", &status);
3037 fits_delete_key(fptr,
"LTM1_2", &status);
3045 if (!fits_read_key_dbl(fptr,
"CRPIX1", &ctemp1, comment, &status) &&
3046 !fits_read_key_dbl(fptr,
"CRPIX2", &ctemp2, comment, &status))
3051 fits_update_key_dbl(fptr,
"CRPIX1", naxis1 - ctemp1, WCS_DECIMALS, comment, &status);
3052 else if (angle == 90)
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);
3057 else if (angle == 180)
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);
3062 else if (angle == 270)
3064 fits_update_key_dbl(fptr,
"CRPIX1", ctemp2, WCS_DECIMALS, comment, &status);
3065 fits_update_key_dbl(fptr,
"CRPIX2", ctemp1, WCS_DECIMALS, comment, &status);
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);
3075 else if (angle == 180)
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);
3080 else if (angle == 270)
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);
3091 if (!fits_read_key_dbl(fptr,
"CDELT1", &ctemp1, comment, &status) &&
3092 !fits_read_key_dbl(fptr,
"CDELT2", &ctemp2, comment, &status))
3097 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3098 else if (angle == 90)
3100 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3101 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3103 else if (angle == 180)
3105 fits_update_key_dbl(fptr,
"CDELT1", ctemp1, WCS_DECIMALS, comment, &status);
3106 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3108 else if (angle == 270)
3110 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3111 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3118 fits_update_key_dbl(fptr,
"CDELT1", -ctemp2, WCS_DECIMALS, comment, &status);
3119 fits_update_key_dbl(fptr,
"CDELT2", ctemp1, WCS_DECIMALS, comment, &status);
3121 else if (angle == 180)
3123 fits_update_key_dbl(fptr,
"CDELT1", -ctemp1, WCS_DECIMALS, comment, &status);
3124 fits_update_key_dbl(fptr,
"CDELT2", -ctemp2, WCS_DECIMALS, comment, &status);
3126 else if (angle == 270)
3128 fits_update_key_dbl(fptr,
"CDELT1", ctemp2, WCS_DECIMALS, comment, &status);
3129 fits_update_key_dbl(fptr,
"CDELT2", -ctemp1, WCS_DECIMALS, comment, &status);
3140 if (!fits_read_key_dbl(fptr,
"CD1_1", &ctemp1, comment, &status))
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);
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);
3153 else if (angle == 90)
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);
3160 else if (angle == 180)
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);
3167 else if (angle == 270)
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);
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);
3184 else if (angle == 180)
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);
3191 else if (angle == 270)
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);
3204 if (!fits_read_key_dbl(fptr,
"CO1_1", &ctemp1, comment, &status))
3209 for (i = 1; i < 13; i++)
3211 sprintf(keyword,
"CO1_%d", i);
3212 fits_delete_key(fptr, keyword, &status);
3214 for (i = 1; i < 13; i++)
3216 sprintf(keyword,
"CO2_%d", i);
3217 fits_delete_key(fptr, keyword, &status);
3225 return m_ImageBuffer;
3230 delete[] m_ImageBuffer;
3231 m_ImageBuffer = buffer;
3234 bool FITSData::checkDebayer()
3237 char bayerPattern[64];
3240 if (fits_read_keyword(fptr,
"BAYERPAT", bayerPattern,
nullptr, &status))
3248 QString pattern(bayerPattern);
3249 pattern = pattern.remove(
'\'').trimmed();
3251 if (pattern ==
"RGGB")
3253 else if (pattern ==
"GBRG")
3255 else if (pattern ==
"GRBG")
3257 else if (pattern ==
"BGGR")
3266 fits_read_key(fptr, TINT,
"XBAYROFF", &debayerParams.
offsetX,
nullptr, &status);
3267 fits_read_key(fptr, TINT,
"YBAYROFF", &debayerParams.
offsetY,
nullptr, &status);
3292 if (bayerBuffer ==
nullptr)
3294 int anynull = 0, status = 0;
3296 bayerBuffer = m_ImageBuffer;
3298 if (fits_read_img(fptr, m_DataType, 1, stats.
samples_per_channel,
nullptr, bayerBuffer, &anynull, &status))
3301 fits_get_errstatus(status, errmsg);
3325 auto * destinationBuffer =
new uint8_t[rgb_size];
3327 if (destinationBuffer ==
nullptr)
3329 KSNotification::error(i18n(
"Unable to allocate memory for temporary bayer buffer."), i18n(
"Debayer error"));
3333 int ds1394_height = stats.
height;
3334 uint8_t * dc1394_source = bayerBuffer;
3336 if (debayerParams.
offsetY == 1)
3338 dc1394_source += stats.
width;
3342 if (debayerParams.
offsetX == 1)
3354 delete[] destinationBuffer;
3358 if (m_Channels == 1)
3360 delete[] m_ImageBuffer;
3361 m_ImageBuffer =
new uint8_t[rgb_size];
3363 if (m_ImageBuffer ==
nullptr)
3365 delete[] destinationBuffer;
3366 KSNotification::error(i18n(
"Unable to allocate memory for temporary bayer buffer."), i18n(
"Debayer error"));
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);
3378 for (
int i = 0; i <= imax; i += 3)
3380 *rBuff++ = destinationBuffer[i];
3381 *gBuff++ = destinationBuffer[i + 1];
3382 *bBuff++ = destinationBuffer[i + 2];
3387 delete[] destinationBuffer;
3388 bayerBuffer =
nullptr;
3397 auto * destinationBuffer =
new uint8_t[rgb_size];
3399 auto * buffer =
reinterpret_cast<uint16_t *
>(bayerBuffer);
3400 auto * dstBuffer =
reinterpret_cast<uint16_t *
>(destinationBuffer);
3402 if (destinationBuffer ==
nullptr)
3404 KSNotification::error(i18n(
"Unable to allocate memory for temporary bayer buffer."), i18n(
"Debayer error"));
3408 int ds1394_height = stats.
height;
3409 uint16_t * dc1394_source = buffer;
3411 if (debayerParams.
offsetY == 1)
3413 dc1394_source += stats.
width;
3417 if (debayerParams.
offsetX == 1)
3423 debayerParams.
method, 16);
3429 delete[] destinationBuffer;
3433 if (m_Channels == 1)
3435 delete[] m_ImageBuffer;
3436 m_ImageBuffer =
new uint8_t[rgb_size];
3438 if (m_ImageBuffer ==
nullptr)
3440 delete[] destinationBuffer;
3441 KSNotification::error(i18n(
"Unable to allocate memory for temporary bayer buffer."), i18n(
"Debayer error"));
3446 buffer =
reinterpret_cast<uint16_t *
>(m_ImageBuffer);
3450 uint16_t * rBuff = buffer;
3451 uint16_t * gBuff = buffer + (stats.
width * stats.
height);
3452 uint16_t * bBuff = buffer + (stats.
width * stats.
height * 2);
3455 for (
int i = 0; i <= imax; i += 3)
3457 *rBuff++ = dstBuffer[i];
3458 *gBuff++ = dstBuffer[i + 1];
3459 *bBuff++ = dstBuffer[i + 2];
3463 delete[] destinationBuffer;
3464 bayerBuffer =
nullptr;
3471 for (
int i = 0; i < m_Channels; i++)
3472 adu += stats.
mean[i];
3474 return (adu / static_cast<double>(m_Channels));
3497 template <
typename T>
3504 for (
int y = 0; y < stats.
height; y++)
3506 size_t yOffset = y * stats.
width;
3507 const T * grayLine =
reinterpret_cast<T *
>(m_ImageBuffer) + yOffset;
3509 const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.
width;
3510 const T * grayLine_p1 = y >= stats.
height - 1 ? grayLine : grayLine + stats.
width;
3512 float * gradientLine = gradient.
data() + yOffset;
3513 float * directionLine = direction.
data() + yOffset;
3515 for (
int x = 0; x < stats.
width; x++)
3517 int x_m1 = x < 1 ? x : x - 1;
3518 int x_p1 = x >= stats.
width - 1 ? x : x + 1;
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];
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];
3526 gradientLine[x] = qAbs(gradX) + qAbs(gradY);
3554 if (gradX == 0 && gradY == 0)
3555 directionLine[x] = 0;
3556 else if (gradX == 0)
3557 directionLine[x] = 3;
3560 qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
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;
3569 directionLine[x] = 3;
3579 for (
int y = 1; y < height - 1; y++)
3581 for (
int x = 1; x < width - 1; x++)
3583 int index = x + y *
width;
3584 float val = gradient[index];
3585 if (val > 0 && ids[index] == 0)
3587 trace(width, height, ++
id, gradient, ids, x, y);
3598 int yOffset = y *
width;
3599 float * cannyLine = image.
data() + yOffset;
3600 int * idLine = ids.
data() + yOffset;
3607 for (
int j = -1; j < 2; j++)
3611 if (nextY < 0 || nextY >= height)
3614 float * cannyLineNext = cannyLine + j *
width;
3616 for (
int i = -1; i < 2; i++)
3620 if (i == j || nextX < 0 || nextX >= width)
3623 if (cannyLineNext[nextX] > 0)
3626 trace(width, height,
id, image, ids, nextX, nextY);
3639 return autoRemoveTemporaryFITS;
3644 autoRemoveTemporaryFITS = value;
3648 template <
typename T>
3649 void FITSData::convertToQImage(
double dataMin,
double dataMax,
double scale,
double zero,
QImage &image)
3651 #pragma GCC diagnostic push
3652 #pragma GCC diagnostic ignored "-Wcast-align"
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();
3660 uint32_t size = w * h;
3666 for (
int j = 0; j < h; j++)
3668 unsigned char * scanLine = image.
scanLine(j);
3670 for (
int i = 0; i < w; i++)
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);
3680 double rval = 0, gval = 0, bval = 0;
3683 for (
int j = 0; j < h; j++)
3685 auto * scanLine =
reinterpret_cast<QRgb *
>((image.
scanLine(j)));
3687 for (
int i = 0; i < w; i++)
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);
3693 value = qRgb(rval * scale + zero, gval * scale + zero, bval * scale + zero);
3695 scanLine[i] = value;
3712 if (future.
result() ==
false)
3719 fitsImage.
fill(Qt::white);
3728 for (
int i = 0; i < 256; i++)
3729 fitsImage.
setColor(i, qRgb(i, i, i));
3736 double dataMin = data.stats.
mean[0] - data.stats.
stddev[0];
3737 double dataMax = data.stats.
mean[0] + data.stats.
stddev[0] * 3;
3739 double bscale = 255. / (dataMax - dataMin);
3740 double bzero = (-dataMin) * (255. / (dataMax - dataMin));
3746 data.convertToQImage<uint8_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3750 data.convertToQImage<int16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3754 data.convertToQImage<uint16_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3758 data.convertToQImage<int32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3762 data.convertToQImage<uint32_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3766 data.convertToQImage<
float>(dataMin, dataMax, bscale, bzero, fitsImage);
3770 data.convertToQImage<int64_t>(dataMin, dataMax, bscale, bzero, fitsImage);
3774 data.convertToQImage<
double>(dataMin, dataMax, bscale, bzero, fitsImage);
3788 qCCritical(KSTARS_FITS) <<
"Failed to convert" << filename <<
"to FITS since format" << format <<
"is not supported in Qt";
3796 qCCritical(KSTARS_FITS) <<
"Failed to open image" <<
filename;
3803 fitsfile *fptr =
nullptr;
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};
3809 if (fits_create_file(&fptr,
QString(
'!' + output).toLatin1().data(), &status))
3811 qCCritical(KSTARS_FITS) <<
"Failed to create FITS file. Error:" << status;
3815 if (fits_create_img(fptr, BYTE_IMG, naxis, naxes, &status))
3817 qCWarning(KSTARS_FITS) <<
"fits_create_img failed:" << error_status;
3819 fits_close_file(fptr, &status);
3824 fits_update_key(fptr, TLONG,
"EXPOSURE", &exposure,
"Total Exposure Time", &status);
3829 nelements = naxes[0] * naxes[1];
3830 if (fits_write_img(fptr, TBYTE, fpixel, nelements, input.
bits(), &status))
3832 fits_get_errstatus(status, error_status);
3833 qCWarning(KSTARS_FITS) <<
"fits_write_img GRAY failed:" << error_status;
3835 fits_close_file(fptr, &status);
3842 nelements = naxes[0] * naxes[1] * 3;
3844 uint8_t *srcBuffer = input.
bits();
3846 uint32_t srcBytes = naxes[0] * naxes[1] * 4 - 4;
3848 uint8_t *rgbBuffer =
new uint8_t[nelements];
3849 if (rgbBuffer ==
nullptr)
3851 qCWarning(KSTARS_FITS) <<
"Not enough memory for RGB buffer";
3852 fits_close_file(fptr, &status);
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)
3861 *subB++ = srcBuffer[i];
3862 *subG++ = srcBuffer[i + 1];
3863 *subR++ = srcBuffer[i + 2];
3866 if (fits_write_img(fptr, TBYTE, fpixel, nelements, rgbBuffer, &status))
3868 fits_get_errstatus(status, error_status);
3869 qCWarning(KSTARS_FITS) <<
"fits_write_img RGB failed:" << error_status;
3871 fits_close_file(fptr, &status);
3872 delete [] rgbBuffer;
3876 delete [] rgbBuffer;
3879 if (fits_flush_file(fptr, &status))
3881 fits_get_errstatus(status, error_status);
3882 qCWarning(KSTARS_FITS) <<
"fits_flush_file failed:" << error_status;
3884 fits_close_file(fptr, &status);
3888 if (fits_close_file(fptr, &status))
3890 fits_get_errstatus(status, error_status);
3891 qCWarning(KSTARS_FITS) <<
"fits_close_file failed:" << error_status;
3901 int status = 0, exttype = 0;
3903 fitsfile * new_fptr;
3906 qCInfo(KSTARS_FITS) <<
"Creating new WCS file:" << newWCSFile <<
"with parameters Orientation:" << orientation
3907 <<
"RA:" << ra <<
"DE:" << dec <<
"Pixel Scale:" << pixscale;
3912 if (fits_create_file(&new_fptr,
QString(
'!' + newWCSFile).toLatin1(), &status))
3914 fits_get_errstatus(status, errMsg);
3916 fits_report_error(stderr, status);
3920 if (fits_movabs_hdu(fptr, 1, &exttype, &status))
3922 fits_get_errstatus(status, errMsg);
3924 fits_report_error(stderr, status);
3928 if (fits_copy_file(fptr, new_fptr, 1, 1, 1, &status))
3930 fits_get_errstatus(status, errMsg);
3932 fits_report_error(stderr, status);
3937 if (fits_close_file(fptr, &status))
3939 fits_get_errstatus(status, errMsg);
3941 fits_report_error(stderr, status);
3947 if (m_isTemporary && autoRemoveTemporaryFITS)
3950 m_isTemporary =
false;
3951 qCDebug(KSTARS_FITS) <<
"Removing FITS File: " << m_Filename;
3954 m_Filename = newWCSFile;
3955 m_isTemporary =
true;
3959 if (fits_movabs_hdu(fptr, 1, &exttype, &status))
3961 fits_get_errstatus(status, errMsg);
3963 fits_report_error(stderr, status);
3968 if (fits_write_img(fptr, m_DataType, 1, nelements, m_ImageBuffer, &status))
3970 fits_get_errstatus(status, errMsg);
3972 fits_report_error(stderr, status);
3979 if (fits_update_key(fptr, TDOUBLE,
"DATAMIN", &(stats.
min),
"Minimum value", &status))
3981 fits_get_errstatus(status, errMsg);
3983 fits_report_error(stderr, status);
3988 if (fits_update_key(fptr, TDOUBLE,
"DATAMAX", &(stats.
max),
"Maximum value", &status))
3990 fits_get_errstatus(status, errMsg);
3992 fits_report_error(stderr, status);
3997 if (fits_update_key(fptr, TUSHORT,
"NAXIS1", &(stats.
width),
"length of data axis 1", &status))
3999 fits_get_errstatus(status, errMsg);
4001 fits_report_error(stderr, status);
4006 if (fits_update_key(fptr, TUSHORT,
"NAXIS2", &(stats.
height),
"length of data axis 2", &status))
4008 fits_get_errstatus(status, errMsg);
4010 fits_report_error(stderr, status);
4014 fits_update_key(fptr, TDOUBLE,
"OBJCTRA", &ra,
"Object RA", &status);
4015 fits_update_key(fptr, TDOUBLE,
"OBJCTDEC", &dec,
"Object DEC", &status);
4019 fits_update_key(fptr, TINT,
"EQUINOX", &epoch,
"Equinox", &status);
4021 fits_update_key(fptr, TDOUBLE,
"CRVAL1", &ra,
"CRVAL1", &status);
4022 fits_update_key(fptr, TDOUBLE,
"CRVAL2", &dec,
"CRVAL1", &status);
4024 char radecsys[8] =
"FK5";
4025 char ctype1[16] =
"RA---TAN";
4026 char ctype2[16] =
"DEC--TAN";
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);
4032 double crpix1 =
width() / 2.0;
4033 double crpix2 =
height() / 2.0;
4035 fits_update_key(fptr, TDOUBLE,
"CRPIX1", &crpix1,
"CRPIX1", &status);
4036 fits_update_key(fptr, TDOUBLE,
"CRPIX2", &crpix2,
"CRPIX2", &status);
4039 double secpix1 = pixscale;
4040 double secpix2 = pixscale;
4042 fits_update_key(fptr, TDOUBLE,
"SECPIX1", &secpix1,
"SECPIX1", &status);
4043 fits_update_key(fptr, TDOUBLE,
"SECPIX2", &secpix2,
"SECPIX2", &status);
4045 double degpix1 = secpix1 / 3600.0;
4046 double degpix2 = secpix2 / 3600.0;
4048 fits_update_key(fptr, TDOUBLE,
"CDELT1", °pix1,
"CDELT1", &status);
4049 fits_update_key(fptr, TDOUBLE,
"CDELT2", °pix2,
"CDELT2", &status);
4052 double rotation = 360 - orientation;
4056 fits_update_key(fptr, TDOUBLE,
"CROTA1", &rotation,
"CROTA1", &status);
4057 fits_update_key(fptr, TDOUBLE,
"CROTA2", &rotation,
"CROTA2", &status);
4060 if (fits_write_date(fptr, &status))
4062 fits_get_errstatus(status, errMsg);
4064 fits_report_error(stderr, status);
4071 if (fits_write_history(fptr, history.
toLatin1(), &status))
4073 fits_get_errstatus(status, errMsg);
4075 fits_report_error(stderr, status);
4079 fits_flush_file(fptr, &status);
4083 qCDebug(KSTARS_FITS) <<
"Finished creating WCS file: " << newWCSFile;
4093 fits_update_key(fptr, TDOUBLE,
"OBJCTRA", &ra,
"Object RA", &status);
4094 fits_update_key(fptr, TDOUBLE,
"OBJCTDEC", &dec,
"Object DEC", &status);
4098 fits_update_key(fptr, TINT,
"EQUINOX", &epoch,
"Equinox", &status);
4100 fits_update_key(fptr, TDOUBLE,
"CRVAL1", &ra,
"CRVAL1", &status);
4101 fits_update_key(fptr, TDOUBLE,
"CRVAL2", &dec,
"CRVAL1", &status);
4103 char radecsys[8] =
"FK5";
4104 char ctype1[16] =
"RA---TAN";
4105 char ctype2[16] =
"DEC--TAN";
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);
4111 double crpix1 =
width() / 2.0;
4112 double crpix2 =
height() / 2.0;
4114 fits_update_key(fptr, TDOUBLE,
"CRPIX1", &crpix1,
"CRPIX1", &status);
4115 fits_update_key(fptr, TDOUBLE,
"CRPIX2", &crpix2,
"CRPIX2", &status);
4118 double secpix1 = pixscale;
4119 double secpix2 = pixscale;
4121 fits_update_key(fptr, TDOUBLE,
"SECPIX1", &secpix1,
"SECPIX1", &status);
4122 fits_update_key(fptr, TDOUBLE,
"SECPIX2", &secpix2,
"SECPIX2", &status);
4124 double degpix1 = secpix1 / 3600.0;
4125 double degpix2 = secpix2 / 3600.0;
4127 fits_update_key(fptr, TDOUBLE,
"CDELT1", °pix1,
"CDELT1", &status);
4128 fits_update_key(fptr, TDOUBLE,
"CDELT2", °pix2,
"CDELT2", &status);
4131 double rotation = 360 - orientation;
4135 fits_update_key(fptr, TDOUBLE,
"CROTA1", &rotation,
"CROTA1", &status);
4136 fits_update_key(fptr, TDOUBLE,
"CROTA2", &rotation,
"CROTA2", &status);
4140 qCDebug(KSTARS_FITS) <<
"Finished update WCS info.";
4147 return (point.
x() >= 0 && point.
y() >= 0 && point.
x() <= stats.
width && point.
y() <= stats.
height);
4152 int x = 0, y = 0, w = stats.
width, h = stats.
height, maxRadius = 50;
4158 w = boundary.
width();
4163 auto * data =
new float[w * h];
4168 getFloatBuffer<uint8_t>(data, x, y, w, h);
4171 getFloatBuffer<int16_t>(data, x, y, w, h);
4174 getFloatBuffer<uint16_t>(data, x, y, w, h);
4177 getFloatBuffer<int32_t>(data, x, y, w, h);
4180 getFloatBuffer<uint32_t>(data, x, y, w, h);
4184 data =
reinterpret_cast<float *
>(m_ImageBuffer);
4187 getFloatBuffer<int64_t>(data, x, y, w, h);
4190 getFloatBuffer<double>(data, x, y, w, h);
4197 float * imback =
nullptr;
4198 double * flux =
nullptr, *fluxerr =
nullptr, *area =
nullptr;
4199 short * flag =
nullptr;
4200 short flux_flag = 0;
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 };
4210 sep_image im = {data,
nullptr,
nullptr,
SEP_TFLOAT, 0, 0, w, h, 0.0,
SEP_NOISE_NONE, 1.0, 0.0};
4214 if (status != 0)
goto exit;
4217 imback = (
float *)malloc((w * h) *
sizeof(float));
4219 if (status != 0)
goto exit;
4223 if (status != 0)
goto exit;
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;
4235 for (
int i = 0; i < catalog->
nobj; i++)
4237 double flux = catalog->
flux[i];
4239 sep_flux_radius(&im, catalog->
x[i], catalog->
y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag);
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];
4246 center->
HFR = center->
width = flux_fractions[0];
4247 if (flux_fractions[1] < maxRadius)
4248 center->
width = flux_fractions[1] * 2;
4253 std::sort(edges.
begin(), edges.
end(), [](
const Edge * edge1,
const Edge * edge2) ->
bool {
return edge1->
width > edge2->width;});
4257 int starCount = qMin(100, edges.
count());
4258 for (
int i = 0; i < starCount; i++)
4259 starCenters.append(edges[i]);
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;
4270 if (stats.
bitpix != FLOAT_IMG)
4282 char errorMessage[512];
4284 qCritical(KSTARS_FITS) << errorMessage;
4288 return starCenters.count();
4291 template <
typename T>
4294 auto * rawBuffer =
reinterpret_cast<T *
>(m_ImageBuffer);
4296 float * floatPtr = buffer;
4301 for (
int y1 = y; y1 < y2; y1++)
4303 int offset = y1 * stats.
width;
4304 for (
int x1 = x; x1 < x2; x1++)
4306 *floatPtr++ = rawBuffer[offset + x1];
static bool autoDebayer()
Get Automatically debayer a FITS image if it is contains a bayer pattern.
FITSSkyObject(SkyObject *object, int xPos, int yPos)
void constructHistogram()
QString & append(QChar ch)
bool load(QIODevice *device, const char *format)
dc1394color_filter_t filter
bool pixelToWCS(const QPointF &wcsPixelPoint, SkyPoint &wcsCoord)
pixelToWCS Convert Pixel coordinates to J2000 world coordinates
const QString & filename() const
bool injectWCS(double orientation, double ra, double dec, double pixscale)
injectWCS Add WCS keywords to file
void sep_bkg_free(sep_bkg *bkg)
void append(const T &value)
void applyFilter(FITSScale type, uint8_t *image=nullptr, QVector< double > *targetMin=nullptr, QVector< double > *targetMax=nullptr)
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.
Stats struct to hold statisical data about the FITS data.
QStringList split(const QString &sep, SplitBehavior behavior, Qt::CaseSensitivity cs) const
void getCenterSelection(int *x, int *y)
void saveStatistics(Statistic &other)
void getFloatBuffer(float *buffer, int x, int y, int w, int h)
bool getAutoRemoveTemporaryFITS() const
static KStarsData * Instance()
const double & Degrees() const
void restoreStatistics(Statistic &other)
void setColorCount(int colorCount)
const CachingDms & dec0() const
int getFlipHCounter() const
void sep_get_errmsg(int status, char *errtext)
QList< SkyObject * > findObjectsInArea(const SkyPoint &p1, const SkyPoint &p2)
iterator erase(iterator pos)
static QString writableLocation(QStandardPaths::StandardLocation type)
void error(const QString &message, const QString &title)
static void rotate(double *a, double *b, double *c)
QString & remove(int position, int n)
int fp_init(fpstate *fpptr)
#define LOW_EDGE_CUTOFF_2
int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype)
static QImage FITSToImage(const QString &filename)
bool copy(const QString &newName)
void setColor(int index, QRgb colorValue)
QList< Key > keys() const
int getRotCounter() const
int getBytesPerPixel() const
int count(const T &value) const
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)
void append(const T &value)
QString fromUtf8(const char *str, int size)
The sky coordinates of a point in the sky.
QVariant property(const char *name) const
int toInt(bool *ok) const
void fill(uint pixelValue)
#define DIFFUSE_THRESHOLD
void setFlipHCounter(int value)
int fp_unpack(char *infits, char *outfits, fpstate fpvar)
QList< Edge * > getStarCentersInSubFrame(QRect subFrame) const
bool startsWith(const QString &s, Qt::CaseSensitivity cs) const
int saveFITS(const QString &newFilename)
static bool auto3DCube()
Get Process 3D FITS Cube (RGB).
void getMinMax(double *min, double *max, uint8_t channel=0) const
QVector< uint32_t > getCumulativeFrequency(int channel=0) const
bool endsWith(const QString &s, Qt::CaseSensitivity cs) const
static double focusThreshold()
Get Relative percentage strength of centroid edge pixel strength to average pixel value...
QFuture< T > run(Function function,...)
void calculateStats(bool refresh=false)
bool contains(const QPoint &point, bool proper) const
QFuture< void > map(Sequence &sequence, MapFunction function)
const CachingDms & ra0() const
QList< FITSSkyObject * > getSkyObjects()
SkyMapComposite is the root object in the object hierarchy of the sky map.
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day...
void setBayerParams(BayerParams *param)
Structure to hold FITS Header records.
int findStars(StarAlgorithm algorithm=ALGORITHM_CENTROID, const QRect &trackingBox=QRect())
SkyMapComposite * skyComposite()
An angle, stored as degrees, but expressible in many ways.
double getHFR(HFRType type=HFR_AVERAGE)
virtual qint64 size() const
int findOneStar(const QRect &boundary)
bool loadFITSFromMemory(const QString &inFilename, void *fits_buffer, size_t fits_buffer_size, bool silent)
loadFITSFromMemory Loading FITS from memory buffer.
QDateTime fromString(const QString &string, Qt::DateFormat format)
#define MINIMUM_EDGE_LIMIT
bool getRecordValue(const QString &key, QVariant &value) const
QFuture< bool > loadFITS(const QString &inFilename, bool silent=true)
loadFITS Loading FITS file asynchronously.
bool contains(const QPointF &point) const
void getBayerParams(BayerParams *param)
bool wcsToPixel(SkyPoint &wcsCoord, QPointF &wcsPixelPoint, QPointF &wcsImagePoint)
wcsToPixel Given J2000 (RA0,DE0) coordinates.
void setImageBuffer(uint8_t *buffer)
QList< FITSSkyObject * > objList
void setRotCounter(int value)
int filterStars(const float innerRadius, const float outerRadius)
dc1394bayer_method_t method
static int findCannyStar(FITSData *data, const QRect &boundary=QRect())
double getBinWidth(int channel=0)
There are several time-dependent values used in position calculations, that are not specific to an ob...
QDateTime currentDateTime()
QByteArray toLatin1() const
QString mid(int position, int n) const
char * toString(const T &value)
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype)
virtual void updateCoordsNow(const KSNumbers *num)
updateCoordsNow Shortcut for updateCoords( const KSNumbers *num, false, nullptr, nullptr, true)
void setFlipVCounter(int value)
bool greaterThan(Edge *s1, Edge *s2)
QList< T > mid(int pos, int length) const
int getFlipVCounter() const
uint8_t * getImageBuffer()
QString getLastError() const
void setMinMax(double newMin, double newMax, uint8_t channel=0)
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(FITSMode fitsMode=FITS_NORMAL)
wcs_point * getWCSCoord()
uint32_t samples_per_channel
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)
void setDec0(dms d)
Sets Dec0, the catalog Declination.
#define LOW_EDGE_CUTOFF_1
void sep_catalog_free(sep_catalog *catalog)
const int MINIMUM_ROWS_PER_CENTER
virtual QString name(void) const
void appendStar(Edge *newCenter)
int sep_background(sep_image *image, int bw, int bh, int fw, int fh, double fthresh, sep_bkg **bkg)
int findSEPStars(const QRect &boundary=QRect())
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.
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
dc1394error_t
Error codes returned by most libdc1394 functions.
static bool ImageToFITS(const QString &filename, const QString &format, QString &output)
ImageToFITS Convert an image file with supported format to a FITS file.
double getJMIndex() const
const T value(const Key &key) const
void setAutoRemoveTemporaryFITS(bool value)