diff --git a/Tests/kstars_ui/CMakeLists.txt b/Tests/kstars_ui/CMakeLists.txt --- a/Tests/kstars_ui/CMakeLists.txt +++ b/Tests/kstars_ui/CMakeLists.txt @@ -2,6 +2,7 @@ SET(KSTARS_UI_TESTS_SRC kstars_ui_tests.cpp test_ekos.cpp + test_ekos_focus.cpp test_ekos_simulator.cpp) include_directories(${CFITSIO_INCLUDE_DIR}) diff --git a/Tests/kstars_ui/kstars_ui_tests.cpp b/Tests/kstars_ui/kstars_ui_tests.cpp --- a/Tests/kstars_ui/kstars_ui_tests.cpp +++ b/Tests/kstars_ui/kstars_ui_tests.cpp @@ -10,6 +10,7 @@ #include "kstars_ui_tests.h" #include "test_ekos.h" #include "test_ekos_simulator.h" +#include "test_ekos_focus.h" #include "auxiliary/kspaths.h" #if defined(HAVE_INDI) @@ -148,12 +149,21 @@ { TestEkos * ek = new TestEkos(); result |= QTest::qExec(ek, argc, argv); + delete ek; } if (!result) { - TestEkosSimulator * eks = new TestEkosSimulator(); - result |= QTest::qExec(eks, argc, argv); + TestEkosSimulator * ek = new TestEkosSimulator(); + result |= QTest::qExec(ek, argc, argv); + delete ek; + } + + if (!result) + { + TestEkosFocus * ek = new TestEkosFocus(); + result |= QTest::qExec(ek, argc, argv); + delete ek; } #endif diff --git a/Tests/kstars_ui/test_ekos_focus.h b/Tests/kstars_ui/test_ekos_focus.h new file mode 100644 --- /dev/null +++ b/Tests/kstars_ui/test_ekos_focus.h @@ -0,0 +1,110 @@ +/* KStars UI tests + Copyright (C) 2020 + Eric Dejouhanet + + This application is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public + License as published by the Free Software Foundation; either + version 2 of the License, or (at your option) any later version. + */ + +#ifndef TESTEKOSFOCUS_H +#define TESTEKOSFOCUS_H + +#include "config-kstars.h" + +#if defined(HAVE_INDI) + +#include + +/** @brief Helper to retrieve a gadget in the Focus tab specifically. + * @param klass is the class of the gadget to look for. + * @param name is the gadget name to look for in the UI configuration. + * @warning Fails the test if the gadget "name" of class "klass" does not exist in the Focus module + */ +#define KTRY_FOCUS_GADGET(klass, name) klass * const name = Ekos::Manager::Instance()->focusModule()->findChild(#name); \ + QVERIFY2(name != nullptr, QString(#klass " '%1' does not exist and cannot be used").arg(#name).toStdString().c_str()) + +/** @brief Helper to click a button in the Focus tab specifically. + * @param button is the gadget name of the button to look for in the UI configuration. + * @warning Fails the test if the button is not currently enabled. + */ +#define KTRY_FOCUS_CLICK(button) do { \ + QTimer::singleShot(200, Ekos::Manager::Instance(), []() { \ + KTRY_FOCUS_GADGET(QPushButton, button); \ + QVERIFY2(button->isEnabled(), QString("QPushButton '%1' is disabled and cannot be clicked").arg(#button).toStdString().c_str()); \ + QTest::mouseClick(button, Qt::LeftButton); }); } while(false) + +/** @brief Helper to set a string text into a QComboBox in the Focus module. + * @param combobox is the gadget name of the QComboBox to look for in the UI configuration. + * @param text is the string text to set in the gadget. + * @note This is a contrived method to set a text into a QComboBox programmatically *and* emit the "activated" message. + * @warning Fails the test if the name does not exist in the Focus UI or if the text cannot be set in the gadget. + */ +#define KTRY_FOCUS_COMBO_SET(combobox, text) do { \ + KTRY_FOCUS_GADGET(QComboBox, combobox); \ + int const cbIndex = combobox->findText(text); \ + QVERIFY(0 <= cbIndex); \ + combobox->setCurrentIndex(cbIndex); \ + combobox->activated(cbIndex); \ + QCOMPARE(combobox->currentText(), QString(text)); } while(false); + +/** @brief Helper for exposure. + * @param exposure is the amount of seconds to expose for. + * @param averaged is the number of frames the procedure should average before computation. + * @note The Focus capture button is disabled during exposure. + * @warning Fails the test if the exposure cannot be entered or if the capture button is + * disabled or does not toggle during exposure or if the stop button is not the opposite of the capture button. + */ +#define KTRY_FOCUS_CAPTURE(exposure, averaged) do { \ + KTRY_FOCUS_GADGET(QDoubleSpinBox, exposureIN); \ + exposureIN->setValue(exposure); \ + KTRY_FOCUS_GADGET(QSpinBox, focusFramesSpin); \ + focusFramesSpin->setValue(averaged); \ + KTRY_FOCUS_GADGET(QPushButton, captureB); \ + KTRY_FOCUS_GADGET(QPushButton, stopFocusB); \ + QTRY_VERIFY_WITH_TIMEOUT(captureB->isEnabled(), 1000); \ + QTRY_VERIFY_WITH_TIMEOUT(!stopFocusB->isEnabled(), 1000); \ + KTRY_FOCUS_CLICK(captureB); \ + QTRY_VERIFY_WITH_TIMEOUT(!captureB->isEnabled(), 1000); \ + QTRY_VERIFY_WITH_TIMEOUT(stopFocusB->isEnabled(), 1000); \ + QTest::qWait(1.2*exposure*averaged); \ + QTRY_VERIFY_WITH_TIMEOUT(captureB->isEnabled(), 5000); \ + QTRY_VERIFY_WITH_TIMEOUT(!stopFocusB->isEnabled(), 5000); } while (false) + +/** brief Helper to configure main star detection parameters. + * @param detection is the name of the star detection method to use. + * @param algorithm is the name of the autofocus algorithm to use. + * @param fieldin is the lower radius of the annulus filtering stars. + * @param fieldout is the upper radius of the annulus filtering stars. + * @warning Fails the test if detection, algorithm, full-field checkbox or annulus fields cannot be used. + */ +#define KTRY_FOCUS_CONFIGURE(detection, algorithm, fieldin, fieldout) do { \ + KTRY_FOCUS_GADGET(QCheckBox, useFullField); \ + useFullField->setCheckState(Qt::CheckState::Checked); \ + KTRY_FOCUS_GADGET(QDoubleSpinBox, fullFieldInnerRing); \ + fullFieldInnerRing->setValue(fieldin); \ + KTRY_FOCUS_GADGET(QDoubleSpinBox, fullFieldOuterRing); \ + fullFieldOuterRing->setValue(fieldout); \ + KTRY_FOCUS_COMBO_SET(focusDetectionCombo, detection); \ + KTRY_FOCUS_COMBO_SET(focusAlgorithmCombo, algorithm); } while (false) + +class TestEkosFocus : public QObject +{ + Q_OBJECT +public: + explicit TestEkosFocus(QObject *parent = nullptr); + +private slots: + void initTestCase(); + void cleanupTestCase(); + + void init(); + void cleanup(); + + void testStarDetection_data(); + void testStarDetection(); +}; + +#endif +#endif // TESTEKOSFOCUS_H diff --git a/Tests/kstars_ui/test_ekos_focus.cpp b/Tests/kstars_ui/test_ekos_focus.cpp new file mode 100644 --- /dev/null +++ b/Tests/kstars_ui/test_ekos_focus.cpp @@ -0,0 +1,152 @@ +/* KStars UI tests + Copyright (C) 2020 + Eric Dejouhanet + + This application is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public + License as published by the Free Software Foundation; either + version 2 of the License, or (at your option) any later version. + */ + +#include "kstars_ui_tests.h" +#include "test_ekos.h" +#include "test_ekos_simulator.h" +#include "test_ekos_focus.h" +#include "ekos/manager.h" +#include "kstars.h" + +#if defined(HAVE_INDI) + +TestEkosFocus::TestEkosFocus(QObject *parent) : QObject(parent) +{ + +} + +void TestEkosFocus::initTestCase() +{ + KVERIFY_EKOS_IS_HIDDEN(); + KTRY_OPEN_EKOS(); + KVERIFY_EKOS_IS_OPENED(); + KTRY_EKOS_START_SIMULATORS(); + +} + +void TestEkosFocus::cleanupTestCase() +{ + KTRY_EKOS_STOP_SIMULATORS(); + KTRY_CLOSE_EKOS(); + KVERIFY_EKOS_IS_HIDDEN(); +} + +void TestEkosFocus::init() +{ + +} + +void TestEkosFocus::cleanup() +{ + +} + +void TestEkosFocus::testStarDetection_data() +{ +#if QT_VERSION < 0x050900 + QSKIP("Skipping fixture-based test on old QT version."); +#else + QTest::addColumn("NAME"); + QTest::addColumn("RA"); + QTest::addColumn("DEC"); + + // Altitude computation taken from SchedulerJob::findAltitude + GeoLocation * const geo = KStarsData::Instance()->geo(); + KStarsDateTime const now(KStarsData::Instance()->lt()); + KSNumbers const numbers(now.djd()); + CachingDms const LST = geo->GSTtoLST(geo->LTtoUT(now).gst()); + + std::list Objects = { "Polaris", "Mizar", "M 51", "M 13", "M 47", "Pherkab", "Dubhe", "Vega", "NGC 2238", "M 81" }; + size_t count = 0; + + foreach (char const *name, Objects) + { + SkyObject const * const so = KStars::Instance()->data()->objectNamed(name); + if (so != nullptr) + { + SkyObject o(*so); + o.updateCoordsNow(&numbers); + o.EquatorialToHorizontal(&LST, geo->lat()); + if (10.0 < so->alt().Degrees()) + { + QTest::addRow("%s", name) + << name + << so->ra().toHMSString() + << so->dec().toDMSString(); + count++; + } + else QWARN(QString("Fixture '%1' altitude is '%2' degrees, discarding.").arg(name).arg(so->alt().Degrees()).toStdString().c_str()); + } + } + + if (!count) + QSKIP("No usable fixture objects, bypassing test."); +#endif +} + +void TestEkosFocus::testStarDetection() +{ +#if QT_VERSION < 0x050900 + QSKIP("Skipping fixture-based test on old QT version."); +#else + Ekos::Manager * const ekos = Ekos::Manager::Instance(); + + QFETCH(QString, NAME); + QFETCH(QString, RA); + QFETCH(QString, DEC); + qDebug("Test focusing on '%s' RA '%s' DEC '%s'", + NAME.toStdString().c_str(), + RA.toStdString().c_str(), + DEC.toStdString().c_str()); + + // Just sync to RA/DEC to make the mount teleport to the object + QWARN("During this test, the mount is not tracking - we leave it as is for the feature in the CCD simulator to trigger a failure."); + QTRY_VERIFY_WITH_TIMEOUT(ekos->mountModule() != nullptr, 5000); + QVERIFY(ekos->mountModule()->sync(RA, DEC)); + + // Wait for Focus to come up, switch to Focus tab + QTRY_VERIFY_WITH_TIMEOUT(ekos->focusModule() != nullptr, 5000); + KTRY_EKOS_GADGET(QTabWidget, toolsWidget); + toolsWidget->setCurrentWidget(ekos->focusModule()); + QTRY_COMPARE_WITH_TIMEOUT(toolsWidget->currentWidget(), ekos->focusModule(), 1000); + + KTRY_FOCUS_GADGET(QLineEdit, starsOut); + + // Run the focus procedure for SEP + KTRY_FOCUS_CONFIGURE("SEP", "Iterative", 0.0, 100.0); + KTRY_FOCUS_CAPTURE(1, 1); + QWARN("No way to know when star detection procedure is fininshed."); + QTest::qWait(1000); + QTRY_VERIFY_WITH_TIMEOUT(starsOut->text().toInt() >= 1, 5000); + + // Run the focus procedure for Centroid + KTRY_FOCUS_CONFIGURE("Centroid", "Iterative", 0.0, 100.0); + KTRY_FOCUS_CAPTURE(1, 1); + QWARN("No way to know when star detection procedure is fininshed."); + QTest::qWait(1000); + QTRY_VERIFY_WITH_TIMEOUT(starsOut->text().toInt() >= 1, 5000); + + // Run the focus procedure for Threshold + KTRY_FOCUS_CONFIGURE("Threshold", "Iterative", 0.0, 100.0); + KTRY_FOCUS_CAPTURE(1, 1); + QWARN("No way to know when procedure is fininshed."); + QTest::qWait(1000); + QTRY_VERIFY_WITH_TIMEOUT(starsOut->text().toInt() >= 1, 5000); + + // Run the focus procedure for Gradient + KTRY_FOCUS_CONFIGURE("Gradient", "Iterative", 0.0, 100.0); + KTRY_FOCUS_CAPTURE(1, 1); + QWARN("No way to know when procedure is fininshed."); + QTest::qWait(1000); + QTRY_VERIFY_WITH_TIMEOUT(starsOut->text().toInt() >= 1, 5000); +#endif +} + +#endif diff --git a/Tests/kstars_ui/test_ekos_simulator.h b/Tests/kstars_ui/test_ekos_simulator.h --- a/Tests/kstars_ui/test_ekos_simulator.h +++ b/Tests/kstars_ui/test_ekos_simulator.h @@ -25,23 +25,26 @@ profileCBox->setCurrentText(p); \ QTRY_COMPARE(profileCBox->currentText(), p); } while(false) +#define KTRY_EKOS_GADGET(klass, name) klass * const name = Ekos::Manager::Instance()->findChild(#name); \ + QVERIFY2(name != nullptr, QString(#klass "'%1' does not exist and cannot be used").arg(#name).toStdString().c_str()) + #define KTRY_EKOS_CLICK(button) do { \ - QPushButton * const b = Ekos::Manager::Instance()->findChild(button); \ - QVERIFY2(b != nullptr, QString("QPushButton '%1' does not exist and cannot be clicked").arg(button).toStdString().c_str()); \ - QVERIFY2(b->isEnabled(), QString("QPushButton '%1' is disabled and cannot be clicked").arg(button).toStdString().c_str()); \ - QTimer::singleShot(200, Ekos::Manager::Instance(), [&]() { QTest::mouseClick(b, Qt::LeftButton); }); } while(false) + QTimer::singleShot(200, Ekos::Manager::Instance(), []() { \ + KTRY_EKOS_GADGET(QPushButton, button); \ + QVERIFY2(button->isEnabled(), QString("QPushButton '%1' is disabled and cannot be clicked").arg(#button).toStdString().c_str()); \ + QTest::mouseClick(button, Qt::LeftButton); }); } while(false) #define KTRY_EKOS_START_SIMULATORS() do { \ KTRY_EKOS_SELECT_PROFILE("Simulators"); \ - KTRY_EKOS_CLICK("processINDIB"); \ + KTRY_EKOS_CLICK(processINDIB); \ QWARN("HACK HACK HACK adding delay here for devices to connect"); \ - QTest::qWait(5000); } while(false) + QTest::qWait(10000); } while(false) #define KTRY_EKOS_STOP_SIMULATORS() do { \ - KTRY_EKOS_CLICK("disconnectB"); \ + KTRY_EKOS_CLICK(disconnectB); \ QWARN("Intentionally leaving a delay here for BZ398192"); \ QTest::qWait(5000); \ - KTRY_EKOS_CLICK("processINDIB"); \ + KTRY_EKOS_CLICK(processINDIB); \ QTest::qWait(1000); } while(false) class TestEkosSimulator : public QObject diff --git a/kstars/CMakeLists.txt b/kstars/CMakeLists.txt --- a/kstars/CMakeLists.txt +++ b/kstars/CMakeLists.txt @@ -261,6 +261,11 @@ fitsviewer/fitshistogram.cpp fitsviewer/fitsview.cpp fitsviewer/fitsdata.cpp + fitsviewer/fitsstardetector.cpp + fitsviewer/fitsthresholddetector.cpp + fitsviewer/fitsgradientdetector.cpp + fitsviewer/fitscentroiddetector.cpp + fitsviewer/fitssepdetector.cpp ) set (fitsui_SRCS fitsviewer/fitsheaderdialog.ui diff --git a/kstars/ekos/auxiliary/darklibrary.cpp b/kstars/ekos/auxiliary/darklibrary.cpp --- a/kstars/ekos/auxiliary/darklibrary.cpp +++ b/kstars/ekos/auxiliary/darklibrary.cpp @@ -282,13 +282,13 @@ FITSData *lightData = lightImage->getImageData(); - T *lightBuffer = reinterpret_cast(lightData->getImageBuffer()); + T *lightBuffer = reinterpret_cast(lightData->getWritableImageBuffer()); int lightW = lightData->width(); int lightH = lightData->height(); int darkW = darkData->width(); int darkoffset = offsetX + offsetY * darkW; - T *darkBuffer = reinterpret_cast(darkData->getImageBuffer()) + darkoffset; + T const *darkBuffer = reinterpret_cast(darkData->getImageBuffer()) + darkoffset; for (int i = 0; i < lightH; i++) { diff --git a/kstars/ekos/focus/focus.h b/kstars/ekos/focus/focus.h --- a/kstars/ekos/focus/focus.h +++ b/kstars/ekos/focus/focus.h @@ -436,6 +436,12 @@ */ void syncTrackingBoxPosition(); + /** @internal Search for stars using the method currently configured. + * @param image_data is the FITS frame to work with. + * @note Updates the frame after search. + */ + double searchForStars(FITSData *image_data); + /// Focuser device needed for focus operation ISD::Focuser *currentFocuser { nullptr }; /// CCD device needed for focus operation diff --git a/kstars/ekos/focus/focus.cpp b/kstars/ekos/focus/focus.cpp --- a/kstars/ekos/focus/focus.cpp +++ b/kstars/ekos/focus/focus.cpp @@ -1049,6 +1049,61 @@ setCaptureComplete(); } +double Focus::searchForStars(FITSData *image_data) +{ + // When we're using FULL field view, we always use either CENTROID algorithm which is the default + // standard algorithm in KStars, or SEP. The other algorithms are too inefficient to run on full frames and require + // a bounding box for them to be effective in near real-time application. + if (Options::focusUseFullField()) + { + Q_ASSERT_X(focusView->getTrackingBox().isNull(), __FUNCTION__, "Tracking box is disabled when detecting in full-field"); + + if (focusDetection != ALGORITHM_CENTROID && focusDetection != ALGORITHM_SEP) + focusView->findStars(ALGORITHM_CENTROID); + else + focusView->findStars(focusDetection); + + focusView->setStarFilterRange(static_cast (fullFieldInnerRing->value() / 100.0), + static_cast (fullFieldOuterRing->value() / 100.0)); + focusView->filterStars(); + focusView->updateFrame(); + + // Get the average HFR of the whole frame + return image_data->getHFR(HFR_AVERAGE); + } + else + { + // If star is already selected then use whatever algorithm currently selected. + if (starSelected) + { + focusView->findStars(focusDetection); + focusView->updateFrame(); + return image_data->getHFR(HFR_MAX); + } + else + { + // Disable tracking box + focusView->setTrackingBoxEnabled(false); + + // If algorithm is set something other than Centeroid or SEP, then force Centroid + // Since it is the most reliable detector when nothing was selected before. + if (focusDetection != ALGORITHM_CENTROID && focusDetection != ALGORITHM_SEP) + focusView->findStars(ALGORITHM_CENTROID); + else + // Otherwise, continue to find use using the selected algorithm + focusView->findStars(focusDetection); + + // Reenable tracking box + focusView->setTrackingBoxEnabled(true); + + focusView->updateFrame(); + + // Get maximum HFR in the frame + return image_data->getHFR(HFR_MAX); + } + } +} + void Focus::setCaptureComplete() { DarkLibrary::Instance()->disconnect(this); @@ -1093,62 +1148,11 @@ // First check that we haven't already search for stars // Since star-searching algorithm are time-consuming, we should only search when necessary if (image_data->areStarsSearched() == false) - { - // Reset current HFR - currentHFR = -1; - - // When we're using FULL field view, we always use either CENTROID algorithm which is the default - // standard algorithm in KStars, or SEP. The other algorithms are too inefficient to run on full frames and require - // a bounding box for them to be effective in near real-time application. - if (Options::focusUseFullField()) - { - if (focusDetection != ALGORITHM_CENTROID && focusDetection != ALGORITHM_SEP) - focusView->findStars(ALGORITHM_CENTROID); - else - focusView->findStars(focusDetection); - focusView->setStarFilterRange(static_cast (fullFieldInnerRing->value() / 100.0), - static_cast (fullFieldOuterRing->value() / 100.0)); - focusView->filterStars(); - focusView->updateFrame(); - - // Get the average HFR of the whole frame - currentHFR = image_data->getHFR(HFR_AVERAGE); - } - else - { - // If star is already selected then use whatever algorithm currently selected. - if (starSelected) - { - focusView->findStars(focusDetection); - focusView->updateFrame(); - currentHFR = image_data->getHFR(HFR_MAX); - } - else - { - // Disable tracking box - focusView->setTrackingBoxEnabled(false); - - // If algorithm is set something other than Centeroid or SEP, then force Centroid - // Since it is the most reliable detector when nothing was selected before. - if (focusDetection != ALGORITHM_CENTROID && focusDetection != ALGORITHM_SEP) - focusView->findStars(ALGORITHM_CENTROID); - else - // Otherwise, continue to find use using the selected algorithm - focusView->findStars(focusDetection); - - // Reenable tracking box - focusView->setTrackingBoxEnabled(true); - - focusView->updateFrame(); - - // Get maximum HFR in the frame - currentHFR = image_data->getHFR(HFR_MAX); - } - } - } + currentHFR = searchForStars(image_data); // Let's now report the current HFR qCDebug(KSTARS_EKOS_FOCUS) << "Focus newFITS #" << HFRFrames.count() + 1 << ": Current HFR " << currentHFR << " Num stars " << (starSelected ? 1 : image_data->getDetectedStars()); + // Add it to existing frames in case we need to take an average HFRFrames.append(currentHFR); @@ -1210,7 +1214,6 @@ return; } - // Let signal the current HFR now depending on whether the focuser is absolute or relative if (canAbsMove) emit newHFR(currentHFR, static_cast(currentPosition)); diff --git a/kstars/ekos/guide/internalguide/gmath.cpp b/kstars/ekos/guide/internalguide/gmath.cpp --- a/kstars/ekos/guide/internalguide/gmath.cpp +++ b/kstars/ekos/guide/internalguide/gmath.cpp @@ -585,63 +585,63 @@ { case TBYTE: { - uint8_t *buffer = imageData->getImageBuffer(); + uint8_t const *buffer = imageData->getImageBuffer(); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TSHORT: { - int16_t *buffer = reinterpret_cast(imageData->getImageBuffer()); + int16_t const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TUSHORT: { - uint16_t *buffer = reinterpret_cast(imageData->getImageBuffer()); + uint16_t const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TLONG: { - int32_t *buffer = reinterpret_cast(imageData->getImageBuffer()); + int32_t const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TULONG: { - uint32_t *buffer = reinterpret_cast(imageData->getImageBuffer()); + uint32_t const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TFLOAT: { - float *buffer = reinterpret_cast(imageData->getImageBuffer()); + float const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TLONGLONG: { - int64_t *buffer = reinterpret_cast(imageData->getImageBuffer()); + int64_t const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } break; case TDOUBLE: { - double *buffer = reinterpret_cast(imageData->getImageBuffer()); + double const *buffer = reinterpret_cast(imageData->getImageBuffer()); for (uint32_t i = 0; i < imgSize; i++) imgFloat[i] = buffer[i]; } @@ -820,9 +820,9 @@ Vector ret; int i, j; double resx, resy, mass, threshold, pval; - T *psrc = nullptr; - T *porigin = nullptr; - T *pptr; + T const *psrc = nullptr; + T const *porigin = nullptr; + T const *pptr; QRect trackingBox = guideView->getTrackingBox(); @@ -856,7 +856,7 @@ return ret; } - T *pdata = reinterpret_cast(imageData->getImageBuffer()); + T const *pdata = reinterpret_cast(imageData->getImageBuffer()); qCDebug(KSTARS_EKOS_GUIDE) << "Tracking Square " << trackingBox; @@ -877,7 +877,7 @@ float i0, i1, i2, i3, i4, i5, i6, i7, i8; int ix = 0, iy = 0; int xM4; - T *p; + T const *p; double average, fit, bestFit = 0; int minx = 0; int maxx = width; diff --git a/kstars/fitsviewer/fitscentroiddetector.h b/kstars/fitsviewer/fitscentroiddetector.h new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitscentroiddetector.h @@ -0,0 +1,78 @@ +/*************************************************************************** + fitscentroiddetector.h - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#ifndef FITSCENTROIDDETECTOR_H +#define FITSCENTROIDDETECTOR_H + +#include +#include "fitsstardetector.h" + +class FITSCentroidDetector: public FITSStarDetector +{ + Q_OBJECT + +public: + explicit FITSCentroidDetector(FITSData *parent): FITSStarDetector(parent) {}; + +public: + /** @brief Find sources in the parent FITS data file. + * @see FITSStarDetector::findSources(). + */ + int findSources(QList &starCenters, QRect const &boundary = QRect()) override; + + /** @brief Configure the detection method. + * @see FITSStarDetector::configure(). + * @note Parameter "initStdDev" defaults to MINIMUM_STDVAR. + * @note Parameter "minEdgeWidth" defaults to MINIMUM_PIXEL_RANGE. + * @todo Provide all constants of this class as parameters, and explain their use. + */ + FITSStarDetector & configure(const QString &setting, const QVariant &value) override; + +public: + /** @group Detection internals + * @{ */ + static constexpr int MINIMUM_STDVAR { 5 }; + static constexpr int MINIMUM_PIXEL_RANGE { 5 }; + static constexpr int MINIMUM_EDGE_LIMIT { 2 }; + static constexpr int MAX_EDGE_LIMIT { 10000 }; + static constexpr double DIFFUSE_THRESHOLD { 0.15 }; + static constexpr int MINIMUM_ROWS_PER_CENTER { 3 }; + static constexpr int LOW_EDGE_CUTOFF_1 { 50 }; + static constexpr int LOW_EDGE_CUTOFF_2 { 10 }; + /** @} */ + +protected: + /** @internal Find sources in the parent FITS data file, dependent of the pixel depth. + * @see FITSGradientDetector::findSources. + */ + template + int findSources(QList &starCenters, const QRect &boundary); + + /** @internal Check whether two sources overlap. + * @param s1, s2 are the two sources to check collision on. + * @return true if the sources collide, else false. + */ + bool checkCollision(Edge * s1, Edge * s2) const; + +protected: + int m_initStdDev { MINIMUM_STDVAR }; + int m_minEdgeWidth { MINIMUM_PIXEL_RANGE }; +}; + +#endif // FITSCENTROIDDETECTOR_H diff --git a/kstars/fitsviewer/fitscentroiddetector.cpp b/kstars/fitsviewer/fitscentroiddetector.cpp new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitscentroiddetector.cpp @@ -0,0 +1,468 @@ +/*************************************************************************** + fitscentroiddetector.cpp - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#include + +#include "fitscentroiddetector.h" +#include "fits_debug.h" + +FITSStarDetector& FITSCentroidDetector::configure(const QString &setting, const QVariant &value) +{ + if (!setting.compare("initStdDev", Qt::CaseInsensitive)) + { + bool ok = false; + double const result = value.toDouble(&ok); + if (ok) + m_initStdDev = result; + } + if (!setting.compare("minEdgeWidth", Qt::CaseInsensitive)) + { + bool ok = false; + int const result = value.toInt(&ok); + if (ok) + m_initStdDev = result; + } + + return *this; +} + +bool FITSCentroidDetector::checkCollision(Edge * s1, Edge * s2) const +{ + int dis; //distance + + int diff_x = s1->x - s2->x; + int diff_y = s1->y - s2->y; + + dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y)); + dis -= s1->width / 2; + dis -= s2->width / 2; + + if (dis <= 0) //collision + return true; + + //no collision + return false; +} + +/*** Find center of stars and calculate Half Flux Radius */ +int FITSCentroidDetector::findSources(QList &starCenters, const QRect &boundary) +{ + switch (parent()->property("dataType").toInt()) + { + case TBYTE: + return findSources(starCenters, boundary); + + case TSHORT: + return findSources(starCenters, boundary); + + case TUSHORT: + return findSources(starCenters, boundary); + + case TLONG: + return findSources(starCenters, boundary); + + case TULONG: + return findSources(starCenters, boundary); + + case TFLOAT: + return findSources(starCenters, boundary); + + case TLONGLONG: + return findSources(starCenters, boundary); + + case TDOUBLE: + return findSources(starCenters, boundary); + + default: + return -1; + } +} + +template +int FITSCentroidDetector::findSources(QList &starCenters, const QRect &boundary) +{ + FITSData const * const image_data = reinterpret_cast(parent()); + + if (image_data == nullptr) + return 0; + + FITSData::Statistic const &stats = image_data->getStatistics(); + FITSMode const m_Mode = static_cast (parent()->property("mode").toInt()); + + double threshold = 0, sum = 0, avg = 0, min = 0; + int starDiameter = 0; + int pixVal = 0; + int minimumEdgeCount = MINIMUM_EDGE_LIMIT; + + auto * buffer = reinterpret_cast(image_data->getImageBuffer()); + + double JMIndex = 100; + +#if 0//ndef KSTARS_LITE + if (histogram) + { + if (!histogram->isConstructed()) + histogram->constructHistogram(); + JMIndex = histogram->getJMIndex(); + } +#endif + + float dispersion_ratio = 1.5; + + QList edges; + + if (JMIndex < DIFFUSE_THRESHOLD) + { + m_minEdgeWidth = JMIndex * 35 + 1; + minimumEdgeCount = m_minEdgeWidth - 1; + } + else + { + m_minEdgeWidth = 6; + minimumEdgeCount = 4; + } + + while (m_initStdDev >= 1) + { + m_minEdgeWidth--; + minimumEdgeCount--; + + m_minEdgeWidth = qMax(3, m_minEdgeWidth); + minimumEdgeCount = qMax(3, minimumEdgeCount); + + if (JMIndex < DIFFUSE_THRESHOLD) + { + // Taking the average out seems to have better result for noisy images + threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - m_initStdDev) * 0.5 + 1); + + min = stats.min[0]; + if (threshold - min < 0) + { + threshold = stats.mean[0] * ((MINIMUM_STDVAR - m_initStdDev) * 0.5 + 1); + min = 0; + } + + dispersion_ratio = 1.4 - (MINIMUM_STDVAR - m_initStdDev) * 0.08; + } + else + { + threshold = stats.mean[0] + stats.stddev[0] * m_initStdDev * (0.3 - (MINIMUM_STDVAR - m_initStdDev) * 0.05); + min = stats.min[0]; + // Ratio between centeroid center and edge + dispersion_ratio = 1.8 - (MINIMUM_STDVAR - m_initStdDev) * 0.2; + } + + qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR; + qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min + << ") minimum edge width" << m_minEdgeWidth << " minimum edge limit " << minimumEdgeCount; + + threshold -= min; + + int subX, subY, subW, subH; + + if (boundary.isNull()) + { + if (m_Mode == FITS_GUIDE || m_Mode == FITS_FOCUS) + { + // Only consider the central 70% + subX = round(stats.width * 0.15); + subY = round(stats.height * 0.15); + subW = stats.width - subX; + subH = stats.height - subY; + } + else + { + // Consider the complete area 100% + subX = 0; + subY = 0; + subW = stats.width; + subH = stats.height; + } + } + else + { + subX = boundary.x(); + subY = boundary.y(); + subW = subX + boundary.width(); + subH = subY + boundary.height(); + } + + // Detect "edges" that are above threshold + for (int i = subY; i < subH; i++) + { + starDiameter = 0; + + for (int j = subX; j < subW; j++) + { + pixVal = buffer[j + (i * stats.width)] - min; + + // If pixel value > threshold, let's get its weighted average + if (pixVal >= threshold) + { + avg += j * pixVal; + sum += pixVal; + starDiameter++; + } + // Value < threshold but avg exists + else if (sum > 0) + { + // We found a potential centroid edge + if (starDiameter >= m_minEdgeWidth) + { + float center = avg / sum + 0.5; + if (center > 0) + { + int i_center = std::floor(center); + + // Check if center is 10% or more brighter than edge, if not skip + if (((buffer[i_center + (i * stats.width)] - min) / + (buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >= + dispersion_ratio) && + ((buffer[i_center + (i * stats.width)] - min) / + (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >= + dispersion_ratio)) + { + qCDebug(KSTARS_FITS) + << "Edge center is " << buffer[i_center + (i * stats.width)] - min + << " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min + << " and ratio is " + << ((buffer[i_center + (i * stats.width)] - min) / + (buffer[i_center + (i * stats.width) - starDiameter / 2] - min)) + << " located at X: " << center << " Y: " << i + 0.5; + + auto * newEdge = new Edge(); + + newEdge->x = center; + newEdge->y = i + 0.5; + newEdge->scanned = 0; + newEdge->val = buffer[i_center + (i * stats.width)] - min; + newEdge->width = starDiameter; + newEdge->HFR = 0; + newEdge->sum = sum; + + edges.append(newEdge); + } + } + } + + // Reset + avg = sum = starDiameter = 0; + } + } + } + + qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count(); + + // In case of hot pixels + if (edges.count() == 1 && m_initStdDev > 1) + { + m_initStdDev--; + continue; + } + + if (edges.count() >= MAX_EDGE_LIMIT) + { + qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count(); + qDeleteAll(edges); + return -1; + } + + if (edges.count() >= minimumEdgeCount) + break; + + qDeleteAll(edges); + edges.clear(); + m_initStdDev--; + } + + int cen_count = 0; + int cen_x = 0; + int cen_y = 0; + int cen_v = 0; + int cen_w = 0; + int width_sum = 0; + + // Let's sort edges, starting with widest + auto const greaterThan = [](Edge const *a, Edge const *b) { return a->sum > b->sum; }; + std::sort(edges.begin(), edges.end(), greaterThan); + + // Now, let's scan the edges and find the maximum centroid vertically + for (int i = 0; i < edges.count(); i++) + { + qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of " + << edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum; + + // If edge scanned already, skip + if (edges[i]->scanned == 1) + { + qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted"; + continue; + } + + qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ..."; + + // Get X, Y, and Val of edge + cen_x = edges[i]->x; + cen_y = edges[i]->y; + cen_v = edges[i]->sum; + cen_w = edges[i]->width; + + float avg_x = 0; + float avg_y = 0; + + sum = 0; + cen_count = 0; + + // Now let's compare to other edges until we hit a maxima + for (int j = 0; j < edges.count(); j++) + { + if (edges[j]->scanned) + continue; + + if (checkCollision(edges[j], edges[i])) + { + if (edges[j]->sum >= cen_v) + { + cen_v = edges[j]->sum; + cen_w = edges[j]->width; + } + + edges[j]->scanned = 1; + cen_count++; + + avg_x += edges[j]->x * edges[j]->val; + avg_y += edges[j]->y * edges[j]->val; + sum += edges[j]->val; + + continue; + } + } + + int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - m_initStdDev)); + + if (edges.count() < LOW_EDGE_CUTOFF_1) + { + if (edges.count() < LOW_EDGE_CUTOFF_2) + cen_limit = 1; + else + cen_limit = 2; + } + + qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << m_initStdDev << " and limit is " + << cen_limit; + + if (cen_limit < 1) + continue; + + // If centroid count is within acceptable range + //if (cen_limit >= 2 && cen_count >= cen_limit) + if (cen_count >= cen_limit) + { + // We detected a centroid, let's init it + auto * rCenter = new Edge(); + + rCenter->x = avg_x / sum; + rCenter->y = avg_y / sum; + width_sum += rCenter->width; + rCenter->width = cen_w; + + qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")"; + + // Calculate Total Flux From Center, Half Flux, Full Summation + double TF = 0; + double HF = 0; + double FSum = 0; + + cen_x = (int)std::floor(rCenter->x); + cen_y = (int)std::floor(rCenter->y); + + if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height) + { + delete rCenter; + continue; + } + + // Complete sum along the radius + //for (int k=0; k < rCenter->width; k++) + for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--) + { + FSum += buffer[cen_x - k + (cen_y * stats.width)] - min; + //qDebug() << image_buffer[cen_x-k+(cen_y*stats.width)] - min; + } + + // Half flux + HF = FSum / 2.0; + + // Total flux starting from center + TF = buffer[cen_y * stats.width + cen_x] - min; + + int pixelCounter = 1; + + // Integrate flux along radius axis until we reach half flux + for (int k = 1; k < rCenter->width / 2; k++) + { + if (TF >= HF) + { + qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels."; + break; + } + + TF += buffer[cen_y * stats.width + cen_x + k] - min; + TF += buffer[cen_y * stats.width + cen_x - k] - min; + + pixelCounter++; + } + + // Calculate weighted Half Flux Radius + rCenter->HFR = pixelCounter * (HF / TF); + // Store full flux + rCenter->val = FSum; + + qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum; + + starCenters.append(rCenter); + } + } + + if (starCenters.count() > 1 && m_Mode != FITS_FOCUS) + { + float width_avg = (float)width_sum / starCenters.count(); + float lsum = 0, sdev = 0; + + for (auto ¢er : starCenters) + lsum += (center->width - width_avg) * (center->width - width_avg); + + sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4; + + // Reject stars > 4 * stddev + foreach (Edge * center, starCenters) + if (center->width > sdev) + starCenters.removeOne(center); + + //foreach(Edge *center, starCenters) + //qDebug() << center->x << "," << center->y << "," << center->width << "," << center->val << endl; + } + + // Release memory + qDeleteAll(edges); + + return starCenters.count(); +} + + diff --git a/kstars/fitsviewer/fitsdata.h b/kstars/fitsviewer/fitsdata.h --- a/kstars/fitsviewer/fitsdata.h +++ b/kstars/fitsviewer/fitsdata.h @@ -23,6 +23,7 @@ #include "bayer.h" #include "fitscommon.h" +#include "fitsstardetector.h" #ifdef WIN32 // This header must be included before fitsio.h to avoid compiler errors with Visual Studio @@ -43,9 +44,6 @@ #endif #endif -#define MINIMUM_PIXEL_RANGE 5 -#define MINIMUM_STDVAR 5 - class QProgressDialog; class SkyObject; @@ -58,17 +56,7 @@ float dec; } wcs_point; -class Edge -{ - public: - float x; - float y; - int val; - int scanned; - float width; - float HFR; - float sum; -}; +class Edge; class FITSSkyObject : public QObject { @@ -172,11 +160,13 @@ // Access functions void clearImageBuffers(); void setImageBuffer(uint8_t *buffer); - uint8_t *getImageBuffer(); + uint8_t const *getImageBuffer() const; + uint8_t *getWritableImageBuffer(); // Statistics void saveStatistics(Statistic &other); void restoreStatistics(Statistic &other); + Statistic const &getStatistics() const { return stats; }; uint16_t width() const { @@ -287,18 +277,10 @@ int findStars(StarAlgorithm algorithm = ALGORITHM_CENTROID, const QRect &trackingBox = QRect()); - void getCenterSelection(int *x, int *y); - int findOneStar(const QRect &boundary); - - // Star Detection - Partially customized Canny edge detection algorithm - static int findCannyStar(FITSData *data, const QRect &boundary = QRect()); - template - static int findCannyStar(FITSData *data, const QRect &boundary); - // Use SEP (Sextractor Library) to find stars template - void getFloatBuffer(float *buffer, int x, int y, int w, int h); - int findSEPStars(const QRect &boundary = QRect()); + void getFloatBuffer(float *buffer, int x, int y, int w, int h) const; + int findSEPStars(QList&, const QRect &boundary = QRect()) const; // Apply ring filter to searched stars int filterStars(const float innerRadius, const float outerRadius); @@ -443,7 +425,6 @@ void loadCommon(const QString &inFilename); bool privateLoad(void *fits_buffer, size_t fits_buffer_size, bool silent); void rotWCSFITS(int angle, int mirror); - bool checkCollision(Edge *s1, Edge *s2); int calculateMinMax(bool refresh = false); bool checkDebayer(); void readWCSKeys(); @@ -462,14 +443,6 @@ // Apply Filter template void applyFilter(FITSScale type, uint8_t *targetImage, QVector * min = nullptr, QVector * max = nullptr); - // Star Detect - Centroid - template - int findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth); - int findCentroid(const QRect &boundary = QRect(), int initStdDev = MINIMUM_STDVAR, - int minEdgeWidth = MINIMUM_PIXEL_RANGE); - // Star Detect - Threshold - template - int findOneStar(const QRect &boundary); template void calculateMinMax(); @@ -483,17 +456,9 @@ template QPair getSquaredSumAndMean(uint32_t start, uint32_t stride); - // Sobel detector by Gonzalo Exequiel Pedone - template - void sobel(QVector &gradient, QVector &direction); - template void convertToQImage(double dataMin, double dataMax, double scale, double zero, QImage &image); - // Give unique IDs to each contiguous region - int partition(int width, int height, QVector &gradient, QVector &ids); - void trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y); - #ifndef KSTARS_LITE FITSHistogram *histogram { nullptr }; // Pointer to the FITS data histogram #endif diff --git a/kstars/fitsviewer/fitsdata.cpp b/kstars/fitsviewer/fitsdata.cpp --- a/kstars/fitsviewer/fitsdata.cpp +++ b/kstars/fitsviewer/fitsdata.cpp @@ -18,8 +18,11 @@ ***************************************************************************/ #include "fitsdata.h" +#include "fitsthresholddetector.h" +#include "fitsgradientdetector.h" +#include "fitscentroiddetector.h" +#include "fitssepdetector.h" -#include "sep/sep.h" #include "fpack.h" #include "kstarsdata.h" @@ -54,21 +57,9 @@ #define ZOOM_LOW_INCR 10 #define ZOOM_HIGH_INCR 50 -const int MINIMUM_ROWS_PER_CENTER = 3; const QString FITSData::m_TemporaryPath = QStandardPaths::writableLocation(QStandardPaths::TempLocation); -#define DIFFUSE_THRESHOLD 0.15 -#define MAX_EDGE_LIMIT 10000 -#define LOW_EDGE_CUTOFF_1 50 -#define LOW_EDGE_CUTOFF_2 10 -#define MINIMUM_EDGE_LIMIT 2 - -bool greaterThan(Edge * s1, Edge * s2) -{ - //return s1->width > s2->width; - return s1->sum > s2->sum; -} FITSData::FITSData(FITSMode fitsMode): m_Mode(fitsMode) { @@ -880,59 +871,6 @@ return false; } -bool FITSData::checkCollision(Edge * s1, Edge * s2) -{ - int dis; //distance - - int diff_x = s1->x - s2->x; - int diff_y = s1->y - s2->y; - - dis = std::abs(sqrt(diff_x * diff_x + diff_y * diff_y)); - dis -= s1->width / 2; - dis -= s2->width / 2; - - if (dis <= 0) //collision - return true; - - //no collision - return false; -} - -int FITSData::findCannyStar(FITSData * data, const QRect &boundary) -{ - switch (data->property("dataType").toInt()) - { - case TBYTE: - return FITSData::findCannyStar(data, boundary); - - case TSHORT: - return FITSData::findCannyStar(data, boundary); - - case TUSHORT: - return FITSData::findCannyStar(data, boundary); - - case TLONG: - return FITSData::findCannyStar(data, boundary); - - case TULONG: - return FITSData::findCannyStar(data, boundary); - - case TFLOAT: - return FITSData::findCannyStar(data, boundary); - - case TLONGLONG: - return FITSData::findCannyStar(data, boundary); - - case TDOUBLE: - return FITSData::findCannyStar(data, boundary); - - default: - break; - } - - return 0; -} - int FITSData::findStars(StarAlgorithm algorithm, const QRect &trackingBox) { int count = 0; @@ -944,19 +882,24 @@ switch (algorithm) { case ALGORITHM_SEP: - count = findSEPStars(trackingBox); + count = FITSSEPDetector(this) + .findSources(starCenters, trackingBox); break; case ALGORITHM_GRADIENT: - count = findCannyStar(this, trackingBox); + count = FITSGradientDetector(this) + .findSources(starCenters, trackingBox); break; case ALGORITHM_CENTROID: - count = findCentroid(trackingBox); + count = FITSCentroidDetector(this) + .findSources(starCenters, trackingBox); break; case ALGORITHM_THRESHOLD: - count = findOneStar(trackingBox); + count = FITSThresholdDetector(this) + .configure("Threshold", Options::focusThreshold()) + .findSources(starCenters, trackingBox); break; } @@ -983,814 +926,6 @@ return starCenters.count(); } -template -int FITSData::findCannyStar(FITSData * data, const QRect &boundary) -{ - int subX = qMax(0, boundary.isNull() ? 0 : boundary.x()); - int subY = qMax(0, boundary.isNull() ? 0 : boundary.y()); - int subW = (boundary.isNull() ? data->width() : boundary.width()); - int subH = (boundary.isNull() ? data->height() : boundary.height()); - - int BBP = data->getBytesPerPixel(); - - uint16_t dataWidth = data->width(); - - // #1 Find offsets - uint32_t size = subW * subH; - uint32_t offset = subX + subY * dataWidth; - - // #2 Create new buffer - auto * buffer = new uint8_t[size * BBP]; - // If there is no offset, copy whole buffer in one go - if (offset == 0) - memcpy(buffer, data->getImageBuffer(), size * BBP); - else - { - uint8_t * dataPtr = buffer; - uint8_t * origDataPtr = data->getImageBuffer(); - uint32_t lineOffset = 0; - // Copy data line by line - for (int height = subY; height < (subY + subH); height++) - { - lineOffset = (subX + height * dataWidth) * BBP; - memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP); - dataPtr += (subW * BBP); - } - } - - // #3 Create new FITSData to hold it - auto * boundedImage = new FITSData(); - boundedImage->stats.width = subW; - boundedImage->stats.height = subH; - boundedImage->stats.bitpix = data->stats.bitpix; - boundedImage->stats.bytesPerPixel = data->stats.bytesPerPixel; - boundedImage->stats.samples_per_channel = size; - boundedImage->stats.ndim = 2; - - boundedImage->setProperty("dataType", data->property("dataType")); - - // #4 Set image buffer and calculate stats. - boundedImage->setImageBuffer(buffer); - - boundedImage->calculateStats(true); - - // #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain - boundedImage->applyFilter(FITS_MEDIAN); - boundedImage->applyFilter(FITS_HIGH_CONTRAST); - - // #6 Perform Sobel to find gradients and their directions - QVector gradients; - QVector directions; - - // TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed - // and discarded whenever necessary. It won't work on noisy images unless this is done. - boundedImage->sobel(gradients, directions); - - QVector ids(gradients.size()); - - int maxID = boundedImage->partition(subW, subH, gradients, ids); - - // Not needed anymore - delete boundedImage; - - if (maxID == 0) - return 0; - - typedef struct - { - float massX = 0; - float massY = 0; - float totalMass = 0; - } massInfo; - - QMap masses; - - // #7 Calculate center of mass for all detected regions - for (int y = 0; y < subH; y++) - { - for (int x = 0; x < subW; x++) - { - int index = x + y * subW; - - int regionID = ids[index]; - if (regionID > 0) - { - float pixel = gradients[index]; - - masses[regionID].totalMass += pixel; - masses[regionID].massX += x * pixel; - masses[regionID].massY += y * pixel; - } - } - } - - // Compare multiple masses, and only select the highest total mass one as the desired star - int maxRegionID = 1; - int maxTotalMass = masses[1].totalMass; - double totalMassRatio = 1e6; - for (auto key : masses.keys()) - { - massInfo oneMass = masses.value(key); - if (oneMass.totalMass > maxTotalMass) - { - totalMassRatio = oneMass.totalMass / maxTotalMass; - maxTotalMass = oneMass.totalMass; - maxRegionID = key; - } - } - - // If image has many regions and there is no significant relative center of mass then it's just noise and no stars - // are probably there above a useful threshold. - if (maxID > 10 && totalMassRatio < 1.5) - return 0; - - auto * center = new Edge; - center->width = -1; - center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5; - center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5; - center->HFR = 1; - - // Maximum Radius - int maxR = qMin(subW - 1, subH - 1) / 2; - - for (int r = maxR; r > 1; r--) - { - int pass = 0; - - for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0) - { - int testX = center->x + std::cos(theta) * r; - int testY = center->y + std::sin(theta) * r; - - // if out of bound, break; - if (testX < 0 || testX >= subW || testY < 0 || testY >= subH) - break; - - if (gradients[testX + testY * subW] > 0) - //if (thresholded[testX + testY * subW] > 0) - { - if (++pass >= 24) - { - center->width = r * 2; - // Break of outer loop - r = 0; - break; - } - } - } - } - - qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width; - - // If no stars were detected - if (center->width == -1) - { - delete center; - return 0; - } - - // 30% fuzzy - //center->width += center->width*0.3 * (running_threshold / threshold); - - double FSum = 0, HF = 0, TF = 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 subPixels; - subPixels.reserve(center->width / resolution); - - const T * origBuffer = reinterpret_cast(data->getImageBuffer()) + offset; - - for (double x = leftEdge; x <= rightEdge; x += resolution) - { - double slice = resolution * (origBuffer[static_cast(floor(x)) + cen_y * dataWidth]); - FSum += slice; - subPixels.append(slice); - } - - // Half flux - HF = FSum / 2.0; - - 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 overpassed HF, let's calculate from last TF how much until we reach HF - - // #1 Accurate calculation, but very sensitive to small variations of flux - center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution; - - // #2 Less accurate calculation, but stable against small variations of flux - //center->HFR = (k - 1) * resolution; - break; - } - - lastTF = TF; - } - - // Correct center for subX and subY - center->x += subX; - center->y += subY; - - data->appendStar(center); - - qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR; - - return 1; -} - -int FITSData::findOneStar(const QRect &boundary) -{ - switch (m_DataType) - { - case TBYTE: - return findOneStar(boundary); - - case TSHORT: - return findOneStar(boundary); - - case TUSHORT: - return findOneStar(boundary); - - case TLONG: - return findOneStar(boundary); - - case TULONG: - return findOneStar(boundary); - - case TFLOAT: - return findOneStar(boundary); - - case TLONGLONG: - return findOneStar(boundary); - - case TDOUBLE: - return findOneStar(boundary); - - default: - break; - } - - return 0; -} - -template -int FITSData::findOneStar(const QRect &boundary) -{ - if (boundary.isEmpty()) - return -1; - - 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(m_ImageBuffer); - - // TODO replace magic number with something more useful to understand - double threshold = stats.mean[0] * Options::focusThreshold() / 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() << "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 0; - } - - // 30% fuzzy - //center->width += center->width*0.3 * (running_threshold / threshold); - - starCenters.append(center); - - 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 subPixels; - subPixels.reserve(center->width / resolution); - - for (double x = leftEdge; x <= rightEdge; x += resolution) - { - //subPixels[x] = resolution * (image_buffer[static_cast(floor(x)) + cen_y * stats.width] - min); - double slice = resolution * (buffer[static_cast(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; - } - - return 1; -} - -/*** Find center of stars and calculate Half Flux Radius */ -int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth) -{ - switch (m_DataType) - { - case TBYTE: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TSHORT: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TUSHORT: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TLONG: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TULONG: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TFLOAT: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TLONGLONG: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - case TDOUBLE: - return findCentroid(boundary, initStdDev, minEdgeWidth); - - default: - return -1; - } -} - -template -int FITSData::findCentroid(const QRect &boundary, int initStdDev, int minEdgeWidth) -{ - double threshold = 0, sum = 0, avg = 0, min = 0; - int starDiameter = 0; - int pixVal = 0; - int minimumEdgeCount = MINIMUM_EDGE_LIMIT; - - auto * buffer = reinterpret_cast(m_ImageBuffer); - - double JMIndex = 100; -#ifndef KSTARS_LITE - if (histogram) - { - if (!histogram->isConstructed()) - histogram->constructHistogram(); - JMIndex = histogram->getJMIndex(); - } -#endif - - float dispersion_ratio = 1.5; - - QList edges; - - if (JMIndex < DIFFUSE_THRESHOLD) - { - minEdgeWidth = JMIndex * 35 + 1; - minimumEdgeCount = minEdgeWidth - 1; - } - else - { - minEdgeWidth = 6; - minimumEdgeCount = 4; - } - - while (initStdDev >= 1) - { - minEdgeWidth--; - minimumEdgeCount--; - - minEdgeWidth = qMax(3, minEdgeWidth); - minimumEdgeCount = qMax(3, minimumEdgeCount); - - if (JMIndex < DIFFUSE_THRESHOLD) - { - // Taking the average out seems to have better result for noisy images - threshold = stats.max[0] - stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); - - min = stats.min[0]; - if (threshold - min < 0) - { - threshold = stats.mean[0] * ((MINIMUM_STDVAR - initStdDev) * 0.5 + 1); - min = 0; - } - - dispersion_ratio = 1.4 - (MINIMUM_STDVAR - initStdDev) * 0.08; - } - else - { - threshold = stats.mean[0] + stats.stddev[0] * initStdDev * (0.3 - (MINIMUM_STDVAR - initStdDev) * 0.05); - min = stats.min[0]; - // Ratio between centeroid center and edge - dispersion_ratio = 1.8 - (MINIMUM_STDVAR - initStdDev) * 0.2; - } - - qCDebug(KSTARS_FITS) << "SNR: " << stats.SNR; - qCDebug(KSTARS_FITS) << "The threshold level is " << threshold << "(actual " << threshold - min - << ") minimum edge width" << minEdgeWidth << " minimum edge limit " << minimumEdgeCount; - - threshold -= min; - - int subX, subY, subW, subH; - - if (boundary.isNull()) - { - if (m_Mode == FITS_GUIDE || m_Mode == FITS_FOCUS) - { - // Only consider the central 70% - subX = round(stats.width * 0.15); - subY = round(stats.height * 0.15); - subW = stats.width - subX; - subH = stats.height - subY; - } - else - { - // Consider the complete area 100% - subX = 0; - subY = 0; - subW = stats.width; - subH = stats.height; - } - } - else - { - subX = boundary.x(); - subY = boundary.y(); - subW = subX + boundary.width(); - subH = subY + boundary.height(); - } - - // Detect "edges" that are above threshold - for (int i = subY; i < subH; i++) - { - starDiameter = 0; - - for (int j = subX; j < subW; j++) - { - pixVal = buffer[j + (i * stats.width)] - min; - - // If pixel value > threshold, let's get its weighted average - if (pixVal >= threshold) - { - avg += j * pixVal; - sum += pixVal; - starDiameter++; - } - // Value < threshold but avg exists - else if (sum > 0) - { - // We found a potential centroid edge - if (starDiameter >= minEdgeWidth) - { - float center = avg / sum + 0.5; - if (center > 0) - { - int i_center = std::floor(center); - - // Check if center is 10% or more brighter than edge, if not skip - if (((buffer[i_center + (i * stats.width)] - min) / - (buffer[i_center + (i * stats.width) - starDiameter / 2] - min) >= - dispersion_ratio) && - ((buffer[i_center + (i * stats.width)] - min) / - (buffer[i_center + (i * stats.width) + starDiameter / 2] - min) >= - dispersion_ratio)) - { - qCDebug(KSTARS_FITS) - << "Edge center is " << buffer[i_center + (i * stats.width)] - min - << " Edge is " << buffer[i_center + (i * stats.width) - starDiameter / 2] - min - << " and ratio is " - << ((buffer[i_center + (i * stats.width)] - min) / - (buffer[i_center + (i * stats.width) - starDiameter / 2] - min)) - << " located at X: " << center << " Y: " << i + 0.5; - - auto * newEdge = new Edge(); - - newEdge->x = center; - newEdge->y = i + 0.5; - newEdge->scanned = 0; - newEdge->val = buffer[i_center + (i * stats.width)] - min; - newEdge->width = starDiameter; - newEdge->HFR = 0; - newEdge->sum = sum; - - edges.append(newEdge); - } - } - } - - // Reset - avg = sum = starDiameter = 0; - } - } - } - - qCDebug(KSTARS_FITS) << "Total number of edges found is: " << edges.count(); - - // In case of hot pixels - if (edges.count() == 1 && initStdDev > 1) - { - initStdDev--; - continue; - } - - if (edges.count() >= MAX_EDGE_LIMIT) - { - qCWarning(KSTARS_FITS) << "Too many edges, aborting... " << edges.count(); - qDeleteAll(edges); - return -1; - } - - if (edges.count() >= minimumEdgeCount) - break; - - qDeleteAll(edges); - edges.clear(); - initStdDev--; - } - - int cen_count = 0; - int cen_x = 0; - int cen_y = 0; - int cen_v = 0; - int cen_w = 0; - int width_sum = 0; - - // Let's sort edges, starting with widest - std::sort(edges.begin(), edges.end(), greaterThan); - - // Now, let's scan the edges and find the maximum centroid vertically - for (int i = 0; i < edges.count(); i++) - { - qCDebug(KSTARS_FITS) << "# " << i << " Edge at (" << edges[i]->x << "," << edges[i]->y << ") With a value of " - << edges[i]->val << " and width of " << edges[i]->width << " pixels. with sum " << edges[i]->sum; - - // If edge scanned already, skip - if (edges[i]->scanned == 1) - { - qCDebug(KSTARS_FITS) << "Skipping check for center " << i << " because it was already counted"; - continue; - } - - qCDebug(KSTARS_FITS) << "Investigating edge # " << i << " now ..."; - - // Get X, Y, and Val of edge - cen_x = edges[i]->x; - cen_y = edges[i]->y; - cen_v = edges[i]->sum; - cen_w = edges[i]->width; - - float avg_x = 0; - float avg_y = 0; - - sum = 0; - cen_count = 0; - - // Now let's compare to other edges until we hit a maxima - for (int j = 0; j < edges.count(); j++) - { - if (edges[j]->scanned) - continue; - - if (checkCollision(edges[j], edges[i])) - { - if (edges[j]->sum >= cen_v) - { - cen_v = edges[j]->sum; - cen_w = edges[j]->width; - } - - edges[j]->scanned = 1; - cen_count++; - - avg_x += edges[j]->x * edges[j]->val; - avg_y += edges[j]->y * edges[j]->val; - sum += edges[j]->val; - - continue; - } - } - - int cen_limit = (MINIMUM_ROWS_PER_CENTER - (MINIMUM_STDVAR - initStdDev)); - - if (edges.count() < LOW_EDGE_CUTOFF_1) - { - if (edges.count() < LOW_EDGE_CUTOFF_2) - cen_limit = 1; - else - cen_limit = 2; - } - - qCDebug(KSTARS_FITS) << "center_count: " << cen_count << " and initstdDev= " << initStdDev << " and limit is " - << cen_limit; - - if (cen_limit < 1) - continue; - - // If centroid count is within acceptable range - //if (cen_limit >= 2 && cen_count >= cen_limit) - if (cen_count >= cen_limit) - { - // We detected a centroid, let's init it - auto * rCenter = new Edge(); - - rCenter->x = avg_x / sum; - rCenter->y = avg_y / sum; - width_sum += rCenter->width; - rCenter->width = cen_w; - - qCDebug(KSTARS_FITS) << "Found a real center with number with (" << rCenter->x << "," << rCenter->y << ")"; - - // Calculate Total Flux From Center, Half Flux, Full Summation - double TF = 0; - double HF = 0; - double FSum = 0; - - cen_x = (int)std::floor(rCenter->x); - cen_y = (int)std::floor(rCenter->y); - - if (cen_x < 0 || cen_x > stats.width || cen_y < 0 || cen_y > stats.height) - { - delete rCenter; - continue; - } - - // Complete sum along the radius - //for (int k=0; k < rCenter->width; k++) - for (int k = rCenter->width / 2; k >= -(rCenter->width / 2); k--) - { - FSum += buffer[cen_x - k + (cen_y * stats.width)] - min; - //qDebug() << image_buffer[cen_x-k+(cen_y*stats.width)] - min; - } - - // Half flux - HF = FSum / 2.0; - - // Total flux starting from center - TF = buffer[cen_y * stats.width + cen_x] - min; - - int pixelCounter = 1; - - // Integrate flux along radius axis until we reach half flux - for (int k = 1; k < rCenter->width / 2; k++) - { - if (TF >= HF) - { - qCDebug(KSTARS_FITS) << "Stopping at TF " << TF << " after #" << k << " pixels."; - break; - } - - TF += buffer[cen_y * stats.width + cen_x + k] - min; - TF += buffer[cen_y * stats.width + cen_x - k] - min; - - pixelCounter++; - } - - // Calculate weighted Half Flux Radius - rCenter->HFR = pixelCounter * (HF / TF); - // Store full flux - rCenter->val = FSum; - - qCDebug(KSTARS_FITS) << "HFR for this center is " << rCenter->HFR << " pixels and the total flux is " << FSum; - - starCenters.append(rCenter); - } - } - - if (starCenters.count() > 1 && m_Mode != FITS_FOCUS) - { - float width_avg = (float)width_sum / starCenters.count(); - float lsum = 0, sdev = 0; - - for (auto ¢er : starCenters) - lsum += (center->width - width_avg) * (center->width - width_avg); - - sdev = (std::sqrt(lsum / (starCenters.count() - 1))) * 4; - - // Reject stars > 4 * stddev - foreach (Edge * center, starCenters) - if (center->width > sdev) - starCenters.removeOne(center); - - //foreach(Edge *center, starCenters) - //qDebug() << center->x << "," << center->y << "," << center->width << "," << center->val << endl; - } - - // Release memory - qDeleteAll(edges); - - return starCenters.count(); -} - double FITSData::getHFR(HFRType type) { // This method is less susceptible to noise @@ -2297,27 +1432,6 @@ return starCentersInSubFrame; } -void FITSData::getCenterSelection(int * x, int * y) -{ - if (starCenters.count() == 0) - return; - - auto * pEdge = new Edge(); - pEdge->x = *x; - pEdge->y = *y; - pEdge->width = 1; - - foreach (Edge * center, starCenters) - if (checkCollision(pEdge, center)) - { - *x = static_cast(center->x); - *y = static_cast(center->y); - break; - } - - delete (pEdge); -} - bool FITSData::checkForWCS() { #ifndef KSTARS_LITE @@ -3254,7 +2368,12 @@ } -uint8_t * FITSData::getImageBuffer() +uint8_t * FITSData::getWritableImageBuffer() +{ + return m_ImageBuffer; +} + +uint8_t const * FITSData::getImageBuffer() const { return m_ImageBuffer; } @@ -3535,161 +2654,6 @@ return (adu / static_cast(m_Channels)); } -/* CannyDetector, Implementation of Canny edge detector in Qt/C++. - * Copyright (C) 2015 Gonzalo Exequiel Pedone - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - * Email : hipersayan DOT x AT gmail DOT com - * Web-Site: https://github.com/hipersayanX/CannyDetector - */ - -template -void FITSData::sobel(QVector &gradient, QVector &direction) -{ - //int size = image.width() * image.height(); - gradient.resize(stats.samples_per_channel); - direction.resize(stats.samples_per_channel); - - for (int y = 0; y < stats.height; y++) - { - size_t yOffset = y * stats.width; - const T * grayLine = reinterpret_cast(m_ImageBuffer) + yOffset; - - const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width; - const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width; - - float * gradientLine = gradient.data() + yOffset; - float * directionLine = direction.data() + yOffset; - - for (int x = 0; x < stats.width; x++) - { - int x_m1 = x < 1 ? x : x - 1; - int x_p1 = x >= stats.width - 1 ? x : x + 1; - - int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] - - 2 * grayLine[x_m1] - grayLine_p1[x_m1]; - - int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] - - 2 * grayLine_p1[x] - grayLine_p1[x_p1]; - - gradientLine[x] = qAbs(gradX) + qAbs(gradY); - - /* Gradient directions are classified in 4 possible cases - * - * dir 0 - * - * x x x - * - - - - * x x x - * - * dir 1 - * - * x x / - * x / x - * / x x - * - * dir 2 - * - * \ x x - * x \ x - * x x \ - * - * dir 3 - * - * x | x - * x | x - * x | x - */ - if (gradX == 0 && gradY == 0) - directionLine[x] = 0; - else if (gradX == 0) - directionLine[x] = 3; - else - { - qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI; - - if (a >= -22.5 && a < 22.5) - directionLine[x] = 0; - else if (a >= 22.5 && a < 67.5) - directionLine[x] = 2; - else if (a >= -67.5 && a < -22.5) - directionLine[x] = 1; - else - directionLine[x] = 3; - } - } - } -} - -int FITSData::partition(int width, int height, QVector &gradient, QVector &ids) -{ - int id = 0; - - for (int y = 1; y < height - 1; y++) - { - for (int x = 1; x < width - 1; x++) - { - int index = x + y * width; - float val = gradient[index]; - if (val > 0 && ids[index] == 0) - { - trace(width, height, ++id, gradient, ids, x, y); - } - } - } - - // Return max id - return id; -} - -void FITSData::trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y) -{ - int yOffset = y * width; - float * cannyLine = image.data() + yOffset; - int * idLine = ids.data() + yOffset; - - if (idLine[x] != 0) - return; - - idLine[x] = id; - - for (int j = -1; j < 2; j++) - { - int nextY = y + j; - - if (nextY < 0 || nextY >= height) - continue; - - float * cannyLineNext = cannyLine + j * width; - - for (int i = -1; i < 2; i++) - { - int nextX = x + i; - - if (i == j || nextX < 0 || nextX >= width) - continue; - - if (cannyLineNext[nextX] > 0) - { - // Trace neighbors. - trace(width, height, id, image, ids, nextX, nextY); - } - } - } -} - QString FITSData::getLastError() const { return lastError; @@ -4213,167 +3177,6 @@ return (point.x() >= 0 && point.y() >= 0 && point.x() <= stats.width && point.y() <= stats.height); } -int FITSData::findSEPStars(const QRect &boundary) -{ - int x = 0, y = 0, w = stats.width, h = stats.height, maxRadius = 50; - - if (!boundary.isNull()) - { - x = boundary.x(); - y = boundary.y(); - w = boundary.width(); - h = boundary.height(); - maxRadius = w; - } - - auto * data = new float[w * h]; - - switch (stats.bitpix) - { - case BYTE_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case SHORT_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case USHORT_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case LONG_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case ULONG_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case FLOAT_IMG: - delete [] data; - data = reinterpret_cast(m_ImageBuffer); - break; - case LONGLONG_IMG: - getFloatBuffer(data, x, y, w, h); - break; - case DOUBLE_IMG: - getFloatBuffer(data, x, y, w, h); - break; - default: - delete [] data; - return -1; - } - - float * imback = nullptr; - double * flux = nullptr, *fluxerr = nullptr, *area = nullptr; - short * flag = nullptr; - short flux_flag = 0; - int status = 0; - sep_bkg * bkg = nullptr; - sep_catalog * catalog = nullptr; - float conv[] = {1, 2, 1, 2, 4, 2, 1, 2, 1}; - double flux_fractions[2] = {0}; - double requested_frac[2] = { 0.5, 0.99 }; - QList edges; - - // #0 Create SEP Image structure - sep_image im = {data, nullptr, nullptr, SEP_TFLOAT, 0, 0, w, h, 0.0, SEP_NOISE_NONE, 1.0, 0.0}; - - // #1 Background estimate - status = sep_background(&im, 64, 64, 3, 3, 0.0, &bkg); - if (status != 0) goto exit; - - // #2 Background evaluation - imback = (float *)malloc((w * h) * sizeof(float)); - status = sep_bkg_array(bkg, imback, SEP_TFLOAT); - if (status != 0) goto exit; - - // #3 Background subtraction - status = sep_bkg_subarray(bkg, im.data, im.dtype); - if (status != 0) goto exit; - - // #4 Source Extraction - // Note that we set deblend_cont = 1.0 to turn off deblending. - status = sep_extract(&im, 2 * bkg->globalrms, SEP_THRESH_ABS, 10, conv, 3, 3, SEP_FILTER_CONV, 32, 1.0, 1, 1.0, &catalog); - if (status != 0) goto exit; - - // TODO - // Must detect edge detection - // Must limit to brightest 100 (by flux) centers - // Should probably use ellipse to draw instead of simple circle? - // Useful for galaxies and also elenogated stars. - for (int i = 0; i < catalog->nobj; i++) - { - double flux = catalog->flux[i]; - // Get HFR - sep_flux_radius(&im, catalog->x[i], catalog->y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag); - - auto * center = new Edge(); - center->x = catalog->x[i] + x + 0.5; - center->y = catalog->y[i] + y + 0.5; - center->val = catalog->peak[i]; - center->sum = flux; - center->HFR = center->width = flux_fractions[0]; - if (flux_fractions[1] < maxRadius) - center->width = flux_fractions[1] * 2; - edges.append(center); - } - - // Let's sort edges, starting with widest - std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->width > edge2->width;}); - - // Take only the first 100 stars - { - int starCount = qMin(100, edges.count()); - for (int i = 0; i < starCount; i++) - starCenters.append(edges[i]); - } - - edges.clear(); - - qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << "#" << "#X" << "#Y" << "#Flux" << "#Width" << "#HFR"; - for (int i = 0; i < starCenters.count(); i++) - qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << i << starCenters[i]->x << starCenters[i]->y - << starCenters[i]->sum << starCenters[i]->width << starCenters[i]->HFR; - -exit: - if (stats.bitpix != FLOAT_IMG) - delete [] data; - sep_bkg_free(bkg); - sep_catalog_free(catalog); - free(imback); - free(flux); - free(fluxerr); - free(area); - free(flag); - - if (status != 0) - { - char errorMessage[512]; - sep_get_errmsg(status, errorMessage); - qCritical(KSTARS_FITS) << errorMessage; - return -1; - } - - return starCenters.count(); -} - -template -void FITSData::getFloatBuffer(float * buffer, int x, int y, int w, int h) -{ - auto * rawBuffer = reinterpret_cast(m_ImageBuffer); - - float * floatPtr = buffer; - - int x2 = x + w; - int y2 = y + h; - - for (int y1 = y; y1 < y2; y1++) - { - int offset = y1 * stats.width; - for (int x1 = x; x1 < x2; x1++) - { - *floatPtr++ = rawBuffer[offset + x1]; - } - } -} - void FITSData::saveStatistics(Statistic &other) { other = stats; diff --git a/kstars/fitsviewer/fitsgradientdetector.h b/kstars/fitsviewer/fitsgradientdetector.h new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsgradientdetector.h @@ -0,0 +1,76 @@ +/*************************************************************************** + fitsgradientdetector.h - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#ifndef FITSGRADIENTDETECTOR_H +#define FITSGRADIENTDETECTOR_H + +#include "fitsstardetector.h" + +class FITSGradientDetector: public FITSStarDetector +{ + Q_OBJECT + +public: + explicit FITSGradientDetector(FITSData *parent): FITSStarDetector(parent) {}; + +public: + /** @brief Find sources in the parent FITS data file. + * @see FITSStarDetector::findSources(). + */ + int findSources(QList &starCenters, QRect const &boundary = QRect()) override; + + /** @brief Configure the detection method. + * @see FITSStarDetector::configure(). + * @note No parameters are currently available for configuration. + */ + FITSStarDetector & configure(const QString &setting, const QVariant &value) override; + +protected: + /** @internal Find sources in the parent FITS data file, dependent of the pixel depth. + * @see FITSGradientDetector::findSources. + */ + template + int findSources(QList &starCenters, const QRect &boundary); + + /** @internal Implementation of the Canny Edge detection (CannyEdgeDetector). + * @copyright 2015 Gonzalo Exequiel Pedone (https://github.com/hipersayanX/CannyDetector). + * @param data is the FITS data to run the detection onto. + * @param gradient is the vector storing the amount of change in pixel sequences. + * @param direction is the vector storing the four directions (horizontal, vertical and two diagonals) the changes stored in 'gradient' are detected in. + */ + template + void sobel(FITSData const * data, QVector &gradient, QVector &direction) const; + + /** @internal Identify gradient connections. + * @param width, height are the dimensions of the frame to work on. + * @param gradient is the vector holding the amount of change in pixel sequences. + * @param ids is the vector storing which gradient was identified for each pixel. + */ + int partition(int width, int height, QVector &gradient, QVector &ids) const; + + /** @internal Trace gradient neighbors. + * @param width, height are the dimensions of the frame to work on. + * @param image is the image to work on, actually gradients extracted using the sobel algorithm. + * @param ids is the vector storing which gradient was identified for each pixel. + * @param x, y locate the pixel to trace from. + */ + void trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y) const; +}; + +#endif // FITSGRADIENTDETECTOR_H diff --git a/kstars/fitsviewer/fitsgradientdetector.cpp b/kstars/fitsviewer/fitsgradientdetector.cpp new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsgradientdetector.cpp @@ -0,0 +1,461 @@ +/*************************************************************************** + fitsgradientdetector.cpp - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#include + +#include "fits_debug.h" +#include "fitsgradientdetector.h" + +FITSStarDetector& FITSGradientDetector::configure(const QString &, const QVariant &) +{ + return *this; +} + +int FITSGradientDetector::findSources(QList &starCenters, const QRect &boundary) +{ + switch (parent()->property("dataType").toInt()) + { + case TBYTE: + return findSources(starCenters, boundary); + + case TSHORT: + return findSources(starCenters, boundary); + + case TUSHORT: + return findSources(starCenters, boundary); + + case TLONG: + return findSources(starCenters, boundary); + + case TULONG: + return findSources(starCenters, boundary); + + case TFLOAT: + return findSources(starCenters, boundary); + + case TLONGLONG: + return findSources(starCenters, boundary); + + case TDOUBLE: + return findSources(starCenters, boundary); + + default: + break; + } + + return 0; +} + +template +int FITSGradientDetector::findSources(QList &starCenters, const QRect &boundary) +{ + FITSData const * const data = reinterpret_cast(parent()); + + if (data == nullptr) + return 0; + + int subX = qMax(0, boundary.isNull() ? 0 : boundary.x()); + int subY = qMax(0, boundary.isNull() ? 0 : boundary.y()); + int subW = (boundary.isNull() ? data->width() : boundary.width()); + int subH = (boundary.isNull() ? data->height() : boundary.height()); + + int BBP = data->getBytesPerPixel(); + + uint16_t dataWidth = data->width(); + + // #1 Find offsets + uint32_t size = subW * subH; + uint32_t offset = subX + subY * dataWidth; + + // #2 Create new buffer + auto * buffer = new uint8_t[size * BBP]; + // If there is no offset, copy whole buffer in one go + if (offset == 0) + memcpy(buffer, data->getImageBuffer(), size * BBP); + else + { + uint8_t * dataPtr = buffer; + uint8_t const * origDataPtr = data->getImageBuffer(); + uint32_t lineOffset = 0; + // Copy data line by line + for (int height = subY; height < (subY + subH); height++) + { + lineOffset = (subX + height * dataWidth) * BBP; + memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP); + dataPtr += (subW * BBP); + } + } + + // #3 Create new FITSData to hold it + auto * boundedImage = new FITSData(); + FITSData::Statistic stats; + stats.width = subW; + stats.height = subH; + stats.bitpix = data->getStatistics().bitpix; + stats.bytesPerPixel = data->getStatistics().bytesPerPixel; + stats.samples_per_channel = size; + stats.ndim = 2; + boundedImage->restoreStatistics(stats); + + boundedImage->setProperty("dataType", parent()->property("dataType")); + + // #4 Set image buffer and calculate stats. + boundedImage->setImageBuffer(buffer); + + boundedImage->calculateStats(true); + + // #5 Apply Median + High Contrast filter to remove noise and move data to non-linear domain + boundedImage->applyFilter(FITS_MEDIAN); + boundedImage->applyFilter(FITS_HIGH_CONTRAST); + + // #6 Perform Sobel to find gradients and their directions + QVector gradients; + QVector directions; + + // TODO Must trace neighbours and assign IDs to each shape so that they can be centered massed + // and discarded whenever necessary. It won't work on noisy images unless this is done. + sobel(boundedImage, gradients, directions); + + QVector ids(gradients.size()); + + int maxID = partition(subW, subH, gradients, ids); + + // Not needed anymore + delete boundedImage; + + if (maxID == 0) + return 0; + + typedef struct + { + float massX = 0; + float massY = 0; + float totalMass = 0; + } massInfo; + + QMap masses; + + // #7 Calculate center of mass for all detected regions + for (int y = 0; y < subH; y++) + { + for (int x = 0; x < subW; x++) + { + int index = x + y * subW; + + int regionID = ids[index]; + if (regionID > 0) + { + float pixel = gradients[index]; + + masses[regionID].totalMass += pixel; + masses[regionID].massX += x * pixel; + masses[regionID].massY += y * pixel; + } + } + } + + // Compare multiple masses, and only select the highest total mass one as the desired star + int maxRegionID = 1; + int maxTotalMass = masses[1].totalMass; + double totalMassRatio = 1e6; + for (auto key : masses.keys()) + { + massInfo oneMass = masses.value(key); + if (oneMass.totalMass > maxTotalMass) + { + totalMassRatio = oneMass.totalMass / maxTotalMass; + maxTotalMass = oneMass.totalMass; + maxRegionID = key; + } + } + + // If image has many regions and there is no significant relative center of mass then it's just noise and no stars + // are probably there above a useful threshold. + if (maxID > 10 && totalMassRatio < 1.5) + return 0; + + auto * center = new Edge; + center->width = -1; + center->x = masses[maxRegionID].massX / masses[maxRegionID].totalMass + 0.5; + center->y = masses[maxRegionID].massY / masses[maxRegionID].totalMass + 0.5; + center->HFR = 1; + + // Maximum Radius + int maxR = qMin(subW - 1, subH - 1) / 2; + + for (int r = maxR; r > 1; r--) + { + int pass = 0; + + for (float theta = 0; theta < 2 * M_PI; theta += (2 * M_PI) / 36.0) + { + int testX = center->x + std::cos(theta) * r; + int testY = center->y + std::sin(theta) * r; + + // if out of bound, break; + if (testX < 0 || testX >= subW || testY < 0 || testY >= subH) + break; + + if (gradients[testX + testY * subW] > 0) + //if (thresholded[testX + testY * subW] > 0) + { + if (++pass >= 24) + { + center->width = r * 2; + // Break of outer loop + r = 0; + break; + } + } + } + } + + qCDebug(KSTARS_FITS) << "FITS: Weighted Center is X: " << center->x << " Y: " << center->y << " Width: " << center->width; + + // If no stars were detected + if (center->width == -1) + { + delete center; + return 0; + } + + // 30% fuzzy + //center->width += center->width*0.3 * (running_threshold / threshold); + + double FSum = 0, HF = 0, TF = 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 subPixels; + subPixels.reserve(center->width / resolution); + + const T * origBuffer = reinterpret_cast(data->getImageBuffer()) + offset; + + for (double x = leftEdge; x <= rightEdge; x += resolution) + { + double slice = resolution * (origBuffer[static_cast(floor(x)) + cen_y * dataWidth]); + FSum += slice; + subPixels.append(slice); + } + + // Half flux + HF = FSum / 2.0; + + 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 overpassed HF, let's calculate from last TF how much until we reach HF + + // #1 Accurate calculation, but very sensitive to small variations of flux + center->HFR = (k - 1 + ((HF - lastTF) / (TF - lastTF)) * 2) * resolution; + + // #2 Less accurate calculation, but stable against small variations of flux + //center->HFR = (k - 1) * resolution; + break; + } + + lastTF = TF; + } + + // Correct center for subX and subY + center->x += subX; + center->y += subY; + + //data->appendStar(center); + starCenters.append(center); + + qCDebug(KSTARS_FITS) << "Flux: " << FSum << " Half-Flux: " << HF << " HFR: " << center->HFR; + + return 1; +} + +/* CannyDetector, Implementation of Canny edge detector in Qt/C++. + * Copyright (C) 2015 Gonzalo Exequiel Pedone + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * Email : hipersayan DOT x AT gmail DOT com + * Web-Site: https://github.com/hipersayanX/CannyDetector + */ + +template +void FITSGradientDetector::sobel(FITSData const *data, QVector &gradient, QVector &direction) const +{ + if (data == nullptr) + return; + + FITSData::Statistic const &stats = data->getStatistics(); + + //int size = image.width() * image.height(); + gradient.resize(stats.samples_per_channel); + direction.resize(stats.samples_per_channel); + + for (int y = 0; y < stats.height; y++) + { + size_t yOffset = y * stats.width; + const T * grayLine = reinterpret_cast(data->getImageBuffer()) + yOffset; + + const T * grayLine_m1 = y < 1 ? grayLine : grayLine - stats.width; + const T * grayLine_p1 = y >= stats.height - 1 ? grayLine : grayLine + stats.width; + + float * gradientLine = gradient.data() + yOffset; + float * directionLine = direction.data() + yOffset; + + for (int x = 0; x < stats.width; x++) + { + int x_m1 = x < 1 ? x : x - 1; + int x_p1 = x >= stats.width - 1 ? x : x + 1; + + int gradX = grayLine_m1[x_p1] + 2 * grayLine[x_p1] + grayLine_p1[x_p1] - grayLine_m1[x_m1] - + 2 * grayLine[x_m1] - grayLine_p1[x_m1]; + + int gradY = grayLine_m1[x_m1] + 2 * grayLine_m1[x] + grayLine_m1[x_p1] - grayLine_p1[x_m1] - + 2 * grayLine_p1[x] - grayLine_p1[x_p1]; + + gradientLine[x] = qAbs(gradX) + qAbs(gradY); + + /* Gradient directions are classified in 4 possible cases + * + * dir 0 + * + * x x x + * - - - + * x x x + * + * dir 1 + * + * x x / + * x / x + * / x x + * + * dir 2 + * + * \ x x + * x \ x + * x x \ + * + * dir 3 + * + * x | x + * x | x + * x | x + */ + if (gradX == 0 && gradY == 0) + directionLine[x] = 0; + else if (gradX == 0) + directionLine[x] = 3; + else + { + qreal a = 180. * atan(qreal(gradY) / gradX) / M_PI; + + if (a >= -22.5 && a < 22.5) + directionLine[x] = 0; + else if (a >= 22.5 && a < 67.5) + directionLine[x] = 2; + else if (a >= -67.5 && a < -22.5) + directionLine[x] = 1; + else + directionLine[x] = 3; + } + } + } +} + +int FITSGradientDetector::partition(int width, int height, QVector &gradient, QVector &ids) const +{ + int id = 0; + + for (int y = 1; y < height - 1; y++) + { + for (int x = 1; x < width - 1; x++) + { + int index = x + y * width; + float val = gradient[index]; + if (val > 0 && ids[index] == 0) + { + trace(width, height, ++id, gradient, ids, x, y); + } + } + } + + // Return max id + return id; +} + +void FITSGradientDetector::trace(int width, int height, int id, QVector &image, QVector &ids, int x, int y) const +{ + int yOffset = y * width; + float * cannyLine = image.data() + yOffset; + int * idLine = ids.data() + yOffset; + + if (idLine[x] != 0) + return; + + idLine[x] = id; + + for (int j = -1; j < 2; j++) + { + int nextY = y + j; + + if (nextY < 0 || nextY >= height) + continue; + + float * cannyLineNext = cannyLine + j * width; + + for (int i = -1; i < 2; i++) + { + int nextX = x + i; + + if (i == j || nextX < 0 || nextX >= width) + continue; + + if (cannyLineNext[nextX] > 0) + { + // Trace neighbors. + trace(width, height, id, image, ids, nextX, nextY); + } + } + } +} diff --git a/kstars/fitsviewer/fitshistogram.cpp b/kstars/fitsviewer/fitshistogram.cpp --- a/kstars/fitsviewer/fitshistogram.cpp +++ b/kstars/fitsviewer/fitshistogram.cpp @@ -188,7 +188,7 @@ uint16_t width = imageData->width(), height = imageData->height(); uint8_t channels = imageData->channels(); - auto * buffer = reinterpret_cast(imageData->getImageBuffer()); + auto * const buffer = reinterpret_cast(imageData->getImageBuffer()); double min, max; for (int i = 0 ; i < 3; i++) @@ -527,7 +527,7 @@ { FITSData * imageData = tab->getView()->getImageData(); - uint8_t * image_buffer = imageData->getImageBuffer(); + uint8_t const * image_buffer = imageData->getImageBuffer(); int totalPixels = imageData->width() * imageData->height() * imageData->channels(); unsigned long totalBytes = totalPixels * imageData->getBytesPerPixel(); @@ -578,7 +578,7 @@ { FITSView * image = tab->getView(); FITSData * imageData = image->getImageData(); - uint8_t * image_buffer = (imageData->getImageBuffer()); + uint8_t const * image_buffer = (imageData->getImageBuffer()); int totalPixels = imageData->width() * imageData->height() * imageData->channels(); @@ -626,7 +626,7 @@ FITSView * image = tab->getView(); FITSData * imageData = image->getImageData(); - uint8_t * image_buffer = imageData->getImageBuffer(); + uint8_t const * image_buffer = imageData->getImageBuffer(); uint8_t * buffer = nullptr; unsigned int size = imageData->width() * imageData->height() * imageData->channels(); @@ -652,7 +652,7 @@ // If it's rotation of flip, no need to calculate delta if (type >= FITS_ROTATE_CW && type <= FITS_FLIP_V) { - imageData->applyFilter(type, image_buffer); + imageData->applyFilter(type); } else { diff --git a/kstars/fitsviewer/fitslabel.cpp b/kstars/fitsviewer/fitslabel.cpp --- a/kstars/fitsviewer/fitslabel.cpp +++ b/kstars/fitsviewer/fitslabel.cpp @@ -91,7 +91,7 @@ double x, y; FITSData *view_data = view->getImageData(); - uint8_t *buffer = view_data->getImageBuffer(); + uint8_t const *buffer = view_data->getImageBuffer(); if (buffer == nullptr) return; @@ -118,31 +118,31 @@ break; case TSHORT: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); break; case TUSHORT: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); break; case TLONG: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); break; case TULONG: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index]); break; case TFLOAT: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index], 'f', 5); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index], 'f', 5); break; case TLONGLONG: - stringValue = QLocale().toString(static_cast((reinterpret_cast(buffer))[index])); + stringValue = QLocale().toString(static_cast((reinterpret_cast(buffer))[index])); break; case TDOUBLE: - stringValue = QLocale().toString((reinterpret_cast(buffer))[index], 'f', 5); + stringValue = QLocale().toString((reinterpret_cast(buffer))[index], 'f', 5); break; diff --git a/kstars/fitsviewer/fitssepdetector.h b/kstars/fitsviewer/fitssepdetector.h new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitssepdetector.h @@ -0,0 +1,57 @@ +/*************************************************************************** + fitssepdetector.h - FITS Image + ------------------- + begin : Sun March 29 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#ifndef FITSSEPDETECTOR_H +#define FITSSEPDETECTOR_H + +#include + +#include "fitsstardetector.h" + +class FITSSEPDetector : public FITSStarDetector +{ + Q_OBJECT + +public: + explicit FITSSEPDetector(FITSData *parent): FITSStarDetector(parent) {}; + +public: + /** @brief Find sources in the parent FITS data file. + * @see FITSStarDetector::findSources(). + */ + int findSources(QList &starCenters, QRect const &boundary = QRect()) override; + + /** @brief Configure the detection method. + * @see FITSStarDetector::configure(). + * @note No parameters are currently available for configuration. + * @todo Provide parameters for detection configuration. + */ + FITSStarDetector & configure(const QString &setting, const QVariant &value) override; + +protected: + /** @internal Consolidate a float data buffer from FITS data. + * @param buffer is the destination float block. + * @param x, y, w, h define a (x,y)-(x+w,y+h) sub-frame to extract from the FITS data out to block 'buffer'. + * @param image_data is the FITS data block to extract from. + */ + template + void getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const * image_data) const; +}; + +#endif // FITSSEPDETECTOR_H diff --git a/kstars/fitsviewer/fitssepdetector.cpp b/kstars/fitsviewer/fitssepdetector.cpp new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitssepdetector.cpp @@ -0,0 +1,201 @@ +/*************************************************************************** + fitssepdetector.cpp - FITS Image + ------------------- + begin : Sun March 29 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#include + +#include "sep/sep.h" +#include "fits_debug.h" +#include "fitssepdetector.h" + +FITSStarDetector& FITSSEPDetector::configure(const QString &, const QVariant &) +{ + return *this; +} + +int FITSSEPDetector::findSources(QList &starCenters, const QRect &boundary) +{ + FITSData const * const image_data = reinterpret_cast(parent()); + + if (image_data == nullptr) + return 0; + + FITSData::Statistic const &stats = image_data->getStatistics(); + + int x = 0, y = 0, w = stats.width, h = stats.height, maxRadius = 50; + + if (!boundary.isNull()) + { + x = boundary.x(); + y = boundary.y(); + w = boundary.width(); + h = boundary.height(); + maxRadius = w; + } + + auto * data = new float[w * h]; + + switch (stats.bitpix) + { + case BYTE_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case SHORT_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case USHORT_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case LONG_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case ULONG_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case FLOAT_IMG: + memcpy(data, image_data->getImageBuffer(), sizeof(float)*w*h); + break; + case LONGLONG_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + case DOUBLE_IMG: + getFloatBuffer(data, x, y, w, h, image_data); + break; + default: + delete [] data; + return -1; + } + + float * imback = nullptr; + double * flux = nullptr, *fluxerr = nullptr, *area = nullptr; + short * flag = nullptr; + short flux_flag = 0; + int status = 0; + sep_bkg * bkg = nullptr; + sep_catalog * catalog = nullptr; + float conv[] = {1, 2, 1, 2, 4, 2, 1, 2, 1}; + double flux_fractions[2] = {0}; + double requested_frac[2] = { 0.5, 0.99 }; + QList edges; + + // #0 Create SEP Image structure + sep_image im = {data, nullptr, nullptr, SEP_TFLOAT, 0, 0, w, h, 0.0, SEP_NOISE_NONE, 1.0, 0.0}; + + // #1 Background estimate + status = sep_background(&im, 64, 64, 3, 3, 0.0, &bkg); + if (status != 0) goto exit; + + // #2 Background evaluation + imback = (float *)malloc((w * h) * sizeof(float)); + status = sep_bkg_array(bkg, imback, SEP_TFLOAT); + if (status != 0) goto exit; + + // #3 Background subtraction + status = sep_bkg_subarray(bkg, im.data, im.dtype); + if (status != 0) goto exit; + + // #4 Source Extraction + // Note that we set deblend_cont = 1.0 to turn off deblending. + status = sep_extract(&im, 2 * bkg->globalrms, SEP_THRESH_ABS, 10, conv, 3, 3, SEP_FILTER_CONV, 32, 1.0, 1, 1.0, &catalog); + if (status != 0) goto exit; + + // TODO + // Must detect edge detection + // Must limit to brightest 100 (by flux) centers + // Should probably use ellipse to draw instead of simple circle? + // Useful for galaxies and also elenogated stars. + for (int i = 0; i < catalog->nobj; i++) + { + double flux = catalog->flux[i]; + // Get HFR + sep_flux_radius(&im, catalog->x[i], catalog->y[i], maxRadius, 5, 0, &flux, requested_frac, 2, flux_fractions, &flux_flag); + + auto * center = new Edge(); + center->x = catalog->x[i] + x + 0.5; + center->y = catalog->y[i] + y + 0.5; + center->val = catalog->peak[i]; + center->sum = flux; + center->HFR = center->width = flux_fractions[0]; + if (flux_fractions[1] < maxRadius) + center->width = flux_fractions[1] * 2; + edges.append(center); + } + + // Let's sort edges, starting with widest + std::sort(edges.begin(), edges.end(), [](const Edge * edge1, const Edge * edge2) -> bool { return edge1->width > edge2->width;}); + + // Take only the first 100 stars + { + int starCount = qMin(100, edges.count()); + for (int i = 0; i < starCount; i++) + starCenters.append(edges[i]); + } + + edges.clear(); + + qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << "#" << "#X" << "#Y" << "#Flux" << "#Width" << "#HFR"; + for (int i = 0; i < starCenters.count(); i++) + qCDebug(KSTARS_FITS) << qSetFieldWidth(10) << i << starCenters[i]->x << starCenters[i]->y + << starCenters[i]->sum << starCenters[i]->width << starCenters[i]->HFR; + +exit: + if (stats.bitpix != FLOAT_IMG) + delete [] data; + sep_bkg_free(bkg); + sep_catalog_free(catalog); + free(imback); + free(flux); + free(fluxerr); + free(area); + free(flag); + + if (status != 0) + { + char errorMessage[512]; + sep_get_errmsg(status, errorMessage); + qCritical(KSTARS_FITS) << errorMessage; + return -1; + } + + return starCenters.count(); +} + +template +void FITSSEPDetector::getFloatBuffer(float * buffer, int x, int y, int w, int h, FITSData const *data) const +{ + auto * rawBuffer = reinterpret_cast(data->getImageBuffer()); + + if (buffer == nullptr) + return; + + float * floatPtr = buffer; + + int x2 = x + w; + int y2 = y + h; + + FITSData::Statistic const &stats = data->getStatistics(); + + for (int y1 = y; y1 < y2; y1++) + { + int offset = y1 * stats.width; + for (int x1 = x; x1 < x2; x1++) + { + *floatPtr++ = rawBuffer[offset + x1]; + } + } +} diff --git a/kstars/fitsviewer/fitsstardetector.h b/kstars/fitsviewer/fitsstardetector.h new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsstardetector.h @@ -0,0 +1,78 @@ +/*************************************************************************** + fitsstardetector.h - FITS Image + ------------------- + begin : Fri March 27 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#ifndef FITSSTARDETECTION_H +#define FITSSTARDETECTION_H + +#include +#include +#include + +#include "fitsdata.h" + +class FITSData; + +class Edge +{ +public: + float x {0}; + float y {0}; + int val {0}; + int scanned {0}; + float width {0}; + float HFR {-1}; + float sum {0}; +}; + +class FITSStarDetector : public QObject +{ + Q_OBJECT + +public: + /** @brief Instantiate a detector for a FITS data file. + */ + explicit FITSStarDetector(FITSData *parent): + QObject(reinterpret_cast(parent)) + {}; + +public: + /** @brief Find sources in the parent FITS data file. + * @param starCenters is the list of sources to append to. + * @param boundary is the rectangle in which to find sources, by default the full frame. + * @return The number of sources detected by the procedure. + */ + virtual int findSources(QList &starCenters, QRect const &boundary = QRect()) = 0; + + /** @brief Configure the detection method. + * @param setting is the name of a detection setting. + * @param value is the value of the detection setting identified by 'setting'. + * @return The detector as a chain to call the overriden findSources. + */ + virtual FITSStarDetector & configure(const QString &setting, const QVariant &value) = 0; + +public: + /** @brief Helper to configure the detection method from a data model. + * @param settings is the list of key/value pairs for the method to use settings from. + * @note Data model 'settings' is considered a key/value list, using column 1 text as case-insensitive keys and column 2 data as values. + * @return The detector as a chain to call the overriden findSources. + */ + FITSStarDetector & configure(QStandardItemModel const &settings); +}; + +#endif // FITSSTARDETECTION_H diff --git a/kstars/fitsviewer/fitsstardetector.cpp b/kstars/fitsviewer/fitsstardetector.cpp new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsstardetector.cpp @@ -0,0 +1,11 @@ +#include "fitsstardetector.h" + +FITSStarDetector& FITSStarDetector::configure(QStandardItemModel const &settings) +{ + Q_ASSERT(2 <= settings.columnCount()); + + for (int row = 0; row < settings.rowCount(); row++) + configure(settings.item(row, 0)->data().toString(), settings.item(row, 1)->data()); + + return *this; +} diff --git a/kstars/fitsviewer/fitsthresholddetector.h b/kstars/fitsviewer/fitsthresholddetector.h new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsthresholddetector.h @@ -0,0 +1,62 @@ +/*************************************************************************** + fitsthresholddetector.h - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#ifndef FITSTHRESHOLDDETECTOR_H +#define FITSTHRESHOLDDETECTOR_H + +#include "fitsstardetector.h" + +class FITSThresholdDetector: public FITSStarDetector +{ + Q_OBJECT + +public: + explicit FITSThresholdDetector(FITSData *parent): FITSStarDetector(parent) {}; + +public: + /** @brief Find sources in the parent FITS data file. + * @see FITSStarDetector::findSources(). + */ + int findSources(QList &starCenters, QRect const &boundary = QRect()) override; + + /** @brief Configure the detection method. + * @see FITSStarDetector::configure(). + * @note Parameter "threshold" defaults to THRESHOLD_PERCENTAGE of the mean pixel value of the frame. + * @todo Provide parameters for detection configuration. + */ + FITSStarDetector & configure(const QString &setting, const QVariant &value) override; + +public: + /** @group Detection internals + * @{ */ + static constexpr int THRESHOLD_PERCENTAGE { 120 }; + /** @} */ + +protected: + /** @internal Find sources in the parent FITS data file, dependent of the pixel depth. + * @see FITSGradientDetector::findSources. + */ + template + int findOneStar(QList &starCenters, const QRect &boundary) const; + +protected: + double focusThreshold { THRESHOLD_PERCENTAGE }; +}; + +#endif // FITSTHRESHOLDDETECTOR_H diff --git a/kstars/fitsviewer/fitsthresholddetector.cpp b/kstars/fitsviewer/fitsthresholddetector.cpp new file mode 100644 --- /dev/null +++ b/kstars/fitsviewer/fitsthresholddetector.cpp @@ -0,0 +1,227 @@ +/*************************************************************************** + fitsthresholddetector.cpp - FITS Image + ------------------- + begin : Sat March 28 2020 + copyright : (C) 2004 by Jasem Mutlaq, (C) 2020 by Eric Dejouhanet + email : eric.dejouhanet@gmail.com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * Some code fragments were adapted from Peter Kirchgessner's FITS plugin* + * See http://members.aol.com/pkirchg for more details. * + ***************************************************************************/ + +#include + +#include "fits_debug.h" +#include "fitsthresholddetector.h" + +FITSStarDetector& FITSThresholdDetector::configure(const QString &setting, const QVariant &value) +{ + bool ok = false; + + if (!setting.compare("threshold", Qt::CaseInsensitive)) + { + double const _focusThreshold = value.toDouble(&ok); + if (ok) + focusThreshold = _focusThreshold; + } + + return *this; +} + +int FITSThresholdDetector::findSources(QList &starCenters, QRect const &boundary) +{ + switch (parent()->property("dataType").toInt()) + { + case TBYTE: + return findOneStar(starCenters, boundary); + + case TSHORT: + return findOneStar(starCenters, boundary); + + case TUSHORT: + return findOneStar(starCenters, boundary); + + case TLONG: + return findOneStar(starCenters, boundary); + + case TULONG: + return findOneStar(starCenters, boundary); + + case TFLOAT: + return findOneStar(starCenters, boundary); + + case TLONGLONG: + return findOneStar(starCenters, boundary); + + case TDOUBLE: + return findOneStar(starCenters, boundary); + + default: + break; + } + + return 0; +} + +template +int FITSThresholdDetector::findOneStar(QList &starCenters, const QRect &boundary) const +{ + FITSData const * const image_data = reinterpret_cast(parent()); + + if (image_data == nullptr) + return 0; + + FITSData::Statistic const &stats = image_data->getStatistics(); + + if (boundary.isEmpty()) + return -1; + + 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(image_data->getImageBuffer()); + + // TODO replace magic number with something more useful to understand + double threshold = stats.mean[0] * focusThreshold / 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() << "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 0; + } + + // 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 subPixels; + subPixels.reserve(center->width / resolution); + + for (double x = leftEdge; x <= rightEdge; x += resolution) + { + //subPixels[x] = resolution * (image_buffer[static_cast(floor(x)) + cen_y * stats.width] - min); + double slice = resolution * (buffer[static_cast(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; + } + + starCenters.append(center); + + return 1; +} diff --git a/kstars/fitsviewer/fitsview.h b/kstars/fitsviewer/fitsview.h --- a/kstars/fitsviewer/fitsview.h +++ b/kstars/fitsviewer/fitsview.h @@ -34,9 +34,6 @@ #include -#define MINIMUM_PIXEL_RANGE 5 -#define MINIMUM_STDVAR 5 - class QAction; class QEvent; class QGestureEvent; diff --git a/kstars/fitsviewer/stretch.h b/kstars/fitsviewer/stretch.h --- a/kstars/fitsviewer/stretch.h +++ b/kstars/fitsviewer/stretch.h @@ -69,7 +69,7 @@ /** * @brief computeParams Automatically generates and sets stretch parameters from the image. */ - StretchParams computeParams(uint8_t *input); + StretchParams computeParams(const uint8_t *input); /** * @brief run run the stretch algorithm according to the params given @@ -81,11 +81,11 @@ * Sampling is applied to the output (that is, with sampling=2, we compute every other output * sample both in width and height, so the output would have about 4X fewer pixels. */ - void run(uint8_t *input, QImage *output_image, int sampling=1); + void run(uint8_t const *input, QImage *output_image, int sampling=1); private: // Adjusts input_range for float and double types. - void recalculateInputRange(uint8_t *input); + void recalculateInputRange(const uint8_t *input); // Inputs. int image_width; diff --git a/kstars/fitsviewer/stretch.cpp b/kstars/fitsviewer/stretch.cpp --- a/kstars/fitsviewer/stretch.cpp +++ b/kstars/fitsviewer/stretch.cpp @@ -26,7 +26,7 @@ // Returns the rough max of the buffer. template -T sampledMax(T *values, int size, int sampleBy) +T sampledMax(T const *values, int size, int sampleBy) { T maxVal = 0; for (int i = 0; i < size; i+= sampleBy) @@ -38,7 +38,7 @@ // Returns the median of the sample values. // The values are not modified. template -T median(T *values, int size, int sampleBy) +T median(T const *values, int size, int sampleBy) { const int downsampled_size = size / sampleBy; std::vector samples(downsampled_size); @@ -223,7 +223,7 @@ // See section 8.5.7 in above link https://pixinsight.com/doc/docs/XISF-1.0-spec/XISF-1.0-spec.html template -void computeParamsOneChannel(T *buffer, StretchParams1Channel *params, +void computeParamsOneChannel(T const *buffer, StretchParams1Channel *params, int inputRange, int height, int width) { // Find the median sample. @@ -316,40 +316,40 @@ input_range = getRange(dataType); } -void Stretch::run(uint8_t *input, QImage *outputImage, int sampling) +void Stretch::run(uint8_t const *input, QImage *outputImage, int sampling) { Q_ASSERT(outputImage->width() == (image_width + sampling - 1) / sampling); Q_ASSERT(outputImage->height() == (image_height + sampling - 1) / sampling); recalculateInputRange(input); switch (dataType) { case TBYTE: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TSHORT: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TUSHORT: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TLONG: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TFLOAT: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TLONGLONG: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; case TDOUBLE: - stretchChannels(reinterpret_cast(input), outputImage, params, + stretchChannels(reinterpret_cast(input), outputImage, params, input_range, image_height, image_width, image_channels, sampling); break; default: @@ -359,20 +359,20 @@ // The input range for float/double is ambiguous, and we can't tell without the buffer, // so we set it to 64K and possibly reduce it when we see the data. -void Stretch::recalculateInputRange(uint8_t *input) +void Stretch::recalculateInputRange(uint8_t const *input) { if (input_range <= 1) return; if (dataType != TFLOAT && dataType != TDOUBLE) return; float mx = 0; if (dataType == TFLOAT) - mx = sampledMax(reinterpret_cast(input), image_height * image_width, 1000); + mx = sampledMax(reinterpret_cast(input), image_height * image_width, 1000); else if (dataType == TDOUBLE) - mx = sampledMax(reinterpret_cast(input), image_height * image_width, 1000); + mx = sampledMax(reinterpret_cast(input), image_height * image_width, 1000); if (mx <= 1.01f) input_range = 1; } -StretchParams Stretch::computeParams(uint8_t *input) +StretchParams Stretch::computeParams(uint8_t const *input) { recalculateInputRange(input); StretchParams result; @@ -385,49 +385,49 @@ { case TBYTE: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TSHORT: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TUSHORT: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TLONG: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TFLOAT: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TLONGLONG: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break; } case TDOUBLE: { - auto buffer = reinterpret_cast(input); + auto buffer = reinterpret_cast(input); computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width); break;