File indexing completed on 2024-04-14 14:10:38

 /*
     SPDX-FileCopyrightText: 2004 Jasem Mutlaq
     SPDX-FileCopyrightText: 2020 Eric Dejouhanet <eric.dejouhanet@gmail.com>
 
     SPDX-License-Identifier: GPL-2.0-or-later
 
     Some code fragments were adapted from Peter Kirchgessner's FITS plugin
     SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg>
 */
 
 #include <math.h>
 #include <cmath>
 #include <QtConcurrent>
 
 #include "fits_debug.h"
 #include "fitsthresholddetector.h"
 #include "fitsdata.h"
 
 //void FITSThresholdDetector::configure(const QString &setting, const QVariant &value)
 //{
 //    if (!setting.compare("THRESHOLD_PERCENTAGE", Qt::CaseInsensitive))
 //        if (value.canConvert<int>())
 //            THRESHOLD_PERCENTAGE = value.value<int>();
 //}
 
 QFuture<bool> FITSThresholdDetector::findSources(QRect const &boundary)
 {
     FITSImage::Statistic const &stats = m_ImageData->getStatistics();
     switch (stats.dataType)
     {
         case TSHORT:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int16_t>, boundary);
 
         case TUSHORT:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint16_t>, boundary);
 
         case TLONG:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int32_t>, boundary);
 
         case TULONG:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint32_t>, boundary);
 
         case TFLOAT:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<float>, boundary);
 
         case TLONGLONG:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<int64_t>, boundary);
 
         case TDOUBLE:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<double>, boundary);
 
         case TBYTE:
         default:
             return QtConcurrent::run(this, &FITSThresholdDetector::findOneStar<uint8_t>, boundary);
     }
 
 }
 
 template <typename T>
 bool FITSThresholdDetector::findOneStar(const QRect &boundary) const
 {
     FITSImage::Statistic const &stats = m_ImageData->getStatistics();
 
     int subX = boundary.x();
     int subY = boundary.y();
     int subW = subX + boundary.width();
     int subH = subY + boundary.height();
 
     float massX = 0, massY = 0, totalMass = 0;
 
     auto * buffer = reinterpret_cast<T const *>(m_ImageData->getImageBuffer());
 
     double THRESHOLD_PERCENTAGE = getValue("THRESHOLD_PERCENTAGE", 120.0).toDouble();
     // TODO replace magic number with something more useful to understand
     double threshold = stats.mean[0] * THRESHOLD_PERCENTAGE / 100.0;
 
     for (int y = subY; y < subH; y++)
     {
         for (int x = subX; x < subW; x++)
         {
             T pixel = buffer[x + y * stats.width];
             if (pixel > threshold)
             {
                 totalMass += pixel;
                 massX += x * pixel;
                 massY += y * pixel;
             }
         }
     }
 
     qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << massX / totalMass << " Y: " << massY / totalMass;
 
     auto * center  = new Edge;
     center->width = -1;
     center->x     = massX / totalMass + 0.5;
     center->y     = massY / totalMass + 0.5;
     center->HFR   = 1;
 
     // Maximum Radius
     int maxR = qMin(subW - 1, subH - 1) / 2;
 
     // Critical threshold
     double critical_threshold = threshold * 0.7;
     double running_threshold  = threshold;
 
     while (running_threshold >= critical_threshold)
     {
         for (int r = maxR; r > 1; r--)
         {
             int pass = 0;
 
             for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 10.0)
             {
                 int testX = center->x + std::cos(theta) * r;
                 int testY = center->y + std::sin(theta) * r;
 
                 // if out of bound, break;
                 if (testX < subX || testX > subW || testY < subY || testY > subH)
                     break;
 
                 if (buffer[testX + testY * stats.width] > running_threshold)
                     pass++;
             }
 
             //qDebug() << Q_FUNC_INFO << "Testing for radius " << r << " passes # " << pass << " @ threshold " << running_threshold;
             //if (pass >= 6)
             if (pass >= 5)
             {
                 center->width = r * 2;
                 break;
             }
         }
 
         if (center->width > 0)
             break;
 
         // Increase threshold fuzziness by 10%
         running_threshold -= running_threshold * 0.1;
     }
 
     // If no stars were detected
     if (center->width == -1)
     {
         delete center;
         return false;
     }
 
     // 30% fuzzy
     //center->width += center->width*0.3 * (running_threshold / threshold);
 
     double FSum = 0, HF = 0, TF = 0, min = stats.min[0];
     const double resolution = 1.0 / 20.0;
 
     int cen_y = qRound(center->y);
 
     double rightEdge = center->x + center->width / 2.0;
     double leftEdge  = center->x - center->width / 2.0;
 
     QVector<double> subPixels;
     subPixels.reserve(center->width / resolution);
 
     for (double x = leftEdge; x <= rightEdge; x += resolution)
     {
         //subPixels[x] = resolution * (image_buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
         double slice = resolution * (buffer[static_cast<int>(floor(x)) + cen_y * stats.width] - min);
         FSum += slice;
         subPixels.append(slice);
     }
 
     // Half flux
     HF = FSum / 2.0;
 
     //double subPixelCenter = center->x - fmod(center->x,resolution);
     int subPixelCenter = (center->width / resolution) / 2;
 
     // Start from center
     TF            = subPixels[subPixelCenter];
     double lastTF = TF;
     // Integrate flux along radius axis until we reach half flux
     //for (double k=resolution; k < (center->width/(2*resolution)); k += resolution)
     for (int k = 1; k < subPixelCenter; k++)
     {
         TF += subPixels[subPixelCenter + k];
         TF += subPixels[subPixelCenter - k];
 
         if (TF >= HF)
         {
             // 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.
             // The second method is not truly HFR but is much more resistant to noise.
 
             // #1 Approximate HFR, accurate and reliable but quite variable to small changes in star flux
             center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution;
 
             // #2 Not exactly HFR, but much more stable
             //center->HFR = (k*resolution) * (HF/TF);
             break;
         }
 
         lastTF = TF;
     }
 
     QList<Edge*> starCenters;
     starCenters.append(center);
     m_ImageData->setStarCenters(starCenters);
 
     return true;
 }