On 9/10/25 11:16, William F Pokorny wrote:
> If you're referring to the sqrt I think you are, you've guessed 
> correctly! 🙂 I have a potential fix for it. The sqrt(-) along with some 
> fine detail in how std::min(), std::max() work when passed a -nan value 
> as the first argument vs the second argument cause the segment blink out 
> issue. Only limited testing thus far. Maybe I'll get to more tonight.

Attached an image of some more testing. It includes the fix for the 
segment blink out issue & that fix looks good thus far. For testing I 
disabled all the point list sanity checking to see what works and not.

---

Not going to go through the image in detail, but in answer to your 
question about on curve points going negative it looks like a no go 
without a deep dive and probable re-work of the current code.

For reasons I don't understand, when I allow on curve points to go 
negative the camera ray direction affects results. This shown in the 
fourth row where the sor slowly disappears then reappears as rotated 
about x.

---

The descending in y on curve points doesn't work at all.

---

Allowing the first and last control points to go negative is fine and I 
plan to add that relaxing to the blank out fix for yuqk's next release.

---

You had mentioned not being able to flip the two control points in y. 
Where everything I'd tried to that point seemed to work. You also 
mentioned this bit of code in sor.cpp :

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

I'm finding it is this run time check which sometimes doesn't allow 
control points flipping in y. For yuqk what will be in the next release 
is some added error text when that check trips so it is easier to 
understand what is wrong.

Error! Problem point in sor segment 0 point set.
fabs(P[2][Y] - P[0][Y] || fabs(P[3][Y] - P[1][Y])
fabs(0.5 - 0.5 || fabs(0.9 - 0) < 4.44089e-08

The center images in the bottom two rows(*) of the attached image trip 
this check and error. So, that code is catching some 'particular' cases 
where given the control point locations relative to other on curve 
points doesn't render correctly.

(*) Those rows are start the control points at one extreme and slowly 
flip them to the opposite extreme in y moving left to right in the lower 
two rows.

---

So with the segment blink out fix, it does look like more could be 
allowed as valid point sets for the sor & yuqk will adopt those 
relaxations. The shapes with negative control points do break up more 
often, but as we found this can happen today with certain point lists.

(I've not really looked at intersections and differences as yet)

Bill P.

Aside: Somewhere I should mentioned too that that blink out fix corrects 
the blink out, but in certain cases it also slightly corrects bounding 
cylinders where due how min/max work there was no blink out. In other 
words some fringe fixes/differences might appear too due the blink out fix.

Post a reply to this message

Attachments:
Download 'sor_testing.png' (114 KB)

Preview of image 'sor_testing.png'

William F Pokorny <ano### [at] anonymousorg> wrote:

> Attached an image of some more testing. It includes the fix for the
> segment blink out issue & that fix looks good thus far.

So what did you do?

in line 1082
x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);

I'd suggest changing that something like

x[n] = sqrt( fabs(r[n] * (r[n] * (r[n] * A + B) + C) + D) );


> Not going to go through the image in detail, but in answer to your
> question about on curve points going negative it looks like a no go
> without a deep dive and probable re-work of the current code.

I'd suggest the above, if you haven't already done that.

Also, it may have something to do with calculating the intersections.

There's a check on Radius2 that I saw, so see what happens when you change:

Lines 1110-1111
Radius1 = xmin;
Radius2 = xmax;

to use fabs()

> For reasons I don't understand, when I allow on curve points to go
> negative the camera ray direction affects results. This shown in the
> fourth row where the sor slowly disappears then reappears as rotated
> about x.


> The descending in y on curve points doesn't work at all.

Not surprising - all the checks for the control point order need to be
rewritten.


> You had mentioned not being able to flip the two control points in y.
> Where everything I'd tried to that point seemed to work.

I'll have to do some more experimenting when I get the chance.

> You also
> mentioned this bit of code in sor.cpp :
>
> if ((fabs(P[i+2][Y] - P[i][Y]) < gkMinIsectDepthReturned) ||
>      (fabs(P[i+3][Y] - P[i+1][Y]) < gkMinIsectDepthReturned))
> {
> throw POV_EXCEPTION_STRING("Incorrect point in surface of revolution.");
> }
>
> I'm finding it is this run time check which sometimes doesn't allow
> control points flipping in y.

That def seems to be related to what you're seeing in those last 2 rows - since
the cp's are nearly coincident in y.


> The shapes with negative control points do break up more
> often, but as we found this can happen today with certain point lists.

Do you think as a quick fix, we could just add -ymin to all the values to
translate everything into the first quadrant, starting at y=0, and then
translate everything back after the math gets done?

Very nice that some of this is getting worked out after 30 years.  :)

- BW

Post a reply to this message

On 9/11/25 10:36, Bald Eagle wrote:
>> Attached an image of some more testing. It includes the fix for the
>> segment blink out issue & that fix looks good thus far.
> So what did you do?

The current fix is:

     if ((r[n] >= y[0]) && (r[n] <= y[1]))
     {
//      x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);
         double tmpVal = (r[n] * (r[n] * (r[n] * A + B) + C) + D);
         if (tmpVal>=0.0)
             x[n] = sqrt(tmpVal);
         else
         {
             tmpVal = std::abs(tmpVal);
             x[n] = sqrt(tmpVal);
             x[n] *= -1.0;
         }
     }

As for the bounding suggestion, I don't know. Of note is that the 
Radius1 value is never used that I see.

Something I did try was faking a user bounding radius multiplier option. 
In some configurations where only a subset of the on curve points were 
negative in x. It often helped with apparent clipping - but there always 
remained some issue with the the ray directions relative to the overall 
sor - where all or part of sor 'segment(s)' would drop away (this an 
issue apart from the entire-segment blink out issue up top).

My wild guess. There are internally calculated A,B,C,D values used for 
the quadratic eqn solved while setting the two bounding radii values. 
Those calculated values involve the x values as given for the sor point 
set. Maybe something could go 'more right' there when negative values 
are present - I don't know.

Bill P.

Post a reply to this message

I'm still working on unraveling the source code for
bool Sor::Intersect, line 239 in sor.cpp

Line 276 has the interesting expression:  r0 = P[X] * D[Z] - P[Z] * D[X];
and I wasn't sure what that exactly did, especially since the comment that
precedes it is
 /* Get distance of the ray from rotation axis (= y axis). */
and that didn't look like any distance formula that I've ever seen.
The minus sign certainly helped obscure its meaning and purpose.

Apparently it's a sort of 2D cross-product, or "perp-product"

https://s.goessner.net/articles/crossProductHarmful.html
Also being covered in Graphics Gems IV:
Hill, F. S., Jr., The Pleasures of `Perp Dot' Products, p. 138-148.

And since I'm speculating that there are other examples this to be found in
source, I'm documenting this fully here for future reference.

immediately following that, I have a comparison of the length of the ray
direction vector D with some variable "a", which due to poor source code
documentation, I have no idea what "a" exactly is, or why if they are equal that
r0 has to get divided by the sqrt of a.

    #if ((a = Dx * D.x + D.z * D.z) > 0.0)
            #local r0 = r0/sqrt(a);
    #end

I also have Entry->A (through D) which I need to figure out.

- BW

Post a reply to this message

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

I have a functional, if not yet rigorous understanding about the ray-radius
intersection test(s).  I am hoping to understand it completely enough to make
one of my explanatory diagrams.

> immediately following that, I have a comparison of the length of the ray
> direction vector D with some variable "a", which due to poor source code
> documentation, I have no idea what "a" exactly is, or why if they are equal that
> r0 has to get divided by the sqrt of a.
>
>     #if ((a = Dx * D.x + D.z * D.z) > 0.0)
>             #local r0 = r0/sqrt(a);
>     #end

I believe I see what's going on here.
in c++, == is a comparison.
= is an assignment.

So the #if line is simultaneously assigning the evaluation of the expression to
a, and then comparing it's value to 0.0.

> I also have Entry->A (through D) which I need to figure out.

This has to do with the arrow operator, or "Class Member Access Operator".

In sor.h, we have

struct Sor_Spline_Entry_Struct final
{
    DBL A, B, C, D;
};

and

struct Sor_Spline_Struct final
{
    int References;
    SOR_SPLINE_ENTRY *Entry;
    BCYL *BCyl;                 /* bounding cylinder.                  */
};

So "Entry" is a pointer to SOR_SPLINE_ENTRY, which has class members A through
D.

"Entry->A" is a whole identifier, which snakes through the 2 structs to identify
and use A.

and

class Sor final : public ObjectBase
{
    public:
        int Number;
        SOR_SPLINE *Spline;      /* List of spline segments     */
.... etc.

So A through D have something to do with the spline segments.
I just don't "see" where they come from or what their exact meaning is.
Perhaps somewhere else farther upstream in the code is a place where the control
points get grouped in sets of 4, and so A through D would be sets of control
points.

- BW

Post a reply to this message

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

> So A through D have something to do with the spline segments.
> I just don't "see" where they come from or what their exact meaning is.
> Perhaps somewhere else farther upstream in the code is a place where the control
> points get grouped in sets of 4, and so A through D would be sets of control
> points.

And this is the challenge of trying to do this piecemeal, in fits and starts
during random available 15 min opportunistic time slots.

Just a few lines above that is:

/* Step through the list of intersections. */

    for (j = 0; j < cnt; j++)
    {
        /* Get current segment. */

        Entry = &Spline->Entry[intervals[j].n];

and I'm kinda guessing that intervals comes from insert_hit in
https://github.com/POV-Ray/povray/blob/master/source/core/bounding/boundingcylinder.cpp

So that's my current working theory.

Post a reply to this message

I also cannot find the following 2 functions in the code repository.
I can't search without being logged in, and username/pswd are scribbled down
somewhere at home.

    MInvTransPoint(P, ray.Origin, Trans);

    MInvTransDirection(D, ray.Direction, Trans);

Anyone?

Post a reply to this message

hi,

"Bald Eagle" <cre### [at] netscapenet> wrote:
> I also cannot find the following 2 functions in the code repository.
> I can't search without being logged in, and username/pswd are scribbled down
> somewhere at home.
>
>     MInvTransPoint(P, ray.Origin, Trans);
>
>     MInvTransDirection(D, ray.Direction, Trans);
>
> Anyone?

core/math/matrix.h has the prototypes, hth.

and thx re other thread.


regards, jr.

Post a reply to this message

"jr" <cre### [at] gmailcom> wrote:
> hi,
>
> "Bald Eagle" <cre### [at] netscapenet> wrote:
> > I also cannot find the following 2 functions in the code repository.
> > I can't search without being logged in, and username/pswd are scribbled down
> > somewhere at home.
> >
> >     MInvTransPoint(P, ray.Origin, Trans);
> >
> >     MInvTransDirection(D, ray.Direction, Trans);
> >
> > Anyone?
>
> core/math/matrix.h has the prototypes, hth.
>
> and thx re other thread.
>
>
> regards, jr.

Thanks, jr!

Ah, I think I (mostly) understand now.
matrix.h declares those inline functions for use in the namespace, and the
MInvTransPoint and MInvTransDirection functions just use the
MTransPoint and MTransDirection functions but somehow using the multiplicative
inverse of the transform matrix as input instead of the regular matrix.

in matrix.cpp, (line 535) we have

void Compute_Scaling_Transform (TRANSFORM *result, const Vector3d& vector)
{
    MIdentity (result->matrix);

    (result->matrix)[0][0] = vector[X];
    (result->matrix)[1][1] = vector[Y];
    (result->matrix)[2][2] = vector[Z];

    MIdentity (result->inverse);

    (result->inverse)[0][0] = 1.0 / vector[X];
    (result->inverse)[1][1] = 1.0 / vector[Y];
    (result->inverse)[2][2] = 1.0 / vector[Z];
}

and at line 155, there is:

void MIdentity (MATRIX result)
{
    int i, j;

    for (i = 0; i < 4; i++)
    {
        for (j = 0; j < 4; j++)
        {
            if (i == j)
            {
                result[i][j] = 1.0;
            }
            else
            {
                result[i][j] = 0.0;
            }
        }
    }
}

I'm still trying to get a feel for all of this, having only ever coded in "c++"
(Arduino) for maybe 6 months, several years back.

(Currently reading www.learncpp.com in small spurts, then I'll probably start
going through things with Stroustrup in hand to more rigorously work things
out.)

I guess perhaps in the same way that we can have operator overloading, where I
have two functions with the same name, but different input types,
a single function can yield one of several different types of outputs if that
class member access operator -> is used ...?

- BW

Post a reply to this message

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

> bool Sor::Intersect, line 239 in sor.cpp
>
> Line 276 has the interesting expression:  r0 = P[X] * D[Z] - P[Z] * D[X];

>  /* Get distance of the ray from rotation axis (= y axis). */
....
> Apparently it's a sort of 2D cross-product, or "perp-product"
....
>     #if ((a = Dx * D.x + D.z * D.z) > 0.0)
>             #local r0 = r0/sqrt(a);
>     #end

OK, so r0 gets calculated by taking the vector cross product of
the PROJECTIONS of the ray origin and the ray direction (after they have been
translated into "object space")
[still hunting down how that Trans gets calculated, but I have an idea]

So, since they are projections onto the xz plane, y=0 for both vectors.
that means that two of the 2x2 determinants drop out due to multiplication by 0,
leaving only the one determinant to multiply by the y basis vector,
giving a y vector that is parallel to the 2 projected vectors.

Thus the shortened form P[X] * D[Z] - P[Z] * D[X]

However, in working this out, I had my vectors in different places in the
matrix, and so I got P[Z] * D[X] - P[X] * D[Z] which should just be the opposite
sign of the calculation in code.  Not really a worry, since it's just a vector
cross product.

The potential problem that I see is the test that uses this value:

    /* Test if ray misses object's bounds. */

    #if (r0 > Radius2)                      // Radius2 is xmax of the sor CP's
        #local Result = false;
    #else

If the sign of the vector cross product is negative, then that rejection test
will always fail, and the ray will always move on to have the intersection
tested.   Which, if I'm interpreting all of this correctly, means that all of
the rays on one side of the sor always get fully evaluated for intersections,
this potentially making rendering slower than it needs to be.

Shouldn't we be using at least abs (r0) to account for P-cross-D AND D-cross-P?



In more detail:

First we translate and then normalize the ray direction vector.

Then we test the ray to see if it is above or below the sor
as well as to the left or right of the sor

starting at line 266
    if (((D[Y] >= 0.0) && (P[Y] >  Height2)) ||
        ((D[Y] <= 0.0) && (P[Y] <  Height1)) ||
        ((D[X] >= 0.0) && (P[X] >  Radius2)) ||
        ((D[X] <= 0.0) && (P[X] < -Radius2)))
    {
        return(false);
    }

Line 266 checks if the ray direction is flat or points up, and the ray origin is
higher than the highest sor control point (Height2 = ymax)
So it will never intersect the rotated sor spline
This is OK

Line 267 checks if the ray direction is flat or points down, and the ray origin
is lower than the lowest sor control point (Height1 = ymin)
So it will never intersect the rotated sor spline
This is OK

Line 268 checks if the ray direction is straight or points right, and the ray
origin is to the right of the rightmost sor control point (Radius2 = xmax)
So it will never intersect the rotated sor spline
This is OK

Line 269 checks if the ray direction is straight or points left, and the ray
origin is to the left of the leftmost (rotated) sor control point (-Radius2 =
-xmax)
So it will never intersect the rotated sor spline
This is OK

BW: since this is before we do any projection, and we're only testing height and
width, can't we have a slanted ray that WOULD intersect the sor spline once it
got rotated around the y-axis?

Then we are doing a calculation that purports to be the DISTANCE between the ray
and the y-axis.

Working that out from first principles, it is the geometric problem of
determining the perpendicular distance (the shortest distance) between a point
and a line.
We simplify this by projecting onto 2D since the vertical angle of the ray makes
no difference at this point.
The vector cross product of 2 vectors gives us the area of the parallelogram
defined by those 2 vectors (just use those 2 vectors to make the remaining 2
sides)
The area is also equal to the base times the height of the parallelogram, and so
we solve for the height (H), which is the perpendicular - and therefore the
shortest distance between the point (y-axis in the xz plane, <0, 0, 0>) and the
line (the ray composed of origin and unit directional vector)
That gives us LENGTH (P-cross-D) = H times LENGTH (D)
Rearranging gives us H = LENGTH (P-cross-D) / LENGTH (D)

Now, in code, we have
r0 = P.x * D.z - P.z * D.x;
and that is only multiplied by the y unit basis vector,
so since there are no other dimensional units, this directly evaluates to the
length of the perpendicular y vector, and we don't have to do the expensive
square root of the sum of the squares.

a is the squared length of the direction vector
which, if non-zero, the vector cross product length gets divided by,
which should be the shortest, perpendicular length.

A few minor problems.
The ray direction vector already gets normalized at the beginning of
Sor_Intersect, so why are we calculating it's length again, and then checking if
it's more than zero, when at this point it should always be ONE.
Which bring up the potential divide-by-zero error when the normalization
calculation is done - and there is NO checking THERE!
And then there is the vector cross product sign issue.

So, after this gets rigorously checked, I'd suggest that
1. the ray direction vector length check gets moved to earlier in the code
2. we can do away with recalculating the length, checking for non-zero length,
and dividing the vector cross product length by a ray length which should
already always be 1.
3. Line 276 r0 = P[X] * D[Z] - P[Z] * D[X];
should be r0 = fabs(P[X] * D[Z] - P[Z] * D[X]);
so that the "length" is always positive.

- BE

Post a reply to this message