3 decimals is not a lot of precision for rotation matrices

[Ben Supnik]

Hi Y'all,

Roland and I have been tossing around part numbers, looking at smoothing results. One thing that surprised me is that we needed a rather large 0.05 LDU error margin to 'snap' the segments of the 6x6 webbed dishes together. (The part is built via 8 sub-sections rotated 45 degrees.) 0.05 seems like a 'big' error because part authors can specify vertices to the 0.001.

Today I cracked open the part and took a look at the numbers, and the large error margin now makes sense.

The fundamental problem is absolute error amplification: if the error on a matrix ratio is 1/1000 (0.001 precision) then the actual error in the position of a rotated part is 0.001 * R where R is the approximate radius of rotation. In other words, a part that extends 60 LDU from its center will have an absolute error of 60/1000 (0.06) after rotation. (In fact, the real results are up to 3x worse because error accumulation is additive and we add the X, Y, and Z components of the matrix).

The actual numbers of the 6x6 dish show that the worst case error could have been much worse than this actual part. The bottom line is: 0.001 just isn't _that_ precise for a rotation component that's going to be used for a moderately big part. :-(

This matters because the rounding tolerance for welding rotated part segments for smoothing defines the smallest 3-d detail an author can make without having their triangles collapsed.

I realize that we have a gajillion parts and my going "we should do X" isn't going to actually change the library, but I wonder if this perspective should result in some kind of updated practices for authoring.

Naively, I would suggest two possible changes...

1. We could allow higher precision for matrices. For example, the part is using 0.707 for a 45-degree rotation. If I go into the text file and replace it with 0.707106781186548 I can cut my error margin down by 4x. (Clearly BrickSmith isn't using all of 0.707106781186548 - I just typed a big long number into the text file...BrickSmith is using 32-bit float.)

2. Authors could avoid angle rotations that aren't 90 degree increments. At 90 degrees, the rotation matrix can be _exactly_ represented - all rotational coefficients are 0, 1, or -1. :-) So the 'lock-up' of vertices around 90 degree rotations should be pretty good - certainly in the same ballpark as the geometry itself.

This would mean 'bigger' parts - 90 degree swathes of dish instead of 45 degrees, but the precision would be a lot better.

Cheers
Ben

[Roland Melkert]

If replacing the values works that well, I'm thinking to do those replacements during type 1 line parsing. At the very least you could do reverse lookups for 3 digit 45 or 15 degree values, it won't effect other processing because it will round down to the old value.

This would prevent editing / bloating the parts, and authors could limit higher precision in the file to weird angles only.

[Philippe Hurbain]

I generally use 5 decimals for matrix coefficient. Defaut setup of Datheader is set that way.

[Tim Gould]

I reckon that should be the standard, and anything less should be considered too low.

[Ben Supnik]

Hi Y'all,

I didn't do a comprehensive (or even remotely detailed) check on what kinds of precisions were in the lib.

We could potentially create an automated process to 'up-sample' matrix transforms that fit known angles (and possibly flag parts where the matrix transforms are 'random', as the odds of the part meshing are pretty slim and it's a weird use of matrices).

Cheers
Ben

[Tim Gould]

Hi Y'all,

I didn't do a comprehensive (or even remotely detailed) check on what kinds of precisions were in the lib.

We could potentially create an automated process to 'up-sample' matrix transforms that fit known angles (and possibly flag parts where the matrix transforms are 'random', as the odds of the part meshing are pretty slim and it's a weird use of matrices).

Cheers
Ben

There's plenty of good reasons to have non-obvious matrix transformations, especially on rects abd other surface prims.

But certainly checking against floor and ceil (to various precision) on sin and cos 15n would be a good start.

Tim

[Ben Supnik]

There's plenty of good reasons to have non-obvious matrix transformations, especially on rects abd other surface prims.

But certainly checking against floor and ceil (to various precision) on sin and cos 15n would be a good start.

Tim

From a 3-d modeling background, building water-tight flat solids out of matrix transforms seems odd, but if the primitive is, well, a primitive, then I guess I can see why you'd do that. :-)

EDIT: it seems odd compared to what most 3-d programs do, but I cannot deny that if we know the location of our solid's corners precisely in cartesian space, the transform on the rect can specify an angled rectangle to fit that shape exactly, because the matrix is composed of scales and offsets per axis. The precision requirements might be higher than for raw vertices, but you guys have already addressed that!


I was thinking more specifically of 'rotational symmetry' - that is, for a rotationally symemtric part, whether the sub-part slice is modeled in 'strange' pie slices. This could be detected by the use of 'weird' rotations on sub-parts rather than primitives...

cheers
Ben

[Tim Gould]

The whole format doesn't make much sense now. But it was invented in the 90s before 3D stuff was really standardised to work on meshes.

Pure rotations can be easily checked as the matrix determinant is one.

Tim

[Roland Melkert]

Side note:

Would it be save to re normalize the 3x3rotation part of type 1 lines (if abs(d)==1) just after loading, it would in theory automatically scale all fractions up.

Or i'm i missing something?

[Tim Gould]

No. Because you can do a scale and rotate at once.

And even if your part is a rotation matrix, simply scaling it will make it no longer a rotation matrix.

There is a safe way to do it, but it's complicated. When I get enough spare time I'll write a bit of a document on safely determining the 3D angle of a rotation matrix, and fixing it. It is not, however, a simple procedure!

Tim

[Ben Supnik]

I think "fixing a rounded rotation" is getting way into the territory of a batch process we should consider running on the library, then fast-tracking if the mechanical transformation can be demonstrated and spot-tested to be good. :-)

What I like about pre-processing is that the code runs once and a human can inspect the results. Every bit of 'mesh repair' I put into BrickSmith is code whose results are somewhat hidden to the part author, and whose results (as data) can only be viewed by a programmer with source code and a debugger. It's an opaque process, so the less complicated that black box, the better.

(Or maybe I am just lazy and want to keep my end of the code simple. ;-)

cheers
ben

[Sergio Reano]

I think "fixing a rounded rotation" is getting way into the territory of a batch process we should consider running on the library, then fast-tracking if the mechanical transformation can be demonstrated and spot-tested to be good. :-)

What I like about pre-processing is that the code runs once and a human can inspect the results. Every bit of 'mesh repair' I put into BrickSmith is code whose results are somewhat hidden to the part author, and whose results (as data) can only be viewed by a programmer with source code and a debugger. It's an opaque process, so the less complicated that black box, the better.

(Or maybe I am just lazy and want to keep my end of the code simple. ;-)

cheers
ben

I Agree!
Furthermore, a lot of code (and time) is spent to fix parts definition problems.

Following a pair of useful (I hope) function to extract the rotation portion of a matrix:

Given the matrix M, you call GetMatrixScale to get scale factor Vs of M along x,y,z axis
vs=GetMatrixScale(M)

then you call the M=PurifyMatrix(M,Vs,vt) to get back M as a pure rotation matrix. The VT parameter is a vector containing original M translation values

Finally you can get the rotation angles along each 3 axles that M was built by:
vRot=GetMatrixRot(M)
vRot is a x,y,z vector containing a value representing the rotation around the specific axle that built M

To recreate M you can multiply the 3 matrix resulting by creating a rotation matrix around any of the 3axles by the amount o vRot:
xM=matrix.rotationaxleX(vrot.x)
yM=matrix.rotationaxleY(vrot.y)
zM=matrix.rotationaxleZ(vrot.z)
M= xM x yM x zM

Hope this could useful for someone...

Sergio

About the code: it is for DirectX matrix notation. For openGL standard you may need to transpose matrix values

Code:

Function GetMatrixScale(ByRef M As Matrix) As Vector3
    Dim V As Vector3
    With M
        V = Vec3(CSng(Math.Sqrt(.M11 * .M11 + .M12 * .M12 + .M13 * .M13)), CSng(Math.Sqrt(.M21 * .M21 + .M22 * .M22 + .M23 * .M23)), CSng(Math.Sqrt(.M31 * .M31 + .M32 * .M32 + .M33 * .M33)))
        V.Y = V.Y * Math.Sign(.M21 + .M22 + .M23)
        Return RoundV(V)
    End With
End Function

Function PurifyMatrix(ByVal M As Matrix, ByRef vS As Vector3, ByRef vT As Vector3) As Matrix
    With M
        .M11 /= vS.X : .M12 /= vS.X : .M13 /= vS.X
        .M21 /= vS.Y : .M22 /= vS.Y : .M23 /= vS.Y
        .M31 /= vS.Z : .M32 /= vS.Z : .M33 /= vS.Z

        vT = Vec3(.M41, .M42, .M43)
        .M41 = 0 : .M42 = 0 : .M43 = 0
    End With

    Return M
End Function

Function GetMatrixRot(ByVal M As Matrix) As Vector3
    Dim vRot As Vector3
    With vRot
        If Abs(M.M13) > 1 Then M.M13 = Sign(M.M13)
        .Y = CSng(Math.Asin(M.M13))
        If Math.Sin(.Y) > 0.99999 Or Math.Sin(.Y) < -0.99999 Then
            .X = 0
            .Z = -CSng(Math.Atan2(M.M21, M.M22))
        Else
            .X = CSng(Math.Atan2(M.M23, M.M33))
            .Z = CSng(Math.Atan2(M.M12, M.M11))
        End If
    End With
    Return vRot
End Function