On 9/10/25 09:27, Bald Eagle wrote:
> I'll make a speculative prediction and say that the algorithm will try to take
> the sqrt of that negative radius and fail - either throwing an error, or just
> blanking out that portion of the curve - or the whole segment.

Thanks for reference! :-)

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.

I see some other not really wrong, but not optimal stuff happening too 
with the cylinder segment bounding in y. Sort of as if it was written to 
support on curve points in ascending or descending order less any 
optimization for the y coordinates (The cylinders get taller and taller 
as you have more segments instead of aligning with each of the segments 
in y).

Bill P.

Post a reply to this message

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