THIS PAGE HAS BEEN MOVED TO SHAREPOINT!
Please refer to this site/make edits here for the most updated information: https://partnershealthcare.sharepoint.com/sites/LCN/SitePages/Morpho-Optimization-Project.aspx
Parent: MorphoOptimizationProject_BetterSerialCode
In mrisurf.c the function mrisComputeDefectMRILogUnlikelihood_wkr does a huge fraction of the mris_fix_topology work, and has had a lot of optimizations applied to it.
Several of these exploit the spatial nature of the data.
One of the last of these has been given a brick [iMin..iMax, jMin..jMax, kMin..kMax] of voxels near a triangular face, and is filling in each of the voxels with the distance to the nearest point in the face. This requires determining whether the nearest point is within the triangle, on an edge of the triangle, or a vertex of the triangle. This determination is complicated, as is the distance calculation once it has been decided.
However the code is determining the minimum of these for several different faces, not just one face. Once the minimum distance to one face is known, the code does a quick calculation to determine whether this face might be closer. It does this by using some geometry and some algebra.
The line from the ijk voxel to the nearest point on the face can be one edge of a triangle, where the other two edges go from the ijk voxel to the center of the face and from the center of the face to the nearest point in the face. This triangle lets us assert that
distance(IJK, Center) <= distance(IJK, nearest_point) + distance(nearest_point, Center)
which can be rewritten as
distance(IJK, Center) - distance(Center, nearest_point) <= distance(IJK,nearest_point)
But the distance from the center to the nearest point must be less than from the center to the furtherest vertex from the center so this can be used to create another inequality
distance(IJK, Center) - distance(Center, furtherest_vertex) <= distance(IJK, nearest_point)
But distance(IJK, nearest_point) will only reduce the minimum when
distance(IJK, nearest_point) <= minimum already discovered
So if
minimum already discovered <= distance(IJK, Center) - distance(Center, furtherest_vertex)
then IJK can not reduce the minimum. This test is used as a filter to discard IJK before deciding where on the face the nearest point is, which results in a significant decrease in execution time.
Calculating the distance(IJK, Center) requires a sqrt, and sqrt is an expensive operation, one that should be eliminated. To do so, first rewrite this inequality as
minimum already discovered + distance(Center, furtherest_vertex) <= distance(IJK, Center)
The LHS is the sum of two positive terms, and hence is positive. The RHS is also positive, and squaring is a monotonically increasing operation, so this inequality is true if and only if the squares are similarly ordered.
square(minimum already discovered + distance(Center, furtherest_vertex)) <= distanceSquared(IJK, Center)
Since the IJK does not appear in the LHS, this term is a loop invariant of the loops over i,j,k, and so the test turns unto three subtracts, three multiples, two adds, and a comparison!
Even if it only eliminated 50% of the IJK from consideration, it would pay for itself, and it does quite a bit better than that.