Differential D28292 Diff 79624 kstars/fitsviewer/fitsdata.cpp

Changeset View

Standalone View

kstars/fitsviewer/fitsdata.cpp

Show All 12 Lines
13	* the Free Software Foundation; either version 2 of the License, or *	13		* the Free Software Foundation; either version 2 of the License, or *
14	* (at your option) any later version. *	14		* (at your option) any later version. *
15	* *	15		* *
16	* Some code fragments were adapted from Peter Kirchgessner's FITS plugin*	16		* Some code fragments were adapted from Peter Kirchgessner's FITS plugin*
17	* See http://members.aol.com/pkirchg for more details. *	17		* See http://members.aol.com/pkirchg for more details. *
18	***************************************************************************/	18		***************************************************************************/
19		19
20	#include "fitsdata.h"	20		#include "fitsdata.h"
		21		#include "fitsthresholddetector.h"
		22		#include "fitsgradientdetector.h"
		23		#include "fitscentroiddetector.h"
		24		#include "fitssepdetector.h"
21		25
22	#include "sep/sep.h"
23	#include "fpack.h"	26		#include "fpack.h"
24		27
25	#include "kstarsdata.h"	28		#include "kstarsdata.h"
26	#include "ksutils.h"	29		#include "ksutils.h"
27	#include "kspaths.h"	30		#include "kspaths.h"
28	#include "Options.h"	31		#include "Options.h"
29	#include "skymapcomposite.h"	32		#include "skymapcomposite.h"
30	#include "auxiliary/ksnotification.h"	33		#include "auxiliary/ksnotification.h"
Show All 19 Lines
50	#include <fits_debug.h>	53		#include <fits_debug.h>
51		54
52	#define ZOOM_DEFAULT 100.0	55		#define ZOOM_DEFAULT 100.0
53	#define ZOOM_MIN 10	56		#define ZOOM_MIN 10
54	#define ZOOM_MAX 400	57		#define ZOOM_MAX 400
55	#define ZOOM_LOW_INCR 10	58		#define ZOOM_LOW_INCR 10
56	#define ZOOM_HIGH_INCR 50	59		#define ZOOM_HIGH_INCR 50
57		60
58	const int MINIMUM_ROWS_PER_CENTER = 3;
59	const QString FITSData::m_TemporaryPath = QStandardPaths::writableLocation(QStandardPaths::TempLocation);	61		const QString FITSData::m_TemporaryPath = QStandardPaths::writableLocation(QStandardPaths::TempLocation);
60		62
61	#define DIFFUSE_THRESHOLD 0.15
62		63
63	#define MAX_EDGE_LIMIT 10000
64	#define LOW_EDGE_CUTOFF_1 50
65	#define LOW_EDGE_CUTOFF_2 10
66	#define MINIMUM_EDGE_LIMIT 2
67
68	bool greaterThan(Edge * s1, Edge * s2)
69	{
70	//return s1->width > s2->width;
71	return s1->sum > s2->sum;
72	}
73		64
74	FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)	65		FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode)
75	{	66		{
76	debayerParams.method = DC1394_BAYER_METHOD_NEAREST;	67		debayerParams.method = DC1394_BAYER_METHOD_NEAREST;
77	debayerParams.filter = DC1394_COLOR_FILTER_RGGB;	68		debayerParams.filter = DC1394_COLOR_FILTER_RGGB;
78	debayerParams.offsetX = debayerParams.offsetY = 0;	69		debayerParams.offsetX = debayerParams.offsetY = 0;
79	}	70		}
80		71
▲ Show 20 Lines • Show All 797 Lines • ▼ Show 20 Line(s)		868	{
878	value = oneRecord->value;	869		value = oneRecord->value;
879	return true;	870		return true;
880	}	871		}
881	}	872		}
882		873
883	return false;	874		return false;
884	}	875		}
885		876
886	bool FITSData::checkCollision(Edge * s1, Edge * s2)
887	{
888	int dis; //distance
889
890	int diff_x = s1->x - s2->x;
891	int diff_y = s1->y - s2->y;
892
893	dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y));
894	dis -= s1->width / 2;
895	dis -= s2->width / 2;
896
897	if (dis <= 0) //collision
898	return true;
899
900	//no collision
901	return false;
902	}
903
904	int FITSData::findCannyStar(FITSData * data, const QRect &boundary)
905	{
906	switch (data->property("dataType").toInt())
907	{
908	case TBYTE:
909	return FITSData::findCannyStar<uint8_t>(data, boundary);
910
911	case TSHORT:
912	return FITSData::findCannyStar<int16_t>(data, boundary);
913
914	case TUSHORT:
915	return FITSData::findCannyStar<uint16_t>(data, boundary);
916
917	case TLONG:
918	return FITSData::findCannyStar<int32_t>(data, boundary);
919
920	case TULONG:
921	return FITSData::findCannyStar<uint16_t>(data, boundary);
922
923	case TFLOAT:
924	return FITSData::findCannyStar<float>(data, boundary);
925
926	case TLONGLONG:
927	return FITSData::findCannyStar<int64_t>(data, boundary);
928
929	case TDOUBLE:
930	return FITSData::findCannyStar<double>(data, boundary);
931
932	default:
933	break;
934	}
935
936	return 0;
937	}
938
939	int FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox)	877		int FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox)
940	{	878		{
941	int count = 0;	879		int count = 0;
942	starAlgorithm = algorithm;	880		starAlgorithm = algorithm;
943		881
944	qDeleteAll(starCenters);	882		qDeleteAll(starCenters);
945	starCenters.clear();	883		starCenters.clear();
946		884
947	switch (algorithm)	885		switch (algorithm)
948	{	886		{
949	case ALGORITHM_SEP:	887		case ALGORITHM_SEP:
950	count = findSEPStars(trackingBox);	888		count = FITSSEPDetector(this)
		889		.findSources(starCenters, trackingBox);
951	break;	890		break;
952		891
953	case ALGORITHM_GRADIENT:	892		case ALGORITHM_GRADIENT:
954	count = findCannyStar(this, trackingBox);	893		count = FITSGradientDetector(this)
		894		.findSources(starCenters, trackingBox);
955	break;	895		break;
956		896
957	case ALGORITHM_CENTROID:	897		case ALGORITHM_CENTROID:
958	count = findCentroid(trackingBox);	898		count = FITSCentroidDetector(this)
		899		.findSources(starCenters, trackingBox);
959	break;	900		break;
960		901
961	case ALGORITHM_THRESHOLD:	902		case ALGORITHM_THRESHOLD:
962	count = findOneStar(trackingBox);	903		count = FITSThresholdDetector(this)
		904		.configure("Threshold", Options::focusThreshold())
		905		.findSources(starCenters, trackingBox);
963	break;	906		break;
964	}	907		}
965		908
966	starsSearched = true;	909		starsSearched = true;
967		910
968	return count;	911		return count;
969	}	912		}
970		913
Show All 10 Lines		922	{
981	long const y = edge->y - this->height() / 2;	924		long const y = edge->y - this->height() / 2;
982	long const sqRadius = x * x + y * y;	925		long const sqRadius = x * x + y * y;
983	return sqRadius < sqInnerRadius \|\| sqOuterRadius < sqRadius;	926		return sqRadius < sqInnerRadius \|\| sqOuterRadius < sqRadius;
984	}), starCenters.end());	927		}), starCenters.end());
985		928
986	return starCenters.count();	929		return starCenters.count();
987	}	930		}
988		931
989	template <typename T>
990	int FITSData::findCannyStar(FITSData * data, const QRect &boundary)
991	{
992	int subX = qMax(0, boundary.isNull() ? 0 : boundary.x());
993	int subY = qMax(0, boundary.isNull() ? 0 : boundary.y());
994	int subW = (boundary.isNull() ? data->width() : boundary.width());
995	int subH = (boundary.isNull() ? data->height() : boundary.height());
996
997	int BBP = data->getBytesPerPixel();
998
999	uint16_t dataWidth = data->width();
1000
1001	// #1 Find offsets
1002	uint32_t size = subW * subH;
1003	uint32_t offset = subX + subY * dataWidth;
1004
1005	// #2 Create new buffer
1006	auto * buffer = new uint8_t[size * BBP];
1007	// If there is no offset, copy whole buffer in one go
1008	if (offset == 0)
1009	memcpy(buffer, data->getImageBuffer(), size * BBP);
1010	else
1011	{
1012	uint8_t * dataPtr = buffer;
1013	uint8_t * origDataPtr = data->getImageBuffer();
1014	uint32_t lineOffset = 0;
1015	// Copy data line by line
1016	for (int height = subY; height < (subY + subH); height++)
1017	{
1018	lineOffset = (subX + height * dataWidth) * BBP;
1019	memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
1020	dataPtr += (subW * BBP);
1021	}
1022	}
1023
1024	// #3 Create new FITSData to hold it
1025	auto * boundedImage = new FITSData();
1026	boundedImage->stats.width = subW;
1027	boundedImage->stats.height = subH;
1028	boundedImage->stats.bitpix = data->stats.bitpix;
1029	boundedImage->stats.bytesPerPixel = data->stats.bytesPerPixel;
1030	boundedImage->stats.samples_per_channel = size;
1031	boundedImage->stats.ndim = 2;
1032
1033	boundedImage->setProperty("dataType", data->property("dataType"));
1034
1035	// #4 Set image buffer and calculate stats.
1036	boundedImage->setImageBuffer(buffer);
1037
1038	boundedImage->calculateStats(true);
1039
1040	// #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain
1041	boundedImage->applyFilter(FITS_MEDIAN);
1042	boundedImage->applyFilter(FITS_HIGH_CONTRAST);
1043
1044	// #6 Perform Sobel to find gradients and their directions
1045	QVector<float> gradients;
1046	QVector<float> directions;
1047
1048	// TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed
1049	// and discarded whenever necessary. It won't work on noisy images unless this is done.
1050	boundedImage->sobel<T>(gradients, directions);
1051
1052	QVector<int> ids(gradients.size());
1053
1054	int maxID = boundedImage->partition(subW, subH, gradients, ids);
1055
1056	// Not needed anymore
1057	delete boundedImage;
1058
1059	if (maxID == 0)
1060	return 0;
1061
1062	typedef struct
1063	{
1064	float massX = 0;
1065	float massY = 0;
1066	float totalMass = 0;
1067	} massInfo;
1068
1069	QMap<int, massInfo> masses;
1070
1071	// #7 Calculate center of mass for all detected regions
1072	for (int y = 0; y < subH; y++)
1073	{
1074	for (int x = 0; x < subW; x++)
1075	{
1076	int index = x + y * subW;
1077
1078	int regionID = ids[index];
1079	if (regionID > 0)
1080	{
1081	float pixel = gradients[index];
1082
1083	masses[regionID].totalMass += pixel;
1084	masses[regionID].massX += x * pixel;
1085	masses[regionID].massY += y * pixel;
1086	}
1087	}
1088	}
1089
1090	// Compare multiple masses, and only select the highest total mass one as the desired star
1091	int maxRegionID = 1;
1092	int maxTotalMass = masses[1].totalMass;
1093	double totalMassRatio = 1e6;
1094	for (auto key : masses.keys())
1095	{
1096	massInfo oneMass = masses.value(key);
1097	if (oneMass.totalMass > maxTotalMass)
1098	{
1099	totalMassRatio = oneMass.totalMass / maxTotalMass;
1100	maxTotalMass = oneMass.totalMass;
1101	maxRegionID = key;
1102	}
1103	}
1104
1105	// If image has many regions and there is no significant relative center of mass then it's just noise and no stars
1106	// are probably there above a useful threshold.
1107	if (maxID > 10 && totalMassRatio < 1.5)
1108	return 0;
1109
1110	auto * center = new Edge;
1111	center->width = -1;
1112	center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
1113	center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
1114	center->HFR = 1;
1115
1116	// Maximum Radius
1117	int maxR = qMin(subW - 1, subH - 1) / 2;
1118
1119	for (int r = maxR; r > 1; r--)
1120	{
1121	int pass = 0;
1122
1123	for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
1124	{
1125	int testX = center->x + std::cos(theta) * r;
1126	int testY = center->y + std::sin(theta) * r;
1127
1128	// if out of bound, break;
1129	if (testX < 0 \|\| testX >= subW \|\| testY < 0 \|\| testY >= subH)
1130	break;
1131
1132	if (gradients[testX + testY * subW] > 0)
1133	//if (thresholded[testX + testY * subW] > 0)
1134	{
1135	if (++pass >= 24)
1136	{
1137	center->width = r * 2;
1138	// Break of outer loop
1139	r = 0;
1140	break;
1141	}
1142	}
1143	}
1144	}
1145
1146	qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width;
1147
1148	// If no stars were detected
1149	if (center->width == -1)
1150	{
1151	delete center;
1152	return 0;
1153	}
1154
1155	// 30% fuzzy
1156	//center->width += center->width0.3 (running_threshold / threshold);
1157
1158	double FSum = 0, HF = 0, TF = 0;
1159	const double resolution = 1.0 / 20.0;
1160
1161	int cen_y = qRound(center->y);
1162
1163	double rightEdge = center->x + center->width / 2.0;
1164	double leftEdge = center->x - center->width / 2.0;
1165
1166	QVector<double> subPixels;
1167	subPixels.reserve(center->width / resolution);
1168
1169	const T * origBuffer = reinterpret_cast<T *>(data->getImageBuffer()) + offset;
1170
1171	for (double x = leftEdge; x <= rightEdge; x += resolution)
1172	{
1173	double slice = resolution * (origBuffer[static_cast<int>(floor(x)) + cen_y * dataWidth]);
1174	FSum += slice;
1175	subPixels.append(slice);
1176	}
1177
1178	// Half flux
1179	HF = FSum / 2.0;
1180
1181	int subPixelCenter = (center->width / resolution) / 2;
1182
1183	// Start from center
1184	TF = subPixels[subPixelCenter];
1185	double lastTF = TF;
1186	// Integrate flux along radius axis until we reach half flux
1187	//for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
1188	for (int k = 1; k < subPixelCenter; k++)
1189	{
1190	TF += subPixels[subPixelCenter + k];
1191	TF += subPixels[subPixelCenter - k];
1192
1193	if (TF >= HF)
1194	{
1195	// We overpassed HF, let's calculate from last TF how much until we reach HF
1196
1197	// #1 Accurate calculation, but very sensitive to small variations of flux
1198	center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
1199
1200	// #2 Less accurate calculation, but stable against small variations of flux
1201	//center->HFR = (k - 1) * resolution;
1202	break;
1203	}
1204
1205	lastTF = TF;
1206	}
1207
1208	// Correct center for subX and subY
1209	center->x += subX;
1210	center->y += subY;
1211
1212	data->appendStar(center);
1213
1214	qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR;
1215
1216	return 1;
1217	}
1218
1219	int FITSData::findOneStar(const QRect &boundary)
1220	{
1221	switch (m_DataType)
1222	{
1223	case TBYTE:
1224	return findOneStar<uint8_t>(boundary);
1225
1226	case TSHORT:
1227	return findOneStar<int16_t>(boundary);
1228
1229	case TUSHORT:
1230	return findOneStar<uint16_t>(boundary);
1231
1232	case TLONG:
1233	return findOneStar<int32_t>(boundary);
1234
1235	case TULONG:
1236	return findOneStar<uint32_t>(boundary);
1237
1238	case TFLOAT:
1239	return findOneStar<float>(boundary);
1240
1241	case TLONGLONG:
1242	return findOneStar<int64_t>(boundary);
1243
1244	case TDOUBLE:
1245	return findOneStar<double>(boundary);
1246
1247	default:
1248	break;
1249	}
1250
1251	return 0;
1252	}
1253
1254	template <typename T>
1255	int FITSData::findOneStar(const QRect &boundary)
1256	{
1257	if (boundary.isEmpty())
1258	return -1;
1259
1260	int subX = boundary.x();
1261	int subY = boundary.y();
1262	int subW = subX + boundary.width();
1263	int subH = subY + boundary.height();
1264
1265	float massX = 0, massY = 0, totalMass = 0;
1266
1267	auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
1268
1269	// TODO replace magic number with something more useful to understand
1270	double threshold = stats.mean[0] * Options::focusThreshold() / 100.0;
1271
1272	for (int y = subY; y < subH; y++)
1273	{
1274	for (int x = subX; x < subW; x++)
1275	{
1276	T pixel = buffer[x + y * stats.width];
1277	if (pixel > threshold)
1278	{
1279	totalMass += pixel;
1280	massX += x * pixel;
1281	massY += y * pixel;
1282	}
1283	}
1284	}
1285
1286	qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << massX / totalMass << " Y: " << massY / totalMass;
1287
1288	auto * center = new Edge;
1289	center->width = -1;
1290	center->x = massX / totalMass + 0.5;
1291	center->y = massY / totalMass + 0.5;
1292	center->HFR = 1;
1293
1294	// Maximum Radius
1295	int maxR = qMin(subW - 1, subH - 1) / 2;
1296
1297	// Critical threshold
1298	double critical_threshold = threshold * 0.7;
1299	double running_threshold = threshold;
1300
1301	while (running_threshold >= critical_threshold)
1302	{
1303	for (int r = maxR; r > 1; r--)
1304	{
1305	int pass = 0;
1306
1307	for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0)
1308	{
1309	int testX = center->x + std::cos(theta) * r;
1310	int testY = center->y + std::sin(theta) * r;
1311
1312	// if out of bound, break;
1313	if (testX < subX \|\| testX > subW \|\| testY < subY \|\| testY > subH)
1314	break;
1315
1316	if (buffer[testX + testY * stats.width] > running_threshold)
1317	pass++;
1318	}
1319
1320	//qDebug() << "Testing for radius " << r << " passes # " << pass << " @ threshold " << running_threshold;
1321	//if (pass >= 6)
1322	if (pass >= 5)
1323	{
1324	center->width = r * 2;
1325	break;
1326	}
1327	}
1328
1329	if (center->width > 0)
1330	break;
1331
1332	// Increase threshold fuzziness by 10%
1333	running_threshold -= running_threshold * 0.1;
1334	}
1335
1336	// If no stars were detected
1337	if (center->width == -1)
1338	{
1339	delete center;
1340	return 0;
1341	}
1342
1343	// 30% fuzzy
1344	//center->width += center->width0.3 (running_threshold / threshold);
1345
1346	starCenters.append(center);
1347
1348	double FSum = 0, HF = 0, TF = 0, min = stats.min[0];
1349	const double resolution = 1.0 / 20.0;
1350
1351	int cen_y = qRound(center->y);
1352
1353	double rightEdge = center->x + center->width / 2.0;
1354	double leftEdge = center->x - center->width / 2.0;
1355
1356	QVector<double> subPixels;
1357	subPixels.reserve(center->width / resolution);
1358
1359	for (double x = leftEdge; x <= rightEdge; x += resolution)
1360	{
1361	//subPixels[x] = resolution * (image_buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
1362	double slice = resolution * (buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
1363	FSum += slice;
1364	subPixels.append(slice);
1365	}
1366
1367	// Half flux
1368	HF = FSum / 2.0;
1369
1370	//double subPixelCenter = center->x - fmod(center->x,resolution);
1371	int subPixelCenter = (center->width / resolution) / 2;
1372
1373	// Start from center
1374	TF = subPixels[subPixelCenter];
1375	double lastTF = TF;
1376	// Integrate flux along radius axis until we reach half flux
1377	//for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
1378	for (int k = 1; k < subPixelCenter; k++)
1379	{
1380	TF += subPixels[subPixelCenter + k];
1381	TF += subPixels[subPixelCenter - k];
1382
1383	if (TF >= HF)
1384	{
1385	// 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.
1386	// The second method is not truly HFR but is much more resistant to noise.
1387
1388	// #1 Approximate HFR, accurate and reliable but quite variable to small changes in star flux
1389	center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
1390
1391	// #2 Not exactly HFR, but much more stable
1392	//center->HFR = (kresolution) (HF/TF);
1393	break;
1394	}
1395
1396	lastTF = TF;
1397	}
1398
1399	return 1;
1400	}
1401
1402	/*** Find center of stars and calculate Half Flux Radius */
1403	int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth)
1404	{
1405	switch (m_DataType)
1406	{
1407	case TBYTE:
1408	return findCentroid<uint8_t>(boundary, initStdDev, minEdgeWidth);
1409
1410	case TSHORT:
1411	return findCentroid<int16_t>(boundary, initStdDev, minEdgeWidth);
1412
1413	case TUSHORT:
1414	return findCentroid<uint16_t>(boundary, initStdDev, minEdgeWidth);
1415
1416	case TLONG:
1417	return findCentroid<int32_t>(boundary, initStdDev, minEdgeWidth);
1418
1419	case TULONG:
1420	return findCentroid<uint32_t>(boundary, initStdDev, minEdgeWidth);
1421
1422	case TFLOAT:
1423	return findCentroid<float>(boundary, initStdDev, minEdgeWidth);
1424
1425	case TLONGLONG:
1426	return findCentroid<int64_t>(boundary, initStdDev, minEdgeWidth);
1427
1428	case TDOUBLE:
1429	return findCentroid<double>(boundary, initStdDev, minEdgeWidth);
1430
1431	default:
1432	return -1;
1433	}
1434	}
1435
1436	template <typename T>
1437	int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth)
1438	{
1439	double threshold = 0, sum = 0, avg = 0, min = 0;
1440	int starDiameter = 0;
1441	int pixVal = 0;
1442	int minimumEdgeCount = MINIMUM_EDGE_LIMIT;
1443
1444	auto * buffer = reinterpret_cast<T *>(m_ImageBuffer);
1445
1446	double JMIndex = 100;
1447	#ifndef KSTARS_LITE
1448	if (histogram)
1449	{
1450	if (!histogram->isConstructed())
1451	histogram->constructHistogram();
1452	JMIndex = histogram->getJMIndex();
1453	}
1454	#endif
1455
1456	float dispersion_ratio = 1.5;
1457
1458	QList<Edge *> edges;
1459
1460	if (JMIndex < DIFFUSE_THRESHOLD)
1461	{
1462	minEdgeWidth = JMIndex * 35 + 1;
1463	minimumEdgeCount = minEdgeWidth - 1;
1464	}
1465	else
1466	{
1467	minEdgeWidth = 6;
1468	minimumEdgeCount = 4;
1469	}
1470
1471	while (initStdDev >= 1)
1472	{
1473	minEdgeWidth--;
1474	minimumEdgeCount--;
1475
1476	minEdgeWidth = qMax(3, minEdgeWidth);
1477	minimumEdgeCount = qMax(3, minimumEdgeCount);
1478
1479	if (JMIndex < DIFFUSE_THRESHOLD)
1480	{
1481	// Taking the average out seems to have better result for noisy images
1482	threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1);
1483
1484	min = stats.min[0];
1485	if (threshold - min < 0)
1486	{
1487	threshold = stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1);
1488	min = 0;
1489	}
1490
1491	dispersion_ratio = 1.4 - (MINIMUM_STDVAR - initStdDev) * 0.08;
1492	}
1493	else
1494	{
1495	threshold = stats.mean[0] + stats.stddev[0] * initStdDev * (0.3 - (MINIMUM_STDVAR - initStdDev) * 0.05);
1496	min = stats.min[0];
1497	// Ratio between centeroid center and edge
1498	dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2;
1499	}
1500
1501	qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR;
1502	qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min
1503	<< ") minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount;
1504
1505	threshold -= min;
1506
1507	int subX, subY, subW, subH;
1508
1509	if (boundary.isNull())
1510	{
1511	if (m_Mode == FITS_GUIDE \|\| m_Mode == FITS_FOCUS)
1512	{
1513	// Only consider the central 70%
1514	subX = round(stats.width * 0.15);
1515	subY = round(stats.height * 0.15);
1516	subW = stats.width - subX;
1517	subH = stats.height - subY;
1518	}
1519	else
1520	{
1521	// Consider the complete area 100%
1522	subX = 0;
1523	subY = 0;
1524	subW = stats.width;
1525	subH = stats.height;
1526	}
1527	}
1528	else
1529	{
1530	subX = boundary.x();
1531	subY = boundary.y();
1532	subW = subX + boundary.width();
1533	subH = subY + boundary.height();
1534	}
1535
1536	// Detect "edges" that are above threshold
1537	for (int i = subY; i < subH; i++)
1538	{
1539	starDiameter = 0;
1540
1541	for (int j = subX; j < subW; j++)
1542	{
1543	pixVal = buffer[j + (i * stats.width)] - min;
1544
1545	// If pixel value > threshold, let's get its weighted average
1546	if (pixVal >= threshold)
1547	{
1548	avg += j * pixVal;
1549	sum += pixVal;
1550	starDiameter++;
1551	}
1552	// Value < threshold but avg exists
1553	else if (sum > 0)
1554	{
1555	// We found a potential centroid edge
1556	if (starDiameter >= minEdgeWidth)
1557	{
1558	float center = avg / sum + 0.5;
1559	if (center > 0)
1560	{
1561	int i_center = std::floor(center);
1562
1563	// Check if center is 10% or more brighter than edge, if not skip
1564	if (((buffer[i_center + (i * stats.width)] - min) /
1565	(buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >=
1566	dispersion_ratio) &&
1567	((buffer[i_center + (i * stats.width)] - min) /
1568	(buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >=
1569	dispersion_ratio))
1570	{
1571	qCDebug(KSTARS_FITS)
1572	<< "Edge center is " << buffer[i_center + (i * stats.width)] - min
1573	<< " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min
1574	<< " and ratio is "
1575	<< ((buffer[i_center + (i * stats.width)] - min) /
1576	(buffer[i_center + (i * stats.width) - starDiameter / 2] - min))
1577	<< " located at X: " << center << " Y: " << i + 0.5;
1578
1579	auto * newEdge = new Edge();
1580
1581	newEdge->x = center;
1582	newEdge->y = i + 0.5;
1583	newEdge->scanned = 0;
1584	newEdge->val = buffer[i_center + (i * stats.width)] - min;
1585	newEdge->width = starDiameter;
1586	newEdge->HFR = 0;
1587	newEdge->sum = sum;
1588
1589	edges.append(newEdge);
1590	}
1591	}
1592	}
1593
1594	// Reset
1595	avg = sum = starDiameter = 0;
1596	}
1597	}
1598	}
1599
1600	qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count();
1601
1602	// In case of hot pixels
1603	if (edges.count() == 1 && initStdDev > 1)
1604	{
1605	initStdDev--;
1606	continue;
1607	}
1608
1609	if (edges.count() >= MAX_EDGE_LIMIT)
1610	{
1611	qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count();
1612	qDeleteAll(edges);
1613	return -1;
1614	}
1615
1616	if (edges.count() >= minimumEdgeCount)
1617	break;
1618
1619	qDeleteAll(edges);
1620	edges.clear();
1621	initStdDev--;
1622	}
1623
1624	int cen_count = 0;
1625	int cen_x = 0;
1626	int cen_y = 0;
1627	int cen_v = 0;
1628	int cen_w = 0;
1629	int width_sum = 0;
1630
1631	// Let's sort edges, starting with widest
1632	std::sort(edges.begin(), edges.end(), greaterThan);
1633
1634	// Now, let's scan the edges and find the maximum centroid vertically
1635	for (int i = 0; i < edges.count(); i++)
1636	{
1637	qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of "
1638	<< edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum;
1639
1640	// If edge scanned already, skip
1641	if (edges[i]->scanned == 1)
1642	{
1643	qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted";
1644	continue;
1645	}
1646
1647	qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ...";
1648
1649	// Get X, Y, and Val of edge
1650	cen_x = edges[i]->x;
1651	cen_y = edges[i]->y;
1652	cen_v = edges[i]->sum;
1653	cen_w = edges[i]->width;
1654
1655	float avg_x = 0;
1656	float avg_y = 0;
1657
1658	sum = 0;
1659	cen_count = 0;
1660
1661	// Now let's compare to other edges until we hit a maxima
1662	for (int j = 0; j < edges.count(); j++)
1663	{
1664	if (edges[j]->scanned)
1665	continue;
1666
1667	if (checkCollision(edges[j], edges[i]))
1668	{
1669	if (edges[j]->sum >= cen_v)
1670	{
1671	cen_v = edges[j]->sum;
1672	cen_w = edges[j]->width;
1673	}
1674
1675	edges[j]->scanned = 1;
1676	cen_count++;
1677
1678	avg_x += edges[j]->x * edges[j]->val;
1679	avg_y += edges[j]->y * edges[j]->val;
1680	sum += edges[j]->val;
1681
1682	continue;
1683	}
1684	}
1685
1686	int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev));
1687
1688	if (edges.count() < LOW_EDGE_CUTOFF_1)
1689	{
1690	if (edges.count() < LOW_EDGE_CUTOFF_2)
1691	cen_limit = 1;
1692	else
1693	cen_limit = 2;
1694	}
1695
1696	qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is "
1697	<< cen_limit;
1698
1699	if (cen_limit < 1)
1700	continue;
1701
1702	// If centroid count is within acceptable range
1703	//if (cen_limit >= 2 && cen_count >= cen_limit)
1704	if (cen_count >= cen_limit)
1705	{
1706	// We detected a centroid, let's init it
1707	auto * rCenter = new Edge();
1708
1709	rCenter->x = avg_x / sum;
1710	rCenter->y = avg_y / sum;
1711	width_sum += rCenter->width;
1712	rCenter->width = cen_w;
1713
1714	qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")";
1715
1716	// Calculate Total Flux From Center, Half Flux, Full Summation
1717	double TF = 0;
1718	double HF = 0;
1719	double FSum = 0;
1720
1721	cen_x = (int)std::floor(rCenter->x);
1722	cen_y = (int)std::floor(rCenter->y);
1723
1724	if (cen_x < 0 \|\| cen_x > stats.width \|\| cen_y < 0 \|\| cen_y > stats.height)
1725	{
1726	delete rCenter;
1727	continue;
1728	}
1729
1730	// Complete sum along the radius
1731	//for (int k=0; k < rCenter->width; k++)
1732	for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--)
1733	{
1734	FSum += buffer[cen_x - k + (cen_y * stats.width)] - min;
1735	//qDebug() << image_buffer[cen_x-k+(cen_y*stats.width)] - min;
1736	}
1737
1738	// Half flux
1739	HF = FSum / 2.0;
1740
1741	// Total flux starting from center
1742	TF = buffer[cen_y * stats.width + cen_x] - min;
1743
1744	int pixelCounter = 1;
1745
1746	// Integrate flux along radius axis until we reach half flux
1747	for (int k = 1; k < rCenter->width / 2; k++)
1748	{
1749	if (TF >= HF)
1750	{
1751	qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels.";
1752	break;
1753	}
1754
1755	TF += buffer[cen_y * stats.width + cen_x + k] - min;
1756	TF += buffer[cen_y * stats.width + cen_x - k] - min;
1757
1758	pixelCounter++;
1759	}
1760
1761	// Calculate weighted Half Flux Radius
1762	rCenter->HFR = pixelCounter * (HF / TF);
1763	// Store full flux
1764	rCenter->val = FSum;
1765
1766	qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum;
1767
1768	starCenters.append(rCenter);
1769	}
1770	}
1771
1772	if (starCenters.count() > 1 && m_Mode != FITS_FOCUS)
1773	{
1774	float width_avg = (float)width_sum / starCenters.count();
1775	float lsum = 0, sdev = 0;
1776
1777	for (auto &center : starCenters)
1778	lsum += (center->width - width_avg) * (center->width - width_avg);
1779
1780	sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4;
1781
1782	// Reject stars > 4 * stddev
1783	foreach (Edge * center, starCenters)
1784	if (center->width > sdev)
1785	starCenters.removeOne(center);
1786
1787	//foreach(Edge *center, starCenters)
1788	//qDebug() << center->x << "," << center->y << "," << center->width << "," << center->val << endl;
1789	}
1790
1791	// Release memory
1792	qDeleteAll(edges);
1793
1794	return starCenters.count();
1795	}
1796
1797	double FITSData::getHFR(HFRType type)	932		double FITSData::getHFR(HFRType type)
1798	{	933		{
1799	// This method is less susceptible to noise	934		// This method is less susceptible to noise
1800	// Get HFR for the brightest star only, instead of averaging all stars	935		// Get HFR for the brightest star only, instead of averaging all stars
1801	// It is more consistent.	936		// It is more consistent.
1802	// TODO: Try to test this under using a real CCD.	937		// TODO: Try to test this under using a real CCD.
1803		938
1804	if (starCenters.empty())	939		if (starCenters.empty())
▲ Show 20 Lines • Show All 493 Lines • ▼ Show 20 Line(s)		1430	{
2298	if(subFrame.contains(x, y))	1433		if(subFrame.contains(x, y))
2299	{	1434		{
2300	starCentersInSubFrame.append(starCenters[i]);	1435		starCentersInSubFrame.append(starCenters[i]);
2301	}	1436		}
2302	}	1437		}
2303	return starCentersInSubFrame;	1438		return starCentersInSubFrame;
2304	}	1439		}
2305		1440
2306	void FITSData::getCenterSelection(int * x, int * y)
2307	{
2308	if (starCenters.count() == 0)
2309	return;
2310
2311	auto * pEdge = new Edge();
2312	pEdge->x = *x;
2313	pEdge->y = *y;
2314	pEdge->width = 1;
2315
2316	foreach (Edge * center, starCenters)
2317	if (checkCollision(pEdge, center))
2318	{
2319	*x = static_cast<int>(center->x);
2320	*y = static_cast<int>(center->y);
2321	break;
2322	}
2323
2324	delete (pEdge);
2325	}
2326
2327	bool FITSData::checkForWCS()	1441		bool FITSData::checkForWCS()
2328	{	1442		{
2329	#ifndef KSTARS_LITE	1443		#ifndef KSTARS_LITE
2330	#ifdef HAVE_WCSLIB	1444		#ifdef HAVE_WCSLIB
2331		1445
2332	int status = 0;	1446		int status = 0;
2333	char * header;	1447		char * header;
2334	int nkeyrec, nreject;	1448		int nkeyrec, nreject;
▲ Show 20 Lines • Show All 920 Lines • ▼ Show 20 Line(s)		2359	{
3255	{	2369		{
3256	sprintf(keyword, "CO2_%d", i);	2370		sprintf(keyword, "CO2_%d", i);
3257	fits_delete_key(fptr, keyword, &status);	2371		fits_delete_key(fptr, keyword, &status);
3258	}	2372		}
3259	}	2373		}
3260		2374
3261	}	2375		}
3262		2376
3263	uint8_t * FITSData::getImageBuffer()	2377		uint8_t * FITSData::getWritableImageBuffer()
		2378		{
		2379		return m_ImageBuffer;
		2380		}
		2381
		2382		uint8_t const * FITSData::getImageBuffer() const
3264	{	2383		{
3265	return m_ImageBuffer;	2384		return m_ImageBuffer;
3266	}	2385		}
3267		2386
3268	void FITSData::setImageBuffer(uint8_t * buffer)	2387		void FITSData::setImageBuffer(uint8_t * buffer)
3269	{	2388		{
3270	delete[] m_ImageBuffer;	2389		delete[] m_ImageBuffer;
3271	m_ImageBuffer = buffer;	2390		m_ImageBuffer = buffer;
▲ Show 20 Lines • Show All 266 Lines • ▼ Show 20 Line(s)
3538	{	2657		{
3539	double adu = 0;	2658		double adu = 0;
3540	for (int i = 0; i < m_Channels; i++)	2659		for (int i = 0; i < m_Channels; i++)
3541	adu += stats.mean[i];	2660		adu += stats.mean[i];
3542		2661
3543	return (adu / static_cast<double>(m_Channels));	2662		return (adu / static_cast<double>(m_Channels));
3544	}	2663		}
3545		2664
3546	/* CannyDetector, Implementation of Canny edge detector in Qt/C++.
3547	* Copyright (C) 2015 Gonzalo Exequiel Pedone
3548	*
3549	* This program is free software: you can redistribute it and/or modify
3550	* it under the terms of the GNU General Public License as published by
3551	* the Free Software Foundation, either version 3 of the License, or
3552	* (at your option) any later version.
3553	*
3554	* This program is distributed in the hope that it will be useful,
3555	* but WITHOUT ANY WARRANTY; without even the implied warranty of
3556	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3557	* GNU General Public License for more details.
3558	*
3559	* You should have received a copy of the GNU General Public License
3560	* along with this program. If not, see <http://www.gnu.org/licenses/>.
3561	*
3562	* Email : hipersayan DOT x AT gmail DOT com
3563	* Web-Site: https://github.com/hipersayanX/CannyDetector
3564	*/
3565
3566	template <typename T>
3567	void FITSData::sobel(QVector<float> &gradient, QVector<float> &direction)
3568	{
3569	//int size = image.width() * image.height();
3570	gradient.resize(stats.samples_per_channel);
3571	direction.resize(stats.samples_per_channel);
3572
3573	for (int y = 0; y < stats.height; y++)
3574	{
3575	size_t yOffset = y * stats.width;
3576	const T * grayLine = reinterpret_cast<T *>(m_ImageBuffer) + yOffset;
3577
3578	const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width;
3579	const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width;
3580
3581	float * gradientLine = gradient.data() + yOffset;
3582	float * directionLine = direction.data() + yOffset;
3583
3584	for (int x = 0; x < stats.width; x++)
3585	{
3586	int x_m1 = x < 1 ? x : x - 1;
3587	int x_p1 = x >= stats.width - 1 ? x : x + 1;
3588
3589	int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] -
3590	2 * grayLine[x_m1] - grayLine_p1[x_m1];
3591
3592	int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] -
3593	2 * grayLine_p1[x] - grayLine_p1[x_p1];
3594
3595	gradientLine[x] = qAbs(gradX) + qAbs(gradY);
3596
3597	/* Gradient directions are classified in 4 possible cases
3598	*
3599	* dir 0
3600	*
3601	* x x x
3602	* - - -
3603	* x x x
3604	*
3605	* dir 1
3606	*
3607	* x x /
3608	* x / x
3609	* / x x
3610	*
3611	* dir 2
3612	*
3613	* \ x x
3614	* x \ x
3615	* x x \
3616	*
3617	* dir 3
3618	*
3619	* x \| x
3620	* x \| x
3621	* x \| x
3622	*/
3623	if (gradX == 0 && gradY == 0)
3624	directionLine[x] = 0;
3625	else if (gradX == 0)
3626	directionLine[x] = 3;
3627	else
3628	{
3629	qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
3630
3631	if (a >= -22.5 && a < 22.5)
3632	directionLine[x] = 0;
3633	else if (a >= 22.5 && a < 67.5)
3634	directionLine[x] = 2;
3635	else if (a >= -67.5 && a < -22.5)
3636	directionLine[x] = 1;
3637	else
3638	directionLine[x] = 3;
3639	}
3640	}
3641	}
3642	}
3643
3644	int FITSData::partition(int width, int height, QVector<float> &gradient, QVector<int> &ids)
3645	{
3646	int id = 0;
3647
3648	for (int y = 1; y < height - 1; y++)
3649	{
3650	for (int x = 1; x < width - 1; x++)
3651	{
3652	int index = x + y * width;
3653	float val = gradient[index];
3654	if (val > 0 && ids[index] == 0)
3655	{
3656	trace(width, height, ++id, gradient, ids, x, y);
3657	}
3658	}
3659	}
3660
3661	// Return max id
3662	return id;
3663	}
3664
3665	void FITSData::trace(int width, int height, int id, QVector<float> &image, QVector<int> &ids, int x, int y)
3666	{
3667	int yOffset = y * width;
3668	float * cannyLine = image.data() + yOffset;
3669	int * idLine = ids.data() + yOffset;
3670
3671	if (idLine[x] != 0)
3672	return;
3673
3674	idLine[x] = id;
3675
3676	for (int j = -1; j < 2; j++)
3677	{
3678	int nextY = y + j;
3679
3680	if (nextY < 0 \|\| nextY >= height)
3681	continue;
3682
3683	float * cannyLineNext = cannyLine + j * width;
3684
3685	for (int i = -1; i < 2; i++)
3686	{
3687	int nextX = x + i;
3688
3689	if (i == j \|\| nextX < 0 \|\| nextX >= width)
3690	continue;
3691
3692	if (cannyLineNext[nextX] > 0)
3693	{
3694	// Trace neighbors.
3695	trace(width, height, id, image, ids, nextX, nextY);
3696	}
3697	}
3698	}
3699	}
3700
3701	QString FITSData::getLastError() const	2665		QString FITSData::getLastError() const
3702	{	2666		{
3703	return lastError;	2667		return lastError;
3704	}	2668		}
3705		2669
3706	bool FITSData::getAutoRemoveTemporaryFITS() const	2670		bool FITSData::getAutoRemoveTemporaryFITS() const
3707	{	2671		{
3708	return autoRemoveTemporaryFITS;	2672		return autoRemoveTemporaryFITS;
▲ Show 20 Lines • Show All 507 Lines • ▼ Show 20 Line(s)		3128	{
4216	return true;	3180		return true;
4217	}	3181		}
4218		3182
4219	bool FITSData::contains(const QPointF &point) const	3183		bool FITSData::contains(const QPointF &point) const
4220	{	3184		{
4221	return (point.x() >= 0 && point.y() >= 0 && point.x() <= stats.width && point.y() <= stats.height);	3185		return (point.x() >= 0 && point.y() >= 0 && point.x() <= stats.width && point.y() <= stats.height);
4222	}	3186		}
4223		3187
4224	int FITSData::findSEPStars(const QRect &boundary)
4225	{
4226	int x = 0, y = 0, w = stats.width, h = stats.height, maxRadius = 50;
4227
4228	if (!boundary.isNull())
4229	{
4230	x = boundary.x();
4231	y = boundary.y();
4232	w = boundary.width();
4233	h = boundary.height();
4234	maxRadius = w;
4235	}
4236
4237	auto * data = new float[w * h];
4238
4239	switch (stats.bitpix)
4240	{
4241	case BYTE_IMG:
4242	getFloatBuffer<uint8_t>(data, x, y, w, h);
4243	break;
4244	case SHORT_IMG:
4245	getFloatBuffer<int16_t>(data, x, y, w, h);
4246	break;
4247	case USHORT_IMG:
4248	getFloatBuffer<uint16_t>(data, x, y, w, h);
4249	break;
4250	case LONG_IMG:
4251	getFloatBuffer<int32_t>(data, x, y, w, h);
4252	break;
4253	case ULONG_IMG:
4254	getFloatBuffer<uint32_t>(data, x, y, w, h);
4255	break;
4256	case FLOAT_IMG:
4257	delete [] data;
4258	data = reinterpret_cast<float *>(m_ImageBuffer);
4259	break;
4260	case LONGLONG_IMG:
4261	getFloatBuffer<int64_t>(data, x, y, w, h);
4262	break;
4263	case DOUBLE_IMG:
4264	getFloatBuffer<double>(data, x, y, w, h);
4265	break;
4266	default:
4267	delete [] data;
4268	return -1;
4269	}
4270
4271	float * imback = nullptr;
4272	double * flux = nullptr, fluxerr = nullptr, area = nullptr;
4273	short * flag = nullptr;
4274	short flux_flag = 0;
4275	int status = 0;
4276	sep_bkg * bkg = nullptr;
4277	sep_catalog * catalog = nullptr;
4278	float conv[] = {1, 2, 1, 2, 4, 2, 1, 2, 1};
4279	double flux_fractions[2] = {0};
4280	double requested_frac[2] = { 0.5, 0.99 };
4281	QList<Edge *> edges;
4282
4283	// #0 Create SEP Image structure
4284	sep_image im = {data, nullptr, nullptr, SEP_TFLOAT, 0, 0, w, h, 0.0, SEP_NOISE_NONE, 1.0, 0.0};
4285
4286	// #1 Background estimate
4287	status = sep_background(&im, 64, 64, 3, 3, 0.0, &bkg);
4288	if (status != 0) goto exit;
4289
4290	// #2 Background evaluation
4291	imback = (float )malloc((w h) * sizeof(float));
4292	status = sep_bkg_array(bkg, imback, SEP_TFLOAT);
4293	if (status != 0) goto exit;
4294
4295	// #3 Background subtraction
4296	status = sep_bkg_subarray(bkg, im.data, im.dtype);
4297	if (status != 0) goto exit;
4298
4299	// #4 Source Extraction
4300	// Note that we set deblend_cont = 1.0 to turn off deblending.
4301	status = sep_extract(&im, 2 * bkg->globalrms, SEP_THRESH_ABS, 10, conv, 3, 3, SEP_FILTER_CONV, 32, 1.0, 1, 1.0, &catalog);
4302	if (status != 0) goto exit;
4303
4304	// TODO
4305	// Must detect edge detection
4306	// Must limit to brightest 100 (by flux) centers
4307	// Should probably use ellipse to draw instead of simple circle?
4308	// Useful for galaxies and also elenogated stars.
4309	for (int i = 0; i < catalog->nobj; i++)
4310	{
4311	double flux = catalog->flux[i];
4312	// Get HFR
4313	sep_flux_radius(&im, catalog->x[i], catalog->y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag);
4314
4315	auto * center = new Edge();
4316	center->x = catalog->x[i] + x + 0.5;
4317	center->y = catalog->y[i] + y + 0.5;
4318	center->val = catalog->peak[i];
4319	center->sum = flux;
4320	center->HFR = center->width = flux_fractions[0];
4321	if (flux_fractions[1] < maxRadius)
4322	center->width = flux_fractions[1] * 2;
4323	edges.append(center);
4324	}
4325
4326	// Let's sort edges, starting with widest
4327	std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->width > edge2->width;});
4328
4329	// Take only the first 100 stars
4330	{
4331	int starCount = qMin(100, edges.count());
4332	for (int i = 0; i < starCount; i++)
4333	starCenters.append(edges[i]);
4334	}
4335
4336	edges.clear();
4337
4338	qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << "#" << "#X" << "#Y" << "#Flux" << "#Width" << "#HFR";
4339	for (int i = 0; i < starCenters.count(); i++)
4340	qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << i << starCenters[i]->x << starCenters[i]->y
4341	<< starCenters[i]->sum << starCenters[i]->width << starCenters[i]->HFR;
4342
4343	exit:
4344	if (stats.bitpix != FLOAT_IMG)
4345	delete [] data;
4346	sep_bkg_free(bkg);
4347	sep_catalog_free(catalog);
4348	free(imback);
4349	free(flux);
4350	free(fluxerr);
4351	free(area);
4352	free(flag);
4353
4354	if (status != 0)
4355	{
4356	char errorMessage[512];
4357	sep_get_errmsg(status, errorMessage);
4358	qCritical(KSTARS_FITS) << errorMessage;
4359	return -1;
4360	}
4361
4362	return starCenters.count();
4363	}
4364
4365	template <typename T>
4366	void FITSData::getFloatBuffer(float * buffer, int x, int y, int w, int h)
4367	{
4368	auto * rawBuffer = reinterpret_cast<T *>(m_ImageBuffer);
4369
4370	float * floatPtr = buffer;
4371
4372	int x2 = x + w;
4373	int y2 = y + h;
4374
4375	for (int y1 = y; y1 < y2; y1++)
4376	{
4377	int offset = y1 * stats.width;
4378	for (int x1 = x; x1 < x2; x1++)
4379	{
4380	*floatPtr++ = rawBuffer[offset + x1];
4381	}
4382	}
4383	}
4384
4385	void FITSData::saveStatistics(Statistic &other)	3188		void FITSData::saveStatistics(Statistic &other)
4386	{	3189		{
4387	other = stats;	3190		other = stats;
4388	}	3191		}
4389		3192
4390	void FITSData::restoreStatistics(Statistic &other)	3193		void FITSData::restoreStatistics(Statistic &other)
4391	{	3194		{
4392	stats = other;	3195		stats = other;
4393	}	3196		}