Arrived at problem 14 of the
requirements: Computing barycentric coordinates given cartesian coordinates.
Barycentric coordinates [b1, b2, b3] of a point in a plane are coordinates defined by a triangle in the plane (defines the plane) where the triangle's first vertex's barycentric coordinates are [1,0,0], second vertex's are [0,1,0], third vertex's are [0,0,1] and one have that b1+b2+b3=1.
Coordinates of any point in the plane can be expressed as the weighted average of the triangle's vertices coordinates, these weights are the barycentric coordinates:
Barycentric Cartesian
(b1, b2, b3)<=>b1v1+b2v2+b3v3
Where
vi is the positional vector for vertex
i.
Barycentric coordinates has several uses and are of theoretical importance outside the scope of this blog entry. Here I give two practical uses:
- Provides a robust test for if a point in triangle's plane with barycentric coordinates [b1, b2, b3] is inside the triangle, on one of the triangle's edge or outside the triangle: If all components are in the range[0,1] point is inside, if bi is 0 and the two others are in the [0,1] range then the point is on the edge oposite of triangle's vertice i, and finaly if at least one component is not in the [0,1] range then point is outside triangle.
- If a function's values is known at the triangle's vertices one can perform linear interpolation inside the triangle at the point: v1*b1+v2*b2+v3*b3. Where vi is the value at vertex i and [b1, b2, b3] is the barycentric coordinates for point to evaluate at. Note that converting from barycentric coordinates to cartesian coordinates is a special case of this.
I choose to convert from cartesian coordinates to barycentric by programming the solution to the equations for the 2D cartesian to barycentric coordinates case:
Q=(y1-y3)*(x2-x3)+(y2-y3)*(x2-x1)
b1=((y-y3)*(x2-x3)+(y2-y3)*(x3-x))/Q
b2=((y-y1)*(x3-x1)+(y3-y1)*(x1-x))/Q
b3=((y-y2)*(x1-x2)+(y1-y2)*(x2-x))/Q
Where (x
i, y
i) are cartesian coordinates for triangle's
i'th vertex and (x, y) is the cartesian coordinates of the point in question. This is the solution to the equation set of three uknown (b1, b2, b3):
x=b1*x1+b2*x2+b3*x3
y=b1*y1+b2*y2+b3*y3
1=b1+b2+b3
To solve 3D cartesian to barycentric I project the triangle and point in question to 2D and solve as a 2D problem.
Implementation is done in two steps:
- Modify the T program so it: A) Also accept a 2D triangle (flag 2 set='2D mode'. B) Project the triangle to the cardinal plane that gives the largest area, this is the triangle used for computing barycentric coordinates. C) Computes the quotient Q in the above equations.
- Implement a program labeled B to accepts a 2D or 3D cartesian tuple and use result from B and C above to compute barycentric coordinates.
Step 1 : The modified T program
Stack Input/Output:2D case: [[x2, y2], [x1, y1], [x0, y0] | L]->XEQ T->[[nx, ny, nz], [x2, y2], [x1, y1], [x0, y0] | L]
3D case: [[x2, y2, z2], [x1, y1, z1], [x0, y0, z0] | L]->XEQ T->[[nx, ny, nz], [x2, y2, z2], [x1, y1, z1], [x0, y0, z0] | L]
Where v0=(x0, y0, z0) is first, v1=(x1, y1, z1) is second, v2=(x2, y2, z2) is the third triangle vertex and [nx, ny, nz] is the cross product (v1-v0)x(v2-v0).
Variables:Writes:
K : First vertex (3D vector).
L : Second vertex (3D vector).
M : Third vertex (3D vector).
Q : Quotient (y1-y3)*(x2-x3)+(y2-y3)*(x2-x1).
U : First vertex of triangle projected to 2D (2D vector).
V : Second vertex of triangle projected to 2D (2D vector).
W : Third vertex of triangle projected to 2D (2D vector).
X : Index of coordinate discarded in 3D->2D projection.
Program:T001 LBL T
T002 XEQ U001
T003 FS? 2
T004 XEQ X068
T005 STO M
T006 RDN
T007 FS? 2
T008 XEQ X068
T009 STO L
T010 RDN
T011 FS? 2
T012 XEQ X068
T013 STO K
T014 R^
T015 X<>Y
T016 -
T017 RDN
T018 RDN
T019 LASTX
T020 -
T021 X<>Y
T022 RDN
T023 XEQ X001
T024 ABS
T025 1E-6
T026 X>Y?
T027 GTO T065
T028 LASTX
T029 XEQ X025
T030 STO X
T031 RCL K
T032 X<>Y
T033 XEQ X052
T034 STO U
T035 RDN
T036 RCL L
T037 X<>Y
T038 XEQ X052
T039 STO V
T040 RDN
T041 RCL M
T042 X<>Y
T043 XEQ X052
T044 STO W
T045 eq ([0,1]xU-[0,1]xW)x([1,0]xV-[1,0]xW)+
([0,1]xV-[0,1]xW)x([1,0]xW-[1,0]xU)
T046 STO Q
T047 LASTX
T048 4
T049 STO I
T050 RDN
T051 RCL(I)
T052 ABS
TO53 RDN
T054 FS? 2
T055 GTO T060
T056 RCL K
T057 RCL L
T058 RCL M
T059 GTO T063
T060 RCL U
T061 RCL V
T062 RCL W
T063 R^
T064 RTN
T065 SF 10
T066 eq DEGENERATED
T067 PSE
T068 CF 10
T069 XEQ U029
T070 RTN
Comments:Terms of use.Mnemonics: T for Triangle.
Uses the
X068 subroutine to place 2D triangle in th XZ plane if flag 2 set and the
X025 routine to find out which coordinate to discard for the 3D to 2D projection and finaly the
X052 routine to actual do the 3D to 2D conversion. In the 2D case this is a rather roundabout way to get end user's input into the variables U, V and W! But I think this worked out quite well in order to have one program for both the 2D and 3D case:
One can compute barycentric coordinates in 3D using four cross product, four dot product and three scalar division operations avoiding projection from 3D to 2D, but involves more scalar math operations than the implemented method. On a vector architecture it may be a better choice, on the HP-35s architecture I think this is a good solution!
The index of the coordinate discarded to do the 3D to 2D projection is remembered in variabe X so the next program to actual compute barycentric coordinates know which coordinate to discard from input.
An unrelated change done to the T program is that it now check if the computed plane normal vector is considered to be of zero length: Then all points are the same or colinear, the message "DEGENERATED" is then shortly displayed and stack registers restored.
Note that one get a cross product calculated normal vector also in the 2D case (direction
[0, 1, 0]), length still twice triangle's area.
Step 2 : Program to convert from Cartesian to Barycentric Coordinates
Stack Input/Output:[[x, y, z], Y, Z, T | L]->XEQ B->[[b1, b2, b3], Y, Z, T | [x, y, z]]
Where (x, y, z) are cartesian coordinates for point in question and (b1, b2, b3) are computed barycentric coordinates.
Variables:Reads:
Q : Quotient (y1-y3)*(x2-x3)+(y2-y3)*(x2-x1). Populated by the P program.
U : First vertex of triangle projected to 2D (2D vector). Populated by the P program.
V : Second vertex of triangle projected to 2D (2D vector). Populated by the P program.
W : Third vertex of triangle projected to 2D (2D vector). Populated by the P program.
X : Index of coordinate discarded in 3D->2D projection. Populated by the P program.
Program:B001 LBL B
B002 XEQ U001
B003 FS? 2
B004 GTO B007
B005 RCL X
B006 XEQ X052
B007 eq [([0,1]xREGX-[0,1]xW)x([1,0]xV-[1,0]xW)+
([0,1]xV-[0,1]xW)x([1,0]xW-[1,0]xREGX),
([0,1]xREGX-[0,1]xU)x([1,0]xW-[1,0]xU)+
([0,1]xW-[0,1]xU)x([1,0]xU-[1,0]xREGX),
([0,1]xREGX-[0,1]xV)x([1,0]xU-[1,0]xV)+
([0,1]xU-[0,1]xV)x([1,0]xV-[1,0]xREGX)]/Q
B008 XEQ U050
B009 RTN
Comments:Terms of use.Mnemonics: B for Barycentric.
If the input is a 2D (flag 2 set) then goes straight to the equation that calculates baycentric coordinates using the data precomputed by the modified T program, else first project the 3D input to 2D using
subroutine X052 and information found in the X variable.
