## Monday, October 15, 2007

### HP-35s Euclid Implementation - XIII : Cartesian to Barycentric coordinates

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:
1. 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.
2. 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))/Qb2=((y-y1)*(x3-x1)+(y3-y1)*(x1-x))/Qb3=((y-y2)*(x1-x2)+(y1-y2)*(x2-x))/Q`
Where (xi, yi) 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*x3y=b1*y1+b2*y2+b3*y31=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:
1. 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.
2. 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 TT002 XEQ U001T003 FS? 2T004 XEQ X068T005 STO MT006 RDNT007 FS? 2T008 XEQ X068T009 STO LT010 RDNT011 FS? 2T012 XEQ X068T013 STO KT014 R^T015 X<>YT016 -T017 RDNT018 RDNT019 LASTXT020 -T021 X<>YT022 RDNT023 XEQ X001T024 ABST025 1E-6T026 X>Y?T027 GTO T065T028 LASTXT029 XEQ X025T030 STO XT031 RCL KT032 X<>YT033 XEQ X052T034 STO UT035 RDNT036 RCL LT037 X<>YT038 XEQ X052T039 STO VT040 RDNT041 RCL MT042 X<>YT043 XEQ X052T044 STO WT045 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 QT047 LASTXT048 4T049 STO IT050 RDNT051 RCL(I)T052 ABSTO53 RDNT054 FS? 2T055 GTO T060T056 RCL KT057 RCL LT058 RCL MT059 GTO T063T060 RCL UT061 RCL VT062 RCL WT063 R^T064 RTNT065 SF 10T066 eq DEGENERATEDT067 PSET068 CF 10T069 XEQ U029T070 RTN`

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:

`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 BB002 XEQ U001B003 FS? 2B004 GTO B007B005 RCL XB006 XEQ X052B007 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)]/QB008 XEQ U050B009 RTN`