I am working on documenting and editing the source, as well as doing some work
in SDL.

Perhaps someone can help me out here.
I have some ideas about using and testing the control points of the spline.

POV-Ray complains when I have a control point with an x component that is
negative.  However, the only place that I can find the error, and the conditions
that trigger it - are here:

Starting at line 1005 of sor.cpp
https://github.com/POV-Ray/povray/blob/master/source/core/shape/sor.cpp

        if ((fabs(P[i+2][Y] - P[i][Y]) < EPSILON) ||
            (fabs(P[i+3][Y] - P[i+1][Y]) < EPSILON))
        {
            throw POV_EXCEPTION_STRING("Incorrect point in surface of
revolution.");
        }

But that code only uses the y-values, and only tests the size of the difference
between the y values of 2 sets of control points.  It certainly does NOT include
the x-value of a control point.

So then how is the error getting triggered if there is no testing the x-value???

- BE


sor {
 6
 <-2.0, 0.5>,
 <-1.0, 1>,
 <-0.1, 2>,
 <-0.2, 3>,
 <-1.0, 4>,
 <-2.0, 5>
 pigment {rgb z}
}

Post a reply to this message

"Bald Eagle" <cre### [at] netscapenet> wrote:

> POV-Ray complains when I have a control point with an x component that is
> negative.

> So then how is the error getting triggered if there is no testing the x-value???

Thanks to jr hunting this down while I was busy with other things!

parser.cpp starting at line 5665
if ((Points[i][X] < 0.0) ||
     ^^^^^^^^^^^^^^^^^^^ and there we go.

            ((i > 1 ) && (i < Object->Number - 1) && (Points[i][Y] <=
Points[i-1][Y])))
        {
            Error("Incorrect point in surface of revolution.");
        }

The more I do this, the more I see exactly what clipka meant about the parser.

All of these tests and guards ought to be isolated in sor.cpp, not smeared all
over the source codebase.

- BE

Post a reply to this message

Here's an overview of what I'm planning on doing with sor {}.

in parser.cpp, I'll be removing
if (P[i][X] < 0.0)

In part because it will allow placing the tangent-determining control points
anywhere,
and because I plan on extending the rendering of the sor to the parts that are
currently
not rendered when they cross or get interpolated across the y-axis into -x

The control points will be used to derive the cubic equation of the spline, and
then the
curve equation to will be evaluated to determine if the curve crosses the
y-axis.  If so,
a warning will be issued.
Then the roots (the points where the curve crosses the y-axis) will be found.
The Catmull-Rom spline will be subdivided by creating new subsegments up to the
roots, between the roots, and beyond the roots.

The cubic equation coefficients get stored as A, B, C, and D.
In regions where the spline goes negative in X will have all of the signs of the
coefficients flipped.
That will yield a valid region of the spline segment.

The whole spline will therefore be visibly rendered.

Guard checks on the x & y values, and the height of the spline segments that
include the tangent-controlling points will be fixed to remove the unnecessary
errors.

I have plans to emit the volume of the sor object, and a way to uv-map the
curvature of the surface. (see attached)

Long-term, it would be great to allow defining the sor profile as a spline {}
object and using that spline in the sor {} block syntax.
That would allow using linear and quadratic splines in the sor as well.
The spline would allow easy graphing the cross-section of the sor if desired.
This would also create a new centripetal Catmull-Rom spline type.

Potentially, the positive and negative portions of the sor spline could be
rendered independently.
Allowing descending y-values for the sor would be nice, to allow for
conceptually modelling a sor in the opposite direction.

Post a reply to this message

Attachments:
Download 'sor_curvature_scene_2026-05-13.png' (53 KB)

Preview of image 'sor_curvature_scene_2026-05-13.png'

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

"boltcapt" <bol### [at] aolcom> wrote:

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

I think that after going through all of that, you may have added some technical
points and terminology/references of which I was not aware.
Those would be good to have in the documentation.

The larger project, is of course, implementing the fixes discovered for:

The unrendered portions of the revolved spline, due to those portions of the
interpolated Catmull-Rom spline being negative

The first control point not being able to be placed at the same Y value as
another control point (possibly raise a warning and maybe auto-fix)

And the false positive bounding tests that call unnecessary root-solving.

Related to root-solving, there would be the issue of possibly implementing
tighter bounding by the "bounding by proxy Bezier spline method" mentioned in
the thread.

That last could possibly solve some related problems with the sphere_sweep {}
object.

The deep-dive into the how and why of the sor {} and reimplementing it all in
SDL along with various fixes too a lot, and I don't know that currently have
enough mental energy and, more importantly, uninterrupted time to scale the
steep learning curve of learning how to competently code in c++, set up a
machine, and do all of the iterative compiling of edited source and bug-hunting.


So,
The proposal/request would be that you could identify the exact files that need
to be edited, and write the fixes based on my SDL versions, and I could follow
along and try to make sure it all jives with my copious notes,

or

I could edit the source, and you could just do the compiling and supervise the
edits "No, you shouldn't do it that way because" or similar dev code coaching,

or

some mix of the two approaches.

The end goal would be to produce new edited source with excellent commentary and
an eye toward allowing any future changes and extensions without having to
discover fire, invent the wheel, and reverse-engineer the sor {} source code all
over again.

It sounds like a lot, and probably is (isn't everything), but I also think that
since most of the theoretical heavy lifting, and implementation in SDL is
already done, that it might also be a bit easier (in the broad-stroke sense)
than it seems at first glance.

Maybe we can start a new thread and we can discuss the details of what needs
doing and go from there.   I can outline things, you can ask questions, and we
can hammer out a game plan before even touching any code.

Thank you,

Bill "Bald Eagle" Walker

Post a reply to this message