vtk-pure-rs 0.2.0

Pure Rust visualization toolkit — data structures, filters, I/O, rendering
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201

use crate::data::{AnyDataArray, DataArray, ImageData, PolyData};

/// Voxelize a triangle mesh into an ImageData volume.
///
/// Voxels inside the mesh are set to 1.0 and outside to 0.0 in the
/// "Voxels" scalar array. Uses ray-casting along the Z axis with parity
/// counting to determine inside/outside.
///
/// `dims` specifies the output grid resolution `[nx, ny, nz]`.
/// `bounds` specifies the spatial extent `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`.
pub fn voxelize(input: &PolyData, dims: [u32; 3], bounds: [[f64; 2]; 3]) -> ImageData {
    let nx: usize = dims[0] as usize;
    let ny: usize = dims[1] as usize;
    let nz: usize = dims[2] as usize;
    let n: usize = nx * ny * nz;

    let sx: f64 = if nx > 1 {
        (bounds[0][1] - bounds[0][0]) / (nx - 1) as f64
    } else {
        1.0
    };
    let sy: f64 = if ny > 1 {
        (bounds[1][1] - bounds[1][0]) / (ny - 1) as f64
    } else {
        1.0
    };
    let sz: f64 = if nz > 1 {
        (bounds[2][1] - bounds[2][0]) / (nz - 1) as f64
    } else {
        1.0
    };

    // Collect triangles
    let mut triangles: Vec<[[f64; 3]; 3]> = Vec::new();
    for cell in input.polys.iter() {
        if cell.len() < 3 {
            continue;
        }
        let v0 = input.points.get(cell[0] as usize);
        for i in 1..cell.len() - 1 {
            let v1 = input.points.get(cell[i] as usize);
            let v2 = input.points.get(cell[i + 1] as usize);
            triangles.push([v0, v1, v2]);
        }
    }

    let mut voxels: Vec<f64> = vec![0.0; n];

    // For each (i, j) column, cast a ray along Z and count crossings
    for j in 0..ny {
        let ry: f64 = bounds[1][0] + j as f64 * sy;
        for i in 0..nx {
            let rx: f64 = bounds[0][0] + i as f64 * sx;

            // Collect Z intersections of ray (rx, ry, z) with all triangles
            let mut z_hits: Vec<f64> = Vec::new();

            for tri in &triangles {
                if let Some(z) = ray_z_intersect(rx, ry, tri) {
                    z_hits.push(z);
                }
            }

            z_hits.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

            // Walk voxels along Z, toggle inside/outside at each crossing
            let mut hit_idx: usize = 0;
            let mut inside: bool = false;
            for k in 0..nz {
                let rz: f64 = bounds[2][0] + k as f64 * sz;

                // Advance past all crossings below this Z
                while hit_idx < z_hits.len() && z_hits[hit_idx] < rz - sz * 0.5 {
                    inside = !inside;
                    hit_idx += 1;
                }

                if inside {
                    voxels[k * ny * nx + j * nx + i] = 1.0;
                }
            }
        }
    }

    let mut output = ImageData::with_dimensions(nx, ny, nz);
    output.set_spacing([sx, sy, sz]);
    output.set_origin([bounds[0][0], bounds[1][0], bounds[2][0]]);
    output
        .point_data_mut()
        .add_array(AnyDataArray::F64(DataArray::from_vec("Voxels", voxels, 1)));
    output
}

/// Test if a Z-axis ray at (rx, ry) intersects a triangle, returning the Z coordinate.
fn ray_z_intersect(rx: f64, ry: f64, tri: &[[f64; 3]; 3]) -> Option<f64> {
    let ax: f64 = tri[0][0];
    let ay: f64 = tri[0][1];
    let bx: f64 = tri[1][0];
    let by: f64 = tri[1][1];
    let cx: f64 = tri[2][0];
    let cy: f64 = tri[2][1];

    // Barycentric coordinates in XY
    let denom: f64 = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy);
    if denom.abs() < 1e-15 {
        return None;
    }
    let inv_denom: f64 = 1.0 / denom;

    let u: f64 = ((by - cy) * (rx - cx) + (cx - bx) * (ry - cy)) * inv_denom;
    let v: f64 = ((cy - ay) * (rx - cx) + (ax - cx) * (ry - cy)) * inv_denom;
    let w: f64 = 1.0 - u - v;

    let eps: f64 = -1e-8;
    if u < eps || v < eps || w < eps {
        return None;
    }

    let z: f64 = u * tri[0][2] + v * tri[1][2] + w * tri[2][2];
    Some(z)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::data::DataSet;

    fn make_box_mesh() -> PolyData {
        // A simple axis-aligned box from (0,0,0) to (1,1,1)
        let pts = vec![
            [0.0, 0.0, 0.0], // 0
            [1.0, 0.0, 0.0], // 1
            [1.0, 1.0, 0.0], // 2
            [0.0, 1.0, 0.0], // 3
            [0.0, 0.0, 1.0], // 4
            [1.0, 0.0, 1.0], // 5
            [1.0, 1.0, 1.0], // 6
            [0.0, 1.0, 1.0], // 7
        ];
        // 12 triangles forming the 6 faces
        let tris = vec![
            [0, 2, 1], [0, 3, 2], // -Z face
            [4, 5, 6], [4, 6, 7], // +Z face
            [0, 1, 5], [0, 5, 4], // -Y face
            [2, 3, 7], [2, 7, 6], // +Y face
            [0, 4, 7], [0, 7, 3], // -X face
            [1, 2, 6], [1, 6, 5], // +X face
        ];
        PolyData::from_triangles(pts, tris)
    }

    #[test]
    fn voxelize_box_produces_correct_dimensions() {
        let mesh = make_box_mesh();
        let result = voxelize(
            &mesh,
            [5, 5, 5],
            [[-0.5, 1.5], [-0.5, 1.5], [-0.5, 1.5]],
        );
        assert_eq!(result.dimensions(), [5, 5, 5]);
        assert_eq!(result.num_points(), 125);
        let arr = result.point_data().get_array("Voxels").unwrap();
        assert_eq!(arr.num_tuples(), 125);
    }

    #[test]
    fn voxelize_box_has_inside_voxels() {
        let mesh = make_box_mesh();
        let result = voxelize(
            &mesh,
            [5, 5, 5],
            [[-0.5, 1.5], [-0.5, 1.5], [-0.5, 1.5]],
        );
        let arr = result.point_data().get_array("Voxels").unwrap();
        let mut buf = [0.0f64];

        // Count how many voxels are inside
        let mut inside_count: usize = 0;
        for i in 0..arr.num_tuples() {
            arr.tuple_as_f64(i, &mut buf);
            if buf[0] > 0.5 {
                inside_count += 1;
            }
        }
        // Should have some inside voxels but not all
        assert!(inside_count > 0);
        assert!(inside_count < 125);
    }

    #[test]
    fn empty_mesh_gives_all_outside() {
        let mesh = PolyData::new();
        let result = voxelize(&mesh, [3, 3, 3], [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]);
        let arr = result.point_data().get_array("Voxels").unwrap();
        let mut buf = [0.0f64];
        for i in 0..arr.num_tuples() {
            arr.tuple_as_f64(i, &mut buf);
            assert!((buf[0]).abs() < 1e-10);
        }
    }
}