quadquality.hpp

Go to the documentation of this file.

 
#include "meshkit/detail/planmodel.hpp"
#ifndef MESHKIT_NATIVE_MESHQUADQUALITY_HPP
#define MESHKIT_NATIVE_MESHQUADQUALITY_HPP
 
#include <algorithm>
#include <cmath>
#include <limits>
 
namespace meshkit::plan
{
 
struct QuadQualityBounds {
    double minAngleDeg = 60.0;
    double maxAngleDeg = 120.0;
    double minScaledJacobian = 0.866;
    double maxAspect = 2.0; 
};
 
struct QuadQuality {
    double minAngleDeg = 0.0;
    double maxAngleDeg = 0.0;
    double scaledJacobian = 0.0; 
    double rectangularity = 0.0; 
    double aspect = 1.0;         
    double skew =
        0.0; 
    double area = 0.0;   
    bool convex = false; 
};
 
inline QuadQuality quadQuality(const Point2 &a, const Point2 &b, const Point2 &c, const Point2 &d)
{
    QuadQuality q;
    const Point2 p[4] = {a, b, c, d};
    Point2 local[4];
    double xScale = 0.0, yScale = 0.0;
    for (int i = 0; i < 4; ++i) {
        if (!std::isfinite(p[i].x()) || !std::isfinite(p[i].y()))
            return {};
        local[i] = p[i] - a;
        if (!std::isfinite(local[i].x()) || !std::isfinite(local[i].y()))
            return {};
        xScale = std::max(xScale, std::abs(local[i].x()));
        yScale = std::max(yScale, std::abs(local[i].y()));
    }
    if (!(xScale > 0.0) || !(yScale > 0.0))
        return {};
    // Local coordinates remove origin cancellation. Independent axis scales
    // keep products finite even for representable highly anisotropic quads.
    const auto areaCross = [&](int i, int j) {
        return (local[i].x() / xScale) * (local[j].y() / yScale) -
               (local[i].y() / yScale) * (local[j].x() / xScale);
    };
    const double area2 = areaCross(1, 2) + areaCross(2, 3);
    if (area2 == 0.0 || !std::isfinite(area2))
        return {};
    const double orient = area2 > 0.0 ? 1.0 : -1.0;
    // Rescale through exponents: neither xScale*yScale nor twice the final
    // area needs to be representable as an intermediate value.
    int areaExponent = 0, xExponent = 0, yExponent = 0;
    const double areaMantissa = std::frexp(std::abs(area2), &areaExponent);
    const double xMantissa = std::frexp(xScale, &xExponent);
    const double yMantissa = std::frexp(yScale, &yExponent);
    q.area =
        std::scalbn(areaMantissa * xMantissa * yMantissa, areaExponent + xExponent + yExponent - 1);
    if (!(q.area > 0.0) || !std::isfinite(q.area))
        return {};
 
    q.minAngleDeg = 360.0;
    q.maxAngleDeg = 0.0;
    q.scaledJacobian = 1.0;
    q.convex = true;
    double maxDev = 0.0;
    double side[4];
    for (int i = 0; i < 4; ++i) {
        const Point2 e1 = p[(i + 3) % 4] - p[i], e2 = p[(i + 1) % 4] - p[i];
        const double l1 = std::hypot(e1.x(), e1.y()), l2 = std::hypot(e2.x(), e2.y());
        if (!(l1 > 0.0) || !(l2 > 0.0) || !std::isfinite(l1) || !std::isfinite(l2))
            return {};
        side[i] = l2;
        const Point2 u1 = e1 / l1, u2 = e2 / l2;
        // Normalize first; raw cross/dot products and length products may
        // overflow or underflow even when the angular metrics are well-defined.
        const double cross = std::clamp(u2.x() * u1.y() - u2.y() * u1.x(), -1.0, 1.0);
        const double dot = std::clamp(u1.x() * u2.x() + u1.y() * u2.y(), -1.0, 1.0);
        const double ang = std::atan2(std::abs(cross), dot) * 180.0 / pi;
        const double sj = cross * orient;
        if (sj <= 0.0)
            q.convex = false;
        q.scaledJacobian = std::min(q.scaledJacobian, sj);
        q.minAngleDeg = std::min(q.minAngleDeg, ang);
        q.maxAngleDeg = std::max(q.maxAngleDeg, ang);
        maxDev = std::max(maxDev, std::abs(ang - 90.0));
    }
    q.rectangularity = std::clamp(1.0 - maxDev / 90.0, 0.0, 1.0);
    const double s02 = 0.5 * side[0] + 0.5 * side[2];
    const double s13 = 0.5 * side[1] + 0.5 * side[3];
    const double lo = std::min(s02, s13), hi = std::max(s02, s13);
    if (!(lo > 0.0))
        return {};
    q.aspect = hi / lo;
    if (!std::isfinite(q.aspect))
        return {};
    const Point2 d1 = c - a, d2 = d - b;
    const double ld1 = std::hypot(d1.x(), d1.y()), ld2 = std::hypot(d2.x(), d2.y());
    if (!(ld1 > 0.0) || !(ld2 > 0.0) || !std::isfinite(ld1) || !std::isfinite(ld2))
        return {};
    const Point2 ud1 = d1 / ld1, ud2 = d2 / ld2;
    q.skew = std::clamp(std::abs(ud1.x() * ud2.x() + ud1.y() * ud2.y()), 0.0, 1.0);
    return q;
}
 
inline QuadQuality quadQuality(const Vector<MeshVertex> &vertices, const MeshTriangle &t)
{
    return quadQuality(vertices[t.v0].xy, vertices[t.v1].xy, vertices[t.v2].xy, vertices[t.v3].xy);
}
 
namespace quadquality_detail
{
inline bool validMetrics(const QuadQuality &q)
{
    // These checks also protect callers constructing QuadQuality manually.
    return q.convex && std::isfinite(q.area) && q.area > 0.0 && std::isfinite(q.minAngleDeg) &&
           std::isfinite(q.maxAngleDeg) && q.minAngleDeg > 0.0 && q.maxAngleDeg < 180.0 &&
           q.minAngleDeg <= q.maxAngleDeg && std::isfinite(q.scaledJacobian) &&
           q.scaledJacobian > 0.0 && q.scaledJacobian <= 1.0 && std::isfinite(q.rectangularity) &&
           q.rectangularity >= 0.0 && q.rectangularity <= 1.0 && std::isfinite(q.aspect) &&
           q.aspect >= 1.0 && std::isfinite(q.skew) && q.skew >= 0.0 && q.skew <= 1.0;
}
inline bool finiteBounds(const QuadQualityBounds &b)
{
    return std::isfinite(b.minAngleDeg) && std::isfinite(b.maxAngleDeg) &&
           std::isfinite(b.minScaledJacobian) && std::isfinite(b.maxAspect);
}
} // namespace quadquality_detail
 
inline bool quadAcceptable(const QuadQuality &q, const QuadQualityBounds &b)
{
    if (!quadquality_detail::validMetrics(q) || !quadquality_detail::finiteBounds(b))
        return false;
    if (q.minAngleDeg < b.minAngleDeg || q.maxAngleDeg > b.maxAngleDeg)
        return false;
    if (q.scaledJacobian < b.minScaledJacobian)
        return false;
    if (b.maxAspect > 0.0 && q.aspect > b.maxAspect)
        return false;
    return true;
}
 
inline double quadScore(const QuadQuality &q, const QuadQualityBounds &b)
{
    if (!quadquality_detail::validMetrics(q) || !quadquality_detail::finiteBounds(b))
        return 0.0;
    const double aspectTerm =
        (b.maxAspect > 0.0 && q.aspect > b.maxAspect) ? b.maxAspect / q.aspect : 1.0;
    return std::clamp(q.scaledJacobian * q.rectangularity * aspectTerm, 0.0, 1.0);
}
 
inline double triangleScaledJacobian(const Point2 &a, const Point2 &b, const Point2 &c)
{
    const Point2 p[3] = {a, b, c};
    double mn = 2.0;
    double area2 = (b.x() - a.x()) * (c.y() - a.y()) - (c.x() - a.x()) * (b.y() - a.y());
    const double orient = area2 >= 0.0 ? 1.0 : -1.0;
    for (int i = 0; i < 3; ++i) {
        const Point2 e1 = p[(i + 2) % 3] - p[i], e2 = p[(i + 1) % 3] - p[i];
        const double l1 = std::hypot(e1.x(), e1.y()), l2 = std::hypot(e2.x(), e2.y());
        if (l1 <= 0.0 || l2 <= 0.0)
            return 0.0;
        mn = std::min(mn, (e2.x() * e1.y() - e2.y() * e1.x()) * orient / (l1 * l2));
    }
    return std::min(mn, 1.0);
}
 
inline double cellCornerAngleDeg(const Vector<MeshVertex> &v, const MeshTriangle &t, int i)
{
    const int n = t.vertexCount();
    const Point2 &prev = v[t.vertex((i + n - 1) % n)].xy;
    const Point2 &cur = v[t.vertex(i)].xy;
    const Point2 &next = v[t.vertex((i + 1) % n)].xy;
    const Point2 e1 = prev - cur, e2 = next - cur;
    return std::atan2(std::abs(e2.x() * e1.y() - e2.y() * e1.x()),
                      e1.x() * e2.x() + e1.y() * e2.y()) *
           180.0 / pi;
}
 
} // namespace meshkit::plan
 
#endif // MESHKIT_NATIVE_MESHQUADQUALITY_HPP