# Thread: [RESOLVED] VB6 and some 3D Linear Algebra

1. ## [RESOLVED] VB6 and some 3D Linear Algebra

Hi Guys,

I've got a linear algebra issue I've been trying to sort out for a while. And, my problem is, I've got linear algebra resources and I've got VB6 resources (i.e., you guys), but I don't have any resources that combine the two. Therefore, you guys are drafted. Also, I bet a couple of you have done some linear algebra at some point in the past.

Basically, I'm trying to build three functions in VB6.

1. Build a 3D segment.
2. Rotate a 3D segment.
3. Extract angles between two 3D segments.

Okay, I need to lay some groundwork. First, what's a 3D segment (not to be confused with a line segment)? Here's a little thing I carry around with me to think about this stuff:

That thing represents a 3D segment. If you just look at the X and Y axes, you can imagine a piece of graph paper with an origin (0,0) in the middle of the paper. The corner of my little device might be at (0,0). However, I've got to think in 3D, so I need that Z axis that sticks straight up. Therefore, the corner of all axes of the device (see picture) can be thought of as (0,0,0).

So, basically, a 3D segment is represented by four points: The origin (i.e., the corner), the X_Axis, the Y_Axis, and the Z_Axis. And, for simplicity, we will assume that each of the axes is exactly one unit long.

Now, let's assume we have a room with a tile floor, and we've designated an origin (0,0,0) in the middle on the floor somewhere. Let's imagine a dot there on the floor in the middle of the room. And also, we've imagined some direction as X+, and another orthogonal (perpendicular) direction as Y+. It's as if we laid a piece of graph paper on the floor of our room.

And obviously (I hope), Z+ will be straight up.

So, if we lay our little wooden 3D segment device on the room's origin (0,0,0) with its X_Axis aligned with X and its Y_Axis aligned with Y, our 3D segment would be defined as:

Origin = (0,0,0)
X_Axis = (1,0,0)
Y_Axis = (0,1,0)
Z_Axis = (0,0,1)

You still with me?

Now, we need to back up just a bit. In addition to these 3D segments, we can define points in space anywhere in our room. For convenience, we will define them as 3D vectors that originate at our origin spot in the room (just like graph paper, but 3D). Because the origin is on the floor, Z will always be positive, but X and Y can be either positive or negative (because it's most convenient to place the origin near the middle of the room).

Therefore, using this system, we can define any point in space as a 3D vector from the room's origin. We'll use a UDT to define this:

Code:
```
Public Type VectorType
x As Double
y As Double
z As Double
End Type
```
Now, let's return to our 3D segment. The 3D segment discussed above (placed at the origin) is really just four of these vectors: One at the origin, one that's one unit along the X_Axis, one that's one unit along the Y_Axis, and one that's one unit along the Z_Axis. So, to define any 3D segment, we can use another UDT:

Code:
```
Public Type SegmentType
Origin As VectorType
AxisX  As VectorType
AxisY  As VectorType
AxisZ  As VectorType
End Type
```
Here's where thing start to get interesting. I've talked about our 3D segment being placed on (and aligned with) the origin, as such:

But that doesn't necessarily have to be the case. In other words, our 3D segment can be both translated and rotated.

Translation is simply the case where our 3D segment's origin isn't positioned on the room's origin. And I'm not really concerned with translation. In the broader scheme of things, translations are easy. They're just vector addition and subtraction of origins.

Rotation is the case where a 3D segment's X_Axis (and/or Y_Axis and/or Z_Axis) is no longer aligned with (or even parallel to) the room's X_Axis (or Y_Axis or Z_Axis, respectively). Here's my attempt to take a picture of my rotated (and translated) little 3D segment device:

It's off the floor (translation), but more importantly, its axes are no longer aligned with the room. However, let's go back to our SegmentType UDT. We can certainly imagine some vector (VectorType) that represents our new 3D segment's origin. Furthermore, we can imagine some vector that represents the 3D points on the ends of each of our little 3D segment's axes. They will no longer be the simple ...
Origin = (0,0,0)
X_Axis = (1,0,0)
Y_Axis = (0,1,0)
Z_Axis = (0,0,1)
... that we started with. But we could imagine measuring them in some way and figuring out all four of those points. And they would "mathematically" represent our rotated and translated 3D segment in space.

I think I'll stop this post here. I haven't gotten to the three functions I'd like, but I'm getting there.

Please feel free to make comments if you'd like. I'll be posting more soon.

Best Regards,
Elroy

Continuing introduction/explanation of this is here.

2. ## Re: VB6 and some 3D Linear Algebra

Segment has start and end co-ordinates in 3D, fex. line segment.

StartPoint.X
StartPoint.Y
StartPoint.Z
EndPoint.X
EndPoint.Y
EndPoint.Z

So translating aircode to VB, one could want to write.
Code:
```Type 3DPoint
X As Double
Y As Double
Z As Double
End Type

Type Face
Vertex(2) As 3DPoint
End Type

Type 3DLine
P1 As 3DPoint
P2 As 3DPoint
End Type

'other types

'Rotating vertices in 3D, using basic sini and cosini functions.
'https://en.wikipedia.org/wiki/Rotation_matrix

Function RotateFace(MyFace As 3DFace, Pivot As 3DPoint, Angle As Double) As 3DFace
RotateFace = MyFace
RotateFace.Vertex(0) = RotatePoint(MyFace.Vertex(0), Pivot, Angle)
RotateFace.Vertex(1) = RotatePoint(MyFace.Vertex(1), Pivot, Angle)
RotateFace.Vertex(2) = RotatePoint(MyFace.Vertex(2), Pivot, Angle)
End Function

Function RotatePoint(MyPoint As 3DPoint, Pivot As 3DPoint, Angle As Double) As 3DPoint
RotatePoint = MyPoint
RotatePoint.X = RotX(MyPoint.X - Pivot.X, MyPoint.Y - Pivot.Y, Angle) + Pivot.X
RotatePoint.Y = RotY(MyPoint.X - Pivot.X, MyPoint.Y - Pivot.Y, Angle) + Pivot.Y
End Function

Function RotX(X1 As Double, Y1 As Double, Angle As Double) As Double
RotX = cHyp(X1, Y1) * Cos((PtAng(X1, Y1) + Angle) * Pi / 180)
End Function

Function RotY(X1 As Double, Y1 As Double, Angle As Double) As Double
RotY = cHyp(X1, Y1) * Sin((PtAng(X1, Y1) + Angle) * Pi / 180)
End Function

Function cHyp(X1 As Double, Y1 As Double) As Double
cHyp = Sqr((X1 * X1) + (Y1 * Y1))
End Function

Function PtAng(X1 As Double, Y1 As Double) As Double
If X1 = 0 Then
If Y1 >= 0 Then
PtAng = 90
Else
PtAng = 270
End If
Exit Function
ElseIf Y1 = 0 Then
If X1 >= 0 Then
PtAng = 0
Else
PtAng = 180
End If
Exit Function
Else
PtAng = Atn(Y1 / X1)
PtAng = PtAng * 180 / Pi
If PtAng < 0 Then PtAng = PtAng + 360
If PtAng > 360 Then PtAng = PtAng - 360
'----------Test for direction-(quadrant check)-------
If X1 < 0 Then PtAng = PtAng + 180
If Y1 < 0 And PtAng < 90 Then PtAng = PtAng + 180
If PtAng < 0 Then PtAng = PtAng + 360
If PtAng > 360 Then PtAng = PtAng - 360
End If
End Function

Function cAngle(MyLine As 3DLine) As Double
If MyLine.P1.X = MyLine.P2.X Then
If MyLine.P1.Y < MyLine.P2.Y Then
cAngle = 90
Else
cAngle = 270
End If
Exit Function
ElseIf MyLine.P1.Y = MyLine.P2.Y Then
If MyLine.P1.X < MyLine.P2.X Then
cAngle = 0
Else
cAngle = 180
End If
Exit Function
Else
cAngle = Atn(CSlope(MyLine))
cAngle = cAngle * 180 / Pi
If cAngle < 0 Then cAngle = cAngle + 360
'----------Test for direction--------
If MyLine.P1.X > MyLine.P2.X And cAngle <> 180 Then cAngle = cAngle + 180
If MyLine.P1.Y > MyLine.P2.Y And cAngle = 90 Then cAngle = cAngle + 180
If cAngle > 360 Then cAngle = cAngle - 360
End If
End Function

Function CSlope(MyLine As 3DLine) As Double
CSlope = (MyLine.P2.Y - MyLine.P1.Y) / (MyLine.P2.X - MyLine.P1.X)
End Function

'Angle between lines.
Function PtPtAngle(P1 As 3DPoint, P2 As 3DPoint) As Double
Dim MyLine As 3DLine
MyLine.P1 = P1
MyLine.P2 = P2
PtPtAngle = cAngle(MyLine)
End Function```

3. ## Re: VB6 and some 3D Linear Algebra

Hi Tech99,

Thanks for playing along. I was worried that I was going to be quite lonely in this thread. It's good to see that at least one other person is interested.

And yes, your conceptualizations are correct. I'm actually using the language of the mocap (motion capture) community, which is a bit different from what you'd hear in a math class (but not all that different).

I think I'll go back and edit my post #1, and always say 3D segment, rather than just segment. I promise, this is what they're called in the mocap community.

And, it's easy to get these 3D segments confused with a line segment, which we'll also be talking about shortly (and which you already introduced).

Also, for those interested, there are basically two "big boys" that sell equipment for the mocap community: Vicon and Motion Analysis Corp. Those two are basically the Intel and AMD of mocap. There are others, but they're way down the list as "also rans". And there is a thriving community that uses this mocap equipment for a variety of purposes. For instance, it's used to make movies (Polar Express, Lord of the Rings, Avatar) and it's used to analyze things such as golf swings or batting swings. However, I use it to better understand how handicapped children walk, so that various therapies may be used to help them walk better. We have databases of "typically developing" children for comparison.

I've got to attend to some other things today, but I'll make another post that completes the definition of a 3D segment (at least the way they're defined in the mocap community).

Again, thanks for your interest.

Best Regards,
Elroy

4. ## Re: VB6 and some 3D Linear Algebra

Okay, to continue with the definition of a 3D segment.

As stated above, we can define points in space with the following UDT:

Code:
```
Public Type VectorType
x As Double
y As Double
z As Double
End Type
```
And then, we can use this to define the UDT for one of our 3D segments, as follows:

Code:
```
Public Type SegmentType
Origin As VectorType
AxisX  As VectorType
AxisY  As VectorType
AxisZ  As VectorType
End Type
```
Now, if our 3D segment is positioned at the origin of the room, things are quite simple. However, when our 3D segment is rotated and/or translated, things aren't quite so simple.

Now, we can imagine our room having an origin, but our 3D segment also has its own origin (the corner_intersection of its three axes). However, this 3D segment's origin isn't typically of importance. It's more often the rotation of our 3D segment that's of more importance. Therefore, when these 3D segments are defined, any translation typically is not included in the definition of the axes:

Code:
```
Public Type SegmentType
Origin As VectorType
AxisX  As VectorType ' <--- Does NOT include any translation offset.
AxisY  As VectorType ' <--- Does NOT include any translation offset.
AxisZ  As VectorType ' <--- Does NOT include any translation offset.
End Type
```
In other words, the axes definitions assume the 3D segment's origin is at the room's origin, even if/when it's not. This simplifies a great many things. For instance, remembering that our axes definitions always define unit axes, we can say:

AxisX.x^2 + AxisX.y^2 + AxisX.z^2 = 1
AxisY.x^2 + AxisY.y^2 + AxisY.z^2 = 1
AxisZ.x^2 + AxisZ.y^2 + AxisZ.z^2 = 1

And this will always be true even if our 3D segment is rotated. Basically, this is just a bit of application of the Pythagorean theorem.

Let me also introduce at this point what simple vector addition and subtraction functions look like:

Code:
```
Public Function AddVectors(v1 As VectorType, v2 As VectorType) As VectorType
AddVectors.x = v1.x + v2.x
AddVectors.y = v1.y + v2.y
AddVectors.z = v1.z + v2.z
End Function

Public Function SubtractVectors(v1 As VectorType, v2 As VectorType) As VectorType
SubtractVectors.x = v1.x - v2.x
SubtractVectors.y = v1.y - v2.y
SubtractVectors.z = v1.z - v2.z
End Function
```
Hopefully, with minimal study, these make sense. If we wanted to include any translation offset into our axes definitions (or possibly remove it), these functions could get this done for us.

Ok, got to work for a bit, but I will continue.

Again, please feel free to jump in. I'd love the company.

Best,
Elroy

5. ## Re: VB6 and some 3D Linear Algebra

Elroy

Sounds interesting.
To me, it begs the question .. what might an example data-set consist of?
• a series of 3-D points with lines drawn to connect them
• a video
• something else

Spoo

6. ## Re: VB6 and some 3D Linear Algebra

Hi Spoo,

Well, this stuff can get complicated really quickly, and I'm doing my best to try and keep it simple.

For instance, in a typical situation, we capture points-in-space dynamically. What "dynamically" means is "through time" at a specified framerate. However, to keep these complexities out of it, we can just imagine static (stationary) points-in-space. (I'll save the "dynamic" aspects for another time.) Basically, it's the difference between a picture and a movie. We'll just deal with pictures.

Just to give you an idea of these points-in-space, here's a link to a mocap picture I previously posted. Those are infra-red (IR) markers that can be "captured" as points-in-space.

With any three (or four) of those IR markers, I can build one of the 3D segments I've been discussing.

Okay, you ask about data files. For this thread, I'll probably just stick with ASCII (Notepad readable) files. However, in production, we get files with a C3D extension. That's an industry standard file type. C3D stands for Coordinate-Three-Dimensional. It's basically a file that has 3D vectors (see the UDT VectorType above) in it. These C3D files can also have angles (i.e. Eulers) in them, but I jump ahead with that. We'll work our way to those angles.

So, without getting too much into the technicalities, we might imagine a "picture" of that guy in that link, with that picture represented by a series of 3D points, with each point having an X, a Y, and a Z measurement to specify precisely where the 3D point is in 3D space.

But the next order of business is to construct those 3D segments, and I'll talk about how we take 3D points and construct 3D segments shortly.

All The Best,
Elroy

EDIT1: Spoo, maybe to give you a visual, here's a screenshot of my C3D reader. You see Point Labels on the left. Those correspond to those IR markers on the picture (linked to above). Once you select a point, you will see the dynamic data for that point. However, for this thread, I'll stay with a static conceptualization. In other words, we'll just deal with the first row of data for each point (i.e., a picture rather than a movie).

And the three columns of data are basically the X-coordinate, the Y-coordinate, and the Z-coordinate for the point being displayed (i.e., the VectorType UDT again), from the room's origin (as discussed in previous posts).

7. ## Re: VB6 and some 3D Linear Algebra

I don't quite understand, why you need to define origin in element types? Origin is kind of 'constantly changing', when rotating objects in 3D space. Vertex point co'ordinates has all info needed, to compute it's whereabouts in space.

8. ## Re: VB6 and some 3D Linear Algebra

Hi Tech99,

It's not me that's defining this stuff. To a large degree, I'm just trying to emulate what some much higher level scripts do in the mocap world. In other words, it's not me that made up defining a 3D segment in the way I've outlined it above.

But I've got to follow their standards, or I'll lose my audience in the mocap community.

Basically, they capture points-in-space, but they rather quickly construct these 3D segments, with all subsequent discussions hinging around these 3D segments. For example, during walking, we typically define the pelvis as the center of the body. We "construct" one of these 3D segments that represents the pelvis moving through space. And then, we construct another of these 3D segments that represents the femur (the bone from the pelvis to the knee). Basically, we can imagine this femur 3D segment attached to the distal (knee) end of the femur.

Then, to proceed down, we can imagine another of these 3D segments attached to the distal (ankle) end of the tibia (bone between knee and ankle).

Carrying this on out, you can construct the entire skeleton. If you're holding a golf club (or a bat), you could even construct that.

Again, this is just the way they do it in the mocap community. Once you have all your 3D segments, you can talk about distances between them as well as rotations (angles) differences between them. Once I lay this groundwork, I'm hoping to get to those things in this thread.

Tech99, I truly do appreciate your interest.

Best Regards,
Elroy

EDIT1: Tech99, I just re-read your post again. Basically, it's a 3D segment's origin that's attached to the distal center of some bone. The direction that the axes are pointing tell us whether that bone is twisted or tilted with respect to the previous 3D segment (such as the pelvis with respect to the femur).

9. ## Re: VB6 and some 3D Linear Algebra

Originally Posted by Elroy
... Once I lay this groundwork, I'm hoping to get to those things in this thread.
It's time to get concrete IMO (because otherwise the interest will be waning, I guess)...

Frankly, I don't really see where the problem is...

If you're familiar with the concept of "Planes" and "Normal-Vectors", you can easily calculate angles between anything.

A Plane in 3D-Space can be determined by any 3 Points.
And a Normal-Vector on such a plane (the direction the Plane points to) is easy to calculate (from these 3 Points which define the Plane).

Now it's only up to you, to define:
- which set of 3 points will make up "a useful plane"
- calculate the Normal on it
- analyze - or compare that Normal-Vector with Normals from "other useful planes"

To give a concrete example...:
Let's say we have "wired" a person with only 3 LEDs (3 Point-Sources), to keep it simple...
- two of them at the right and left hip
- one of them at the nape of the neck

The Triangle between these 3 points is sufficient to define "a useful plane" (let's name it: "upper-body-plane").
If the person is standing still, erect and untwisted (and we assume the idealized case, that the neck-point was mounted "straight above the hip-points" somehow),
then the Normal-Vector of that given plane will point "straight ahead" (parallel to the ground, in the same direction the person is gazing).

If the person bends its upper-body by 90 degrees now, then the Normal-Vector of the upper-body-plane will point straight down.

If the person remains "ideally erect", but is twisting its upper-body (its hips) now, then the Normal-Vector will veer in the appropriate direction(s),
whilst remaining parallel to the ground...

And so on...

HTH

Olaf

10. ## Re: VB6 and some 3D Linear Algebra

Originally Posted by Elroy
EDIT1: Tech99, I just re-read your post again. Basically, it's a 3D segment's origin that's attached to the distal center of some bone. The direction that the axes are pointing tell us whether that bone is twisted or tilted with respect to the previous 3D segment (such as the pelvis with respect to the femur).
Sure, it is so - but that joint origin moves along, when upper joint and connected vertices below it rotates - so there is no use of 'origin' value in vertex types. Origin must be calculated, but you don't need that info (where the origin position is in space) anyway.

11. ## Re: VB6 and some 3D Linear Algebra

Originally Posted by Schmidt
It's time to get concrete IMO (because otherwise the interest will be waning, I guess)...

Hi Olaf,

I'm glad to see you reply. Sorry I haven't been more productive on this thread. I've just got a few different balls (or are they knives) in the air this week. Things will slow down dramatically next week (and even this weekend). I promise I'll get there.

Also, I've thought of you repeatedly as I was composing this thread. I don't have anything concrete to refer to, but I just had an idea that you were versed in linear algebra.

More to come very shortly.

@Tech99: I'm certainly not saying you're wrong. You're not, and I can see precisely what you're saying. However, this is just the way they do it.

Also, please remember (as outlined in post #4) that origin offsets are not stored in the three axes that define the segment. Therefore, the axes only give us the rotation. We must also have the 3D segment's origin if we are to know its location in space. Again, you're absolutely right that there's a bit of redundancy built into these 3D segments. You could eliminate any one of the four vectors in one of these 3D segments and still have a complete definition (well, that is, if the offsets were built into the axes) (But you can still "lose" any axis, and just do a right-hand-rule Cross Product to recover it).

12. ## Re: VB6 and some 3D Linear Algebra

Okay, Olaf has me feeling like I need to get to a question. Therefore, I'll jump ahead.

Above, I've outlined the concept of what I (and the mocap community) call a 3D Segment. What I need is a way to extract the angles that would rotate one 3D Segment to another 3D Segment, and then a way to use those angles to rotate the first 3D Segment to precisely match the second 3D Segment. Basically, that's #3 and #2 from the OP.

I've got the #3 piece, and it's working correctly. However, I can't quite get the #2 piece working correctly. I was going to talk about how these 3D segments are constructed (#1 in the OP). But I'll forego that for now. Let's assume that we've got two different 3D segments (keep imagining my little wooden figure as these 3D segments).

We can imagine a set-of-angles that tell us how to rotate one 3D segment so that it matches the rotation of a second 3D segment. There are several ways to represent these angles. One way is through quaternions. Yet another is through something called Euler angles. Yet a third is an approach that extracts the angles using a fixed 3D segment. (And there are others.) We'll be focusing on the Euler and fixed segment approaches (leaving quaternions for later). There is actually a precise relationship between fixed segment angles and Euler angles, and we'll talk about that.

Also, when dealing with Euler angles, we must attend to the order in which we calculate the angles. To simplify this discussion, I'll just always assume that we're extracting angles in an XYZ order.

Okay, I've started talking in a very high-level way. I hope some are still following.

Now, before we extract a set of Euler angles between two 3D segments, we need a UDT for storing them. However, it turns out we already have one. A set of Euler angles is just three angles, one for X, one for Y, and one for Z. Again, we can just think of rotating the X, Y, & Z axes of our 3D segments. We've already got a VectorType UDT. This UDT will serve perfectly for storing Euler angles. This time, we won't be using the UDT to represent a point-in-space, but it is still a type of vector (just not a 3D coordinate vector).

Therefore, we're set. Here's code to extract Euler angles between two 3D segments:

Code:
```Option Explicit

Public Type VectorType
x As Double
y As Double
z As Double
End Type
'
Public Type SegmentType
Origin As VectorType
AxisX  As VectorType
AxisY  As VectorType
AxisZ  As VectorType
End Type
'
Public Type LineType
Point1 As VectorType
Point2 As VectorType
End Type
'
Public Type ColType
c(1 To 3) As Double
End Type
'
Public Type MatrixType
r(1 To 3) As ColType
End Type

Public Function AnglesBetweenSegments(Child As SegmentType, Parent As SegmentType) As VectorType
'
' This function is designed to replicate the line in BodyBuilder that allows you to create a Euler Angle function (but just for one frame):
'       AnglesBetweenSegments = <Child, Parent, RotationOrder>
'
' Use the NegateVector function if you wish to do the -<Parent, Child, Rotation> method.
'
' There are two approaches to calculating angles such as these:  1) Fixed axis, 2) Moving axis.
'   The "moving axis" approach is what is typically thought of as Euler.
'   The "fixed axis" approach always returns to the parent to determine the next rotation.
'   By default, this function uses the "fixed axis" approach, however, just reverse Parent and Child, and negate, and it'll be "moving axis" (Euler).
'
Dim ChildMatrix     As MatrixType
Dim ParentMatrix    As MatrixType
Dim NewMatrix       As MatrixType
Dim AngleX  As Double
Dim AngleY  As Double
Dim AngleZ  As Double
'
ChildMatrix = MatrixFromSegment(Child)        ' Ignores the origin.
ParentMatrix = MatrixFromSegment(Parent)      ' Ignores the origin.
'
NewMatrix = TransposeMatrix(MultMatrices(TransposeMatrix(ChildMatrix), ParentMatrix))
'
AngleY = ArcSin(NewMatrix.r(1).c(3)): AngleZ = ArcTan2yx(-NewMatrix.r(1).c(2), NewMatrix.r(1).c(1)): AngleX = ArcTan2yx(-NewMatrix.r(2).c(3), NewMatrix.r(3).c(3))
'
'
AnglesBetweenSegments = MakeVector(AngleX, AngleY, AngleZ)
'
End Function

Public Function MatrixFromSegment(Segment As SegmentType) As MatrixType
' This just copies the axes, not the origin.
' The XYZ of each axis make up the columns.
' The Separate axes make up the rows.
'
MatrixFromSegment.r(1).c(1) = Segment.AxisX.x
MatrixFromSegment.r(1).c(2) = Segment.AxisX.y
MatrixFromSegment.r(1).c(3) = Segment.AxisX.z
MatrixFromSegment.r(2).c(1) = Segment.AxisY.x
MatrixFromSegment.r(2).c(2) = Segment.AxisY.y
MatrixFromSegment.r(2).c(3) = Segment.AxisY.z
MatrixFromSegment.r(3).c(1) = Segment.AxisZ.x
MatrixFromSegment.r(3).c(2) = Segment.AxisZ.y
MatrixFromSegment.r(3).c(3) = Segment.AxisZ.z
End Function

Public Function TransposeMatrix(m As MatrixType) As MatrixType
TransposeMatrix.r(1).c(1) = m.r(1).c(1)
TransposeMatrix.r(1).c(2) = m.r(2).c(1)
TransposeMatrix.r(1).c(3) = m.r(3).c(1)
TransposeMatrix.r(2).c(1) = m.r(1).c(2)
TransposeMatrix.r(2).c(2) = m.r(2).c(2)
TransposeMatrix.r(2).c(3) = m.r(3).c(2)
TransposeMatrix.r(3).c(1) = m.r(1).c(3)
TransposeMatrix.r(3).c(2) = m.r(2).c(3)
TransposeMatrix.r(3).c(3) = m.r(3).c(3)
End Function

Public Function MultMatrices(m1 As MatrixType, m2 As MatrixType) As MatrixType
Dim r1Ptr As Long
Dim c1Ptr As Long
Dim c2Ptr As Long
'
For r1Ptr = 1 To 3
For c2Ptr = 1 To 3
For c1Ptr = 1 To 3
' This uses column ptr for row, but that's okay.  It works this way.
MultMatrices.r(r1Ptr).c(c2Ptr) = MultMatrices.r(r1Ptr).c(c2Ptr) + m1.r(r1Ptr).c(c1Ptr) * m2.r(c1Ptr).c(c2Ptr)
Next c1Ptr
Next c2Ptr
Next r1Ptr
End Function

Public Function ArcSin(d As Double) As Double
' d is forced between 1 and -1.  This prevents IEEE rounding from causing any problems.
Dim dRoot As Double
Dim d2 As Double
'
Select Case True
Case d > 1:     d2 = 1
Case d < -1:    d2 = -1
Case Else:      d2 = d
End Select
'
dRoot = Sqr(1 - d2 * d2)
If Abs(dRoot) <> 0 Then
ArcSin = Atn(d2 / dRoot)
Else
ArcSin = (pi / 2) * Sgn(d2)
End If
End Function

Public Function ArcTan2yx(y As Double, x As Double) As Double
If x = 0 Then
ArcTan2yx = 0
Exit Function
End If
'
ArcTan2yx = (pi / 2) * Sgn(y)
'
On Error Resume Next
ArcTan2yx = Atn(y / x)
On Error GoTo 0
'
Exit Function
End Function

Public Function Rad2Deg(rad As Double) As Double
Rad2Deg = rad / (pi / 180)
End Function

Public Function MakeVector(x As Double, y As Double, z As Double) As VectorType
MakeVector.x = x
MakeVector.y = y
MakeVector.z = z
End Function

Public Function pi() As Double
Static dPI As Double
If dPI = 0 Then dPI = 4 * Atn(1)
pi = dPI
End Function```
Now, I hope that everyone isn't overwhelmed. It's that AnglesBetweenSegments that I'm introducing. The rest is just support. Now, the way it's written, it returns fixed segment angles. However, if we swap Parent and Child, and then negate all the returned angles, it returns Euler angles. And this is the way it's typically used (i.e., Euler angles).

Now, that function works. It's been thoroughly tested, and it matches the results from other non-VB6 code.

However, here's the catch. With the returned Euler angle vector, it would seem to be somewhat straightforward to use those angles and to rotate the the Parent 3D segment such that its rotation matches the Child 3D segment. That's what I'm struggling with.

I've got some attempts, but they're not working out. I'm going to get some dinner, and maybe I'll post that code later. To relate to the OP, the code above is essentially point #3 (Extract angles between two 3D segments), and I'm struggling to develop code that does point #2 (Rotate a 3D segment).

Again, if I can find the time, I'll post my attempt later this evening.

All The Best, and sorry for going so high-level with this.
Elroy

13. ## Re: VB6 and some 3D Linear Algebra

Okay, just to conclude this thread, I found my problem. I had a bug in my ATan2 function that crept up under specific circumstances.

I'm now good to go. For anyone interested, here's a pretty good ATan2 function:

Code:
```
Option Explicit

Public Function ArcTan2yx(y As Double, x As Double) As Double
'
On Error GoTo Atan2Error
ArcTan2yx = Atn(y / x)
If (x < 0) Then If (y < 0) Then ArcTan2yx = ArcTan2yx - pi Else ArcTan2yx = ArcTan2yx + pi
Exit Function

Atan2Error:
If Abs(y) > Abs(x) Then     ' Must be an overflow.
If y > 0 Then ArcTan2yx = pi / 2 Else ArcTan2yx = -pi / 2
Else
ArcTan2yx = 0           ' Must be an underflow.
End If
Resume Next
End Function

Public Function pi() As Double
Static dPI As Double
If dPI = 0 Then dPI = 4 * Atn(1)
pi = dPI
End Function
```
All is good with me now, at least on this issue. Resolved.

14. ## Re: VB6 and some 3D Linear Algebra

Code:
```Private Const PI As Double = 3.14159265358979
Private Const PI_2 As Double = 1.5707963267949

Function ArcTan2yx(Y As Double, X As Double) As Double
Select Case X
Case Is > 0
ArcTan2yx = Atn(Y / X)
Case Is < 0
ArcTan2yx = Atn(Y / X) + PI * Sgn(Y)
If Y = 0 Then ArcTan2yx = ArcTan2yx + PI
Case Is = 0
ArcTan2yx = PI_2 * Sgn(Y)
End Select
End Function```

15. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Hi Tech99,

The only problem I see with yours is that X might be very small, but not quite zero, and still cause an error. The following code illustrates the problem:

Code:
```
Option Explicit

Dim x As Double
x = 9.88131291682493E-324

Debug.Print x > 0   ' Reports "True".

Debug.Print 5 / x   ' Generates an error.

End Sub
```
I'm just not sure there's a bullet-proof way to avoid error trapping.

Best Regards,
Elroy

16. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

I don't think that, you can get that small co-ordinate value via input or computation, to cause a division by zero error. I haven't faced that kind of error either.

17. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Code:
```Private Const PI As Double = 3.14159265358979
Private Const PIh As Double = 1.5707963267949
Function ArcTan2yx(Y As Double, X As Double) As Double
If X Then
ArcTan2yx= -PI + Atn(Y / X) - (X > 0#) * PI
Else
ArcTan2yx= -PIh - (Y > 0#) * PI
End If
End Function```

18. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

To everyone who shared an Atan2 function so far - you will get better performance by passing X and Y as ByVal. This is true for all value types, including Doubles. (And since you aren't modifying X/Y inside the function, why wouldn't you use ByVal?)

@Elroy: how critical is performance? There are many Atan2() implementations that don't require explicit error-handling. If you have a specific degree/radian threshold you need to maintain, I imagine we could find you a faster implementation. (Also, why aren't you hard-coding the value of PI as a constant? I think you've set a new record for the slowest conceivable way to use PI at run-time! )

19. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Yeah, Tanner, you're probably right about PI. Ok ok, I'll re-code it as a constant.

Regarding performance, I don't think any of this math stuff is where my performance issues will be, as I'll also be doing file I/O while this stuff is being performed. I suspect my file I/O is orders of magnitude slower than any math I'll be doing.

But I don't just want to throw away performance for the heck of it.

Now, regarding an ATan2() function, I'm certainly willing to listen/consider others, so long as they address the "near zero" problem. I certainly believe Tech99 when he says he's never seen it, but that doesn't mean it won't happen. As suggested in the earlier posts of this thread, I'm rotating these 3D segments all over the place (and then calculating Euler angles between them), and there's absolutely nothing that says I might not hit a "near zero" situation. I suppose I could test for an epsilon value, but error trapping just seemed easier. I'd be curious how other bullet-proof routines are doing it (possibly in other languages).

Best Regards,
Elroy

20. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Thanks for the reply, Elroy. I've actually seen that 4 * Atn(1) hack before, but only in languages that support building against different architectures. If you can't guarantee the resolution of floating-point values in advance, that hack can give you a run-time approximation of the "most accurate" version of pi on your current architecture. So it's not without merit - but VB6 only targets one architecture, so floating-point accuracy is guaranteed in advance. This means we can use hard-coded constants for maximum precision and much better performance.

As for Atan2() implementations, here's a very fast 4-quadrant approximation I use frequently. Credit due to Jim Shima for the original derivation (source here).

Code:
```Private Const PI_14 As Double = 0.785398163397448
Private Const PI_34 As Double = 2.35619449019234

Private Function Atan2_Fastest(ByVal y As Double, ByVal x As Double) As Double

'Cheap non-branching workaround for cases where y approaches 0.0
Dim absY As Double
absY = Abs(y) + 0.0000000001

If (x >= 0#) Then
Atan2_Fastest = PI_14 - PI_14 * (x - absY) / (x + absY)
Else
Atan2_Fastest = PI_34 - PI_14 * (x + absY) / (absY - x)
End If

If (y < 0#) Then Atan2_Fastest = -Atan2_Fastest

End Function```
That implementation has a maximum error of 0.07 radians. If that's too high, the source link above provides slightly modified alternatives for better accuracy, at some cost to performance.

That's the fastest implementation I have. If that max error is too high, a classic alternative is Hastings approximation. My source is the old Apple dev lists:

Code:
```Private Const PI_HALF As Double = 1.5707963267949
Private Const PI As Double = 3.14159265358979
Private Function Atan2_Faster(ByVal y As Double, ByVal x As Double) As Double

If (x = 0#) Then
If (y > 0#) Then
Atan2_Faster = PI_HALF
ElseIf (y = 0#) Then
Atan2_Faster = 0#
Else
Atan2_Faster = -PI_HALF
End If
Else
Dim z As Double
z = y / x
If (Abs(z) < 1#) Then
Atan2_Faster = z / (1# + 0.28 * z * z)
If (x < 0#) Then
If (y < 0#) Then
Atan2_Faster = Atan2_Faster - PI
Else
Atan2_Faster = Atan2_Faster + PI
End If
End If
Else
Atan2_Faster = PI_HALF - z / (z * z + 0.28)
If (y < 0#) Then Atan2_Faster = Atan2_Faster - PI
End If
End If

End Function```
Max error is quite a bit lower, <= 0.005 radians, and if passed values tend to be similar, branch-prediction in modern CPUs eliminates most of the performance penalties associated with the nested Ifs.

One last question (sorry for the post length). Where are your x/y values coming from? If they're coming from a hardware source, there's likely a limit to their accuracy. Hardware can't typically approach anything near the 1e-324 "tiny" number you demonstrated above, which is why most solutions don't deal with potential "near zero but not zero" errors. Fast Atan2 workarounds are very common in game and graphics development, and in the same way, near-zero issues simply don't exist because the resolution of passed (x, y) values is always guaranteed to be much lower than what a Single or Double actually supports.

(And if your hardware can return values as small as 1e-324, you can just forcibly add something like 1e-10 to incoming x/y values to ensure a minimum resolution.)

21. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Hi Tanner,

Thanks a million for the ATan2 algorithms. I'll definitely take a look at them.

Regarding where the X & Y are coming from, it's a combination of hardware inputs and internal calculations. I'll just briefly outline it.

Basically, everything starts with 3D (x,y,z) points moving in space. These are actually measured and reported in IEEE single precision millimeters. However, once these are in a file (through motion capture), they're almost immediately read, with many "virtual" points created (such as taking the mid-point between two points to find things like a joint center). Then, these "real" and "virtual" points are used to create 3D segments (as outlined above). So far, all we've done is some 3D vector math and some cross-products.

However, once these 3D segments are assembled, Euler angles are calculated between them. That's where ATan2() comes in. For instance, we may want to model the rotational relationship between the top of the femur (the hip joint) and the top of the tibia (the knee joint). If done dynamically (through time), this would basically tell us how the knee joint is moving on all three planes (flexion-extension, varus-valgus, axial-rotation), with respect to the femur.

Again, Thanks,
Elroy

22. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Very cool project, Elroy. I'm jealous!

Working in mm space should guarantee no problem with OOB errors. Single-precision floats bottom-out at 10^-38, so unless your equipment is in violation of all known laws of physics (e.g. with precision smaller than the Planck length), you won't get (x, y) values anywhere near the upper/lower limits of floats. And, even if it did produce results that precise, performing your own calculations with Doubles provides an added safeguard. You should be good to go!

Best of luck with the rest of your project.

23. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Hmmm, I forgot to say it in post #21, but yes, I do all internal math with Doubles.

Also, I've been rather exclusively working on this project since 1999, and it'll never be finished. As soon as I get one "chunk" done, they immediately put me to work on something else. I've done specific pieces of linear algebra projects on the past for this VB6 project. However, this time, I'm attempting to build more general purpose functions.

I'm modeling all of my work on an extremely high-level language called BodyBuilder by Vicon. Vicon is the "big boy" in the motion capture (mocap) industry. They've even assisted with several movies (Polar Express, Avatar. etc.). However, all of my stuff is medical/clinical.

In BodyBuilder, everything is IEEE Single, but there are basically five types of variables:
• A scalar (which is just a float number).
• A 3D point (x,y,z, millimeters).
• A line (which is effectively defined by two 3D points).
• A 3D segment (which I tried to outline in the first part of this thread).
• A set of Euler angles (three of them, one for each plane).

And, beyond all the standard math operators (which are just like VB6), and a relatively short list of functions (such as Sin(), Cos(), etc.), there are effectively three high-level operations:

• Construct (and return) one of those 3D segments from a point (the origin) and two lines. Basically, if you recognize that a line is just two points, it constructs a 3D segment from five points.
• Rotate (and return) a 3D segment, given the starting segment, an axis of rotation, and the degrees to rotate. For those who recognize it, this is basically the same as a quaternion rotation.
• Return the Euler angles between two 3D segments.

Also, with the "construct 3D segment" and the "return Euler angles", you must specify an "order constant". There are six of them: xyz, xzy, yzx, yxz, zxy, & zyx. If you've ever dealt with Euler angles, this will be familiar. Fortunately, quaternion rotations are the same regardless of order.

Outlined like that, it all seems pretty simple. However, it gets complicated pretty quickly.

Currently, I'm working on what we call a Spine Protocol for kids with idiopathic scoliosis. Basically, we put a bunch of infra-red (IR) markers on kids, and then have them go through a series of motion (think ... forward bending, back bending, left rotation, side bending, etc). Then, we take those IR markers, treating them as points in space, calculate those 3D segments, and figure out what their motion is. Nuthin' to it.

Actually, using Vicon hardware, that particular protocol is all done. However, I'm currently working on bringing the linear algebra into VB6 (and out of the Vicon BodyBuilder code) so that other non-Vicon hardware systems can use the protocol. I'm also close to done with that as well, but it'll still be a couple of weeks.

Hey, thanks for being interested. I love talking about this stuff.

Tanner, you take care, and have a GREAT weekend.
Elroy

24. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

@reexre...

Your ATan2 in post #17 has a problem similar to my original problem. For example, try x=-6, y=5.

The return should be = 2.446854, but your function returns -3.8363309. In other words, it's 180 degrees out.

@Tech99...

I tried yours under all possible positive, negative, and zero conditions, and all seems copacetic. Since Tanner has convinced me that I don't need to worry about the "near zero" conditions, that's the one I'm going with. I did tweak it just a bit, but that's just me. Here's my final version:

Code:
```
Option Explicit

Public Const pi = 3.14159265358979
Public Const pihalf = 1.5707963267949

Public Function ArcTan2yx(Y As Double, X As Double) As Double
'
Select Case X
Case Is > 0#: ArcTan2yx = Atn(Y / X)
Case Is < 0#: ArcTan2yx = Atn(Y / X) + pi * Sgn(Y): If Y = 0# Then ArcTan2yx = ArcTan2yx + pi
Case Else:    ArcTan2yx = pihalf * Sgn(Y)
End Select
End Function
```

@Tanner...

I truly do appreciate your approximations, and I must admit that I do use them for things such as a Gaussian curve or F curve, but the exact method (even though maybe a tad slower) seems easy in this case.

Best Regards To All,
Elroy

EDIT1: Actually, reexe, after staring at it some more, I've realized that yours is okay. It's just returning values greater than pi, which gives me problems.

EDIT2: I just can't leave this thing alone. I tweaked it to the following:

Code:
```
Function ArcTan2yx(Y As Double, X As Double) As Double
Const epsilon As Double = 1.4E-45
'
Select Case X
Case Is > epsilon:  ArcTan2yx = Atn(Y / X)
Case Is < -epsilon: ArcTan2yx = Atn(Y / X) + pi * Sgn(Y): If Y = 0# Then ArcTan2yx = ArcTan2yx + pi
Case Else:          ArcTan2yx = pihalf * Sgn(Y)
End Select
End Function
```
Unless Y is VERY large coming in (which it never will be in my cases), this version should protect us from the "near zero" cases.

25. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

@reexre...

Your ATan2 in post #17 has a problem similar to my original problem. For example, try x=-6, y=5.

The return should be = 2.446854, but your function returns -3.8363309. In other words, it's 180 degrees out.
2.446854-atn(1)*8 = -3.83633130717959
atn(1)*8 = 2PI = 360

EDIT1: Actually, reexe, after staring at it some more, I've realized that yours is okay. It's just returning values greater than pi, which gives me problems.
??? did not understand

@Tanner
Thank you for your approximation function. By the way , why byVal is faster than byRef ? , I always thought the opposite...

26. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Originally Posted by reexre
2.446854-atn(1)*8 = -3.83633130717959
atn(1)*8 = 2PI = 360

??? did not understand
Some people use 0 to 360, i.e. 0 to 2PI, others require -180 to 180, so -PI to PI. Its just a convention requirement.

27. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Originally Posted by reexre
@Tanner
Thank you for your approximation function. By the way , why byVal is faster than byRef ? , I always thought the opposite...
Hi reexre. ByRef is only faster if the source data is very large. For value data types, like Doubles, passing objects ByRef requires the CPU to use some form of indirect addressing which carries a performance penalty. If compiling to native code, you can somewhat reduce the penalty with the "Assume No Aliasing" optimization, but that carries its own risks. (This is true in all compiled languages.)

Because passing ByRef is also riskier (you can accidentally modify input parameters, aliasing issues, etc), there is generally not much reason to use it... unless you actually need to modify those parameters, or if VB doesn't physically allow you to pass an object ByVal (e.g. arrays). IMO it's unfortunate that ByRef is the default mode in VB... one of those historical oddities they thankfully rectified in VB.Net.

(Other languages have the benefit of a "const" keyword for reference parameters, which mitigates some of these issues, but I guess we'll need to wait for VB7 for that... )

28. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Elroy, your 'epsilon constant' function at #24 post fails, fex. try with very large Y and very small X input values.

Code:
`Debug.Print ArcTan2yx(1E+46, 1E-46)`
With these (imaginary) input values, it falls to case else.

29. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

Originally Posted by Tech99
Elroy, your 'epsilon constant' function at #24 post fails, fex. try with very large Y and very small X input values..
Hi Tech99,

Yeah, I knew that there would still be situations where it would fail. In fact, I mentioned that in my post.

Unless Y is VERY large coming in (which it never will be in my cases), this version should protect us from the "near zero" cases.
However, it's really only the "very small X" values that I'm worried about. My X and Y are millimeters measured in a room that would be 20 meters at its largest point (and even that is WAY over the limit). Therefore, my largest 3D value would be 20000, which is far less than 1E+46.

However, through calculations and mathematical massaging, I did feel it was possible to wind up with a number that was very near zero (but not quite zero). I think Tanner is correct in that it would happen very seldom, but I still saw it as a possibility.

Therefore, for my purposes, I think my epsilon value keeps me out of trouble, and allows me to proceed without error trapping. But I think your caution is warranted for those who may use that ATan2 function for other purposes.

Best Regards,
Elroy

30. ## Re: [RESOLVED] VB6 and some 3D Linear Algebra

@Tanner
Thank you, I'll try to use more often byVal, following your advice.
even if I was used to use byRef. (I almost never used byVal) . Thanks again.

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•

Featured

Click Here to Expand Forum to Full Width