probabilitycontours2d.cpp

Go to the documentation of this file.

/**
 * @file probabilitycontours2d.cpp
 * @brief Implements the probabilitycontours2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/probabilitycontours2d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "core/quantum/probabilitycontours2d.h"
#include "core/math/complexfield2d.h"
 
#include <algorithm>
#include <cmath>
 
namespace qpe {
namespace {
 
struct P2 { float x; float y; };
 
/**
 * @brief Returns probability.
 * @param field Field data to read or modify.
 * @param x X coordinate/value.
 * @param y Y coordinate/value.
 */
float probability(const ComplexField2D &field, std::size_t x, std::size_t y)
{
    return static_cast<float>(std::norm(field(x, y)));
}
 
/**
 * @brief Performs x coordinate.
 * @param field Field data to read or modify.
 * @param x X coordinate/value.
 */
float xCoordinate(const ComplexField2D &field, std::size_t x)
{
    const double centre = 0.5 * static_cast<double>(field.width() - 1u);
    return static_cast<float>((static_cast<double>(x) - centre) * static_cast<double>(field.dx()));
}
 
/**
 * @brief Performs y coordinate.
 * @param field Field data to read or modify.
 * @param y Y coordinate/value.
 */
float yCoordinate(const ComplexField2D &field, std::size_t y)
{
    const double centre = 0.5 * static_cast<double>(field.height() - 1u);
    return static_cast<float>((static_cast<double>(y) - centre) * static_cast<double>(field.dy()));
}
 
/**
 * @brief Performs interpolate.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param va The va value used by the operation.
 * @param vb The vb value used by the operation.
 * @param level The level value used by the operation.
 */
P2 interpolate(const P2 &a, const P2 &b, float va, float vb, float level)
{
    const float dv = vb - va;
    float t = 0.5f;
    if (std::fabs(dv) > 1.0e-30f) {
        t = (level - va) / dv;
    }
    t = std::max(0.0f, std::min(1.0f, t));
    P2 p;
    p.x = a.x + t * (b.x - a.x);
    p.y = a.y + t * (b.y - a.y);
    return p;
}
 
/**
 * @brief Adds segment.
 * @param out The out value used by the operation.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param level The level value used by the operation.
 * @param normalizedLevel The normalized level value used by the operation.
 */
void addSegment(ProbabilityContours2DResult &out,
                const P2 &a,
                const P2 &b,
                float level,
                float normalizedLevel)
{
    ProbabilityContourSegment2D s;
    s.x1 = a.x;
    s.y1 = a.y;
    s.x2 = b.x;
    s.y2 = b.y;
    s.level = level;
    s.normalizedLevel = normalizedLevel;
    out.segments.push_back(s);
}
 
std::size_t automaticStride(std::size_t width,
                            std::size_t height,
                            std::size_t targetSamplesPerAxis)
{
    const std::size_t target = std::max<std::size_t>(32u, targetSamplesPerAxis);
    const std::size_t largest = std::max(width, height);
    if (largest <= target) {
        return 1u;
    }
    return std::max<std::size_t>(1u, (largest + target - 1u) / target);
}
 
} // namespace
 
ProbabilityContours2DResult ProbabilityContours2D::compute(
    const ComplexField2D &field,
    std::size_t requestedLevels,
    ProbabilityContourSpacing2D spacing,
    float logDecades,
    std::size_t targetSamplesPerAxis)
{
    ProbabilityContours2DResult out;
    out.spacing = spacing;
    out.logDecades = std::max(0.25f, std::min(12.0f, logDecades));
 
    if (field.empty() || field.width() < 2u || field.height() < 2u || requestedLevels == 0u) {
        return out;
    }
 
    float maxProbability = 0.0f;
    for (std::size_t y = 0u; y < field.height(); ++y) {
        for (std::size_t x = 0u; x < field.width(); ++x) {
            maxProbability = std::max(maxProbability, probability(field, x, y));
        }
    }
    out.maximumProbability = maxProbability;
    if (!(maxProbability > 1.0e-30f) || !std::isfinite(static_cast<double>(maxProbability))) {
        return out;
    }
 
    const std::size_t count = std::max<std::size_t>(1u, requestedLevels);
    out.levels.reserve(count);
    for (std::size_t i = 0u; i < count; ++i) {
        const float f = static_cast<float>(i + 1u) / static_cast<float>(count + 1u);
        float level = 0.0f;
        if (spacing == ProbabilityContourSpacing2D::Logarithmic) {
            // f runs from faint to bright.  End points stay strictly inside
            // the data range so each level has a chance to produce a contour.
            const float exponent = -out.logDecades * (1.0f - f);
            level = maxProbability * std::pow(10.0f, exponent);
        } else {
            level = maxProbability * f;
        }
        out.levels.push_back(level);
    }
 
    out.samplingStride = automaticStride(field.width(), field.height(), targetSamplesPerAxis);
    const std::size_t stride = out.samplingStride;
 
    for (std::size_t li = 0u; li < out.levels.size(); ++li) {
        const float level = out.levels[li];
        // Colour rank is based on contour order rather than raw probability.
        // With logarithmic levels the latter would compress nearly every faint
        // contour into the same blue colour.
        const float normalizedLevel =
            static_cast<float>(li + 1u) / static_cast<float>(out.levels.size() + 1u);
 
        for (std::size_t y0 = 0u; y0 + 1u < field.height(); y0 += stride) {
            const std::size_t y1 = std::min(field.height() - 1u, y0 + stride);
            if (y1 == y0) {
                continue;
            }
            for (std::size_t x0 = 0u; x0 + 1u < field.width(); x0 += stride) {
                const std::size_t x1 = std::min(field.width() - 1u, x0 + stride);
                if (x1 == x0) {
                    continue;
                }
 
                const float v0 = probability(field, x0, y0); // bottom-left
                const float v1 = probability(field, x1, y0); // bottom-right
                const float v2 = probability(field, x1, y1); // top-right
                const float v3 = probability(field, x0, y1); // top-left
 
                int mask = 0;
                if (v0 > level) mask |= 1;
                if (v1 > level) mask |= 2;
                if (v2 > level) mask |= 4;
                if (v3 > level) mask |= 8;
                if (mask == 0 || mask == 15) {
                    continue;
                }
 
                const P2 p0 = { xCoordinate(field, x0), yCoordinate(field, y0) };
                const P2 p1 = { xCoordinate(field, x1), yCoordinate(field, y0) };
                const P2 p2 = { xCoordinate(field, x1), yCoordinate(field, y1) };
                const P2 p3 = { xCoordinate(field, x0), yCoordinate(field, y1) };
 
                P2 edge[4];
                edge[0] = interpolate(p0, p1, v0, v1, level);
                edge[1] = interpolate(p1, p2, v1, v2, level);
                edge[2] = interpolate(p2, p3, v2, v3, level);
                edge[3] = interpolate(p3, p0, v3, v0, level);
 
                switch (mask) {
                case 1: case 14: addSegment(out, edge[3], edge[0], level, normalizedLevel); break;
                case 2: case 13: addSegment(out, edge[0], edge[1], level, normalizedLevel); break;
                case 3: case 12: addSegment(out, edge[3], edge[1], level, normalizedLevel); break;
                case 4: case 11: addSegment(out, edge[1], edge[2], level, normalizedLevel); break;
                case 6: case 9:  addSegment(out, edge[0], edge[2], level, normalizedLevel); break;
                case 7: case 8:  addSegment(out, edge[3], edge[2], level, normalizedLevel); break;
                case 5: {
                    const float centre = 0.25f * (v0 + v1 + v2 + v3);
                    if (centre > level) {
                        addSegment(out, edge[3], edge[2], level, normalizedLevel);
                        addSegment(out, edge[0], edge[1], level, normalizedLevel);
                    } else {
                        addSegment(out, edge[3], edge[0], level, normalizedLevel);
                        addSegment(out, edge[2], edge[1], level, normalizedLevel);
                    }
                    break;
                }
                case 10: {
                    const float centre = 0.25f * (v0 + v1 + v2 + v3);
                    if (centre > level) {
                        addSegment(out, edge[3], edge[0], level, normalizedLevel);
                        addSegment(out, edge[2], edge[1], level, normalizedLevel);
                    } else {
                        addSegment(out, edge[3], edge[2], level, normalizedLevel);
                        addSegment(out, edge[0], edge[1], level, normalizedLevel);
                    }
                    break;
                }
                default:
                    break;
                }
            }
        }
    }
 
    return out;
}
 
} // namespace qpe