Here is what I have. It is based on this source:
- `povray-3.7.0.10/source/backend/shape/sor.cpp` - `Sor::Compute_Sor()` (the
spline fit), `Sor::Intersect()`, `Sor::Inside()`, `Sor::Normal()`
- `povray-3.7.0.10/source/backend/parser/parse.cpp` - `Parser::Parse_Sor()`

The upstream `master/source/core/shape/sor.cpp` is the same algorithm,
only renamed/namespaced.

---

The SOR curve is a **single-valued function of height**: radius-squared as a
function of height,

r^2 = f(h) = A*h^3 + B*h^2 + C*h + D        (one cubic per segment)

(`sor.cpp:47-49`). It stores `r^2`, not `r`, because the ray-surface equation
`(Pu + k Du)^2 + (Pv + k Dv)^2 = f(Pw + k Dw)` is then a **cubic in k**
(`sor.cpp:51-56`, set up at `sor.cpp:370-381`). That is the "solving a cubic
polynomial" claim on the wiki - it is true *for the SOR primitive specifically*
because the SOR spline is always this height-cubic. (See section 6 on the lathe
comparison you flagged.)

Because it is a function of height, **interior control points must strictly
increase in height** - enforced at `parse.cpp:5816-5820`. You cannot make a
curve that doubles back in h with a SOR; that is what `lathe` is for.

You asked "which control points get used?"

`Parse_Sor()` (`parse.cpp:5780-5843`):

- Requires **>= 4 points** (`parse.cpp:5799`).
- Reads `Number` points, each a 2D vector `<radius, height>` (x = radius,
  y = height; z unused) (`parse.cpp:5808-5821`).
- Then **`Object->Number -= 3`** (`parse.cpp:5839`) - comment: *"There are
  Number-3 segments!"*

So for **N input points there are N-3 segments.** Indexing the points P[0..N-1]:

- The **drawn surface runs from P[1] (base) to P[N-2] (cap)**.
- **P[0] and P[N-1] are "phantom" points** - they are never on the surface.
  They exist only to supply the end tangents (`sor.cpp:75-77`).

Segment `i` (`i = 0 .. N-4`, loop at `sor.cpp:1028`) is bounded by **P[i+1]
(lower) and P[i+2] (upper)**, and additionally reads **P[i] and P[i+3]** as the
outboard neighbors used purely for slope estimation. So **four control points
feed every segment**: the two endpoints plus one neighbor on each side. That is
the Catmull-Rom stencil, and it is the direct answer to "which control points
get used."

The interior tangents are **central finite differences** over the two
neighbors. In `Compute_Sor` (`sor.cpp:1038-1044`):

k[0] = P[i+1].x^2;                                  // r1^2  = f(h1)
k[1] = P[i+2].x^2;                                  // r2^2  = f(h2)
k[2] = (P[i+2].x - P[i].x)   / (P[i+2].y - P[i].y); // dr/dh at lower point
(secant of its neighbors)
k[3] = (P[i+3].x - P[i+1].x) / (P[i+3].y - P[i+1].y);// dr/dh at upper point
k[2] *= 2.0 * P[i+1].x;                             // chain rule: d(r^2)/dh =
2r * dr/dh
k[3] *= 2.0 * P[i+2].x;

Using P1,P2 for the segment endpoints and P0,P3 for the outboard neighbors
(heights h0..h3, radii r0..r3):

- tangent at the lower point = `(r2 - r0)/(h2 - h0)`
- tangent at the upper point = `(r3 - r1)/(h3 - h1)`

That is the **non-uniform Catmull-Rom tangent** (the cardinal-spline tension-0
case), written in finite-difference form for unequally spaced knots. Two
nuances worth putting in the docs:

1. The tangent is computed in **radius** space, then converted to **radius^2**
   space by the chain rule `d(r^2)/dh = 2r*dr/dh` (the `2.0 * P[..].x` factors).
   So the fit is a cubic Hermite in `(h, r^2)` whose slopes come from a
   Catmull-Rom estimate in `(h, r)`.
2. Tangents use the actual heights (`h2-h0`, `h3-h1`), i.e. it is the
   **non-uniform / centripetal-style parameterization by h**, not the uniform
   `(P_{i+1}-P_{i-1})/2` form.

This is exactly the Andrew-Hung "Catmull-Rom in plain English" construction you
found - applied to r^2 vs h.

I think this gets at what you are looking for:

For each segment we have a cubic with 4 unknowns `[A, B, C, D]` and 4 Hermite
constraints (value + slope at each endpoint):

f(h1)  = r1^2
f(h2)  = r2^2
f'(h1) = 2*r1*(r2 - r0)/(h2 - h0)
f'(h2) = 2*r2*(r3 - r1)/(h3 - h1)

with `f(h)  = A h^3 + B h^2 + C h + D`
and  `f'(h) = 3A h^2 + 2B h + C`.

Writing the four constraints as `M x = b`:

        | h1^3   h1^2   h1   1 |   | A |   | r1^2          |   | k0 |
  M x = | h2^3   h2^2   h2   1 | * | B | = | r2^2          | = | k1 | = b
        | 3h1^2  2h1    1    0 |   | C |   | f'(h1)        |   | k2 |
        | 3h2^2  2h2    1    0 |   | D |   | f'(h2)        |   | k3 |

This is **built verbatim** at `sor.cpp:1046-1068`:

| code                                      | row              | meaning
     |
|-------------------------------------------|------------------|----------------------|
| `Mat[0]` with `w=P[i+1].y` (`:1048-1051`) | `[w^3 w^2 w 1]`  | value at lower
point |
| `Mat[1]` with `w=P[i+2].y` (`:1060-1063`) | `[w^3 w^2 w 1]`  | value at upper
point |
| `Mat[2]` with `w=P[i+1].y` (`:1053-1056`) | `[3w^2 2w 1 0]`  | slope at lower
point |
| `Mat[3]` with `w=P[i+2].y` (`:1065-1068`) | `[3w^2 2w 1 0]`  | slope at upper
point |

**The key point the wiki omits:** in `b = M x`,
- **`b` is the *known* right-hand side** - the radius^2 values and the slopes,
  computed straight from the control points (the `k[]` array).
- **`M` is *known*** - just powers of the two endpoint heights.
- **`x = [A,B,C,D]` is the *unknown*** we are solving for.

So you solve it the only way you can: **`x = M^-1 * b`.** your instinct
("x = b/M") is correct - for a matrix that "division" is inversion. The code
does precisely this:

MInvers(Mat, Mat);                                   // sor.cpp:1070  -> Mat =
M^-1
A = k0*Mat[0][0] + k1*Mat[0][1] + k2*Mat[0][2] + k3*Mat[0][3];  // :1074  (row 0
of M^-1 * b)
B = ... Mat[1] ...                                   // :1075
C = ... Mat[2] ...                                   // :1076
D = ... Mat[3] ...                                   // :1077

So the wiki's "calculate b, then b = M x" reads backwards because it never says
*x is the unknown and you invert M*. The doc should state: build M from the two
heights, build b from the control-point radii and the Catmull-Rom slopes, then
`[A B C D]^T = M^-1 b`.

You had a concern about division by zero. The zeros in rows 3 and 4
(`Mat[2][3]=0`, `Mat[3][3]=0`, `sor.cpp:1056,1068`)
are simply the **D-column of the derivative rows**: differentiating
`A h^3+B h^2+C h+D` kills the constant term, so the coefficient on D in `f'(h)`
is
0. They are *entries of the matrix*, not *divisors*. Inversion (`MInvers`) is
Gaussian elimination; it never divides by those structural zeros.

The matrix is non-singular as long as the **two endpoint heights differ**,
`h1 != h2`, which is guaranteed by the strict-monotonic-height rule
(`parse.cpp:5816-5820`) and re-checked against degenerate neighbor spacing at
`sor.cpp:1030-1034` (`fabs(P[i+2].y - P[i].y) < EPSILON` -> error). So there is
no division by zero. (The only divisions in the whole construction are the two
finite-difference denominators `h2-h0` and `h3-h1`, both protected by that same
check.)

After solving, coefficients with magnitude `< EPSILON` are snapped to 0
(`sor.cpp:1079-1082`) purely to speed later polynomial solves.

You are right about the two meanings of h. The doc reuses one symbol for
two different things:

- In **M and b**, `h` (the code's `w = P[i+1].y`, `P[i+2].y`) is a **fixed
  number**: the height of a specific control point. These are constants of the
  linear system.
- In **`f(h) = A h^3 + B h^2 + C h + D`**, `h` is the **free running variable**:
  any height inside the segment `[h1, h2]`. Once A,B,C,D are known you evaluate
  the cubic at any such h to get r^2 there. That evaluation **is** the
  interpolation you said you "couldn't see."

Doc fix: use distinct symbols, e.g. `h_j , h_{j+1}` for the endpoint heights in
M/b, and a separate `h` (or `w`) for the argument of `f`.

Where the running form is used in the engine, for cross-reference:
- `Inside()` evaluates `f(h)` by Horner: `r = h*(h*(h*A+B)+C)+D`
(`sor.cpp:479`).
- `Normal()` uses `f'(h)/2`: `0.5*(h*(3A h+2B)+C)` -> normal `(x, -f'/2, z)`
  (`sor.cpp:547-551`).
- Cap/base radii are `f(h)` at the top/bottom heights (`sor.cpp:1143-1167`).

You are right that the wiki sentence about the sor vs lathe polynomial order.
It is imprecise. The SOR primitive is *always* the
the height-cubic of section 1, so its ray test is *always* a cubic (degree 3,
`Solve_Polynomial(3, ...)` at `sor.cpp:381`). The `lathe` object, by contrast,
supports `linear_spline`, `quadratic_spline`, `cubic_spline`, and
`bezier_spline` (`parse.cpp:2950-2956`), and its splines are parametric `(r(t),
h(t))` rather than `r^2=f(h)`. The "6th-order polynomial" figure applies only to
lathe's **cubic** spline case, not to lathe in general. The doc should say "a
cubic lathe segment requires a 6th-degree solve" rather than implying every
lathe does.

I hope this captures what you were looking for. If not, let me know so I can
refocus.

Post a reply to this message