VOXEL-BASED SURFACE FLATTENING

1. VOXEL-BASED SURFACE FLATTENING Nahum Kiryati* Dept. of Electrical Engineering – Systems Tel Aviv University * Joint work with Ruth Grossmann and Ron Kimmel.

2. SURFACE FLATTENING Transform a surface from 3D space to 2D plane. • Some distortion is usually unavoidable • Need to define the most important property to preserve (area, angles,...) • Use cuts (interruptions) to selectively reduce distortion Peters Projection: area accuarte Cuts reduce distortion

3. Applications: Medical • Cortical surface flattening for visualization, registration, etc. • Virtual colonoscopy Gray-white matter interface Courtesy of Gabriele Lohmann, MPG Leipzig

4. Applications: Computer Graphics Texture Mapping: Projecting a planar pattern onto a 3D surface • Flatten the surface • Map the 2D pattern onto the flattened surface • Map the flattened surface back to 3D Texture Mapping From www.okino.com

5. Some previous works • E.L. Schwartz, A. Shaw and E. Wolfson, “A Numerical Solution to the Generalized Mapmaker’s Problem: Flattening Nonconvex Polyhedral Surfaces”, IEEE Trans. Pattern Anal. Mach. Intell., Vol. 11, pp. 1005-1008, 1989. • C. Bennis, J.M. Vezien and G. Iglesias, “Piecewise Surface Flattening for Non-Distorted Texture Mapping”, Computer Graphics, Vol. 25, pp. 237-247, 1991. • H.A. Drury, D.C. Van Essen, C.A. Anderson, C.W. Lee, T.A. Coogan and J.W. Lewis, “Computerized Mappings of the Cerebral Cortex: A Multiresolution Flattening Method and Surface-Based Coordinate System”, J. of Cognitive Neuroscience, Vol. 8, pp. 1-28, 1996. • S. Angenent, S. Haker, A. Tanenbaum and R. Kikinis, “Conformal Geometry and Brain Flattening”, Proc. 2nd Intl. Conf. on Medical Image Computing and Computer Assisted Intervention (MICCAI), Cambridge, England, pp. 269-278, 1999. • A.M. Dale, B. Fischl and M.I. Sereno, “Cortical Surface-Based Analysis. II: Inflation, Flattening and a Surface-Based Coordinate System”, NeuroImage, Vol. 9, pp. 195-207, 1999. • S. Haker, S. Angenent, A. Tanenbaum, R. Kikinis, G. Sapiro and M. Halle, “Conformal Surface Parameterization for Texture Mapping”, IEEE Trans. on Visualization and Computer Graphics, Vol. 6, No. 2, 2000.

6. New Method: Unique Features • Operates on voxel data: Triangulation not necessary • Optimal (global preservation of geodesic distances) • Solid mathematical basis. Analytic solution (in principle) • Computationally efficient Voxel data From www.chez.com

7. Our Approach to Surface Flattening • Mark N (very many) points on the 3D surface. • Estimate the geodesic distances between each point and all others (how?) • Optimally map the N points onto the 2D plane, such that the 2D Euclidean interpoint distances will be as close as possible to the corresponding 3D geodesic distances (how?) • Regard the N points as control points, and map the rest by interpolation.

8. Crucial problem #1: Estimate the minimal (geodesic) distance between points on a 3D surface. Kiryati & Szekely, 1993: Accurate and efficient solution for surface represented by voxels. • Represent the surface as a weighted graph: • The weight of an edge depends on the link type • Minimal geodesic distance ~ shortest path in graph • Graph is sparse: very fast algorithm direct, minor diagonal, major diagonal Question: what weights should we use?

9. Crucial problem #1: Estimate the minimal (geodesic) distance between points on a 3D surface. Kiryati & Székely, 1993: Accurate and efficient solution for surface represented by voxels. • Represent the surface as a weighted graph: • The weight of an edge depends on the link type • Minimal geodesic distance ~ shortest path in graph • Graph is sparse: very fast algorithm direct, minor diagonal, major diagonal Question: what weights should we use? Wrong answer:

10. Bias in length estimation (2D example) • Continuous curve (length L) Digital curve • Digital curve Chain code (8-connected) • Estimate length: Unbiased, minimum MSE* length estimation: *Unbiased and minimum MSE for straight lines. Unbiased for lines  wonderful for curves (errors cancel out).

11. Kiryati & Kübler, 1992: Unbiased, minimum MSE length estimation for 3D curves • RMS error for straight lines: 2.88% • Straight lines are the worst case • (much better with curves) Use the weights (0.9016, 1.289, 1.615) in the algorithm for finding shortest paths (geodesics) and minimal distances on digitized surfaces. Verwer, Beckers & Smeulders

12. Crucial problem #1: Estimate the minimal distance between points on a digitized 3D surface. • Solved (Kiryati & Székely, 1993)

13. Reminder: Our Approach to Surface Flattening Mark N (very many) points on the 3D surface. Estimate the geodesic distances between each point and all others (how?) Optimally map the N points onto the 2D plane, such that the 2D Euclidean interpoint distances will be as close as possible to the corresponding 3D geodesic distances (how?) Regard the N points as control points, and map the rest by interpolation.

14. Crucial problem #2: Given N points on a surface in 3D, map them onto the 2D plane, such that the 2D Euclidean interpoint distances will be as close as possible to the corresponding geodesic distances on the 3D surface.

15. Crucial problem #2: Given N points on a surface in 3D, map them onto the 2D plane, such that the 2D Euclidean interpoint distances will be as close as possible to the corresponding geodesic distances on the 3D surface. Special case of Multi-Dimensional Scaling (MDS)

16. Goal: Place N points in the 2D plane, such that the sets and are as similar as possible. Objective functions The analytic solution for minimal strain is called classical scaling (Borg & Groenen, 1997).

17. Classical Scaling XCoordinate Matrix:(N x 2) matrix with the (x,y) coordinates of N points in the plane1(unknown) BScalar Product Matrix:(N x N)matrix (unknown) d* Squared Distances Matrix:(N x N)matrix containing the squared Euclidean distances between points in the plane(unknown) Suppose that d* is known. It can be shown that: where and 1 is a vector of 1’s. Having B, X can be found via eigendecomposition: holds exactly two (positive) non-zero eigenvalues. So where and are (2x 2) and (N x 2) submatrices. 1Column centered: coordinate system origin at the point-set centroid.

18. Classical Scaling (cont’d) But: we don’t have the Squared Distances Matrix d*, since interpoint distances in the plane are unknown. We do have the interpoint geodesic distances on the 3D surface. Algorithm • Create D* (Squared Geodesic Distance Matrix) • Compute the estimated Scalar Product Matrix: • Compute the eigendecomposition of up to rank 2: • The flattened Coordinate Matrix is: • Classical scaling minimizes the error metric (strain). Zero for a true Scalar Product Matrix obtained from a true Squared Distance Matrix d* containing squared planar Euclidean distances. • Partial eigendecomposition (up to rank 2) is easy and cheap (power method).

19. Reminder: Our Approach to Surface Flattening Mark N (very many) points on the 3D surface. Estimate the geodesic distances between each point and all others. Optimally map the N points onto the 2D plane, such that the 2D Euclidean interpoint distances will be as close as possible to the corresponding 3D geodesic distances. Regard the N points as control points, and map the rest by interpolation. We used N ~ 1000 control points, and mapped the rest using radial basis function interpolation. Any reasonable interpolation method can be used.

20. Unfolded Curled planar surface (synthetic) Euclidean distance (2D) vs. geodesic distance (3D)

21. Unfolded Curled planar surface (synthetic) Can you figure out why errors over small distances are larger than errors over large distances in this case? Euclidean distance (2D) vs. geodesic distance (3D)

22. Unfolded Curled planar surface (synthetic) Can you figure out why errors over small distances are larger than errors over large distances in this case? On this surface, many short geodesic paths are straight, or nearly straight. Straight lines are the worst case for the geodesic distance estimator. Over long distances, paths are not straight, and errors cancel out. Euclidean distance (2D) vs. geodesic distance (3D)

23. A depth image of a human face. Euclidean distance (2D) vs. geodesic distance (3D)

24. A depth image of a human face. Euclidean distance (2D) vs. geodesic distance (3D) The flattened surface is not shown, because the depth image has no texture.

25. A depth image of a human face. Euclidean distance (2D) vs. geodesic distance (3D) Errors are larger than in the previous example, because this surface is not developable – it cannot be flattened without some distortion.

26. Textures

27. Texture Mapping via Surface Flattening

28. Texture Mapping via Surface Flattening

29. Highlights • New surface flattening algorithm • Operates on voxel data. Triangulation is not necessary. • Distance measurements on 3D surface + classical scaling. • Global distance preservation • Essentially analytic solution: Global optimum guaranteed • Surface flattening, texture mapping

30. Appeared

31. Companion paper • Flattening with triangulated data • Geodesic distance measurements via Fast Marching on Triangulated Domains + classical scaling • State of the art texture mapping quality