Converting a short piece of C/C++ code

[IdeasVacuum]

The definition of three float values in C/C++

Code: Select all

float    ax = (u.x >= 0 ? u.x : -u.x);
float    ay = (u.y >= 0 ? u.y : -u.y);
float    az = (u.z >= 0 ? u.z : -u.z);

I am uncertain how to convert to PB. I tried this, where u is a structured Vector (3D geometry):

Code: Select all

Structure vector3d
x.d               
y.d               
z.d               
EndStructure

Define u.vector3d
Define.d ax, ay, az

If(u\x >= 0) : ax = u\x : Else : ax = -(u\x) : EndIf
If(u\y >= 0) : ay = u\y : Else : ay = -(u\y) : EndIf
If(u\z >= 0) : az = u\z : Else : az = -(u\z) : EndIf

.... but it doesn't work

[wilbert]

Does this work ?

Code: Select all

Structure vector3d
x.d               
y.d               
z.d               
EndStructure

Define u.vector3d
Define.d ax, ay, az

ax = Abs(u\x)
ay = Abs(u\y)
az = Abs(u\z)

[Marc56us]

.... but it doesn't work

No, it's works, but result will be alway positive (-+-=+) so wilbert is right to abreviate this with Abs()

Code: Select all

;u\x = 1
;u\x = 0
u\x = -1

If (u\x >= 0) 
    Debug ">= 0: no change"
    ax = u\x 
Else 
    Debug "< 0 negative but - + - = + so +"
    ax = -(u\x) 
EndIf

[IdeasVacuum]

Hi Marc

Edit: Interestingly, Wilbert's snippet doesn't work but the result is similar....

However, I'm wondering if I did get it right (90%) but the correct result is not as expected. I know the result is partially wrong as it's the value that I can expect (except that, in my test case, it should be positive.

So I'm just going to try it with a tolerance #DBL_EPSILON

[TI-994A]

...doesn't work but the result is similar...

Code: Select all

#include <stdio.h>

struct vector3d { 
    float x;
    float y;
    float z;
};

int main() {
    struct vector3d u;
    
    u.x = 1;
    u.y = -1;
    u.z = 2;
    
    float  ax = (u.x >= 0 ? u.x : -u.x);
    float  ay = (u.y >= 0 ? u.y : -u.y);
    float  az = (u.z >= 0 ? u.z : -u.z);

    printf("%f \n", ax);   // 1.0
    printf("%f \n", ay);   // 1.0
    printf("%f \n", az);   // 2.0

    return 0;
}

Code: Select all

Structure vector3d
  x.d               
  y.d               
  z.d               
EndStructure

Define u.vector3d
Define.d ax, ay, az

u\x = 1
u\y = -1
u\z = 2

If (u\x >= 0) : ax = u\x : Else : ax = -(u\x) : EndIf
If (u\y >= 0) : ay = u\y : Else : ay = -(u\y) : EndIf
If (u\z >= 0) : az = u\z : Else : az = -(u\z) : EndIf

Debug ax   ; 1.0
Debug ay   ; 1.0
Debug az   ; 2.0

Code: Select all

Structure vector3d
  x.d               
  y.d               
  z.d               
EndStructure

Define u.vector3d
Define.d ax, ay, az

u\x = 1
u\y = -1
u\z = 2

ax = Abs(u\x)
ay = Abs(u\y)
az = Abs(u\z)

Debug ax   ; 1.0
Debug ay   ; 1.0
Debug az   ; 2.0

What's different, exactly?

[IdeasVacuum]

I can verify what the results should be in a CAD program.
The X value is always too small (by nearly 200) and signed incorrectly, the Y value is always the correct number but is always negative when it should always be positive (in my test cases). Also, the resulting line length is minuscule, which will cause problems further on......

The full C Code

Code: Select all

// intersect3D_2Planes(): find the 3D intersection of two planes
//    Input:  two planes Pn1 and Pn2
//    Output: *L = the intersection line (when it exists)
//    Return: 0 = disjoint (no intersection)
//            1 = the two  planes coincide
//            2 =  intersection in the unique line *L
int
intersect3D_2Planes( Plane Pn1, Plane Pn2, Line* L )
{
    Vector   u = Pn1.n * Pn2.n;          // cross product
    float    ax = (u.x >= 0 ? u.x : -u.x);
    float    ay = (u.y >= 0 ? u.y : -u.y);
    float    az = (u.z >= 0 ? u.z : -u.z);

    // test if the two planes are parallel
    if ((ax+ay+az) < SMALL_NUM) {        // Pn1 and Pn2 are near parallel
        // test if disjoint or coincide
        Vector   v = Pn2.V0 -  Pn1.V0;
        if (dot(Pn1.n, v) == 0)          // Pn2.V0 lies in Pn1
            return 1;                    // Pn1 and Pn2 coincide
        else
            return 0;                    // Pn1 and Pn2 are disjoint
    }

    // Pn1 and Pn2 intersect in a line
    // first determine max abs coordinate of cross product
    int      maxc;                       // max coordinate
    if (ax > ay) {
        if (ax > az)
             maxc =  1;
        else maxc = 3;
    }
    else {
        if (ay > az)
             maxc =  2;
        else maxc = 3;
    }

    // next, to get a point on the intersect line
    // zero the max coord, and solve for the other two
    Point    iP;                // intersect point
    float    d1, d2;            // the constants in the 2 plane equations
    d1 = -dot(Pn1.n, Pn1.V0);  // note: could be pre-stored  with plane
    d2 = -dot(Pn2.n, Pn2.V0);  // ditto

    switch (maxc) {             // select max coordinate
    case 1:                     // intersect with x=0
        iP.x = 0;
        iP.y = (d2*Pn1.n.z - d1*Pn2.n.z) /  u.x;
        iP.z = (d1*Pn2.n.y - d2*Pn1.n.y) /  u.x;
        break;
    case 2:                     // intersect with y=0
        iP.x = (d1*Pn2.n.z - d2*Pn1.n.z) /  u.y;
        iP.y = 0;
        iP.z = (d2*Pn1.n.x - d1*Pn2.n.x) /  u.y;
        break;
    case 3:                     // intersect with z=0
        iP.x = (d2*Pn1.n.y - d1*Pn2.n.y) /  u.z;
        iP.y = (d1*Pn2.n.x - d2*Pn1.n.x) /  u.z;
        iP.z = 0;
    }
    L->P0 = iP;
    L->P1 = iP + u;
    return 2;
}

As you can see, most of it will easily convert to PB. It is not necessary to test for parallel or coincident planes in my project, there simply aren't any.

My PB Version:

Code: Select all

EnableExplicit

#DBL_EPSILON = 0.0001

Structure pt3d
x.d
y.d
z.d
EndStructure

Structure ln3d
pt1.pt3d
pt2.pt3d
EndStructure

Structure vector3d
x.d
y.d
z.d
EndStructure

Structure plane3d
VecNorm.vector3d
PtV0.pt3d
EndStructure

Procedure PlanePlaneIntersection(*IntersectionLineReturn.ln3d, *PlaneA.plane3d, *PlaneB.plane3d)
;#----------------------------------------------------------------------------------------------
; find the 3D intersection of two planes, each represented by a Normal Vector and a point V0

Protected VecU.vector3d, VecV.vector3d
Protected dX.d, dY.d, dZ.d
Protected iMaxC.i     ;max coordinate
Protected Pt.pt3d     ;intersect point
Protected d1.d, d2.d  ;the constants in the 2 plane equations

               VecU\x = (*PlaneA\VecNorm\y * *PlaneB\VecNorm\z) - (*PlaneB\VecNorm\y * *PlaneA\VecNorm\z)
               VecU\y = (*PlaneA\VecNorm\z * *PlaneB\VecNorm\x) - (*PlaneB\VecNorm\z * *PlaneA\VecNorm\x)
               VecU\z = (*PlaneA\VecNorm\x * *PlaneB\VecNorm\y) - (*PlaneB\VecNorm\x * *PlaneA\VecNorm\y)

               If(VecU\x >= #DBL_EPSILON) : dX = VecU\x : Else : dX = -(VecU\x) : EndIf
               If(VecU\y >= #DBL_EPSILON) : dY = VecU\y : Else : dY = -(VecU\y) : EndIf
               If(VecU\z >= #DBL_EPSILON) : dZ = VecU\z : Else : dZ = -(VecU\z) : EndIf


               ;=== Test to see if the two planes are parallel or they coincide is not required - the planes are guaranteed not to be

               ;Determine max abs coordinate of cross product

               If (dX > dY)

                         If (dX > dZ)

                                   iMaxC = 1
                         Else
                                   iMaxC = 3
                         EndIf

               Else
                         If (dY > dZ)

                                   iMaxC = 2
                         Else
                                   iMaxC = 3
                         EndIf
               EndIf

               ;next, to get a point on the intersect line
               ;zero the max coord, and solve for the other two

               d1 = -(*PlaneA\VecNorm\x * *PlaneA\PtV0\x) + (*PlaneA\VecNorm\y * *PlaneA\PtV0\y) + (*PlaneA\VecNorm\z * *PlaneA\PtV0\z)
               d2 = -(*PlaneB\VecNorm\x * *PlaneB\PtV0\x) + (*PlaneB\VecNorm\y * *PlaneB\PtV0\y) + (*PlaneB\VecNorm\z * *PlaneB\PtV0\z)

               Select iMaxC        ;select max coordinate

                         Case 1    ;intersect with x=0

                                   Pt\x = 0
                                   Pt\y = (d2 * *PlaneA\VecNorm\z - d1 * *PlaneB\VecNorm\z) / VecU\x
                                   Pt\z = (d1 * *PlaneB\VecNorm\y - d2 * *PlaneA\VecNorm\y) / VecU\x

                         Case 2    ;intersect with y=0

                                   Pt\x = (d1 * *PlaneB\VecNorm\z - d2 * *PlaneA\VecNorm\z) / VecU\y
                                   Pt\y = 0
                                   Pt\z = (d2 * *PlaneA\VecNorm\x - d1 * *PlaneB\VecNorm\x) / VecU\y

                         Case 3    ;intersect with z=0

                                   Pt\x = (d2 * *PlaneA\VecNorm\y - d1 * *PlaneB\VecNorm\y) / VecU\z
                                   Pt\y = (d1 * *PlaneB\VecNorm\x - d2 * *PlaneA\VecNorm\x) / VecU\z
                                   Pt\z = 0
               EndSelect

               *IntersectionLineReturn\pt1\x = Pt\x
               *IntersectionLineReturn\pt1\y = Pt\y
               *IntersectionLineReturn\pt1\z = Pt\z

               *IntersectionLineReturn\pt2\x = Pt\x + VecU\x
               *IntersectionLineReturn\pt2\y = Pt\y + VecU\y
               *IntersectionLineReturn\pt2\z = Pt\z + VecU\z
EndProcedure

[/size]

[netmaestro]

I don't think it's any more complicated than this, if so someone correct me:

Code: Select all

; float    ax = (u.x >= 0 ? u.x : -u.x);
; float    ay = (u.y >= 0 ? u.y : -u.y);
; float    az = (u.z >= 0 ? u.z : -u.z);

Define ax.d, ay.d, az.d, ux.d, uy.d, uz.d

; test values
ux =  1.0
uy = -2.0
uz =  3.0

If Bool(ux > 0) : ax = ux : Else : ax = -ux : EndIf
If Bool(uy > 0) : ay = uy : Else : ay = -uy : EndIf
If Bool(uz > 0) : az = uz : Else : az = -uz : EndIf

Debug ax
Debug ay
Debug az

Absolute value, right? Removes the negative sign. Or just Abs() should work.

[alter Mann]

I would calculate the intersection line like this:

Code: Select all

EnableExplicit

#DBL_EPSILON = 0.0001

Structure pt3d
x.d
y.d
z.d
EndStructure

Structure ln3d
pt1.pt3d
pt2.pt3d
EndStructure

Structure vector3d
x.d
y.d
z.d
EndStructure

Structure plane3d
VecNorm.vector3d
PtV0.pt3d
EndStructure

Procedure PlanePlaneIntersection(*IntersectionLineReturn.ln3d, *PlaneA.plane3d, *PlaneB.plane3d)
;#----------------------------------------------------------------------------------------------
; find the 3D intersection of two planes, each represented by a Normal Vector and a point V0

Protected VecU.vector3d, VecV.vector3d, VecH.vector3d
Protected dA.d, dH.d

               VecU\x = (*PlaneA\VecNorm\y * *PlaneB\VecNorm\z) - (*PlaneB\VecNorm\y * *PlaneA\VecNorm\z) ; direction of intersection line
               VecU\y = (*PlaneA\VecNorm\z * *PlaneB\VecNorm\x) - (*PlaneB\VecNorm\z * *PlaneA\VecNorm\x)
               VecU\z = (*PlaneA\VecNorm\x * *PlaneB\VecNorm\y) - (*PlaneB\VecNorm\x * *PlaneA\VecNorm\y)
               
               VecV\x = (*PlaneA\VecNorm\y * VecU\z) - (VecU\y * *PlaneA\VecNorm\z) ; normal direction to plane normal and direction of intersection line
               VecV\y = (*PlaneA\VecNorm\z * VecU\x) - (VecU\z * *PlaneA\VecNorm\x)
               VecV\z = (*PlaneA\VecNorm\x * VecU\y) - (VecU\x * *PlaneA\VecNorm\y)

               VecH\x = *PlaneB\PtV0\y - *PlaneA\PtV0\z                                                             ;
               VecH\y = *PlaneB\PtV0\z - *PlaneA\PtV0\x                                                             ;
               VecH\z = *PlaneB\PtV0\x - *PlaneA\PtV0\y                                                             ;
                                                                                                                    ;
               dH = *PlaneB\VecNorm\x * VecV\x + *PlaneB\VecNorm\y * VecV\y + *PlaneB\VecNorm\z * VecV\z            ;
               dA = (*PlaneB\VecNorm\x * VecH\x + *PlaneB\VecNorm\y * VecH\y + *PlaneB\VecNorm\z * VecH\z) / dH     ;
                                                                                                                    ;
               *IntersectionLineReturn\pt1\x = *PlaneA\PtV0\x + VecV\x * dA                                         ;   intersection point between       
               *IntersectionLineReturn\pt1\y = *PlaneA\PtV0\y + VecV\y * dA                                         ;   line (from point of PlaneA and VecV)
               *IntersectionLineReturn\pt1\z = *PlaneA\PtV0\z + VecV\z * dA                                         ;   and PlaneB                       

               *IntersectionLineReturn\pt2\x = *IntersectionLineReturn\pt1\x + VecU\x
               *IntersectionLineReturn\pt2\y = *IntersectionLineReturn\pt1\y + VecU\y
               *IntersectionLineReturn\pt2\z = *IntersectionLineReturn\pt1\z + VecU\z
EndProcedure

[IdeasVacuum]

netmaestro!

I don't think it's to ensure an absolute value, given that C\CPP have an abs function, I read it as "if the value of u is zero or more, a is a positive value, else a is a negative value.

The snag is, unless I have misunderstood something else further down in the code, this snippet seems to matter a lot to the final value, but having tested more I find I must have that "something else" wrong too.

I'm going to make a sample set of known values.

[wilbert]

@IdeasVacuum, you translate

Code: Select all

    d1 = -dot(Pn1.n, Pn1.V0);  // note: could be pre-stored  with plane
    d2 = -dot(Pn2.n, Pn2.V0);  // ditto

into

Code: Select all

               d1 = -(*PlaneA\VecNorm\x * *PlaneA\PtV0\x) + (*PlaneA\VecNorm\y * *PlaneA\PtV0\y) + (*PlaneA\VecNorm\z * *PlaneA\PtV0\z)
               d2 = -(*PlaneB\VecNorm\x * *PlaneB\PtV0\x) + (*PlaneB\VecNorm\y * *PlaneB\PtV0\y) + (*PlaneB\VecNorm\z * *PlaneB\PtV0\z)

Shouldn't it be

Code: Select all

               d1 = -((*PlaneA\VecNorm\x * *PlaneA\PtV0\x) + (*PlaneA\VecNorm\y * *PlaneA\PtV0\y) + (*PlaneA\VecNorm\z * *PlaneA\PtV0\z))
               d2 = -((*PlaneB\VecNorm\x * *PlaneB\PtV0\x) + (*PlaneB\VecNorm\y * *PlaneB\PtV0\y) + (*PlaneB\VecNorm\z * *PlaneB\PtV0\z))

[IdeasVacuum]

Hi alter Mann, I'm going to give that a go

[IdeasVacuum]

Aha. The results with alter Mann's code are correct in terms of values. Mostly the wrong sign shows up but I can change that in-process no problem.

Much thanks to alter Mann for the enhanced code and thanks to all for leaping to my assistance!

I'm pretty sure this is going to lead into another problem a bit further down the line

[IdeasVacuum]

Shouldn't it be

You are right wilbert. Actually I have found using the "impose sign" method failed when applied to an equation and now I see why.

[IdeasVacuum]

Well, I'm actually not really done and dusted. This is a drawing to illustrate what I'm doing:



So I'm creating intersection lines between planes (just planes, not faces, they are not graphically represented in the project).

The blue square represents the plane used each time by the procedure to create an intersection line with one of the other planes. In most cases, the planes are not orthogonal, i.e. the green planes are not perpendicular to the blue plane and the blue plane is at a slight angle too. So, nearly square to each other but quite.

The goal is actually the intersection points at the intersections of the intersection lines (again, no graphical representation in the project, they are a means to an end).

Example set of planes:
Blue Plane VecNorm\x = 0.0055
Blue Plane VecNorm\y = -0.9998
Blue Plane VecNorm\z = -0.0166

Blue Plane PtV0\x = 20.8553
Blue Plane PtV0\y = 797.4753
Blue Plane PtV0\z = -584.7176

Plane Aft VecNorm\x = 0.0007
Plane Aft VecNorm\y = -0.0632
Plane Aft VecNorm\z = -0.9980

Plane Aft PtV0\x = 120.6621
Plane Aft PtV0\y = 805.9420
Plane Aft PtV0\z = -692.0133

Plane LHS VecNorm\x = -0.0729
Plane LHS VecNorm\y = 0.9967
Plane LHS VecNorm\z = 0.0360

Plane LHS PtV0\x = -451.9908
Plane LHS PtV0\y = 799.1656
Plane LHS PtV0\z = -541.5027

Plane RHS VecNorm\x = 0.9997
Plane RHS VecNorm\y = 0.0240
Plane RHS VecNorm\z = 0.0025

Plane RHS PtV0\x = 447.0984
Plane RHS PtV0\y = 800.9244
Plane RHS PtV0\z = -525.6143

Plane FWD VecNorm\x = -0.0011
Plane FWD VecNorm\y = 0.9978
Plane FWD VecNorm\z = 0.0659

Plane FWD PtV0\x = 134.1713
Plane FWD PtV0\y = 786.6722
Plane FWD PtV0\z = -478.3818

Firstly, the Lines are not right, not planar - they should of course all be on the Blue Plane. Second problem is, when the lines are right, I need to find the Line-Line intersections, which are virtual since they do not "physically" cross each other.
Anyone interested in a virtual beer?

[IdeasVacuum]

For the line intersections I found this by Paul Bourke but I don't know - seems this method needs the lines to physically overlap.

Code: Select all

Procedure.i LineIntersect3d(*ptIntersect.pt3d, *ln1.ln3d, *ln2.ln3d)
;#--------------------------------------------------------------------
;3d intersection of two lines
;http://paulbourke.net/geometry/pointlineplane/lineline.c
;XYZ p1,XYZ p2,XYZ p3,XYZ p4,XYZ *pa,XYZ *pb, double *mua, double *mub)
;p1,p2 = ln1; p3,p4 = ln2

Protected p13.pt3d, p43.pt3d, p21.pt3d ;XYZ p13,p43,p21
Protected d1343.d, d4321.d, d1321.d, d4343.d, d2121.d
Protected dNum.d, dNom.d, dMua.d

               p13\x = (*ln1\pt1\x - *ln2\pt1\x)
               p13\y = (*ln1\pt1\y - *ln2\pt1\y)
               p13\z = (*ln1\pt1\z - *ln2\pt1\z)

               p43\x = (*ln2\pt2\x - *ln2\pt1\x)
               p43\y = (*ln2\pt2\y - *ln2\pt1\y)
               p43\z = (*ln2\pt2\z - *ln2\pt1\z)

               If( (Abs(p43\x) < #DBL_EPSILON) And (Abs(p43\y) < #DBL_EPSILON) And (Abs(p43\z) < #DBL_EPSILON) ) : ProcedureReturn(#False) : EndIf

               p21\x = (*ln1\pt2\x - *ln1\pt1\x)
               p21\y = (*ln1\pt2\y - *ln1\pt1\y)
               p21\z = (*ln1\pt2\z - *ln1\pt1\z)

               If( (Abs(p21\x) < #DBL_EPSILON) And (Abs(p21\y) < #DBL_EPSILON) And (Abs(p21\z) < #DBL_EPSILON) ) : ProcedureReturn(#False) : EndIf

               d1343 = p13\x * p43\x + p13\y * p43\y + p13\z * p43\z
               d4321 = p43\x * p21\x + p43\y * p21\y + p43\z * p21\z
               d1321 = p13\x * p21\x + p13\y * p21\y + p13\z * p21\z
               d4343 = p43\x * p43\x + p43\y * p43\y + p43\z * p43\z
               d2121 = p21\x * p21\x + p21\y * p21\y + p21\z * p21\z

               dNom = (d2121 * d4343 - d4321 * d4321)

               If( (Abs(dNom) < #DBL_EPSILON) ) : ProcedureReturn(#False) : EndIf

               dNum = (d1343 * d4321 - d1321 * d4343)

               dMua = (dNum / dNom)

               *ptIntersect\x = *ln1\pt1\x + dMua * p21\x
               *ptIntersect\y = *ln1\pt1\y + dMua * p21\y
               *ptIntersect\z = *ln1\pt1\z + dMua * p21\z

               ProcedureReturn(#True)
EndProcedure

[/size]