Jungfraujoch/viewer/charts/JFJochDatasetInfoChartView.cpp

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <QMenu>
#include <QApplication>
#include <QClipboard>
#include <QCategoryAxis>

#include "JFJochDatasetInfoChartView.h"

JFJochDatasetInfoChartView::JFJochDatasetInfoChartView(QWidget *parent)
    : QChartView(new QChart(), parent) {
    chart()->legend()->hide();
    setRenderHint(QPainter::Antialiasing);
   // setRubberBand(QChartView::RubberBand::RectangleRubberBand);
    setMouseTracking(true);

    m_hoverLoadTimer = new QTimer(this);
    m_hoverLoadTimer->setSingleShot(true);
    connect(m_hoverLoadTimer, &QTimer::timeout, this, &JFJochDatasetInfoChartView::onHoverLoadTimeout);
}

void JFJochDatasetInfoChartView::setImage(int64_t val) {
    if (!currentSeries)
        return;

    curr_image = val;

    if (val < values.size() && val >= 0) {
        currentSeries->clear();
        if (std::isfinite(values[curr_image])) {
            const double disp = values[curr_image];
            if (std::isfinite(disp))
                currentSeries->append(curr_image, disp);
        }
    }
}

void JFJochDatasetInfoChartView::mousePressEvent(QMouseEvent *event) {
    if (event->button() == Qt::LeftButton) {
        QPointF clickedPoint = event->pos();
        QPointF chartCoord = chart()->mapToValue(clickedPoint);
        int64_t val = std::lround(chartCoord.x());
        if (val >= 0 && val < values.size())
            emit imageSelected(val);
    }
    QChartView::mousePressEvent(event); // Call the base implementation
}

void JFJochDatasetInfoChartView::resetZoom() {
    chart()->zoomReset();
}

void JFJochDatasetInfoChartView::updateChart() {
    chart()->removeAllSeries();
    if (m_hoverLine) {
        chart()->scene()->removeItem(m_hoverLine);
        delete m_hoverLine;
        m_hoverLine = nullptr;
    }

#ifdef JFJOCH_USE_FFTW
    if (m_showFFT) {
        buildFFTChart();
        return;
    }
#endif

    buildTimeDomainChart();
}

void JFJochDatasetInfoChartView::buildTimeDomainChart() {
    if (values.size() >= binning) {
        // At least one full point

        series = new QLineSeries(this);
        currentSeries = new QScatterSeries(this);

        double dispMin = std::numeric_limits<double>::infinity();
        double dispMax = -std::numeric_limits<double>::infinity();

        if (binning == 1) {
            for (int i = 0; i < values.size(); i++) {
                if (!std::isfinite(values[i])) continue;
                const double disp = values[i];
                if (!std::isfinite(disp)) continue;
                series->append(i, disp);
                if (disp < dispMin) dispMin = disp;
                if (disp > dispMax) dispMax = disp;
            }
        } else {
            for (int i = 0; i < static_cast<int>(values.size() / binning); i++) {
                double tmp = 0.0;
                int64_t count = 0;
                for (int b = 0; b < binning; b++) {
                    const double v = values[i * binning + b];
                    if (std::isfinite(v)) {
                        tmp += v;
                        count++;
                    }
                }
                if (count > 0) {
                    const double mean = tmp / static_cast<double>(count);
                    const double disp = mean;
                    if (std::isfinite(disp)) {
                        series->append((i + 0.5) * binning, disp);
                        if (disp < dispMin) dispMin = disp;
                        if (disp > dispMax) dispMax = disp;
                    }
                }
            }
        }

        if (curr_image < values.size() && curr_image >= 0 && std::isfinite(values[curr_image])) {
            currentSeries->append(curr_image, values[curr_image]);
        }

        chart()->addSeries(series);
        chart()->addSeries(currentSeries);
        chart()->createDefaultAxes();

        // ----- X axis handling -----
        QValueAxis *axisX = qobject_cast<QValueAxis *>(chart()->axisX(series));
        if (axisX) {
            if (goniometer_axis.has_value() && m_xUseGoniometerAxis) {
                // Hide labels on numeric axis and move it to the top

                axisX->setTitleText(QString(""));
                axisX->setLabelsVisible(false);

                // Re-attach numeric axis on top side (default axis on other side)
                chart()->removeAxis(axisX);
                chart()->addAxis(axisX, Qt::AlignTop);
                series->attachAxis(axisX);
                currentSeries->attachAxis(axisX);

                // Build a visible category axis on the bottom with goniometer angles
                auto *axXcat = new QCategoryAxis();
                axXcat->setLabelsPosition(QCategoryAxis::AxisLabelsPositionOnValue);
                axXcat->setGridLineVisible(false);
                axXcat->setMinorGridLineVisible(false);
                axXcat->setTitleText(QStringLiteral("Rotation angle (deg.)"));

                const int tickCountX = std::max(2, axisX->tickCount());
                const double xmin = axisX->min();
                const double xmax = axisX->max();
                const double xstep = (tickCountX > 1) ? (xmax - xmin) / (tickCountX - 1) : 0.0;
                const int64_t lastIdx = static_cast<int64_t>(values.empty() ? 0 : values.size() - 1);

                for (int i = 0; i < tickCountX; ++i) {
                    const double xv = (i == tickCountX - 1) ? xmax : (xmin + i * xstep);
                    // Map tick position to closest image index
                    int64_t imgIdx = static_cast<int64_t>(std::llround(xv));
                    if (imgIdx < 0) imgIdx = 0;
                    if (imgIdx > lastIdx) imgIdx = lastIdx;

                    double angleDeg = 0.0;
                    if (lastIdx >= 0) {
                        angleDeg = goniometer_axis->GetAngle_deg(imgIdx);
                    }

                    QString lab = QString::number(angleDeg, 'f', 2);
                    axXcat->append(lab, xv);
                }

                chart()->addAxis(axXcat, Qt::AlignBottom);
                series->attachAxis(axXcat);
                currentSeries->attachAxis(axXcat);
            } else {
                axisX->setLabelsVisible(true);
                axisX->setTitleText(QStringLiteral("Image number"));
            }
        }

        // ----- Y-axis handling -----
        QValueAxis *axisY = qobject_cast<QValueAxis *>(chart()->axisY(series));
        if (axisY) {
            if (std::isfinite(dispMin) && std::isfinite(dispMax)) {
                if (m_minYZeroEnabled) {
                    const double minY = 0.0;
                    const double maxY = (dispMax > minY) ? dispMax : (minY + 1.0);
                    axisY->setRange(minY, maxY);
                } else {
                    // Default: tight range to data
                    if (!(dispMax > dispMin)) {
                        // Avoid zero-height range
                        dispMax = dispMin + 1.0;
                    }
                    axisY->setRange(dispMin, dispMax);
                }
            }

            if (m_yOneOverD) {
                // Keep value axis for numeric range + grid, but move it to the RIGHT
                axisY->setLabelsVisible(false);
                chart()->removeAxis(axisY);
                chart()->addAxis(axisY, Qt::AlignRight);
                series->attachAxis(axisY);
                currentSeries->attachAxis(axisY);

                // Build a mirrored visible axis with labels in d (Å) on the LEFT
                auto *axYcat = new QCategoryAxis();
                axYcat->setLabelsPosition(QCategoryAxis::AxisLabelsPositionOnValue);
                axYcat->setGridLineVisible(false);
                axYcat->setMinorGridLineVisible(false);
                axYcat->setTitleText(QStringLiteral("d (Å)"));

                const int tickCountY = std::max(2, axisY->tickCount());
                const double ymin = axisY->min();
                const double ymax = axisY->max();
                const double ystep = (tickCountY > 1) ? (ymax - ymin) / (tickCountY - 1) : 0.0;

                for (int i = 0; i < tickCountY; ++i) {
                    const double yv = (i == tickCountY - 1) ? ymax : (ymin + i * ystep);
                    QString lab;
                    if (!(yv > 0.0)) {
                        lab = QStringLiteral("—"); // invalid for d
                    } else if (std::abs(yv) < 1e-300) {
                        lab = QStringLiteral("∞");
                    } else {
                        const double d = 1.0 / std::sqrt(yv);
                        lab = QString("%1 Å").arg(d, 0, 'f', 2);
                    }
                    axYcat->append(lab, yv);
                }

                chart()->addAxis(axYcat, Qt::AlignLeft);
                series->attachAxis(axYcat);
                currentSeries->attachAxis(axYcat);

                // Give a bit more room on the left so labels are not clipped
                QMargins m = chart()->margins();
                if (m.left() < 12) {
                    m.setLeft(12);
                    chart()->setMargins(m);
                }
            } else {
                // Normal numeric labels, axis on the LEFT
                chart()->removeAxis(axisY);
                chart()->addAxis(axisY, Qt::AlignLeft);
                series->attachAxis(axisY);
                currentSeries->attachAxis(axisY);

                axisY->setTitleText(QString());
                axisY->setLabelsVisible(true);
            }
        }
    }
}


void JFJochDatasetInfoChartView::setBinning(int64_t val) {
    if (val >= 1) {
        binning = val;
        updateChart();
    }
}

void JFJochDatasetInfoChartView::contextMenuEvent(QContextMenuEvent *event) {
    QMenu menu(this);
    QAction *copyXY = menu.addAction("Copy (x y) points");
    copyXY->setEnabled(!values.empty());

    QAction *sep1 = menu.addSeparator();
    Q_UNUSED(sep1);

    QAction *actMinYZero = menu.addAction("Y min at 0");
    actMinYZero->setCheckable(true);
    actMinYZero->setChecked(m_minYZeroEnabled);

    QAction *actXGoniometer = menu.addAction("Use goniometer X-axis");
    actXGoniometer->setCheckable(true);
    actXGoniometer->setChecked(m_xUseGoniometerAxis);
    actXGoniometer->setEnabled(goniometer_axis.has_value());

#ifdef JFJOCH_USE_FFTW
    QAction *actShowFFT = menu.addAction("Show FFT (amplitude vs Hz)");
    actShowFFT->setCheckable(true);
    actShowFFT->setChecked(m_showFFT);
    // Require valid sampling interval
    actShowFFT->setEnabled(!values.empty() && image_time_us > 0.0);
#endif

    QAction *chosen = menu.exec(event->globalPos());
    if (chosen == copyXY) {
        QString out;
        out.reserve(static_cast<int>(values.size() * 16)); // rough prealloc
        for (size_t i = 0; i < values.size(); ++i) {
            out.append(QString::number(i));
            out.append(' ');
            out.append(QString::number(values[i], 'g', 10));
            if (i + 1 < values.size()) out.append('\n');
        }
        QClipboard *cb = QApplication::clipboard();
        cb->setText(out);
    } else if (chosen == actMinYZero) {
        m_minYZeroEnabled = !m_minYZeroEnabled;
        updateChart();
    } else if (chosen == actXGoniometer) {
        m_xUseGoniometerAxis = !m_xUseGoniometerAxis;
        updateChart();
#ifdef JFJOCH_USE_FFTW
    } else if (chosen == actShowFFT) {
        m_showFFT = !m_showFFT;
        updateChart();
#endif
    }

}

void JFJochDatasetInfoChartView::mouseMoveEvent(QMouseEvent *event) {
    QChartView::mouseMoveEvent(event);
    if (!series || values.empty())
        return;

#ifdef JFJOCH_USE_FFTW
    if (m_showFFT && !m_fftFrequenciesHz.empty()) {
        // FFT mode: x is frequency in Hz
        const QPointF chartPos = chart()->mapToValue(event->pos(), series);
        double f = chartPos.x();
        if (!std::isfinite(f))
            return;

        // If we only have DC, nothing meaningful to show
        if (m_fftFrequenciesHz.size() <= 1)
            return;

        // Find nearest FFT bin, excluding k = 0 (DC component)
        int64_t bestIdx = -1;
        double bestDiff = std::numeric_limits<double>::infinity();
        for (size_t i = 1; i < m_fftFrequenciesHz.size(); ++i) {
            const double diff = std::abs(m_fftFrequenciesHz[i] - f);
            if (diff < bestDiff) {
                bestDiff = diff;
                bestIdx = static_cast<int64_t>(i);
            }
        }
        if (bestIdx < 1)
            return;

        const double fBin = m_fftFrequenciesHz[static_cast<size_t>(bestIdx)];
        const double amp = m_fftMagnitudes[static_cast<size_t>(bestIdx)];

        // Map x position of that bin to scene coords for vertical line
        const QRectF plotArea = chart()->plotArea();
        const QPointF ptOnChart = chart()->mapToPosition(QPointF(fBin, 0.0), series);

        if (!m_hoverLine) {
            m_hoverLine = new QGraphicsLineItem;
            m_hoverLine->setPen(QPen(QColor(200, 0, 0, 150), 1.0));
            chart()->scene()->addItem(m_hoverLine);
        }

        m_hoverLine->setLine(QLineF(ptOnChart.x(), plotArea.top(),
                                    ptOnChart.x(), plotArea.bottom()));

        QString text = QString("f = %1 Hz, amplitude = %2")
                .arg(fBin, 0, 'g', 6)
                .arg(amp, 0, 'g', 6);
        emit writeStatusBar(text, 6000);

        // No image loading in FFT mode
        m_hoverLoadTimer->stop();
        m_hoverPendingIdx = -1;
        return;
    }
#endif

    if (values.empty())
        return;

    // Map mouse position to chart coordinates
    const QPointF chartPos = chart()->mapToValue(event->pos(), series);
    double xVal = chartPos.x();

    // Convert x to closest image index
    int64_t idx = std::lround(xVal);
    if (idx < 0 || idx >= static_cast<int64_t>(values.size()))
        return;

    // Map that x position to scene coords for the vertical line
    const QRectF plotArea = chart()->plotArea();
    const QPointF ptOnChart = chart()->mapToPosition(QPointF(static_cast<double>(idx), 0.0), series);

    if (!m_hoverLine) {
        m_hoverLine = new QGraphicsLineItem;
        m_hoverLine->setPen(QPen(QColor(200, 0, 0, 150), 1.0));
        chart()->scene()->addItem(m_hoverLine);
    }

    m_hoverLine->setLine(QLineF(ptOnChart.x(), plotArea.top(),
                                ptOnChart.x(), plotArea.bottom()));

    // Status bar text
    const double yv = values[static_cast<size_t>(idx)];
    QString text;

    if (m_yOneOverD) {
        if (std::isfinite(yv) && yv > 0.0) {
            const double d = 1.0 / std::sqrt(yv);
            text = QString("image = %1 d = %2 Å")
                    .arg(idx)
                    .arg(d, 0, 'f', 2);
        } else {
            text = QString("image = %1, no resolution estimate").arg(idx);
        }
    } else {
        text = QString("image = %1 value = %2")
                .arg(idx)
                .arg(yv, 0, 'g', 6);
    }
    emit writeStatusBar(text, 6000);

    // Debounced image load on hover when Shift is pressed
    if (event->modifiers() & Qt::ShiftModifier) {
        if (!m_hoverLoadTimer->isActive()) {
            m_hoverPendingIdx = -1;
            if (idx != curr_image)
                emit imageSelected(idx);
            m_hoverLoadTimer->start(500); // debounce
        } else
            m_hoverPendingIdx = idx;
    } else {
        m_hoverLoadTimer->stop();
        m_hoverPendingIdx = -1;
    }
}

void JFJochDatasetInfoChartView::leaveEvent(QEvent *event) {
    QChartView::leaveEvent(event);
    if (m_hoverLine) {
        chart()->scene()->removeItem(m_hoverLine);
        delete m_hoverLine;
        m_hoverLine = nullptr;
    }
    m_hoverLoadTimer->stop();
    m_hoverPendingIdx = -1;
    emit writeStatusBar(QString(), 0);
}

void JFJochDatasetInfoChartView::onHoverLoadTimeout() {
    if (!(QApplication::keyboardModifiers() & Qt::ShiftModifier))
        return;

    if (m_hoverPendingIdx >= 0 &&
        m_hoverPendingIdx < static_cast<int64_t>(values.size())) {
        if (m_hoverPendingIdx != curr_image) {
            emit imageSelected(m_hoverPendingIdx);
        }
    }
}

#ifdef JFJOCH_USE_FFTW
void JFJochDatasetInfoChartView::buildFFTChart() {
    const size_t N = values.size();

    if (N == 0 || !image_time_us.has_value() || image_time_us <= 0.0) {
        return;
    }

    // Prepare input buffer (single precision, NaN/inf treated as 0)
    std::vector<float> in(N, 0.0f);
    for (size_t i = 0; i < N; ++i) {
        const double v = values[i];
        in[i] = std::isfinite(v) ? static_cast<float>(v) : 0.0f;
    }

    const int n = static_cast<int>(N);
    const int nComplex = n / 2 + 1;

    std::vector<fftwf_complex> out(static_cast<size_t>(nComplex));

    fftwf_plan plan = fftwf_plan_dft_r2c_1d(
            n,
            in.data(),
            out.data(),
            FFTW_ESTIMATE);

    if (!plan) {
        return;
    }

    fftwf_execute(plan);
    fftwf_destroy_plan(plan);

    // Compute amplitude spectrum and frequencies (0 .. Nyquist)
    m_fftMagnitudes.resize(static_cast<size_t>(nComplex));
    m_fftFrequenciesHz.resize(static_cast<size_t>(nComplex));

    const double dt = image_time_us.value() * 1e-6;             // seconds per sample
    const double fs = 1.0 / dt;                         // sampling frequency
    const double df = fs / static_cast<double>(n);      // frequency resolution

    for (int k = 0; k < nComplex; ++k) {
        const double re = out[static_cast<size_t>(k)][0];
        const double im = out[static_cast<size_t>(k)][1];
        const double mag = std::hypot(re, im);          // amplitude

        m_fftMagnitudes[static_cast<size_t>(k)] = mag;
        m_fftFrequenciesHz[static_cast<size_t>(k)] = static_cast<double>(k) * df;
    }

    // Build chart series: X = frequency (Hz), Y = amplitude
    series = new QLineSeries(this);
    currentSeries = nullptr; // no "current image" marker in FFT mode

    double magMin = std::numeric_limits<double>::infinity();
    double magMax = -std::numeric_limits<double>::infinity();

    for (int k = 1; k < nComplex; ++k) {
        const double f = m_fftFrequenciesHz[static_cast<size_t>(k)];
        const double mag = m_fftMagnitudes[static_cast<size_t>(k)];
        series->append(f, mag);
        if (mag < magMin) magMin = mag;
        if (mag > magMax) magMax = mag;
    }

    chart()->addSeries(series);
    chart()->createDefaultAxes();

    QValueAxis *axisX = qobject_cast<QValueAxis *>(chart()->axisX(series));
    QValueAxis *axisY = qobject_cast<QValueAxis *>(chart()->axisY(series));

    if (axisX) {
        axisX->setTitleText(QStringLiteral("Frequency (Hz)"));
        axisX->setLabelsVisible(true);
    }

    if (axisY) {
        if (std::isfinite(magMin) && std::isfinite(magMax)) {
            if (!(magMax > magMin)) {
                magMax = magMin + 1.0;
            }
            axisY->setRange(magMin, magMax);
        }
        axisY->setTitleText(QStringLiteral("Amplitude"));
        axisY->setLabelsVisible(true);
    }
}
#endif