meshquadquality.h

Go to the documentation of this file.

 
#ifndef OPENSWMMVIS_MESH_MESHQUADQUALITY_H
#define OPENSWMMVIS_MESH_MESHQUADQUALITY_H
 
#include "mesh/meshresult.h"
 
#include <QPointF>
#include <QVector>
#include <algorithm>
#include <cmath>
#include <limits>
 
namespace mesh {
 
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 QPointF &a, const QPointF &b,
                               const QPointF &c, const QPointF &d)
{
    QuadQuality q;
    const QPointF p[4] = {a, b, c, d};
    double area2 = 0.0;
    for (int i = 0; i < 4; ++i)
    {
        const QPointF &u = p[i], &v = p[(i + 1) % 4];
        area2 += u.x() * v.y() - v.x() * u.y();
    }
    q.area = 0.5 * std::abs(area2);
    const double orient = area2 >= 0.0 ? 1.0 : -1.0;
 
    q.minAngleDeg = 360.0; q.maxAngleDeg = 0.0;
    q.scaledJacobian = 2.0;
    q.convex = true;
    double maxDev = 0.0;
    double side[4];
    for (int i = 0; i < 4; ++i)
    {
        const QPointF &prev = p[(i + 3) % 4], &cur = p[i], &next = p[(i + 1) % 4];
        const QPointF e1 = prev - cur, e2 = next - cur;
        const double l1 = std::hypot(e1.x(), e1.y()), l2 = std::hypot(e2.x(), e2.y());
        side[i] = l2;
        if (l1 <= 0.0 || l2 <= 0.0) { q.convex = false; q.scaledJacobian = 0.0; continue; }
        const double cross = e2.x() * e1.y() - e2.y() * e1.x();   // sign: + when turning consistently with `orient`
        const double dot   = e1.x() * e2.x() + e1.y() * e2.y();
        const double ang   = std::atan2(std::abs(cross), dot) * 180.0 / M_PI;
        const double sj    = cross * orient / (l1 * l2);           // signed sine of the interior angle
        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.scaledJacobian = std::min(q.scaledJacobian, 1.0);   // rounding can give sin θ = 1 + ε for a rectangle
    q.rectangularity = std::max(0.0, 1.0 - maxDev / 90.0);
    const double s02 = 0.5 * (side[0] + side[2]), s13 = 0.5 * (side[1] + side[3]);
    const double lo = std::min(s02, s13), hi = std::max(s02, s13);
    q.aspect = lo > 0.0 ? hi / lo : std::numeric_limits<double>::infinity();
    const QPointF d1 = c - a, d2 = d - b;
    const double ld1 = std::hypot(d1.x(), d1.y()), ld2 = std::hypot(d2.x(), d2.y());
    q.skew = (ld1 > 0.0 && ld2 > 0.0)
               ? std::abs((d1.x() * d2.x() + d1.y() * d2.y()) / (ld1 * ld2)) : 1.0;
    return q;
}
 
inline QuadQuality quadQuality(const QVector<MeshVertex> &vertices, const MeshTriangle &t)
{
    return quadQuality(vertices[t.v0].xy, vertices[t.v1].xy,
                       vertices[t.v2].xy, vertices[t.v3].xy);
}
 
inline bool quadAcceptable(const QuadQuality &q, const QuadQualityBounds &b)
{
    if (!q.convex) 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 (!q.convex || q.scaledJacobian <= 0.0) return 0.0;
    const double aspectTerm = (b.maxAspect > 0.0 && q.aspect > b.maxAspect)
                                  ? b.maxAspect / q.aspect : 1.0;
    return q.scaledJacobian * aspectTerm * std::sqrt(std::max(0.0, 1.0 - q.skew));
}
 
inline double triangleScaledJacobian(const QPointF &a, const QPointF &b, const QPointF &c)
{
    const QPointF 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 QPointF 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 QVector<MeshVertex> &v, const MeshTriangle &t, int i)
{
    const int n = t.vertexCount();
    const QPointF &prev = v[t.vertex((i + n - 1) % n)].xy;
    const QPointF &cur  = v[t.vertex(i)].xy;
    const QPointF &next = v[t.vertex((i + 1) % n)].xy;
    const QPointF 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 / M_PI;
}
 
} // namespace mesh
 
#endif // OPENSWMMVIS_MESH_MESHQUADQUALITY_H