Triangle skew comparison

2025-10-15 22:06 UTC gpt-5-thinking Open in ChatGPT ↗

Trying to check how skewed a triangle is, by comparing its vertices to the vertices of a best fit equilateral triangle.

im checking the differences of differences, as is obvious in the python code. review:

import Rhino.Geometry as rhg

dists = [] side = x[0].DistanceTo(x[1])

for i,p in enumerate(x): if i > 2: break

c = y.ClosestPoint(p)
d = p.DistanceTo(c)
dists.append(d)

diffs = [dists[0]-dists[1], dists[1]-dists[2], dists[2]-dists[0]] diffs = [abs(d) for d in diffs] a = diffs

Fun idea! A few quick notes on your snippet first, then a sturdier way to measure “how equilateral” a triangle is.

Quick review of your code

side = x[0].DistanceTo(x[1]) is never used.
for i,p in enumerate(x): if i > 2: break is redundant if x already has 3 points.
Using ClosestPoint onto y (a curve/polygon?) measures perpendicular distances to the edges of your reference, not vertex-to-vertex correspondence. That can give misleading “skew” if the triangle is just translated/rotated.
Taking “differences of distances” throws away sign and couples errors oddly. A simple, stable scalar score (RMS error, max error, angle variance, etc.) tends to be more interpretable.

A more robust approach

Best-align a template equilateral triangle to your input by a similarity transform (uniform scale + rotation + translation).
Measure the residuals (RMS / max) between the aligned template vertices and your triangle vertices.
(Optional) Also report side-length variance and angle variance—both are zero only for equilateral.

Below is a Grasshopper-friendly Python (RhinoCommon) snippet that:

Fits the triangle’s plane,
Projects to 2D,
Tries all 3 cyclic vertex matchings (so you don’t have to pre-order vertices),
Finds the best uniform-scale Procrustes fit,
Returns RMS error, max error, side/angle deviations, and the aligned equilateral as 3D points.

import Rhino.Geometry as rg
import math

# x: list[rg.Point3d] — your triangle vertices (3 pts)
# returns: dict with metrics + aligned equilateral vertices in 3D
def triangle_equilateral_fit_metrics(x):
    if len(x) != 3:
        raise ValueError("Need exactly 3 points.")

    # --- 1) Fit plane and get 2D coordinates (UV) ---
    success, plane = rg.Plane.FitPlaneToPoints(x)
    if not success:
        raise ValueError("Could not fit a plane to points.")

    to_plane = rg.Transform.PlaneToPlane(plane, rg.Plane.WorldXY)
    from_plane = rg.Transform.PlaneToPlane(rg.Plane.WorldXY, plane)

    pts2d = []
    for p in x:
        q = rg.Point3d(p)
        q.Transform(to_plane)
        pts2d.append(rg.Point2d(q.X, q.Y))

    # Convert to tuples for math
    def p2t(p): return (p.X, p.Y)
    P = [p2t(p) for p in pts2d]  # target triangle (2D)

    # --- 2) Make a centered, unit equilateral template in 2D ---
    # vertices at angles 0, 120, 240; then center to centroid
    T = [(math.cos(a), math.sin(a)) for a in (0.0, 2*math.pi/3, 4*math.pi/3)]
    cx = sum(t[0] for t in T)/3.0
    cy = sum(t[1] for t in T)/3.0
    T = [(t[0]-cx, t[1]-cy) for t in T]

    # helpers
    def centroid(A):
        return (sum(a[0] for a in A)/3.0, sum(a[1] for a in A)/3.0)

    def center(A, c):
        return [(a[0]-c[0], a[1]-c[1]) for a in A]

    def dot(u, v): return u[0]*v[0] + u[1]*v[1]

    def procrustes_similarity(X, Y):
        """
        Find s,R,t to map X -> Y with uniform scale.
        X,Y: lists of 2D tuples, len=3, already corresponded.
        Returns (s, cosθ, sinθ, t=(tx,ty)), and residuals.
        """
        cx = centroid(X); cy = centroid(Y)
        X0 = center(X, cx); Y0 = center(Y, cy)

        # Covariance terms
        a = sum(dot(x, y) for x, y in zip(X0, Y0))
        b = sum(x[0]*y[1] - x[1]*y[0] for x, y in zip(X0, Y0))  # cross term
        normX2 = sum(dot(x, x) for x in X0)
        if normX2 == 0:
            return None

        # rotation (cosθ, sinθ), scale s
        r = math.hypot(a, b)
        if r == 0:
            cos_t, sin_t = 1.0, 0.0
        else:
            cos_t, sin_t = a/r, b/r

        s = r / normX2

        # translation t: Y ≈ s R X + t
        # t = cy - s R cx
        Rc = (cos_t*cx[0] - sin_t*cx[1], sin_t*cx[0] + cos_t*cx[1])
        t = (cy[0] - s*Rc[0], cy[1] - s*Rc[1])

        # Map X -> Yhat and compute residuals
        Yhat = []
        for x in X:
            Rx = (cos_t*x[0] - sin_t*x[1], sin_t*x[0] + cos_t*x[1])
            yhat = (s*Rx[0] + t[0], s*Rx[1] + t[1])
            Yhat.append(yhat)

        errs = [math.hypot(y[0]-yh[0], y[1]-yh[1]) for y, yh in zip(Y, Yhat)]
        rms = math.sqrt(sum(e*e for e in errs)/3.0)
        return (s, cos_t, sin_t, t, Yhat, errs, rms)

    # Try all 3 cyclic matchings between template T and target P
    best = None
    for k in range(3):
        Tk = [T[(i+k)%3] for i in range(3)]
        fit = procrustes_similarity(Tk, P)
        if fit is None:
            continue
        if (best is None) or (fit[-1] < best[-1]):
            best = (k, ) + fit
    if best is None:
        raise ValueError("Alignment failed.")

    k, s, c, s_, t, Yhat, errs, rms = best
    max_err = max(errs)

    # --- 3) Side/angle diagnostics on the input triangle ---
    def sides(A):
        return [
            math.hypot(A[1][0]-A[0][0], A[1][1]-A[0][1]),
            math.hypot(A[2][0]-A[1][0], A[2][1]-A[1][1]),
            math.hypot(A[0][0]-A[2][0], A[0][1]-A[2][1])
        ]

    def angles(A):
        # return internal angles in radians
        def ang(a,b,c):
            # angle at b of triangle abc
            ux, uy = a[0]-b[0], a[1]-b[1]
            vx, vy = c[0]-b[0], c[1]-b[1]
            du = math.hypot(ux, uy); dv = math.hypot(vx, vy)
            if du*dv == 0: return 0.0
            cosv = max(-1.0, min(1.0, (ux*vx + uy*vy)/(du*dv)))
            return math.acos(cosv)
        return [ang(P[2],P[0],P[1]), ang(P[0],P[1],P[2]), ang(P[1],P[2],P[0])]

    S = sides(P)
    meanS = sum(S)/3.0
    side_rel_std = (math.sqrt(sum((s-meanS)**2 for s in S)/3.0) / meanS) if meanS > 0 else 0.0

    Ang = angles(P)
    ang_dev = math.sqrt(sum((a - (math.pi/3))**2 for a in Ang)/3.0)  # rms deviation from 60°

    # --- 4) Build the aligned equilateral in 3D for visualization ---
    # Inverse transform from 2D back to plane, using the Yhat positions:
    eq3d = []
    for (u,v) in Yhat:
        q2d = rg.Point3d(u, v, 0)
        q3d = rg.Point3d(q2d)
        q3d.Transform(from_plane)
        eq3d.append(q3d)

    return {
        "rms_err": rms,          # distance units
        "max_err": max_err,      # distance units
        "side_rel_std": side_rel_std,  # 0 for equilateral
        "angle_rms_dev_rad": ang_dev,  # 0 for equilateral (in radians)
        "equilateral_3d": eq3d   # fitted equilateral vertices (match order of x)
    }

# Example use:
metrics = triangle_equilateral_fit_metrics(x)
a = [metrics["rms_err"], metrics["max_err"], metrics["side_rel_std"], metrics["angle_rms_dev_rad"]]
eq = metrics["equilateral_3d"]  # optional: output for preview

Why this is nicer

Translation/rotation/scale invariant (thanks to Procrustes alignment).
Vertex-to-vertex comparison (no “closest edge” ambiguity).
A single scalar like rms_err or side_rel_std can act as your “skew” score:
- rms_err == 0, side_rel_std == 0, and angle_rms_dev_rad == 0 only for a perfect equilateral.
Works in 3D input because we fit a plane and go to 2D only for the math, then map back.

If you’d like, I can trim outputs to exactly match your a = diffs shape, or plug this into your existing GH component IO.

Here’s what each output means and how to use them.

What each metric means

rms_err (units of length)
Root-mean-square distance between your triangle’s vertices and the fitted equilateral’s vertices (after optimal uniform scale + rotation + translation).
- 0 only for a perfect equilateral that lies exactly on the fitted one.
- Sensitive to overall scale (a larger triangle with the same shape gives larger rms_err).
max_err (units of length)
The largest of the three vertex-to-vertex distances after alignment.
- Highlights the worst local deviation.
- Also scale-dependent.
side_rel_std (dimensionless)
The coefficient of variation of the three side lengths:
$ \text{side_rel_std}=\frac{\sqrt{\frac{1}{3}\sum_i (s_i-\bar s)^2}}{\bar s} $
- 0 for equilateral.
- Scale-invariant and easy to compare across triangles.
- Responds only to side-length imbalance, not to angle distortions that keep sides nearly equal.
angle_rms_dev_rad (radians; convert to degrees if you prefer)
RMS deviation of the three internal angles from $60^\circ$ (i.e., $\pi/3$ rad):
$ \sqrt{\frac{1}{3}\sum_i (\alpha_i-\pi/3)^2} $
- 0 for equilateral.
- Scale-invariant, focuses on angular distortion.
equilateral_3d
The three vertices of the best-fit equilateral in 3D (so you can draw it next to/on top of your input).

Which one to use?

For a single scalar “how equilateral?” score:
Use side_rel_std (dimensionless, intuitive: 0 is perfect, >0 means side imbalance).
If you care about angle shape too, combine it with the angle metric, e.g. $ \text{shape_score} = \sqrt{\text{side_rel_std}^2 + \left(\frac{\text{angle_rms_dev_deg}}{60}\right)^2} $ (now both terms are unitless and roughly comparable).
For visual debugging / viewport overlays:
Use the per-vertex residuals (the three errors that produced rms_err/max_err).
- Draw the fitted equilateral and connect corresponding vertices.
- Color each input vertex by its residual magnitude, or draw short error vectors—this immediately shows where the skew lives.
  Then put side_rel_std as a label/legend number.
When absolute tolerance matters (fabrication, meshing):
Use rms_err or max_err directly (they’re in your model units).

Are they normalized?

side_rel_std: yes (dimensionless, scale-invariant).
angle_rms_dev_rad: yes in the sense it’s scale-invariant, but it’s in radians; you can convert to degrees to be more interpretable.
rms_err / max_err: not normalized (scale-dependent). If you want a unitless version, divide by a characteristic length, e.g. mean side length $\bar s$:
- $ \text{rms_err_rel} = \frac{\text{rms_err}}{\bar s} $
- $ \text{max_err_rel} = \frac{\text{max_err}}{\bar s} $

Practical tips

Show equilateral_3d in the viewport plus error lines from each input vertex → its fitted counterpart; color by residual.
For lists/fields of triangles, sort or filter by side_rel_std (or rms_err_rel) to quickly find the outliers.
Report angles in degrees for readability (e.g., angle_rms_dev_deg = angle_rms_dev_rad * 180/pi).
If two triangles have nearly equal sides but feel “twisted,” angle_rms_dev will reveal that even when side_rel_std is small—use both when you need nuance.

If you want, I can add a couple of lines to the function to return rms_err_rel, max_err_rel, and angle_rms_dev_deg so you have both absolute and normalized metrics ready to plot.

Here you go—compact and Grasshopper-friendly.

import math

# If you already have the function triangle_equilateral_fit_metrics(x) from before,
# this just wraps it to compute the unitless "shape_score".
# Definition recap:
# shape_score = sqrt( side_rel_std^2 + (angle_rms_dev_deg / 60)^2 )

def triangle_shape_score(x, angle_norm_deg=60.0):
    """
    x: list[Rhino.Geometry.Point3d] for the triangle (3 pts)
    angle_norm_deg: normalization for the angular term (default 60°)

    Returns:
        score (float): unitless "how non-equilateral" number (0 = perfect equilateral)
        info (dict): components used to build the score
    """
    metrics = triangle_equilateral_fit_metrics(x)  # from the previous snippet
    side_rel_std = metrics["side_rel_std"]
    angle_rms_dev_deg = metrics["angle_rms_dev_rad"] * 180.0 / math.pi

    score = math.hypot(side_rel_std, angle_rms_dev_deg / angle_norm_deg)

    info = {
        "side_rel_std": side_rel_std,              # dimensionless
        "angle_rms_dev_deg": angle_rms_dev_deg,    # degrees
        "angle_norm_deg": angle_norm_deg,          # normalization used
        "rms_err": metrics["rms_err"],             # absolute residuals if you need them
        "max_err": metrics["max_err"],
        "equilateral_3d": metrics["equilateral_3d"]
    }
    return score, info

# Example GH output:
score, info = triangle_shape_score(x)
a = score                 # single scalar you can visualize/sort by
eq = info["equilateral_3d"]  # optional: draw the fitted equilateral

Notes:

The score is unitless and scale-invariant.
Keep angle_norm_deg=60 for a balanced contribution between side imbalance and angular distortion; tweak it if you want to weight one term more.

Got it—Procrustes (often “Procrustes alignment/analysis”) in our context is just the optimal way to line up two point sets (your triangle vs. an equilateral template) using a similarity transform: one uniform scale + one rotation + one translation. No shear, no non-uniform stretch.

Intuition

First, put both triangles’ centroids at the origin (remove translation).
Then find the rotation that best turns one onto the other.
Also find the single scale that best matches sizes.
After this “best pose,” the remaining distances between corresponding points are the residuals—your “how non-equilateral” signal.

What “best” means

“Best” = least squares: minimize the sum of squared distances between corresponding vertices after the transform.

The tiny bit of math (2D)

Let $X={x_i}$ be the template points (equilateral) and $Y={y_i}$ your triangle points, $i=1..3$.

Center: $ \tilde X_i = x_i - \bar x,\quad \tilde Y_i = y_i - \bar y $

Compute two scalars: $ a=\sum_i \tilde x_i\cdot \tilde y_i,\qquad b=\sum_i \tilde x_i \times \tilde y_i $ (where $\cdot$ is dot product and $\times$ is the 2D scalar “cross”: $x_x y_y - x_y y_x$).
Rotation: $ r=\sqrt{a^2+b^2},\quad \cos\theta=\frac{a}{r},\ \sin\theta=\frac{b}{r} $
Uniform scale: $ s=\frac{r}{\sum_i \lVert \tilde x_i\rVert^2} $
Translation: $ t=\bar y - s,R_\theta,\bar x $ with $R_\theta=\begin{bmatrix}\cos\theta&-\sin\theta\ \sin\theta&\cos\theta\end{bmatrix}$.

Then the fitted points are: $ \hat y_i = s,R_\theta,x_i + t $ Residuals $e_i = \lVert y_i-\hat y_i\rVert$. RMS of ${e_i}$ is your rms_err; max is max_err.

Why this is good

Invariant to where/which way/what size your triangle is—only shape matters.
Closed form, fast, stable (no iterative solver needed).
For triangles, you also try the 3 cyclic correspondences (ABC↔ABC, BCA, CAB) and keep the smallest error so vertex ordering doesn’t matter.

In our code

We center both sets, compute $a,b$, get $\theta$ and $s$, then $t$.
We map the equilateral to your triangle, compute per-vertex residuals → from these we derive rms_err/max_err, and also build the aligned equilateral for display.

That’s Procrustes in a nutshell: “best pose” match under (scale + rotate + translate), leaving only shape discrepancy in the residuals.