Differential D28292 Diff 78929 kstars/fitsviewer/fitssepdetector.cpp

Changeset View

Standalone View

View Options

kstars/fitsviewer/fitssepdetector.cpp

  • This file was added.
1​/***************************************************************************
2​ fitssepdetector.cpp - FITS Image
3​ -------------------
4​ begin : Sun March 29 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 "sep/sep.h"
23​#include "fits_debug.h"
24​#include "fitssepdetector.h"
25​
26​FITSStarDetector& FITSSEPDetector::configure(const QString &, const QVariant &)
27​{
28​ return *this;
29​}
30​
31​int FITSSEPDetector::findSources(QList<Edge*> &starCenters, const QRect &boundary)
32​{
33​ FITSData const * const image_data = reinterpret_cast<FITSData const *>(parent());
34​
35​ if (image_data == nullptr)
36​ return 0;
37​
38​ FITSData::Statistic const &stats = image_data->getStatistics();
39​
40​ int x = 0, y = 0, w = stats.width, h = stats.height, maxRadius = 50;
41​
42​ if (!boundary.isNull())
43​ {
44​ x = boundary.x();
45​ y = boundary.y();
46​ w = boundary.width();
47​ h = boundary.height();
48​ maxRadius = w;
49​ }
50​
51​ auto * data = new float[w * h];
52​
53​ switch (stats.bitpix)
54​ {
55​ case BYTE_IMG:
56​ getFloatBuffer<uint8_t>(data, x, y, w, h, image_data);
57​ break;
58​ case SHORT_IMG:
59​ getFloatBuffer<int16_t>(data, x, y, w, h, image_data);
60​ break;
61​ case USHORT_IMG:
62​ getFloatBuffer<uint16_t>(data, x, y, w, h, image_data);
63​ break;
64​ case LONG_IMG:
65​ getFloatBuffer<int32_t>(data, x, y, w, h, image_data);
66​ break;
67​ case ULONG_IMG:
68​ getFloatBuffer<uint32_t>(data, x, y, w, h, image_data);
69​ break;
70​ case FLOAT_IMG:
71​ memcpy(data, image_data->getImageBuffer(), sizeof(float)*w*h);
72​ break;
73​ case LONGLONG_IMG:
74​ getFloatBuffer<int64_t>(data, x, y, w, h, image_data);
75​ break;
76​ case DOUBLE_IMG:
77​ getFloatBuffer<double>(data, x, y, w, h, image_data);
78​ break;
79​ default:
80​ delete [] data;
81​ return -1;
82​ }
83​
84​ float * imback = nullptr;
85​ double * flux = nullptr, *fluxerr = nullptr, *area = nullptr;
86​ short * flag = nullptr;
87​ short flux_flag = 0;
88​ int status = 0;
89​ sep_bkg * bkg = nullptr;
90​ sep_catalog * catalog = nullptr;
91​ float conv[] = {1, 2, 1, 2, 4, 2, 1, 2, 1};
92​ double flux_fractions[2] = {0};
93​ double requested_frac[2] = { 0.5, 0.99 };
94​ QList<Edge *> edges;
95​
96​ // #0 Create SEP Image structure
97​ sep_image im = {data, nullptr, nullptr, SEP_TFLOAT, 0, 0, w, h, 0.0, SEP_NOISE_NONE, 1.0, 0.0};
98​
99​ // #1 Background estimate
100​ status = sep_background(&im, 64, 64, 3, 3, 0.0, &bkg);
101​ if (status != 0) goto exit;
102​
103​ // #2 Background evaluation
104​ imback = (float *)malloc((w * h) * sizeof(float));
105​ status = sep_bkg_array(bkg, imback, SEP_TFLOAT);
106​ if (status != 0) goto exit;
107​
108​ // #3 Background subtraction
109​ status = sep_bkg_subarray(bkg, im.data, im.dtype);
110​ if (status != 0) goto exit;
111​
112​ // #4 Source Extraction
113​ // Note that we set deblend_cont = 1.0 to turn off deblending.
114​ status = sep_extract(&im, 2 * bkg->globalrms, SEP_THRESH_ABS, 10, conv, 3, 3, SEP_FILTER_CONV, 32, 1.0, 1, 1.0, &catalog);
115​ if (status != 0) goto exit;
116​
117​ // TODO
118​ // Must detect edge detection
119​ // Must limit to brightest 100 (by flux) centers
120​ // Should probably use ellipse to draw instead of simple circle?
121​ // Useful for galaxies and also elenogated stars.
122​ for (int i = 0; i < catalog->nobj; i++)
123​ {
124​ double flux = catalog->flux[i];
125​ // Get HFR
126​ sep_flux_radius(&im, catalog->x[i], catalog->y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag);
127​
128​ auto * center = new Edge();
129​ center->x = catalog->x[i] + x + 0.5;
130​ center->y = catalog->y[i] + y + 0.5;
131​ center->val = catalog->peak[i];
132​ center->sum = flux;
133​ center->HFR = center->width = flux_fractions[0];
134​ if (flux_fractions[1] < maxRadius)
135​ center->width = flux_fractions[1] * 2;
136​ edges.append(center);
137​ }
138​
139​ // Let's sort edges, starting with widest
140​ std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->width > edge2->width;});
141​
142​ // Take only the first 100 stars
143​ {
144​ int starCount = qMin(100, edges.count());
145​ for (int i = 0; i < starCount; i++)
146​ starCenters.append(edges[i]);
147​ }
148​
149​ edges.clear();
150​
151​ qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << "#" << "#X" << "#Y" << "#Flux" << "#Width" << "#HFR";
152​ for (int i = 0; i < starCenters.count(); i++)
153​ qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << i << starCenters[i]->x << starCenters[i]->y
154​ << starCenters[i]->sum << starCenters[i]->width << starCenters[i]->HFR;
155​
156​exit:
157​ if (stats.bitpix != FLOAT_IMG)
158​ delete [] data;
159​ sep_bkg_free(bkg);
160​ sep_catalog_free(catalog);
161​ free(imback);
162​ free(flux);
163​ free(fluxerr);
164​ free(area);
165​ free(flag);
166​
167​ if (status != 0)
168​ {
169​ char errorMessage[512];
170​ sep_get_errmsg(status, errorMessage);
171​ qCritical(KSTARS_FITS) << errorMessage;
172​ return -1;
173​ }
174​
175​ return starCenters.count();
176​}
177​
178​template <typename T>
179​void FITSSEPDetector::getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const *data) const
180​{
181​ auto * rawBuffer = reinterpret_cast<T const *>(data->getImageBuffer());
182​
183​ if (buffer == nullptr)
184​ return;
185​
186​ float * floatPtr = buffer;
187​
188​ int x2 = x + w;
189​ int y2 = y + h;
190​
191​ FITSData::Statistic const &stats = data->getStatistics();
192​
193​ for (int y1 = y; y1 < y2; y1++)
194​ {
195​ int offset = y1 * stats.width;
196​ for (int x1 = x; x1 < x2; x1++)
197​ {
198​ *floatPtr++ = rawBuffer[offset + x1];
199​ }
200​ }
201​}