Differential D28292 Diff 79624 kstars/fitsviewer/fitsgradientdetector.cpp

Changeset View

Standalone View

kstars/fitsviewer/fitsgradientdetector.cpp

This file was added.

		1		/***************************************************************************
		2		fitsgradientdetector.cpp - FITS Image
		3		-------------------
		4		begin : Sat March 28 2020
		5		copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet
		6		email : eric.dejouhanet@gmail.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
		22		#include "fits_debug.h"
		23		#include "fitsgradientdetector.h"
		24
		25		FITSStarDetector& FITSGradientDetector::configure(const QString &, const QVariant &)
		26		{
		27		return *this;
		28		}
		29
		30		int FITSGradientDetector::findSources(QList<Edge*> &starCenters, const QRect &boundary)
		31		{
		32		switch (parent()->property("dataType").toInt())
		33		{
		34		case TBYTE:
		35		return findSources<uint8_t>(starCenters, boundary);
		36
		37		case TSHORT:
		38		return findSources<int16_t>(starCenters, boundary);
		39
		40		case TUSHORT:
		41		return findSources<uint16_t>(starCenters, boundary);
		42
		43		case TLONG:
		44		return findSources<int32_t>(starCenters, boundary);
		45
		46		case TULONG:
		47		return findSources<uint16_t>(starCenters, boundary);
		48
		49		case TFLOAT:
		50		return findSources<float>(starCenters, boundary);
		51
		52		case TLONGLONG:
		53		return findSources<int64_t>(starCenters, boundary);
		54
		55		case TDOUBLE:
		56		return findSources<double>(starCenters, boundary);
		57
		58		default:
		59		break;
		60		}
		61
		62		return 0;
		63		}
		64
		65		template <typename T>
		66		int FITSGradientDetector::findSources(QList<Edge*> &starCenters, const QRect &boundary)
		67		{
		68		FITSData const * const data = reinterpret_cast<FITSData const *>(parent());
		69
		70		if (data == nullptr)
		71		return 0;
		72
		73		int subX = qMax(0, boundary.isNull() ? 0 : boundary.x());
		74		int subY = qMax(0, boundary.isNull() ? 0 : boundary.y());
		75		int subW = (boundary.isNull() ? data->width() : boundary.width());
		76		int subH = (boundary.isNull() ? data->height() : boundary.height());
		77
		78		int BBP = data->getBytesPerPixel();
		79
		80		uint16_t dataWidth = data->width();
		81
		82		// #1 Find offsets
		83		uint32_t size = subW * subH;
		84		uint32_t offset = subX + subY * dataWidth;
		85
		86		// #2 Create new buffer
		87		auto * buffer = new uint8_t[size * BBP];
		88		// If there is no offset, copy whole buffer in one go
		89		if (offset == 0)
		90		memcpy(buffer, data->getImageBuffer(), size * BBP);
		91		else
		92		{
		93		uint8_t * dataPtr = buffer;
		94		uint8_t const * origDataPtr = data->getImageBuffer();
		95		uint32_t lineOffset = 0;
		96		// Copy data line by line
		97		for (int height = subY; height < (subY + subH); height++)
		98		{
		99		lineOffset = (subX + height * dataWidth) * BBP;
		100		memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP);
		101		dataPtr += (subW * BBP);
		102		}
		103		}
		104
		105		// #3 Create new FITSData to hold it
		106		auto * boundedImage = new FITSData();
		107		FITSData::Statistic stats;
		108		stats.width = subW;
		109		stats.height = subH;
		110		stats.bitpix = data->getStatistics().bitpix;
		111		stats.bytesPerPixel = data->getStatistics().bytesPerPixel;
		112		stats.samples_per_channel = size;
		113		stats.ndim = 2;
		114		boundedImage->restoreStatistics(stats);
		115
		116		boundedImage->setProperty("dataType", parent()->property("dataType"));
		117
		118		// #4 Set image buffer and calculate stats.
		119		boundedImage->setImageBuffer(buffer);
		120
		121		boundedImage->calculateStats(true);
		122
		123		// #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain
		124		boundedImage->applyFilter(FITS_MEDIAN);
		125		boundedImage->applyFilter(FITS_HIGH_CONTRAST);
		126
		127		// #6 Perform Sobel to find gradients and their directions
		128		QVector<float> gradients;
		129		QVector<float> directions;
		130
		131		// TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed
		132		// and discarded whenever necessary. It won't work on noisy images unless this is done.
		133		sobel<T>(boundedImage, gradients, directions);
		134
		135		QVector<int> ids(gradients.size());
		136
		137		int maxID = partition(subW, subH, gradients, ids);
		138
		139		// Not needed anymore
		140		delete boundedImage;
		141
		142		if (maxID == 0)
		143		return 0;
		144
		145		typedef struct
		146		{
		147		float massX = 0;
		148		float massY = 0;
		149		float totalMass = 0;
		150		} massInfo;
		151
		152		QMap<int, massInfo> masses;
		153
		154		// #7 Calculate center of mass for all detected regions
		155		for (int y = 0; y < subH; y++)
		156		{
		157		for (int x = 0; x < subW; x++)
		158		{
		159		int index = x + y * subW;
		160
		161		int regionID = ids[index];
		162		if (regionID > 0)
		163		{
		164		float pixel = gradients[index];
		165
		166		masses[regionID].totalMass += pixel;
		167		masses[regionID].massX += x * pixel;
		168		masses[regionID].massY += y * pixel;
		169		}
		170		}
		171		}
		172
		173		// Compare multiple masses, and only select the highest total mass one as the desired star
		174		int maxRegionID = 1;
		175		int maxTotalMass = masses[1].totalMass;
		176		double totalMassRatio = 1e6;
		177		for (auto key : masses.keys())
		178		{
		179		massInfo oneMass = masses.value(key);
		180		if (oneMass.totalMass > maxTotalMass)
		181		{
		182		totalMassRatio = oneMass.totalMass / maxTotalMass;
		183		maxTotalMass = oneMass.totalMass;
		184		maxRegionID = key;
		185		}
		186		}
		187
		188		// If image has many regions and there is no significant relative center of mass then it's just noise and no stars
		189		// are probably there above a useful threshold.
		190		if (maxID > 10 && totalMassRatio < 1.5)
		191		return 0;
		192
		193		auto * center = new Edge;
		194		center->width = -1;
		195		center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5;
		196		center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5;
		197		center->HFR = 1;
		198
		199		// Maximum Radius
		200		int maxR = qMin(subW - 1, subH - 1) / 2;
		201
		202		for (int r = maxR; r > 1; r--)
		203		{
		204		int pass = 0;
		205
		206		for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0)
		207		{
		208		int testX = center->x + std::cos(theta) * r;
		209		int testY = center->y + std::sin(theta) * r;
		210
		211		// if out of bound, break;
		212		if (testX < 0 \|\| testX >= subW \|\| testY < 0 \|\| testY >= subH)
		213		break;
		214
		215		if (gradients[testX + testY * subW] > 0)
		216		//if (thresholded[testX + testY * subW] > 0)
		217		{
		218		if (++pass >= 24)
		219		{
		220		center->width = r * 2;
		221		// Break of outer loop
		222		r = 0;
		223		break;
		224		}
		225		}
		226		}
		227		}
		228
		229		qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width;
		230
		231		// If no stars were detected
		232		if (center->width == -1)
		233		{
		234		delete center;
		235		return 0;
		236		}
		237
		238		// 30% fuzzy
		239		//center->width += center->width0.3 (running_threshold / threshold);
		240
		241		double FSum = 0, HF = 0, TF = 0;
		242		const double resolution = 1.0 / 20.0;
		243
		244		int cen_y = qRound(center->y);
		245
		246		double rightEdge = center->x + center->width / 2.0;
		247		double leftEdge = center->x - center->width / 2.0;
		248
		249		QVector<double> subPixels;
		250		subPixels.reserve(center->width / resolution);
		251
		252		const T * origBuffer = reinterpret_cast<T const *>(data->getImageBuffer()) + offset;
		253
		254		for (double x = leftEdge; x <= rightEdge; x += resolution)
		255		{
		256		double slice = resolution * (origBuffer[static_cast<int>(floor(x)) + cen_y * dataWidth]);
		257		FSum += slice;
		258		subPixels.append(slice);
		259		}
		260
		261		// Half flux
		262		HF = FSum / 2.0;
		263
		264		int subPixelCenter = (center->width / resolution) / 2;
		265
		266		// Start from center
		267		TF = subPixels[subPixelCenter];
		268		double lastTF = TF;
		269		// Integrate flux along radius axis until we reach half flux
		270		//for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
		271		for (int k = 1; k < subPixelCenter; k++)
		272		{
		273		TF += subPixels[subPixelCenter + k];
		274		TF += subPixels[subPixelCenter - k];
		275
		276		if (TF >= HF)
		277		{
		278		// We overpassed HF, let's calculate from last TF how much until we reach HF
		279
		280		// #1 Accurate calculation, but very sensitive to small variations of flux
		281		center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
		282
		283		// #2 Less accurate calculation, but stable against small variations of flux
		284		//center->HFR = (k - 1) * resolution;
		285		break;
		286		}
		287
		288		lastTF = TF;
		289		}
		290
		291		// Correct center for subX and subY
		292		center->x += subX;
		293		center->y += subY;
		294
		295		//data->appendStar(center);
		296		starCenters.append(center);
		297
		298		qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR;
		299
		300		return 1;
		301		}
		302
		303		/* CannyDetector, Implementation of Canny edge detector in Qt/C++.
		304		* Copyright (C) 2015 Gonzalo Exequiel Pedone
		305		*
		306		* This program is free software: you can redistribute it and/or modify
		307		* it under the terms of the GNU General Public License as published by
		308		* the Free Software Foundation, either version 3 of the License, or
		309		* (at your option) any later version.
		310		*
		311		* This program is distributed in the hope that it will be useful,
		312		* but WITHOUT ANY WARRANTY; without even the implied warranty of
		313		* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
		314		* GNU General Public License for more details.
		315		*
		316		* You should have received a copy of the GNU General Public License
		317		* along with this program. If not, see <http://www.gnu.org/licenses/>.
		318		*
		319		* Email : hipersayan DOT x AT gmail DOT com
		320		* Web-Site: https://github.com/hipersayanX/CannyDetector
		321		*/
		322
		323		template <typename T>
		324		void FITSGradientDetector::sobel(FITSData const *data, QVector<float> &gradient, QVector<float> &direction) const
		325		{
		326		if (data == nullptr)
		327		return;
		328
		329		FITSData::Statistic const &stats = data->getStatistics();
		330
		331		//int size = image.width() * image.height();
		332		gradient.resize(stats.samples_per_channel);
		333		direction.resize(stats.samples_per_channel);
		334
		335		for (int y = 0; y < stats.height; y++)
		336		{
		337		size_t yOffset = y * stats.width;
		338		const T * grayLine = reinterpret_cast<T const *>(data->getImageBuffer()) + yOffset;
		339
		340		const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width;
		341		const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width;
		342
		343		float * gradientLine = gradient.data() + yOffset;
		344		float * directionLine = direction.data() + yOffset;
		345
		346		for (int x = 0; x < stats.width; x++)
		347		{
		348		int x_m1 = x < 1 ? x : x - 1;
		349		int x_p1 = x >= stats.width - 1 ? x : x + 1;
		350
		351		int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] -
		352		2 * grayLine[x_m1] - grayLine_p1[x_m1];
		353
		354		int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] -
		355		2 * grayLine_p1[x] - grayLine_p1[x_p1];
		356
		357		gradientLine[x] = qAbs(gradX) + qAbs(gradY);
		358
		359		/* Gradient directions are classified in 4 possible cases
		360		*
		361		* dir 0
		362		*
		363		* x x x
		364		* - - -
		365		* x x x
		366		*
		367		* dir 1
		368		*
		369		* x x /
		370		* x / x
		371		* / x x
		372		*
		373		* dir 2
		374		*
		375		* \ x x
		376		* x \ x
		377		* x x \
		378		*
		379		* dir 3
		380		*
		381		* x \| x
		382		* x \| x
		383		* x \| x
		384		*/
		385		if (gradX == 0 && gradY == 0)
		386		directionLine[x] = 0;
		387		else if (gradX == 0)
		388		directionLine[x] = 3;
		389		else
		390		{
		391		qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI;
		392
		393		if (a >= -22.5 && a < 22.5)
		394		directionLine[x] = 0;
		395		else if (a >= 22.5 && a < 67.5)
		396		directionLine[x] = 2;
		397		else if (a >= -67.5 && a < -22.5)
		398		directionLine[x] = 1;
		399		else
		400		directionLine[x] = 3;
		401		}
		402		}
		403		}
		404		}
		405
		406		int FITSGradientDetector::partition(int width, int height, QVector<float> &gradient, QVector<int> &ids) const
		407		{
		408		int id = 0;
		409
		410		for (int y = 1; y < height - 1; y++)
		411		{
		412		for (int x = 1; x < width - 1; x++)
		413		{
		414		int index = x + y * width;
		415		float val = gradient[index];
		416		if (val > 0 && ids[index] == 0)
		417		{
		418		trace(width, height, ++id, gradient, ids, x, y);
		419		}
		420		}
		421		}
		422
		423		// Return max id
		424		return id;
		425		}
		426
		427		void FITSGradientDetector::trace(int width, int height, int id, QVector<float> &image, QVector<int> &ids, int x, int y) const
		428		{
		429		int yOffset = y * width;
		430		float * cannyLine = image.data() + yOffset;
		431		int * idLine = ids.data() + yOffset;
		432
		433		if (idLine[x] != 0)
		434		return;
		435
		436		idLine[x] = id;
		437
		438		for (int j = -1; j < 2; j++)
		439		{
		440		int nextY = y + j;
		441
		442		if (nextY < 0 \|\| nextY >= height)
		443		continue;
		444
		445		float * cannyLineNext = cannyLine + j * width;
		446
		447		for (int i = -1; i < 2; i++)
		448		{
		449		int nextX = x + i;
		450
		451		if (i == j \|\| nextX < 0 \|\| nextX >= width)
		452		continue;
		453
		454		if (cannyLineNext[nextX] > 0)
		455		{
		456		// Trace neighbors.
		457		trace(width, height, id, image, ids, nextX, nextY);
		458		}
		459		}
		460		}
		461		}