Differential D29364 Diff 81746 kstars/fitsviewer/fitsbahtinovdetector.cpp

Changeset View

Standalone View

kstars/fitsviewer/fitsbahtinovdetector.cpp

This file was added.

		1		/***************************************************************************
		2		fitsbahtinovdetector.cpp - FITS Image
		3		-------------------
		4		begin : Wed April 15 2020
		5		copyright : (C) 2020 by Patrick Molenaar
		6		email : pr_molenaar@hotmail.com
		7		***************************************************************************/
		8
		9		/***************************************************************************
		10		* *
		11		* This program is free software; you can redistribute it and/or modify *
		12		* it under the terms of the GNU General Public License as published by *
		13		* the Free Software Foundation; either version 2 of the License, or *
		14		* (at your option) any later version. *
		15		* *
		16		* Some code fragments were adapted from Peter Kirchgessner's FITS plugin*
		17		* See http://members.aol.com/pkirchg for more details. *
		18		***************************************************************************/
		19
		20		#include <math.h>
		21		#include <cmath>
		22
		23		#include "fits_debug.h"
		24		#include "fitsbahtinovdetector.h"
		25		#include "hough/houghline.h"
		26
		27		#include <QElapsedTimer>
		28
		29		FITSStarDetector& FITSBahtinovDetector::configure(const QString &setting, const QVariant &value)
		30		{
		31		if (!setting.compare("NUMBER_OF_AVERAGE_ROWS", Qt::CaseInsensitive))
		32		{
		33		if (value.canConvert <typeof (NUMBER_OF_AVERAGE_ROWS)> ())
		34		{
		35		NUMBER_OF_AVERAGE_ROWS = value.value <typeof (NUMBER_OF_AVERAGE_ROWS)> ();
		36
		37		// Validate number of average rows value
		38		if (NUMBER_OF_AVERAGE_ROWS % 2 == 0)
		39		{
		40		NUMBER_OF_AVERAGE_ROWS--;
		41		qCInfo(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to "
		42		<< NUMBER_OF_AVERAGE_ROWS;
		43		}
		44		// Rows must be a positive number!
		45		if (NUMBER_OF_AVERAGE_ROWS < 1)
		46		{
		47		NUMBER_OF_AVERAGE_ROWS = 1;
		48		qCInfo(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to "
		49		<< NUMBER_OF_AVERAGE_ROWS;
		50		}
		51		}
		52		}
		53		return *this;
		54		}
		55
		56		int FITSBahtinovDetector::findSources(QList<Edge*> &starCenters, QRect const &boundary)
		57		{
		58		switch (parent()->property("dataType").toInt())
		59		{
		60		case TBYTE:
		61		return findBahtinovStar<uint8_t>(starCenters, boundary);
		62
		63		case TSHORT:
		64		return findBahtinovStar<int16_t>(starCenters, boundary);
		65
		66		case TUSHORT:
		67		return findBahtinovStar<uint16_t>(starCenters, boundary);
		68
		69		case TLONG:
		70		return findBahtinovStar<int32_t>(starCenters, boundary);
		71
		72		case TULONG:
		73		return findBahtinovStar<uint32_t>(starCenters, boundary);
		74
		75		case TFLOAT:
		76		return findBahtinovStar<float>(starCenters, boundary);
		77
		78		case TLONGLONG:
		79		return findBahtinovStar<int64_t>(starCenters, boundary);
		80
		81		case TDOUBLE:
		82		return findBahtinovStar<double>(starCenters, boundary);
		83
		84		default:
		85		break;
		86		}
		87
		88		return 0;
		89		}
		90
		91		template <typename T>
		92		int FITSBahtinovDetector::findBahtinovStar(QList<Edge*> &starCenters, const QRect &boundary)
		93		{
		94		FITSData const * const image_data = reinterpret_cast<FITSData const *>(parent());
		95
		96		if (image_data == nullptr)
		97		return 0;
		98
		99		if (boundary.isEmpty())
		100		return -1;
		101
		102		int subX = qMax(0, boundary.isNull() ? 0 : boundary.x());
		103		int subY = qMax(0, boundary.isNull() ? 0 : boundary.y());
		104		int subW = (boundary.isNull() ? image_data->width() : boundary.width());
		105		int subH = (boundary.isNull() ? image_data->height() : boundary.height());
		106
		107		int BBP = image_data->getBytesPerPixel();
		108		uint16_t dataWidth = image_data->width();
		109
		110		// #1 Find offsets
		111		uint32_t size = subW * subH;
		112		uint32_t offset = subX + subY * dataWidth;
		113
		114		// #2 Create new buffer
		115		auto * buffer = new uint8_t[size * BBP];
		116		// If there is no offset, copy whole buffer in one go
		117		if (offset == 0)
		118		{
		119		memcpy(buffer, image_data->getImageBuffer(), size * BBP);
		120		}
		121		else
		122		{
		123		uint8_t * dataPtr = buffer;
		124		const uint8_t * origDataPtr = image_data->getImageBuffer();
		125		uint32_t lineOffset = 0;
		126		// Copy data line by line
		127		for (int height = subY; height < (subY + subH); height++)
		128		{
		129		lineOffset = (subX + height * dataWidth) * BBP;
		130		memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
		131		dataPtr += (subW * BBP);
		132		}
		133		}
		134
		135		// #3 Create new FITSData to hold it
		136		FITSData::Statistic stats;
		137		stats.width = subW;
		138		stats.height = subH;
		139		stats.bitpix = image_data->getStatistics().bitpix;
		140		stats.bytesPerPixel = BBP;
		141		stats.samples_per_channel = size;
		142		FITSData* boundedImage = new FITSData();
		143		boundedImage->restoreStatistics(stats);
		144		boundedImage->setProperty("dataType", image_data->property("dataType"));
		145
		146		// #4 Set image buffer
		147		boundedImage->setImageBuffer(buffer);
		148
		149		boundedImage->calculateStats(true);
		150
		151		QElapsedTimer timer1;
		152
		153		// Rotate image 180 degrees in steps of 1 degree
		154		QMap<int, BahtinovLineAverage> lineAveragesPerAngle;
		155		int steps = 180;
		156		double radPerStep = M_PI / (double)steps;
		157
		158		timer1.start();
		159
		160		for (int angle = 0; angle < steps; angle++)
		161		{
		162		// TODO Apply multi threading to speed up calculation
		163		BahtinovLineAverage lineAverage = calculateMaxAverage<T>(boundedImage, angle);
		164		// Store line average in map
		165		lineAveragesPerAngle.insert(angle, lineAverage);
		166		}
		167
		168		qCDebug(KSTARS_FITS) << "Getting max average for all 180 rotations took" << timer1.elapsed() << "milliseconds";
		169
		170		// Not needed anymore
		171		delete boundedImage;
		172
		173		// Calculate Bahtinov angles
		174		QVector<HoughLine*> bahtinov_angles;
		175
		176		// For all three Bahtinov angles
		177		for (int index1 = 0; index1 < 3; index1++)
		178		{
		179		double maxAverage = 0.0;
		180		double maxAngle = 0.0;
		181		int maxAverageOffset = 0;
		182		for (int angle = 0; angle < steps; angle++)
		183		{
		184		BahtinovLineAverage lineAverage = lineAveragesPerAngle[angle];
		185		if (lineAverage.average > maxAverage)
		186		{
		187		maxAverage = lineAverage.average;
		188		maxAverageOffset = lineAverage.offset;
		189		maxAngle = angle;
		190		}
		191		}
		192		HoughLine* pHoughLine = new HoughLine(maxAngle * radPerStep, maxAverageOffset, subW, subH, (int)maxAverage);
		193		if (pHoughLine != nullptr)
		194		{
		195		bahtinov_angles.append(pHoughLine);
		196		}
		197
		198		// Remove data around peak to prevent it from being detected again
		199		int minBahtinovAngleOffset = 18;
		200		for (int subAngle = maxAngle - minBahtinovAngleOffset; subAngle < maxAngle + minBahtinovAngleOffset; subAngle++)
		201		{
		202		int angleInRange = subAngle;
		203		if (angleInRange < 0)
		204		{
		205		angleInRange += 180;
		206		}
		207		if (angleInRange >= 180)
		208		{
		209		angleInRange -= 180;
		210		}
		211		lineAveragesPerAngle.remove(angleInRange);
		212		}
		213		}
		214
		215		// Proceed with focus offset calculation, but only when at least 3 lines have been detected
		216		QVector<HoughLine*> top3Lines;
		217		if (bahtinov_angles.size() >= 3)
		218		{
		219		HoughLine::getSortedTopThreeLines(bahtinov_angles, top3Lines);
		220
		221		// Debug output
		222		qCDebug(KSTARS_FITS) << "Sorted bahtinov angles:";
		223		foreach (HoughLine* ln, top3Lines)
		224		{
		225		ln->printHoughLine();
		226		}
		227
		228		// Determine intersection between outer lines
		229		HoughLine* oneLine = top3Lines[0];
		230		HoughLine* otherLine = top3Lines[2];
		231		QPointF intersection;
		232		HoughLine::IntersectResult result = oneLine->Intersect(*otherLine, intersection);
		233		if (result == HoughLine::INTERESECTING) {
		234
		235		qCDebug(KSTARS_FITS) << "Intersection: " << intersection.x() << ", " << intersection.y();
		236
		237		// Determine offset between intersection and middle line
		238		HoughLine* midLine = top3Lines[1];
		239		QPointF intersectionOnMidLine;
		240		double distance;
		241		if (midLine->DistancePointLine(intersection, intersectionOnMidLine, distance)) {
		242		qCDebug(KSTARS_FITS) << "Distance between intersection and midline is " << distance
		243		<< " at mid line point " << intersectionOnMidLine.x() << ", "
		244		<< intersectionOnMidLine.y();
		245
		246		// Add star center to selected stars
		247		// Maximum Radius
		248		int maxR = qMin(subW - 1, subH - 1) / 2;
		249		BahtinovEdge* center = new BahtinovEdge();
		250		center->width = maxR / 3;
		251		center->x = subX + intersection.x();
		252		center->y = subY + intersection.y();
		253		// Set distance value in HFR
		254		center->HFR = distance;
		255
		256		center->offset.setX(subX + intersectionOnMidLine.x());
		257		center->offset.setY(subY + intersectionOnMidLine.y());
		258		oneLine->Offset(subX, subY);
		259		midLine->Offset(subX, subY);
		260		otherLine->Offset(subX, subY);
		261		center->line.append(*oneLine);
		262		center->line.append(*midLine);
		263		center->line.append(*otherLine);
		264		starCenters.append(center);
		265		}
		266		else
		267		{
		268		qCWarning(KSTARS_FITS) << "Closest point does not fall within the line segment.";
		269		}
		270		}
		271		else
		272		{
		273		qCWarning(KSTARS_FITS) << "Lines are not intersecting (result: " << result << ")";
		274		}
		275		}
		276
		277		// Clean up Bahtinov line array (of pointers) as they are no longer needed
		278		for (int index = 0; index < bahtinov_angles.size(); index++)
		279		{
		280		HoughLine* pLineAverage = bahtinov_angles[index];
		281		if (pLineAverage != nullptr) {
		282		delete pLineAverage;
		283		}
		284		}
		285		bahtinov_angles.clear();
		286
		287		// Clean up line averages array as they are no longer needed
		288		lineAveragesPerAngle.clear();
		289
		290		top3Lines.clear();
		291
		292		return 1;
		293		}
		294
		295		template <typename T>
		296		BahtinovLineAverage FITSBahtinovDetector::calculateMaxAverage(const FITSData *data, int angle)
		297		{
		298		int BBP = data->getBytesPerPixel();
		299		int size = data->getStatistics().samples_per_channel;
		300		int width = data->width();
		301		int height = data->height();
		302		int numChannels = data->channels();
		303
		304		BahtinovLineAverage lineAverage;
		305		auto * rotimage = new T[size * BBP];
		306
		307		rotateImage(data, angle, rotimage);
		308
		309		// Calculate average pixel value for each row
		310		auto * rotBuffer = reinterpret_cast<T *>(rotimage);
		311
		312		// printf("Angle;%d;Width;%d;Height;%d;Rows;%d;;RowSum;", angle, width, height, NUMBER_OF_AVERAGE_ROWS);
		313
		314		for (int y = 0; y < height; y++)
		315		{
		316		int yMin = y - ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
		317		int yMax = y + ((NUMBER_OF_AVERAGE_ROWS - 1) / 2);
		318
		319		unsigned long multiRowSum = 0;
		320		// Calculate average over multiple rows
		321		for (int y1 = yMin; y1 <= yMax; y1++)
		322		{
		323		int y2 = y1;
		324		if (y2 < 0) { y2 += height; }
		325		if (y2 >= height) { y2 -= height; }
		326		if (y2 < 0 \|\| y2 >= height) {
		327		qCWarning(KSTARS_FITS) << "Y still out of bounds: 0 <=" << y2 << "<" << height;
		328		}
		329
		330		for (int x1 = 0; x1 < width; x1++)
		331		{
		332		int index = y2 * width + x1;
		333		unsigned long channelAverage = 0;
		334		for (int i = 0; i < numChannels; i++)
		335		{
		336		int offset = size * i;
		337		channelAverage += rotBuffer[index + offset];
		338		}
		339		multiRowSum += qRound(channelAverage / (double)numChannels);
		340		}
		341		}
		342		// printf("%lu;", multiRowSum);
		343
		344		double average = multiRowSum / (double)(width * NUMBER_OF_AVERAGE_ROWS);
		345		if (average > lineAverage.average)
		346		{
		347		lineAverage.average = average;
		348		lineAverage.offset = y;
		349		}
		350		}
		351		// printf(";;MaxAverage;%.3f;MaxAverageIndex;%lu\r\n", lineAverage.average, lineAverage.offset);
		352		// fflush(stdout);
		353
		354		rotBuffer = nullptr;
		355		delete[] rotimage;
		356		rotimage = nullptr;
		357
		358		return lineAverage;
		359		}
		360
		361		/** Rotate an image by angle degrees.
		362		* verbose generates extra info on stdout.
		363		* return true if successful and rotated image.
		364		* @param angle The angle over which the image needs to be rotated
		365		* @param rotImage The image that needs to be rotated, also the rotated image return value
		366		*/
		367		template <typename T>
		368		bool FITSBahtinovDetector::rotateImage(const FITSData data, int angle, T rotImage)
		369		{
		370		int BBP = data->getBytesPerPixel();
		371		int size = data->getStatistics().samples_per_channel;
		372		int width = data->width();
		373		int height = data->height();
		374		int numChannels = data->channels();
		375
		376		int hx, hy;
		377		int offset = 0;
		378		size_t bufferSize;
		379
		380		/* Check allocation buffer for rotated image */
		381		if (rotImage == nullptr)
		382		{
		383		qWarning() << "No memory allocated for rotated image buffer!";
		384		return false;
		385		}
		386
		387		while (angle < 0)
		388		{
		389		angle = angle + 360;
		390		}
		391		while (angle >= 360)
		392		{
		393		angle = angle - 360;
		394		}
		395
		396		hx = qFloor((width + 1) / 2.0);
		397		hy = qFloor((height + 1) / 2.0);
		398
		399		bufferSize = size * numChannels * BBP;
		400		memset(rotImage, 0, bufferSize);
		401
		402		auto * rotBuffer = reinterpret_cast<T *>(rotImage);
		403		auto * buffer = reinterpret_cast<const T *>(data->getImageBuffer());
		404
		405		double innerCircleRadius = (0.5 * qSqrt(2.0) * qMin(hx, hy));
		406		double angleInRad = angle * M_PI / 180.0;
		407		double sinAngle = qSin(angleInRad);
		408		double cosAngle = qCos(angleInRad);
		409		int leftEdge = qCeil(hx - innerCircleRadius);
		410		int rightEdge = qFloor(hx + innerCircleRadius);
		411		int topEdge = qCeil(hy - innerCircleRadius);
		412		int bottomEdge = qFloor(hy + innerCircleRadius);
		413
		414		for (int i = 0; i < numChannels; i++)
		415		{
		416		offset = size * i;
		417		for (int x1 = leftEdge; x1 < rightEdge; x1++)
		418		{
		419		for (int y1 = topEdge; y1 < bottomEdge; y1++)
		420		{
		421		// translate point back to origin:
		422		double x2 = x1 - hx;
		423		double y2 = y1 - hy;
		424
		425		// rotate point
		426		double xnew = x2 * cosAngle - y2 * sinAngle;
		427		double ynew = x2 * sinAngle + y2 * cosAngle;
		428
		429		// translate point back:
		430		x2 = xnew + hx;
		431		y2 = ynew + hy;
		432
		433		int orgIndex = y1 * height + x1;
		434		int newIndex = qRound(y2) * height + qRound(x2);
		435
		436		if (newIndex >= 0 && newIndex < (int)(bufferSize))
		437		{
		438		rotBuffer[newIndex + offset] = buffer[orgIndex + offset];
		439		} // else index out of bounds, do not update pixel
		440		}
		441		}
		442		}
		443
		444		return true;
		445		}