File indexing completed on 2024-03-24 03:46:32
0001 /* 0002 SPDX-FileCopyrightText: 2004 Jasem Mutlaq 0003 SPDX-FileCopyrightText: 2020 Eric Dejouhanet <eric.dejouhanet@gmail.com> 0004 0005 SPDX-License-Identifier: GPL-2.0-or-later 0006 0007 Some code fragments were adapted from Peter Kirchgessner's FITS plugin 0008 SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg> 0009 */ 0010 0011 #include <math.h> 0012 #include <cmath> 0013 #include <QtConcurrent> 0014 0015 #include "fitscentroiddetector.h" 0016 #include "fits_debug.h" 0017 #include "fitsdata.h" 0018 0019 //void FITSCentroidDetector::configure(const QString &setting, const QVariant &value) 0020 //{ 0021 // if (!setting.compare("MINIMUM_STDVAR", Qt::CaseInsensitive)) 0022 // if (value.canConvert <int> ()) 0023 // MINIMUM_STDVAR = value.value <int> (); 0024 0025 // if (!setting.compare("MINIMUM_PIXEL_RANGE", Qt::CaseInsensitive)) 0026 // if (value.canConvert <int> ()) 0027 // MINIMUM_PIXEL_RANGE = value.value <int> (); 0028 0029 // if (!setting.compare("JMINDEX", Qt::CaseInsensitive)) 0030 // if (value.canConvert <double> ()) 0031 // JMINDEX = value.value <double> (); 0032 //} 0033 0034 bool FITSCentroidDetector::checkCollision(Edge * s1, Edge * s2) const 0035 { 0036 int dis; //distance 0037 0038 int diff_x = s1->x - s2->x; 0039 int diff_y = s1->y - s2->y; 0040 0041 dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y)); 0042 dis -= s1->width / 2; 0043 dis -= s2->width / 2; 0044 0045 if (dis <= 0) //collision 0046 return true; 0047 0048 //no collision 0049 return false; 0050 } 0051 0052 /*** Find center of stars and calculate Half Flux Radius */ 0053 QFuture<bool> FITSCentroidDetector::findSources(const QRect &boundary) 0054 { 0055 FITSImage::Statistic const &stats = m_ImageData->getStatistics(); 0056 switch (stats.dataType) 0057 { 0058 case TBYTE: 0059 default: 0060 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint8_t const>, boundary); 0061 0062 case TSHORT: 0063 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int16_t const>, boundary); 0064 0065 case TUSHORT: 0066 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint16_t const>, boundary); 0067 0068 case TLONG: 0069 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int32_t const>, boundary); 0070 0071 case TULONG: 0072 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<uint32_t const>, boundary); 0073 0074 case TFLOAT: 0075 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<float const>, boundary); 0076 0077 case TLONGLONG: 0078 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<int64_t const>, boundary); 0079 0080 case TDOUBLE: 0081 return QtConcurrent::run(this, &FITSCentroidDetector::findSources<double const>, boundary); 0082 0083 } 0084 } 0085 0086 template <typename T> 0087 bool FITSCentroidDetector::findSources(const QRect &boundary) 0088 { 0089 FITSImage::Statistic const &stats = m_ImageData->getStatistics(); 0090 FITSMode const m_Mode = static_cast<FITSMode>(m_ImageData->property("mode").toInt()); 0091 0092 int MINIMUM_STDVAR = getValue("MINIMUM_STDVAR", 5).toInt(); 0093 int minEdgeWidth = getValue("MINIMUM_PIXEL_RANGE", 5).toInt(); 0094 double JMIndex = getValue("JMINDEX", 100.0).toDouble(); 0095 0096 int initStdDev = MINIMUM_STDVAR; 0097 double threshold = 0, sum = 0, avg = 0, min = 0; 0098 int starDiameter = 0; 0099 int pixVal = 0; 0100 int minimumEdgeCount = MINIMUM_EDGE_LIMIT; 0101 0102 auto * buffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer()); 0103 0104 float dispersion_ratio = 1.5; 0105 0106 QList<Edge *> edges; 0107 0108 if (JMIndex < DIFFUSE_THRESHOLD) 0109 { 0110 minEdgeWidth = JMIndex * 35 + 1; 0111 minimumEdgeCount = minEdgeWidth - 1; 0112 } 0113 else 0114 { 0115 minEdgeWidth = 6; 0116 minimumEdgeCount = 4; 0117 } 0118 0119 while (initStdDev >= 1) 0120 { 0121 minEdgeWidth--; 0122 minimumEdgeCount--; 0123 0124 minEdgeWidth = qMax(3, minEdgeWidth); 0125 minimumEdgeCount = qMax(3, minimumEdgeCount); 0126 0127 if (JMIndex < DIFFUSE_THRESHOLD) 0128 { 0129 // Taking the average out seems to have better result for noisy images 0130 threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); 0131 0132 min = stats.min[0]; 0133 if (threshold - min < 0) 0134 { 0135 threshold = stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); 0136 min = 0; 0137 } 0138 0139 dispersion_ratio = 1.4 - (MINIMUM_STDVAR - initStdDev) * 0.08; 0140 } 0141 else 0142 { 0143 threshold = stats.mean[0] + stats.stddev[0] * initStdDev * (0.3 - (MINIMUM_STDVAR - initStdDev) * 0.05); 0144 min = stats.min[0]; 0145 // Ratio between centeroid center and edge 0146 dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2; 0147 } 0148 0149 qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR; 0150 qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min 0151 << ") minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount; 0152 0153 threshold -= min; 0154 0155 int subX, subY, subW, subH; 0156 0157 if (boundary.isNull()) 0158 { 0159 if (m_Mode == FITS_GUIDE || m_Mode == FITS_FOCUS) 0160 { 0161 // Only consider the central 70% 0162 subX = round(stats.width * 0.15); 0163 subY = round(stats.height * 0.15); 0164 subW = stats.width - subX; 0165 subH = stats.height - subY; 0166 } 0167 else 0168 { 0169 // Consider the complete area 100% 0170 subX = 0; 0171 subY = 0; 0172 subW = stats.width; 0173 subH = stats.height; 0174 } 0175 } 0176 else 0177 { 0178 subX = boundary.x(); 0179 subY = boundary.y(); 0180 subW = subX + boundary.width(); 0181 subH = subY + boundary.height(); 0182 } 0183 0184 // Detect "edges" that are above threshold 0185 for (int i = subY; i < subH; i++) 0186 { 0187 starDiameter = 0; 0188 0189 for (int j = subX; j < subW; j++) 0190 { 0191 pixVal = buffer[j + (i * stats.width)] - min; 0192 0193 // If pixel value > threshold, let's get its weighted average 0194 if (pixVal >= threshold) 0195 { 0196 avg += j * pixVal; 0197 sum += pixVal; 0198 starDiameter++; 0199 } 0200 // Value < threshold but avg exists 0201 else if (sum > 0) 0202 { 0203 // We found a potential centroid edge 0204 if (starDiameter >= minEdgeWidth) 0205 { 0206 float center = avg / sum + 0.5; 0207 if (center > 0) 0208 { 0209 int i_center = std::floor(center); 0210 0211 // Check if center is 10% or more brighter than edge, if not skip 0212 if (((buffer[i_center + (i * stats.width)] - min) / 0213 (buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >= 0214 dispersion_ratio) && 0215 ((buffer[i_center + (i * stats.width)] - min) / 0216 (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >= 0217 dispersion_ratio)) 0218 { 0219 qCDebug(KSTARS_FITS) 0220 << "Edge center is " << buffer[i_center + (i * stats.width)] - min 0221 << " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min 0222 << " and ratio is " 0223 << ((buffer[i_center + (i * stats.width)] - min) / 0224 (buffer[i_center + (i * stats.width) - starDiameter / 2] - min)) 0225 << " located at X: " << center << " Y: " << i + 0.5; 0226 0227 auto * newEdge = new Edge(); 0228 0229 newEdge->x = center; 0230 newEdge->y = i + 0.5; 0231 newEdge->scanned = 0; 0232 newEdge->val = buffer[i_center + (i * stats.width)] - min; 0233 newEdge->width = starDiameter; 0234 newEdge->HFR = 0; 0235 newEdge->sum = sum; 0236 0237 edges.append(newEdge); 0238 } 0239 } 0240 } 0241 0242 // Reset 0243 avg = sum = starDiameter = 0; 0244 } 0245 } 0246 } 0247 0248 qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count(); 0249 0250 // In case of hot pixels 0251 if (edges.count() == 1 && initStdDev > 1) 0252 { 0253 initStdDev--; 0254 continue; 0255 } 0256 0257 if (edges.count() >= MAX_EDGE_LIMIT) 0258 { 0259 qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count(); 0260 qDeleteAll(edges); 0261 return -1; 0262 } 0263 0264 if (edges.count() >= minimumEdgeCount) 0265 break; 0266 0267 qDeleteAll(edges); 0268 edges.clear(); 0269 initStdDev--; 0270 } 0271 0272 int cen_count = 0; 0273 int cen_x = 0; 0274 int cen_y = 0; 0275 int cen_v = 0; 0276 int cen_w = 0; 0277 int width_sum = 0; 0278 0279 // Let's sort edges, starting with widest 0280 auto const greaterThan = [](Edge const * a, Edge const * b) 0281 { 0282 return a->sum > b->sum; 0283 }; 0284 std::sort(edges.begin(), edges.end(), greaterThan); 0285 0286 QList<Edge*> starCenters; 0287 // Now, let's scan the edges and find the maximum centroid vertically 0288 for (int i = 0; i < edges.count(); i++) 0289 { 0290 qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of " 0291 << edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum; 0292 0293 // If edge scanned already, skip 0294 if (edges[i]->scanned == 1) 0295 { 0296 qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted"; 0297 continue; 0298 } 0299 0300 qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ..."; 0301 0302 // Get X, Y, and Val of edge 0303 cen_x = edges[i]->x; 0304 cen_y = edges[i]->y; 0305 cen_v = edges[i]->sum; 0306 cen_w = edges[i]->width; 0307 0308 float avg_x = 0; 0309 float avg_y = 0; 0310 0311 sum = 0; 0312 cen_count = 0; 0313 0314 // Now let's compare to other edges until we hit a maxima 0315 for (int j = 0; j < edges.count(); j++) 0316 { 0317 if (edges[j]->scanned) 0318 continue; 0319 0320 if (checkCollision(edges[j], edges[i])) 0321 { 0322 if (edges[j]->sum >= cen_v) 0323 { 0324 cen_v = edges[j]->sum; 0325 cen_w = edges[j]->width; 0326 } 0327 0328 edges[j]->scanned = 1; 0329 cen_count++; 0330 0331 avg_x += edges[j]->x * edges[j]->val; 0332 avg_y += edges[j]->y * edges[j]->val; 0333 sum += edges[j]->val; 0334 0335 continue; 0336 } 0337 } 0338 0339 int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev)); 0340 0341 if (edges.count() < LOW_EDGE_CUTOFF_1) 0342 { 0343 if (edges.count() < LOW_EDGE_CUTOFF_2) 0344 cen_limit = 1; 0345 else 0346 cen_limit = 2; 0347 } 0348 0349 qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is " 0350 << cen_limit; 0351 0352 if (cen_limit < 1) 0353 continue; 0354 0355 // If centroid count is within acceptable range 0356 //if (cen_limit >= 2 && cen_count >= cen_limit) 0357 if (cen_count >= cen_limit) 0358 { 0359 // We detected a centroid, let's init it 0360 auto * rCenter = new Edge(); 0361 0362 rCenter->x = avg_x / sum; 0363 rCenter->y = avg_y / sum; 0364 width_sum += rCenter->width; 0365 rCenter->width = cen_w; 0366 0367 qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")"; 0368 0369 // Calculate Total Flux From Center, Half Flux, Full Summation 0370 double TF = 0; 0371 double HF = 0; 0372 double FSum = 0; 0373 0374 cen_x = (int)std::floor(rCenter->x); 0375 cen_y = (int)std::floor(rCenter->y); 0376 0377 if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height) 0378 { 0379 delete rCenter; 0380 continue; 0381 } 0382 0383 // Complete sum along the radius 0384 //for (int k=0; k < rCenter->width; k++) 0385 for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--) 0386 { 0387 FSum += buffer[cen_x - k + (cen_y * stats.width)] - min; 0388 //qDebug() << Q_FUNC_INFO << image_buffer[cen_x-k+(cen_y*stats.width)] - min; 0389 } 0390 0391 // Half flux 0392 HF = FSum / 2.0; 0393 0394 // Total flux starting from center 0395 TF = buffer[cen_y * stats.width + cen_x] - min; 0396 0397 int pixelCounter = 1; 0398 0399 // Integrate flux along radius axis until we reach half flux 0400 for (int k = 1; k < rCenter->width / 2; k++) 0401 { 0402 if (TF >= HF) 0403 { 0404 qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels."; 0405 break; 0406 } 0407 0408 TF += buffer[cen_y * stats.width + cen_x + k] - min; 0409 TF += buffer[cen_y * stats.width + cen_x - k] - min; 0410 0411 pixelCounter++; 0412 } 0413 0414 // Calculate weighted Half Flux Radius 0415 rCenter->HFR = pixelCounter * (HF / TF); 0416 // Store full flux 0417 rCenter->val = FSum; 0418 0419 qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum; 0420 0421 starCenters.append(rCenter); 0422 } 0423 } 0424 0425 if (starCenters.count() > 1 && m_Mode != FITS_FOCUS) 0426 { 0427 float width_avg = (float)width_sum / starCenters.count(); 0428 float lsum = 0, sdev = 0; 0429 0430 for (auto ¢er : starCenters) 0431 lsum += (center->width - width_avg) * (center->width - width_avg); 0432 0433 sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4; 0434 0435 // Reject stars > 4 * stddev 0436 foreach (Edge * center, starCenters) 0437 if (center->width > sdev) 0438 starCenters.removeOne(center); 0439 0440 //foreach(Edge *center, starCenters) 0441 //qDebug() << Q_FUNC_INFO << center->x << "," << center->y << "," << center->width << "," << center->val << Qt::endl; 0442 } 0443 0444 m_ImageData->setStarCenters(starCenters); 0445 // Release memory 0446 qDeleteAll(edges); 0447 0448 return true; 0449 } 0450 0451