zigsdof/zigsdof/src/stl.zig

366 lines
12 KiB
Zig

const std = @import("std");
const vec3 = @import("./vec3.zig");
fn vecToFixedPoint(vec: vec3.Vec3F64, digits: u5) [3]u64 {
return .{
@round(vec.x() * std.math.pow(f64, 10.0, digits)),
@round(vec.y() * std.math.pow(f64, 10.0, digits)),
@round(vec.z() * std.math.pow(f64, 10.0, digits)),
};
}
pub const Mesh = struct {
vertices: std.ArrayList(vec3.Vec3F64),
triangle_vertices: std.ArrayList([3]usize),
triangle_normals: std.ArrayList(vec3.Vec3F64),
pub fn fromFile(
file_path: []const u8,
allocator: std.mem.Allocator,
io: std.Io,
) !Mesh {
const cwd = std.Io.Dir.cwd();
const file = try cwd.openFile(io, file_path, .{ .mode = .read_only });
defer file.close(io);
var file_reader = file.reader(io, &.{});
const bytes = try file_reader.interface.allocRemaining(allocator, .unlimited);
defer allocator.free(bytes);
return fromByteArray(bytes, allocator);
}
pub fn fromByteArray(
bytes: []const u8,
allocator: std.mem.Allocator,
) !Mesh {
// Skip header
var offset: usize = 80;
const triangle_num: u32 = std.mem.readInt(u32, @ptrCast(bytes[offset .. offset + 4].ptr), .little);
offset += 4;
var triangle_vertices: std.ArrayList([3]usize) = try .initCapacity(allocator, triangle_num);
var triangle_normals: std.ArrayList(vec3.Vec3F64) = try .initCapacity(allocator, triangle_num);
var vertices: std.ArrayList(vec3.Vec3F64) = try .initCapacity(allocator, triangle_num * 3);
var vertex_map: std.AutoHashMap([3]u64, usize) = .init(allocator);
defer vertex_map.deinit();
for (0..triangle_num) |_| {
var normal = vec3.Vec3F64.init(0, 0, 0);
inline for (0..3) |i| {
const bits = std.mem.readInt(
u32,
@ptrCast(bytes[offset .. offset + 4].ptr),
.little,
);
offset += 4;
const val: f32 = @bitCast(bits);
normal.data[i] = val;
}
var indices: [3]usize = .{ 0, 0, 0 };
inline for (0..3) |i| {
var vertex = vec3.Vec3F64.init(0, 0, 0);
inline for (0..3) |j| {
const bits = std.mem.readInt(
u32,
@ptrCast(bytes[offset .. offset + 4].ptr),
.little,
);
offset += 4;
const val: f32 = @bitCast(bits);
vertex.data[j] = val;
}
const gop = try vertex_map.getOrPut(vecToFixedPoint(vertex, 5));
if (gop.found_existing) {
indices[i] = gop.value_ptr.*;
} else {
const new_idx = vertices.items.len;
try vertices.append(allocator, vertex);
gop.value_ptr.* = new_idx;
indices[i] = new_idx;
}
}
// skip attribute bytes
offset += 2;
try triangle_vertices.append(allocator, indices);
try triangle_normals.append(allocator, normal);
}
return .{
.triangle_vertices = triangle_vertices,
.triangle_normals = triangle_normals,
.vertices = vertices,
};
}
pub fn deinit(mesh: *Mesh, allocator: std.mem.Allocator) void {
mesh.triangle_vertices.deinit(allocator);
mesh.triangle_normals.deinit(allocator);
mesh.vertices.deinit(allocator);
}
};
pub fn rotateMeshAroundPoint(mesh: *Mesh, p: vec3.Vec3F64, q: vec3.QuatF64) void {
for (mesh.vertices.items) |*vertex| {
const d = vec3.vec3Sub(vertex.*, p);
const dprime = vec3.quatApply(q, d);
vertex.* = vec3.vec3Add(dprime, p);
}
return;
}
pub fn translateMesh(mesh: *Mesh, p: vec3.Vec3F64) void {
for (mesh.vertices.items) |*vertex| {
vertex.* = vec3.vec3Add(vertex.*, p);
}
return;
}
pub fn writeVtk(allocator: std.mem.Allocator, io: std.Io, mesh: *Mesh, file_path: []const u8) !void {
const out_string = try createVtkString(allocator, mesh);
defer allocator.free(out_string);
const cwd = std.Io.Dir.cwd();
const file = try cwd.openFile(io, file_path, .{ .mode = .read_write });
defer file.close(io);
const buf: [266]u8 = undefined;
const writer = file.writer(io, buf);
try writer.interface.writeAll(out_string);
try writer.interface.flush();
return;
}
// Create VTK file contents from a mesh.
// NOTE: caller is responsible for cleaning up memory.
fn createVtkString(allocator: std.mem.Allocator, mesh: *Mesh) ![]const u8 {
var output = try std.ArrayList(u8).initCapacity(
allocator,
6 * 1000 + mesh.vertices.items.len * 16 + mesh.triangle_vertices.items.len * 16,
);
const header =
\\# vtk DataFile Version 2.0
\\
\\ASCII
\\DATASET POLYDATA
\\
;
try output.appendSlice(allocator, header);
var buf: [256]u8 = undefined;
const points_header = try std.fmt.bufPrint(
&buf,
"POINTS {} float\n",
.{mesh.vertices.items.len},
);
try output.appendSlice(allocator, points_header);
for (mesh.vertices.items) |vertex| {
const vtx = try std.fmt.bufPrint(&buf, "{} {} {}\n", .{
vertex.x(),
vertex.y(),
vertex.z(),
});
try output.appendSlice(allocator, vtx);
}
const polydata_header = try std.fmt.bufPrint(
&buf,
"POLYGONS {} {}\n",
.{
mesh.triangle_vertices.items.len,
mesh.triangle_vertices.items.len * 4,
},
);
try output.appendSlice(allocator, polydata_header);
for (mesh.triangle_vertices.items) |tri| {
const triangle = try std.fmt.bufPrint(&buf, "3 {} {} {}\n", .{
tri[0],
tri[1],
tri[2],
});
try output.appendSlice(allocator, triangle);
}
return output.toOwnedSlice(allocator);
}
test "Verify fixed point conversion works" {
const vec = vec3.Vec3F64.init(1.0, 2.134, 1010.11);
const fixed_point = vecToFixedPoint(vec, 2);
try std.testing.expect(fixed_point[0] == 100);
try std.testing.expect(fixed_point[1] == 213);
try std.testing.expect(fixed_point[2] == 101011);
}
test "Test binary stl is parsed correctly" {
const binary_stl = "\x42\x69\x6e\x61\x72\x79\x20\x53\x54\x4c\x3a\x20\x32\x20\x74\x72\x69\x61" ++
"\x6e\x67\x6c\x65\x73\x3b\x20\x6e\x6f\x72\x6d\x61\x6c\x73\x20\x66\x72\x6f" ++
"\x6d\x20\x72\x69\x67\x68\x74\x2d\x68\x61\x6e\x64\x2d\x72\x75\x6c\x65\x20" ++
"\x76\x65\x72\x74\x65\x78\x20\x77\x69\x6e\x64\x69\x6e\x67\x00\x00\x00\x00" ++
"\x00\x00\x00\x00\x00\x00\x00\x00\x02\x00\x00\x00\x00\x00\x80\x3f\x00\x00" ++
"\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00" ++
"\x00\x00\x00\x00\x00\x00\x80\x3f\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00" ++
"\x80\x3f\x00\x00\x80\x3f\x00\x00\x00\x00\x00\x00\xf3\x04\x35\xbf\xf3\x04" ++
"\x35\x3f\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x80\x3f" ++
"\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x80\x3f\x00\x00\x80\x3f\x00\x00" ++
"\x80\x3f\x00\x00";
const allocator = std.testing.allocator;
var mesh = try Mesh.fromByteArray(binary_stl, allocator);
defer mesh.deinit(allocator);
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[0],
vec3.Vec3F64.init(0, 0, 0),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[1],
vec3.Vec3F64.init(0, 1, 0),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[2],
vec3.Vec3F64.init(0, 1, 1),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[3],
vec3.Vec3F64.init(1, 0, 0),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[4],
vec3.Vec3F64.init(1, 1, 1),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.triangle_normals.items[0],
vec3.Vec3F64.init(1, 0, 0),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.triangle_normals.items[1],
vec3.Vec3F64.init(0, -0.7071067690849304, 0.7071067690849304),
1e-12,
));
const expected_first_tri = .{ 0, 1, 2 };
try std.testing.expect(std.mem.eql(usize, &expected_first_tri, &mesh.triangle_vertices.items[0]));
const expected_second_tri = .{ 0, 3, 4 };
try std.testing.expect(std.mem.eql(usize, &expected_second_tri, &mesh.triangle_vertices.items[1]));
}
fn createDummyMesh(allocator: std.mem.Allocator) !Mesh {
var vertices: std.ArrayList(vec3.Vec3F64) = .empty;
try vertices.appendSlice(allocator, &.{
vec3.Vec3F64.init(0, 0, 0),
vec3.Vec3F64.init(0.5, std.math.sqrt(3.0) / 2.0, 0),
vec3.Vec3F64.init(1, 0, 0),
});
var triangle_normals: std.ArrayList(vec3.Vec3F64) = .empty;
try triangle_normals.appendSlice(allocator, &.{
vec3.Vec3F64.init(1, 0, 0),
});
var tri_indices: std.ArrayList([3]usize) = .empty;
try tri_indices.appendSlice(allocator, &.{.{ 0, 1, 2 }});
return .{
.vertices = vertices,
.triangle_normals = triangle_normals,
.triangle_vertices = tri_indices,
};
}
test "mesh is rotated correctly" {
const allocator = std.testing.allocator;
var mesh = try createDummyMesh(allocator);
defer mesh.deinit(allocator);
var orig_mesh: Mesh = .{
.vertices = try mesh.vertices.clone(allocator),
.triangle_normals = try mesh.triangle_normals.clone(allocator),
.triangle_vertices = try mesh.triangle_vertices.clone(allocator),
};
defer orig_mesh.deinit(allocator);
rotateMeshAroundPoint(
&mesh,
vec3.Vec3F64.init(0.5, std.math.sqrt(3.0) / 6.0, 0),
vec3.yawPitchRollToQuat(2.0 * std.math.pi / 3.0, 0, 0),
);
try std.testing.expect(vec3.vec3ApproxEq(
orig_mesh.vertices.items[0],
mesh.vertices.items[1],
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
orig_mesh.vertices.items[1],
mesh.vertices.items[2],
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
orig_mesh.vertices.items[2],
mesh.vertices.items[0],
1e-12,
));
}
test "mesh is translated correctly" {
const allocator = std.testing.allocator;
var mesh = try createDummyMesh(allocator);
defer mesh.deinit(allocator);
translateMesh(&mesh, vec3.Vec3F64.init(1, 1, 1));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[0],
vec3.Vec3F64.init(1, 1, 1),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[1],
vec3.Vec3F64.init(1.5, 1 + std.math.sqrt(3.0) / 2.0, 1),
1e-12,
));
try std.testing.expect(vec3.vec3ApproxEq(
mesh.vertices.items[2],
vec3.Vec3F64.init(2, 1, 1),
1e-12,
));
}
test "vtk string is written correctly" {
const allocator = std.testing.allocator;
var mesh = try createDummyMesh(allocator);
defer mesh.deinit(allocator);
const vtk_string = try createVtkString(allocator, &mesh);
defer allocator.free(vtk_string);
const expected =
\\# vtk DataFile Version 2.0
\\
\\ASCII
\\DATASET POLYDATA
\\POINTS 3 float
\\0 0 0
\\0.5 0.8660254037844386 0
\\1 0 0
\\POLYGONS 1 4
\\3 0 1 2
\\
;
try std.testing.expect(std.mem.eql(u8, vtk_string, expected));
}