potentialcontours2d.cpp

Go to the documentation of this file.

/**
 * @file potentialcontours2d.cpp
 * @brief Implements the potentialcontours2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/potentialcontours2d.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/potentialcontours2d.h"
#include "core/quantum/potentialfield2d.h"
 
#include <algorithm>
#include <cmath>
 
namespace qpe {
namespace {
 
struct P2 { float x; float y; };
 
/**
 * @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-20f) {
        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.
 */
void addSegment(PotentialContours2DResult &out, const P2 &a, const P2 &b, float level)
{
    PotentialContourSegment2D s;
    s.x1 = a.x; s.y1 = a.y;
    s.x2 = b.x; s.y2 = b.y;
    s.level = level;
    out.segments.push_back(s);
}
 
} // namespace
 
PotentialContours2DResult PotentialContours2D::compute(const PotentialField2D &potential,
                                                        std::size_t requestedLevels)
{
    PotentialContours2DResult out;
    if (potential.empty() || potential.width() < 2u || potential.height() < 2u || requestedLevels == 0u) {
        return out;
    }
 
    out.minimum = potential.minimum();
    out.maximum = potential.maximum();
    const float range = out.maximum - out.minimum;
    if (!(range > 1.0e-12f)) {
        return out;
    }
 
    bool onlyExtremes = true;
    const float eps = std::max(1.0e-7f, std::fabs(range) * 1.0e-6f);
    for (std::size_t y = 0; y < potential.height() && onlyExtremes; ++y) {
        for (std::size_t x = 0; x < potential.width(); ++x) {
            const float v = potential(x, y);
            if (std::fabs(v - out.minimum) > eps && std::fabs(v - out.maximum) > eps) {
                onlyExtremes = false;
                break;
            }
        }
    }
    out.binaryField = onlyExtremes;
 
    if (onlyExtremes) {
        out.levels.push_back(out.minimum + 0.5f * range);
    } else {
        const std::size_t count = std::max<std::size_t>(1u, requestedLevels);
        for (std::size_t i = 0; i < count; ++i) {
            const float f = static_cast<float>(i + 1u) / static_cast<float>(count + 1u);
            out.levels.push_back(out.minimum + f * range);
        }
    }
 
    for (std::size_t li = 0; li < out.levels.size(); ++li) {
        const float level = out.levels[li];
        for (std::size_t y = 0; y + 1u < potential.height(); ++y) {
            for (std::size_t x = 0; x + 1u < potential.width(); ++x) {
                const float v0 = potential(x, y);         // bottom-left
                const float v1 = potential(x + 1u, y);    // bottom-right
                const float v2 = potential(x + 1u, y + 1u); // top-right
                const float v3 = potential(x, y + 1u);    // 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 = { static_cast<float>(potential.xCoordinate(x)),
                                static_cast<float>(potential.yCoordinate(y)) };
                const P2 p1 = { static_cast<float>(potential.xCoordinate(x + 1u)),
                                static_cast<float>(potential.yCoordinate(y)) };
                const P2 p2 = { static_cast<float>(potential.xCoordinate(x + 1u)),
                                static_cast<float>(potential.yCoordinate(y + 1u)) };
                const P2 p3 = { static_cast<float>(potential.xCoordinate(x)),
                                static_cast<float>(potential.yCoordinate(y + 1u)) };
 
                P2 edge[4];
                edge[0] = interpolate(p0, p1, v0, v1, level); // bottom
                edge[1] = interpolate(p1, p2, v1, v2, level); // right
                edge[2] = interpolate(p2, p3, v2, v3, level); // top
                edge[3] = interpolate(p3, p0, v3, v0, level); // left
 
                switch (mask) {
                case 1: case 14: addSegment(out, edge[3], edge[0], level); break;
                case 2: case 13: addSegment(out, edge[0], edge[1], level); break;
                case 3: case 12: addSegment(out, edge[3], edge[1], level); break;
                case 4: case 11: addSegment(out, edge[1], edge[2], level); break;
                case 6: case 9:  addSegment(out, edge[0], edge[2], level); break;
                case 7: case 8:  addSegment(out, edge[3], edge[2], level); break;
                case 5: {
                    const float center = 0.25f * (v0 + v1 + v2 + v3);
                    if (center > level) {
                        addSegment(out, edge[3], edge[2], level);
                        addSegment(out, edge[0], edge[1], level);
                    } else {
                        addSegment(out, edge[3], edge[0], level);
                        addSegment(out, edge[2], edge[1], level);
                    }
                    break;
                }
                case 10: {
                    const float center = 0.25f * (v0 + v1 + v2 + v3);
                    if (center > level) {
                        addSegment(out, edge[3], edge[0], level);
                        addSegment(out, edge[2], edge[1], level);
                    } else {
                        addSegment(out, edge[3], edge[2], level);
                        addSegment(out, edge[0], edge[1], level);
                    }
                    break;
                }
                default: break;
                }
            }
        }
    }
 
    return out;
}
 
} // namespace qpe