IsoContouring and LevelSets Roger Crawfis Contributors Roger Crawfis

Iso-Contouring and Level-Sets Roger Crawfis Contributors: Roger Crawfis, Han-Wei Shen, Raghu Machiraju, Torsten Moeller, Huang Zhiyong, Fan Ding, and Charles Dyer R. Crawfis, Ohio State Univ.

Iso-contour/surface Extractions 2 D Iso-contour 2/27/2021 3 D Iso-surface R. Crawfis, Ohio State Univ. 2

Contouring - The Problem z. Extracting an iso- surface from an implicit function, that is, z. Extracting a surface from volume data (discrete implicit function), f (x , y , z ) = T 2/27/2021 R. Crawfis, Ohio State Univ. 3

More Formally z. A scalar visualization technique that creates curves (in 2 D) or surfaces (in 3 D) representing a constant scalar value across a scalar field. z. Contour lines are called isovalue lines or isolines. z. Contour surfaces are called isovalue surfaces or isosurfaces 2/27/2021 R. Crawfis, Ohio State Univ. 4

Contouring a 2 D structured grid with contour line value = 5 0 1 1 3 2 1 3 6 6 3 3 7 9 7 3 2 7 8 6 2 1 2 3 4 3 2/27/2021 1. Using interpolation to generate points along edges with the constant value 2. Connect these points into contours using a few different approaches. One of the approaches: . Detects edge intersection. Tracks this contour as it moves across cell boundary. Repeat for all contours R. Crawfis, Ohio State Univ. 5

2 D Contouring – Not so easy z Annotations z Smooth Curves Topographic Map of Jerusalem (Contour interval 10 meters) North is at the top of the map. The Mount of Olives is on the far right, Mount Zion on the left. Mount Moriah rises as a long ridge at the south end of the City of David and continues on past the present Temple Mount, and reaches its highest point outside the Northern walls of the Old City, at the top of the map. http: //www. skullandcrossbones. org/articles/solomontemple 2. htm 2/27/2021 R. Crawfis, Ohio State Univ. 6

Quantitative and Qualitative z In this figure, we can actually read off any value to about three significant digits. z Hard and tedious problem to determine where to place labels. z Usually done as vector graphics, rather than raster graphics for precision. http: //omnimap. com/catalog/cats/fish-contour. htm 2/27/2021 R. Crawfis, Ohio State Univ. 7

Divide and Conquer z. Partition space into rectangular subregions. 2/27/2021 R. Crawfis, Ohio State Univ. 8

Dividing Cubes/Squares z Generates Points and renders them z Rendering points is faster than polygons z Principle - Divide a square until contour passes thru the cell 0 1 1 3 2 1 3 6 6 3 3 7 9 7 3 2 7 8 6 2 1 2 3 4 3 Iso-value=5 2/27/2021 R. Crawfis, Ohio State Univ. 9

Dividing Cubes 1. Create a cell 2. Determine if the cell is a surface cell 3. Subdivide surface cells to image resolution 4. Determine intensity or color of refined cell 5. Output a point for each remaining surface cell 2/27/2021 R. Crawfis, Ohio State Univ. 10

Step 3: Subdivide z Subdivide each cell such that the resolution is higher than the image resolution. For example, when a 128 2 gridded field is rendered to a 5122 image, cells containing the contour will be divided to 4 x 4 smaller squares. classify 2/27/2021 subdivide R. Crawfis, Ohio State Univ. output points 11

Surface Tracking 0 z Given function S(t, x, y, z) that returns 1 1 if (x, y, z) is on the surface t and 0 otherwise 3 z Needed: a seed cell p that is on the surface. 2 1 1 3 2 3 6 6 3 7 9 7 3 7 8 6 2 2 3 4 3 1 Seed cell 2/27/2021 R. Crawfis, Ohio State Univ. Iso-value=5 12

Surface Tracking Let Q be a queue of cells. Push p onto Q Flag p as “visited” While Q is not empty do Pop q from Q and output it. For each cell v in the neighborhood of q do If S(t, v) = 1 then If v was not visited then Push v onto Q Flag it as “visited” end {while} 2/27/2021 R. Crawfis, Ohio State Univ. 13

Tracking and Smooth Curves z. Recall that the gradient to a scalar field is normal to the isolines. z. Therefore, the curve will be towards the direction perpendicular to the gradient. z. Can also use the gradient to calculate tangents for each point on the curve to generate smooth curves. 2/27/2021 R. Crawfis, Ohio State Univ. 14

Restricted Topological Cases • A closed curve can not be contained in the cell. • The curve can only enter and exit the cell once per edge. Filled circle vertices indicate scalar value is above the contour value Unfilled ones which scalar values are below the contour values 2/27/2021 R. Crawfis, Ohio State Univ. 15

Not Allowed 2/27/2021 R. Crawfis, Ohio State Univ. 16

Linear Topology z The restricted set is topologically equivalent to linear segments. z Topologically equivalent if we restrict the edge intersections to lie at the mid-points of each edge. 2/27/2021 R. Crawfis, Ohio State Univ. 17

Using Symmetry to Reduce to 5 Contour ambiguity: either solid or dash lines look OK 2/27/2021 R. Crawfis, Ohio State Univ. 18

2 D Iso-contour (0) Remember bi-linear interpolation p 2 p 3 P =? p 4 p 0 2/27/2021 p 5 To know the value of P, we can first compute p 4 and P 5 and then linearly interpolate P p 1 R. Crawfis, Ohio State Univ. 19

2 D Iso-contour (1) Consider a simple case: one cell data set The problem of extracting an iso-contour is an inverse of value interpolation. That is: p 2 p 3 Given f(p 0)=v 0, f(p 1)=v 1, f(p 2)=v 2, f(p 3)=v 3 Find the point(s) P within the cell that have values F(p) = C p 0 2/27/2021 p 1 R. Crawfis, Ohio State Univ. 20

2 D Iso-contour (2) We can solve the problem based on linear interpolation p 2 p 3 (1) Identify edges that contain points P that have value f(P) = C (2) Calculate the positions of P p 0 2/27/2021 p 1 (3) Connect the points with lines R. Crawfis, Ohio State Univ. 21

2 D Iso-contouring – Step 1 (1) Identify edges that contain points P that have value f(P) = C If v 1 < C < v 2 then the edge contains such a point v 1 2/27/2021 v 2 R. Crawfis, Ohio State Univ. 22

2 D Iso-contouring – Step 2 (2) Calculate the position of P p 1 P p 2 v 1 C v 2 2/27/2021 Use linear interpolation: P = P 1 + (C-v 1)/(v 2 -v 1) * (P 2 – P 1) R. Crawfis, Ohio State Univ. 23

2 D Iso-contouring – Step 3 p 2 p 0 2/27/2021 p 3 p 1 Connect the points with line(s) Based on the principle of linear variation, all the points on the line have values equal C R. Crawfis, Ohio State Univ. 24

Inside or Outside? Just a naming convention 1. If a value is smaller than the iso-value, we call it “Outside” 2. If a value is greater than the iso-value, we call it “Inside” - + p 2 p 3 p 0 p 1 inside cell 2/27/2021 - p 2 p 3 p 0 p 1 outside cell R. Crawfis, Ohio State Univ. 25

Contouring in 3 D z. Treat volume as a set of 2 D slices y. Apply 2 D Contouring algorithm on each slice. y. Or given as a set of hand-drawn contours z. Stitch the slices together. 2/27/2021 R. Crawfis, Ohio State Univ. 26

Contour Stitching • Problem: Given: 2 two-dimensional closed curves Curve #1 has m points Curve #2 has n points Which point(s) does vertex i on curve one correspond to on curve two? 2/27/2021 R. Crawfis, Ohio State Univ. i ? ? 27

A Solution z. Fuchs, et. al. y. Optimization problem y 1 stitch consists of: x 2 spans between curves x 1 contour segment Qj Pi Pi+1 y. Triangles of {Pi, Qj, Pi+1} or {Qj+1, Pi, Qj} x. Consistent normal directions 2/27/2021 R. Crawfis, Ohio State Univ. 28

Fuchs, et. al. z. Left span y. Pi. Qj => go up z. Right span (either) y. Pi+1 Qj => go down y. Pi. Qj+1 => go down 2/27/2021 R. Crawfis, Ohio State Univ. 29

Fuchs, et. al. z. Constraints y. Each contour segment is used once and only once. y. If a span appears as a left span, then it must also appear as a right span. y. If a span appears as a right span, then it must also appear as a left span. 2/27/2021 R. Crawfis, Ohio State Univ. 30

Fuchs, et. al. z. This produces an acceptable surface (from a topological point of view) y. No holes z. We would like an optimal one in some sense. 2/27/2021 R. Crawfis, Ohio State Univ. 31

Fuchs et. al. z. Graph problem z. Vertices Vij = span between Pi and Qj z. Edges are constructed from a left span to a right span. Qj y. Only two valid right spans for a left span. Pi 2/27/2021 R. Crawfis, Ohio State Univ. Qj+1 Pi+1 32

Fuchs, et. al. z. Organize these edges as a grid or matrix. j Q i P Pi. Qj. Pi+1 2/27/2021 R. Crawfis, Ohio State Univ. Qj. Pi. Qj+1 33

Fuchs, et. al. z. Acceptable graphs y. Exactly one vertical arc between Pi and Pi+1 y. Exactly one horiz. arc between Qj and Qj+1 y. Either xindegree(Vij) = outdegree(Vij) = 0 xboth > 0 • (if a right, also has to be a left) 2/27/2021 R. Crawfis, Ohio State Univ. 34

Fuchs, et. al. z. Claim: y. An acceptable graph, S, is weakly connected. z. Lemma 2 y. Only 0 or 1 vertex of S has an indegree = 2. x. E. g. , Two cones touching in the center. y. All other vertices have indegree=1 y(recall indegree = outdegree) 2/27/2021 R. Crawfis, Ohio State Univ. 35

Fuchs, et. al. z. Thereom 1: S is an acceptable surface if and only if: y. S has one and only one horizontal arc between adjacent columns. y. S has one and only one vertical arc between adjacent rows. y. S is Eulerian (closed walk with every arc only once). 2/27/2021 R. Crawfis, Ohio State Univ. 36

Fuchs, et. al. z. Number of arcs is thus m+n. z. Many possible solutions still!!! z. Associate costs with each edge y. Area of resulting triangle y. Aspect ratio of resulting triangle 2/27/2021 R. Crawfis, Ohio State Univ. 37

Fuchs, et. al. z. Note that Vi 0 is in S for some i. z. Note that V 0 j is in S for some j. z. With the weights (costs), we can compute the minimum path from a starting node Vi 0. y. Since we do not know which Vi 0, we compute the paths for all of them and take the minimum. x. Some cost saving are achievable. 2/27/2021 R. Crawfis, Ohio State Univ. 38

Marching Cubes z. Predominant method used today. z. Efficient and simple z. It was independently reported by Wyvill and Mc. Peeters in 1986, Lorenson and Cline in 1987 2/27/2021 R. Crawfis, Ohio State Univ. 39

Marching Cubes z. Treat each cube individual y. No 2 D contour curves z. Allow intersections only on the edges or at vertices. z. Pre-calculate all of the necessary information to construct a surface. 2/27/2021 R. Crawfis, Ohio State Univ. 40

Marching Cubes z. Consider a single cube y. All vertices above the contour threshold y. All vertices below y. Mixed above and below + + - + + 2/27/2021 R. Crawfis, Ohio State Univ. + 41

Marching Cubes z Binary label each node => (above/below) z Examine all possible cases of above or below for each vertex. z 8 vertices implies 256 possible cases. + R. Crawfis, Ohio State Univ. - + + 2/27/2021 + + 42

How many cases? Now we have 8 vertices So it is: 2 8= 256 How many unique topological cases? 2/27/2021 R. Crawfis, Ohio State Univ. 43

Case Reduction Value Symmetry _ _ _ + _ 2/27/2021 _ _ R. Crawfis, Ohio State Univ. + + _ _ + + + 44

Case Reduction Rotation Symmetry _ _ _ + 2/27/2021 _ _ _ + R. Crawfis, Ohio State Univ. 45

Case Reduction Mirror Symmetry _ + _ _ _ + + + _ _ + + By inspection, we can reduce 256 2/27/2021 _ + R. Crawfis, Ohio State Univ. 15 46

Marching Cube - Summary z. Create a cube z. Classify each voxel z. Build an index z. Lookup edge list z. Interpolate triangle vertices z. Calculate and interpolate normals 2/27/2021 R. Crawfis, Ohio State Univ. 47

Step 1: Create a Cube z. Consider a cube defined by eight data values: four from slice K, and four from slice K+ 1 2/27/2021 R. Crawfis, Ohio State Univ. 48

Step 2: Classify Each Voxel z. Binary classification of each vertex of the cube as to whether it lies youtside the surface (voxel value < isosurface value) yor inside the surface (voxel value <= isosurface value). 2/27/2021 R. Crawfis, Ohio State Univ. 49

Step 3: Build an Index z. Use the binary labeling of each voxel to create an 8 -bit index: y 8 vertices - 256 cases y. Requires a canonical ordering of the cube. 2/27/2021 R. Crawfis, Ohio State Univ. 50

Step 4: Look-up the Topology z Given the index for each cell, a table lookup is performed to identify the edges that intersect the isosurface. Index intersection edges z Also requires a canonical ordering 0 e 1, e 3, e 5 of the edges: e 1, e 2, …, e 12 … 1 2 e 7 e 11 e 8 e 3 e 4 e 9 3 e 12 e 6 e 5 e 2 14 e 10 e 1 2/27/2021 R. Crawfis, Ohio State Univ. 51

Step 4: Look-up the Topology z all 256 cases could be derived from 15 base cases. z Difficult to determine rotational symmetry and map it correctly to the table. 2/27/2021 R. Crawfis, Ohio State Univ. 52

Step 4: Example z Index = 10001101 z triangle 1 = e 4, e 7, e 12 z triangle 2 = e 1, e 7, e 4 z triangle 3 = e 1, e 6, e 7 z triangle 4 = e 1, e 10, e 6 e 7 e 12 v 7 e 6 e 4 v 1 2/27/2021 R. Crawfis, Ohio State Univ. e 1 v 2 e 10 e 2 53

Step 5: Interpolate Triangle Vertices z For each edge, find the vertex location along the edge using linear interpolation of the voxel values. 2/27/2021 R. Crawfis, Ohio State Univ. 54

Step 6: Compute Normals z. Calculate the normal at each cube vertex z. Use linear interpolation to interpolate the polygon vertex normal 2/27/2021 R. Crawfis, Ohio State Univ. 55

Why is it called marching cubes? Linear search through cells • Row by row, layer by layer • Reuse the interpolated points for adjacent cells No neighborhood information is required!!! 2/27/2021 R. Crawfis, Ohio State Univ. 56

Problems z. Ambiguities: Are the curves below right or wrong? 2/27/2021 R. Crawfis, Ohio State Univ. 57

Problems z. Inconsistencies: Does the result produce the correct topology? z. An iso-contour surface should always be a manifold (with boundary). y. We know whether these cases are right or wrong 2/27/2021 R. Crawfis, Ohio State Univ. 58

2 D Contouring Ambiguity Either choice is acceptable because contour lines are continuous and closed or Independent of other choices Break contour 2/27/2021 R. Crawfis, Ohio State Univ. Join contour 59

Marching Cubes z. Topological inconsistencies in the 15 cases y. Turns out positive and negative are not symmetric. + + + 2/27/2021 + - - + + + R. Crawfis, Ohio State Univ. - - - 60

Topological Inconsistencies z. The problem is that the edge intersections match, but the face intersections do not. z. In 2 D we only had edge intersections, so this was not a problem. z. In 4 D and higher the problem gets even worse. 2/27/2021 R. Crawfis, Ohio State Univ. 61

A Bad Example Resolution: Favor positive pathways thru the voxel (or negative, but be consistent). Case 3 2/27/2021 R. Crawfis, Ohio State Univ. Case 3 Case 6 (compliment symmetry) 62

Ambiguous Cases z Ambiguous cases: 3, 4, 6, 7, 10, 12, 13 z Adjacent vertices: different states z Diagonal vertices: same state z Case 4 is not usually considered ambiguous, but we could favor a tunnel through the long diagonal. 2/27/2021 R. Crawfis, Ohio State Univ. or or 63

Ambiguous Cases z The ambiguous cases lead to a continuous manifold, but due to improper sampling of the volume, a (somewhat arbitrary) decision can be made on whether to close two curves or connect them. z Asymptotic Decider - by Nielson and Hamann (IEEE Vis’ 91) examines neighboring cells to determine this choice. 2/27/2021 R. Crawfis, Ohio State Univ. 64

Asymptotic Decider Based on bilinear interpolation over faces B 11 B 01 (s, t) B 00 B(s, t) = (1 -s, s) B 10 B 01 B 10 B 11 1 -t t The contour curves of B: {(s, t) | B(s, t) = } are hyperbolas 2/27/2021 R. Crawfis, Ohio State Univ. 65

Asymptotic Decider (2) (1, 1) Where the hyperbolas go through the cell depends on the values at the corners, I. e. , B 00, B 01, B 10, B 11 (0, 0) 2/27/2021 R. Crawfis, Ohio State Univ. 66

Asymptotic Decider (3) (1, 1) (S , T ) If < B(S , T ) (0, 0) Asymptote 2/27/2021 R. Crawfis, Ohio State Univ. 67

Asymptotic Decider (4) (1, 1) (S , T ) If > B(S , T ) (0, 0) Asymptote 2/27/2021 R. Crawfis, Ohio State Univ. 68

Asymptotic Decider (5) (1, 1) (S , T ) S = B 00 - B 01 B 00 + B 11 – B 01 – B 10 T = B 00 – B 10 B 00 + B 11 – B 01 – B 10 (0, 0) B(S , T ) = B 00 B 11 + B 10 B 01 B 00 + B 11 – B 01 – B 10 2/27/2021 R. Crawfis, Ohio State Univ. 69

Asymptotic Decider (6) Based on the result of asymptotic decider, we expand the marching cube case 3, 6, 12, 10, 7, 13 (These are the cases with at least one ambiguous faces). 2/27/2021 R. Crawfis, Ohio State Univ. 70

Marching Cubes Demo z Animating the contour value z Special functions for contouring z Varying speeds and numbers of triangles 2/27/2021 R. Crawfis, Ohio State Univ. 71

Marching Cubes z. Data Structures/Tables static int const Hexa. Edges[12][2] = { {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6}, {6, 7}, {7, 4}, {0, 4}, {1, 5}, {3, 7}, {2, 6}}; typedef struct { EDGE_LIST Hexa. Edges[16]; } HEXA_TRIANGLE_CASES; /* Edges to intersect. static const {-1, -1, { 0, 8, 3, { 0, 1, 9, { 1, 8, 3, Three at a time form a triangle. */ HEXA_TRIANGLE_CASES -1, -1, -1, -1, 9, 8, 1, -1, Hexa. Tri. Cases[] = { -1, -1, -1, -1, -1, -1, -1, -1}, /* /* 0 1 2 3 */ */ . . . 2/27/2021 R. Crawfis, Ohio State Univ. 72

Marching Cubes - How simple /* Determine the marching cubes index */ for ( i=0, index = 0; i < 8; i++) if (val 1[nodes[i]] >= thresh) /* If the nodal value is above the */ index |= CASE_MASK[i]; /* threshold, set the appropriate bit. */ tri. Case = Hexa. Tri. Cases[index]; /* tri. Case indexes into the MC table. */ edge = tri. Case->Hexa. Edges; /* edge points to the list of intersected edges */ } for ( ; edge[0] > -1; edge += 3 ) /* stop if we hit the -1 flag */ { for (i=0; i<3; i++) /* Calculate and store three edge intersections */ { vert = Hexa. Edges[edge[i]]; n 0 = nodes[vert[0]]; n 1 = nodes[vert[1]]; t = (thresh - val 1[n 0]) / (val 1[n 1] - val 1[n 0]); tri_ptr[i] = add_intersection( n 0, n 1, t ); /* Save an index to the pt. */ } add_triangle( tri_ptr[0], tri_ptr[1], tri_ptr[2], zone. ID); /* Store the triangle */ } 2/27/2021 R. Crawfis, Ohio State Univ. 73

Iso-Contouring – Advanced Issues 1. Efficiently determining which cells to examine. 2. Using iso-contouring as a slicing mechanism 3. Iso-contouring in higher dimensions 4. Texturing and coloring of iso-contours 5. Polygonal simplification of contours. 6. Choosing a good iso-value 2/27/2021 R. Crawfis, Ohio State Univ. 74

Iso-surface cell search z. Iso-surface cells: cells that contain isosurface. min < iso-value < max z. Marching cubes algorithm performs a linear search to locate the iso-surface cells – not very efficient for large-scale data sets. 2/27/2021 R. Crawfis, Ohio State Univ. 75

Iso-surface Cells z For a given iso-value, only a smaller portion of the cells contain part of the iso-surface. z For a volume with n x n cells, the n average number of the iso-surface cells is O(n x n) (ratio of surface v. s. volume) 2/27/2021 R. Crawfis, Ohio State Univ. n n 76

Efficient Searching z. With < 10% of the voxels contributing to the surface, it is a waste to look at every voxel. z. A voxel can be specified in terms of its interval, its minimum and maximum values. 2/27/2021 R. Crawfis, Ohio State Univ. 77

Efficient iso-surface cell search z. Problem statement: y. Given a scalar field with N cells: xc 1, c 2, …, cn y. With min-max range: x(a 1, b 1), (a 2, b 2), …, (an, bn) Find {Ck | ak < C < bk; C=iso-value} 2/27/2021 R. Crawfis, Ohio State Univ. 78

Efficient search methods 1. Spatial subdivision (domain search) 2. Value subdivision (range search) 3. Contour propagation 2/27/2021 R. Crawfis, Ohio State Univ. 79

Domain search • • Subdivide the space into several sub-domains, check the min/max values for each sub-domain If the min/max values (extreme values) do not contain the iso-value, we skip the entire region Min/max Complexity = O( k log(n/k) ) 2/27/2021 R. Crawfis, Ohio State Univ. 80

Span Space High Gradient Low Gradient Threshold 2/27/2021 Minimum R. Crawfis, Ohio State Univ. 81

Span Space max z Given: iso-value v z Candidate cells: y Minimum cell value less than the iso-value. y Maximum cell value greater than the iso-value. z Different algorithms query this span-space with different efficiencies. Possible candidates v v 2/27/2021 R. Crawfis, Ohio State Univ. min 82

Span Space max z Create two lists: y Store (in order) the minima of each cell into the first list y Store (in order) the maxima of each cell into the second list. z Can easily get the set of cells with minimum less than the iso-value. v 2/27/2021 R. Crawfis, Ohio State Univ. min 83

Span Space z Can easily get the set of cells with maximum greater than the iso-value. z Need the intersection of the two sets. z Not very efficient. 2/27/2021 max v R. Crawfis, Ohio State Univ. v min 84

Span Space - Active List z Given: iso-value v and maximum cell range, r z gracefully update to new value nv max nv+r nv v v-r 2/27/2021 R. Crawfis, Ohio State Univ. nv-r v nv min 85

Span Space – NOISE Algorithm • K-d Trees Split Min. axis Split Max. axis Minimum 2/27/2021 R. Crawfis, Ohio State Univ. 86

Span Space Search Method (2): ISSUE, discretized span space O(log(N/L)) O(1) Complexity = ? ? O(log(N/L)) 2/27/2021 R. Crawfis, Ohio State Univ. 87

Range Search: Interval Tree: Id I left … 2/27/2021 I right … Sort all the data points (x 1, x 2, x 3, x 4, …. , xn) Let d = x n/2 (mid point) We use d to divide the cells into three sets Id, I left, and I right Id : cells that have I left: cells that have I right: cells that have R. Crawfis, Ohio State Univ. min < d max < d min > < max d 88

Contour Propagation Basic Idea: Given an initial cell that contains iso-surface, the remainder of the iso-surface can be found by propagation Initial cell: A C E A B 2/27/2021 D FIFO Queue A Enqueue: B, C BC Dequeue: B C Enqueue: D CD … …. R. Crawfis, Ohio State Univ. Breadth-First Search 89

Challenges Need to know the initial cells! For any given iso-value C, finding the initial cells to start the propagation is almost as hard as finding the iso-surface cells. You could do a global search, but … 2/27/2021 R. Crawfis, Ohio State Univ. 90

Solutions (1) Extrema Graph (Itoh vis’ 95) (2) Seed Sets (Bajaj volvis’ 96) Problem Statement: Given a scalar field with a cell set G, find a subset S G, such that for any given iso-value C, the set S contains initial cells to start the propagation. We need search through S, but S is usually (hopefully) much smaller than G. We will only talk about extrema graph due to time constraint 2/27/2021 R. Crawfis, Ohio State Univ. 91

Extrema Graph (1) Basic Idea: If we find all the local minimum and maximum points (Extrema), and connect them together by straight lines (Arcs), then any closed Iso-contour is intersected by at least one of the arcs. 2/27/2021 R. Crawfis, Ohio State Univ. 92

Extrema Graph (2) 2/27/2021 R. Crawfis, Ohio State Univ. 93

Extrema Graph (3) E 2 E 1 a 2 a 3 E 4 a 5 a 4 E 5 2/27/2021 Extreme Graph: E 7 E 6 a 7 { E, A: E: extrema points A: Arcs conneccts E } An ‘arc’ consists of cells that connect extrema points (we only store min/max of the arc though) E 8 R. Crawfis, Ohio State Univ. 94

Extrema Graph (4) Algorithm: Given an iso-value 1) Search the arcs of the extrema graph (to find the arcs that have min/max values that contain the iso-value 2) Walk through the cells along each of these arcs to find the seed cells 3) Start the iso-contour propagation from the seed cells 4) …. There is something more that needs to be done… 2/27/2021 R. Crawfis, Ohio State Univ. 95

We are not done yet … What ? ! We just mentioned that all the closed iso-contours will intersect with the arcs connecting the extrema points How about non-closed iso-contours? (or called open isocontours) 2/27/2021 R. Crawfis, Ohio State Univ. 96

Extrema Graph (5) Contours missed These open iso-contours will intersect with ? ? cells Boundary Cells!! 2/27/2021 R. Crawfis, Ohio State Univ. 97

Extrema Graph (6) Algorithm (continued) Given an iso-value: 1) Search the arcs of the extrema graph (to find the arcs that have min/max values that contain the iso-value. 2) Walk through the cells along each of the arcs to find the seed cells. 3) Start the iso-contour propagation from the seed cells. 4) Search the cells along the boundary and find additional seed cells. 5) Propagate from these new seed cells for the open iso-contours. 2/27/2021 R. Crawfis, Ohio State Univ. 98

Extrema Graph z Efficiency - Number of cells visited: y extrema graph - N 0. 33 y boundary - N 0. 66 y Iso-surface - N 0. 66 z based on tetrahedra - will create more surface triangles. . . z should extract the same number of cells/ triangulation as Marching Cubes 2/27/2021 R. Crawfis, Ohio State Univ. 99

Extrema Graph z. Storage Costs: y. Extrema graph – very small amount of memory for most datasets. y. For unstructured grids, I need to be able to access my neighbors. This requires a very expensive data structure, more than the span -space data structures. 2/27/2021 R. Crawfis, Ohio State Univ. 100

Selecting iso-values z. Which iso-values should I examine to best comprehend my dataset? z. Some data has very specific values: y. The aluminum structure in my simulation will start to fail at a pressure of 5672 psi. y. The molecular state changes from a liquid to a gas at 100° C. 2/27/2021 R. Crawfis, Ohio State Univ. 101

Contour Spectrum z Basic Idea: Calculate and present to the user several properties of an iso-contour. z Do this for all iso-contours. z This leads to several functions in terms of the iso-value, . z Present these functions to the user as an aid in picking contour values. These slides are from Baja, Pascucci and Schikore (IEEE Visualization 1997) 2/27/2021 R. Crawfis, Ohio State Univ. 102

Outline z. User Interface z Signature Computation z Real Time Quantitative Queries z Rule-based Contouring z Topological Information z Future Directions 2/27/2021 R. Crawfis, Ohio State Univ. 103

Graphical User Interface for Static Data • • The horizontal axis spans the scalar values Plot of a set of signatures (length, area, gradient. . . ) as functions of the scalar value . Vertical axis spans normalized ranges of each signature. White vertical bars mark current selected isovalues. 2/27/2021 R. Crawfis, Ohio State Univ. 104

Graphical User Interface for time varying data high ( , t ) --> c t • • The horizontal axis spans the scalar value dimension . The vertical axis spans the time dimension t. 2/27/2021 R. Crawfis, Ohio State Univ. The color, c, is mapped to the magnitude of a signature function of time, t, and isovalue low magnitude c 105

User Interface - MRI of a human torso • • In real time the exact value of each signature is displayed. The isocontour that bounds the region of interest is obtained by selecting the maximum of the gradient signature. 2/27/2021 R. Crawfis, Ohio State Univ. 106

Signature Computation z Consider a terrain of which you want to compute the length of each isocontour and the area contained inside each isocontour 2/27/2021 R. Crawfis, Ohio State Univ. 107

Signature Computation • The length of each contour is a C 0 spline function. • The area inside/outside each isocontour is a C 1 spline function. 2/27/2021 R. Crawfis, Ohio State Univ. 108

Signature Computation z In general, the size (surface area) of each iso-contour of a scalar field of dimension d is a spline function of d-2 continuity. z The size (volume) of the region inside/ outside is given by a spline function of d-1 continuity 2/27/2021 R. Crawfis, Ohio State Univ. 109

Real Time Quantitative Queries (agricultural yield data) y size and position of the region with unsatisfactory production y size and position of the region where wrong data acquisition occurred Spectrum space is a useful space to visualize by itself. 2/27/2021 R. Crawfis, Ohio State Univ. 110

Rule-based Contouring (CT scan of an engine) z 2/27/2021 R. Crawfis, Ohio State Univ. The contour spectrum allows the development of an adaptive ability to separate interesting isovalues from the others. 111

Rule-based Contouring (foot of the Visible Human) z 2/27/2021 R. Crawfis, Ohio State Univ. The contour spectrum allows the development of an adaptive ability to separate interesting isovalues from the others. 112

Topological information. • number of components per isocontour • which isocontours merge together or split while modifying the isovalue. an isocontour with three connected components two of which are about to merge 2/27/2021 } R. Crawfis, Ohio State Univ. 113

Other Topics z. Interval Volumes z 4 D Iso-contouring z. Smooth surfaces z. Polygon Reduction z. Triangle Strip generation z. Time-varying iso-contours z. Texture map parameterization 2/27/2021 R. Crawfis, Ohio State Univ. Penny Rheingans 114

Level Sets Roger Crawfis Slides collected from: Fan Ding, Charles Dyer, Donald Tanguay, Theodore Kim and Roger Crawfis R. Crawfis, Ohio State Univ.

Motivation: 2/27/2021 R. Crawfis, Ohio State Univ. 116

Easy Case – Use Marching Cubes 2/27/2021 R. Crawfis, Ohio State Univ. 117

Input Data Noisy 2/27/2021 R. Crawfis, Ohio State Univ. 118

Non-uniform Exposure 2/27/2021 R. Crawfis, Ohio State Univ. 119

Intensity Varies 2/27/2021 R. Crawfis, Ohio State Univ. 120

Multiple problems 2/27/2021 R. Crawfis, Ohio State Univ. 121

Applying Marching Cubes (thresholding) 2/27/2021 R. Crawfis, Ohio State Univ. 122

Applying A Threshold 2/27/2021 R. Crawfis, Ohio State Univ. 123

Applying A Threshold 2/27/2021 R. Crawfis, Ohio State Univ. 124

Four (contour) Levels 2/27/2021 R. Crawfis, Ohio State Univ. 125

What To Do? z User Intervention!!! y We see it!! It’s right there!!! y Well, the edges get pretty fuzzy . z Two step process: 1. Draw an initial curve (or surface) within the desired region. 2. Expand that interface outward towards the edge of our desired region. 2/27/2021 R. Crawfis, Ohio State Univ. 126

Interfaces An interface (or front) is a boundary between two regions: “inside” and “outside. ” In 2 -D, an interface is a simple closed curve: 2/27/2021 R. Crawfis, Ohio State Univ. 127

Propagating Interfaces How does an interface evolve over time? At a specific moment, the speed function F (L, G, I) describes the motion of the interface in the normal direction. 2/27/2021 R. Crawfis, Ohio State Univ. 128

Propagating Interfaces Speed F(L, G, I) is a function of 3 types of properties: • Local – depend on local geometric information (e. g. , curvature and normal direction) • Global – depend on the shape and position of the front (e. g. , integrals along the front, heat diffusion) • Independent – do not depend on the shape of the front (e. g. , an underlying fluid velocity that passively transports the front) 2/27/2021 R. Crawfis, Ohio State Univ. 129

Motion Under Curvature Example: Motion by curvature. Each piece moves perpendicular to the curve with speed proportional to the local curvature. small negative motion large positive motion 2/27/2021 R. Crawfis, Ohio State Univ. 130

Motion Under Curvature κ is the inverse of the radius r of the osculating circle. 2/27/2021 R. Crawfis, Ohio State Univ. 131

Motion Under Curvature http: //math. berkeley. edu/~sethian/Applets/java_curve_flow. html 2/27/2021 R. Crawfis, Ohio State Univ. 132

Functional Representation Eulerian framework: define fixed coordinate system on the world. For every world point x, there is (at most) one value y = ft(x). Falling snow example: 2/27/2021 R. Crawfis, Ohio State Univ. 133

Functional Representation However, many simple shapes are multivalued; they are not functions regardless of the orientation of the coordinate system. 2/27/2021 R. Crawfis, Ohio State Univ. 134

Parametric Representation Spatially parameterize the curve x by s so that at time t the curve is xt(s), where 0 s S and the curve is closed: xt(0) = xt(S). Points on initial curve. Gradient (wrt time) is the speed in normal direction. Normal is perpendicular to curve, as is curvature. 2/27/2021 R. Crawfis, Ohio State Univ. 135

Parametric Representation For motion under curvature, speed F depends only on local curvature κ – the equation of motion is thus: where curvature is and the normal is 2/27/2021 R. Crawfis, Ohio State Univ. 136

Particle Methods In order to compute, discretize the parameterization into moving particles which reconstruct the front. Known under a variety of names: marker particle techniques, string methods, nodal methods. = time step = parameterization step = # mesh particles Δs = location of point iΔs at time nΔt 2/27/2021 R. Crawfis, Ohio State Univ. 137

Problem Statement z. Generally, given: y. An initial front y. Equations that govern its evolution z. How do we simulate the front’s evolution? z. Called an ‘initial value problem’ y. Given the initial position y. Solve for a position at a future time 2/27/2021 R. Crawfis, Ohio State Univ. 138

More formally: z. Given some initial front Г: 2/27/2021 R. Crawfis, Ohio State Univ. 139

More formally: (2) z. And a function F that specifies the velocity of the front in the normal direction: 2/27/2021 R. Crawfis, Ohio State Univ. 140

More formally: (3) z. Solve for Г at some future time z. Level set methods are used to track an interface z. Water/air interface, for example 2/27/2021 R. Crawfis, Ohio State Univ. 141

What’s Wrong with the Obvious Solution? z. Why is a level set method necessary? z. There seems to be a more intuitive way to solve this problem 2/27/2021 R. Crawfis, Ohio State Univ. 142

Marker/String Methods z. Why not just connect some control points (in 3 D, triangulate): 2/27/2021 R. Crawfis, Ohio State Univ. 143

Marker/String Methods (2) z. And run the simulation on the points? 2/27/2021 R. Crawfis, Ohio State Univ. 144

Ocean Waves z. Think of an air/water interface with two waves racing towards each other: 2/27/2021 R. Crawfis, Ohio State Univ. 145

Ocean Waves (2) z. What happens to the control points when the waves collide? 2/27/2021 R. Crawfis, Ohio State Univ. 146

Shocks z. Event known as a ‘shock’ z. Below formation called a ‘swallowtail’ 2/27/2021 R. Crawfis, Ohio State Univ. 147

Shocks (2) z. How to fix up the control points? z. Fixing swallowtails known as ‘de-looping’ y. Very difficult y. Some methods exist in 2 D y. No robust 3 D methods so far 2/27/2021 R. Crawfis, Ohio State Univ. 148

Changing Topology Example: two fires merge into a single fire. 2/27/2021 R. Crawfis, Ohio State Univ. 149

Changing Topology In particle methods: • Difficult (and expensive) to detect and change the particle chains • Much more difficult as dimensionality increases 2/27/2021 Buoys! R. Crawfis, Ohio State Univ. 150

Difficulties With Particle Methods z. Instability z. Local singularities z. Management of particles: remove, redistribute, connect 2/27/2021 R. Crawfis, Ohio State Univ. 151

Level Set Formulation Recast problem with one additional dimension – the distance from the interface. And then use Marching Cubes to extract the surface. 2/27/2021 R. Crawfis, Ohio State Univ. 152

Level Set Methods z Contour evolution method due to J. Sethian and S. Osher, 1988 z www. math. berkeley. edu/~sethian/level_set. html z Difficulties with snake-type methods y. Hard to keep track of contour if it self-intersects during its evolution y. Hard to deal with changes in topology 2/27/2021 R. Crawfis, Ohio State Univ. 153

z The level set approach: y. Define problem in 1 higher dimension y. Define level set function z = (x, y, t=0) where the (x, y) plane contains the contour, and z = signed Euclidean distance transform value (negative means inside closed contour, positive means outside contour) 2/27/2021 R. Crawfis, Ohio State Univ. 154

How to Move the Contour? z Move the level set function, (x, y, t), so that it rises, falls, expands, etc. z Contour = cross section at z = 0 2/27/2021 R. Crawfis, Ohio State Univ. 155

Level Set Surface z The zero level set (in blue) at one point in time is a slice of the level set surface (in red) 2/27/2021 R. Crawfis, Ohio State Univ. 156

Level Set Surface z Later in time the level set surface (red) has moved and the new zero level set (blue) defines the new contour 2/27/2021 R. Crawfis, Ohio State Univ. 157

Level Set Surface 2/27/2021 R. Crawfis, Ohio State Univ. 158

Level Set Formulation The interface always lies at the zeroth level set of the function , i. e. , the interface is defined by the implicit equation t (x, y) = 0. 2/27/2021 R. Crawfis, Ohio State Univ. 159

How to Move the Level Set Surface? 1. Define a velocity field, F, that specifies how contour points move in time y Based on application-specific physics such as time, position, normal, curvature, image gradient magnitude 2. Build an initial value for the level set function, (x, y, t=0), based on the initial contour position 3. Adjust over time; current contour defined by (x(t), y(t), t) = 0 2/27/2021 R. Crawfis, Ohio State Univ. 160

Speed Function 2/27/2021 R. Crawfis, Ohio State Univ. 161

Example: Shape Simplification z F = 1 – 0. 1 where is the curvature at each contour point 2/27/2021 R. Crawfis, Ohio State Univ. 162

Example: Segmentation z Digital Subtraction Angiogram z F based on image gradient and contour curvature 2/27/2021 R. Crawfis, Ohio State Univ. 163

Example (cont. ) z Initial contour specified manually 2/27/2021 R. Crawfis, Ohio State Univ. 164

Result – segmentation using Fast marching 2/27/2021 No level set tuning R. Crawfis, Ohio State Univ. 165

Results – vein segmentation No level set tuning 2/27/2021 With level set tuning R. Crawfis, Ohio State Univ. 166

Results – vein segmentation continued Original 2/27/2021 Our result (Fast marching + Level set tuning) R. Crawfis, Ohio State Univ. Sethian’s result (Level set only) 167

Result – segmentation using Fast marching No level set tuning 2/27/2021 R. Crawfis, Ohio State Univ. 168

Results – brain segmentation continued No level set tuning 2/27/2021 With level set tuning R. Crawfis, Ohio State Univ. 169

Results – brain image segmentation # of iterations = 9000 # of iterations = 12000 Fast marching only, no level set tuning 2/27/2021 R. Crawfis, Ohio State Univ. 170