diff --git a/Project.toml b/Project.toml index edb9afd..448eb5c 100644 --- a/Project.toml +++ b/Project.toml @@ -10,12 +10,13 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" +ImplicitBVH = "932a18dc-bb55-4cd5-bdd6-1368ec9cea29" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] Adapt = "4" @@ -25,13 +26,14 @@ GeometryBasics = "0.5" ImageCore = "0.10" ImageIO = "0.6" ImageMagick = "1" +ImplicitBVH = "0.7.0" JSON = "1" LinearAlgebra = "1" MeshIO = "0.5" NearestNeighbors = "0.4" StaticArrays = "1" -StructTypes = "1.11.0" Test = "1" +TimerOutputs = "0.5.29" julia = "1" [extras] diff --git a/dev_script.jl b/dev_script.jl index 562e997..8785b30 100644 --- a/dev_script.jl +++ b/dev_script.jl @@ -1,12 +1,43 @@ using GeometryBasics +using ImplicitBVH using Slicey using StaticArrays +using TimerOutputs +# geometry = read_geometry("test/stl_files/Stanford_Bunny.stl") geometry = read_geometry("test/stl_files/cube.stl") x_min = 200.0 * (1920 / 1080) -build_stage = Slicey.DLPBuildStage((0.0, x_min), (0.0, 200.0), (1920, 1080)) -printer = Slicey.DLPPrinter(build_stage) -slicer = Slicey.DLPSlicer(geometry, printer) +build_stage = BuildStage((0.0, x_min), (0.0, 200.0)) +projector = Slicey.DLPProjector(build_stage, (1920, 1080)) +printer = Slicey.DLPPrinter(build_stage, projector) +slicer = Slicey.DLPSlicer(geometry, printer, TimerOutput()) Slicey.center_geometry!(slicer) -Slicey.slice(slicer, 0.1) \ No newline at end of file + +bvh = Slicey.build_triangle_bvh(slicer.geometry) +pixel_centroids = Slicey.build_pixel_centroids(slicer.printer, 20.0) +dirs = zeros(eltype(pixel_centroids), size(pixel_centroids)...) +# dirs[1, :] .= 1f0 +for n in axes(dirs, 2) + dirs[:, n] .= Slicey.RAY_DIR +end +# dirs + +trav = traverse_rays(bvh, pixel_centroids, dirs) + +counts = zeros(Int, size(dirs, 2)) +vs = slicer.geometry.vertices + +for (ray_idx, leaf_idx) in trav.contacts + @show leaf_idx + tri = slicer.geometry.faces[leaf_idx] + p0 = SVector{3, Float32}(@views pixel_centroids[:, ray_idx]...) + v0 = SVector{3, Float32}(vs[tri[1]]...) + v1 = SVector{3, Float32}(vs[tri[2]]...) + v2 = SVector{3, Float32}(vs[tri[3]]...) + if Slicey.ray_intersects_triangle(p0, Slicey.RAY_DIR, v0, v1, v2) + counts[ray_idx] += 1 + end +end +# Slicey.slice(slicer, 0.1; sample_method = :bvh) +# slicer.timer \ No newline at end of file diff --git a/src/Geometry.jl b/src/Geometry.jl index ad7e5c9..eddb2a9 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -76,6 +76,20 @@ function build_face_centroid_kdtree(geometry::STLGeometry) return kdtree end +function build_triangle_bvh( + geometry::STLGeometry; + volume_type = ImplicitBVH.BBox{Float32} +) + vs = geometry.vertices + bs = Vector{ImplicitBVH.BoundingVolume{volume_type, Int32, UInt32}}(undef, length(geometry.faces)) + for (i, face) in enumerate(geometry.faces) + # @show face + bs[i] = ImplicitBVH.BoundingVolume(volume_type(vs[face[1]], vs[face[2]], vs[face[3]]), Int32(i), UInt32(0)) + # @show sphere + end + return BVH(bs, ImplicitBVH.BBox{Float32}) +end + function centroid(geometry::STLGeometry) bb = bounding_box(geometry) return (bb.origin + (bb.origin + bb.widths)) / 2 @@ -91,12 +105,46 @@ end function point_in_mesh( geometry::STLGeometry, - tree, + bvh::BVH, + p::SVector{3, T} +) where T <: Number + face_pts = face_points(geometry) + # Fixed ray direction + dir = Vec3{Float32}(1f0, 0f0, 0f0) + invdir = Vec3{Float32}(1f0 / dir[1], Inf32, Inf32) + hits = 0 + traverse(bvh) do node + if node.isleaf + for i in node.indices + # tri = triangles[idx] + v0 = SVector{3, T}(face_pts[1, i]...) + v1 = SVector{3, T}(face_pts[2, i]...) + v2 = SVector{3, T}(face_pts[3, i]...) + + if ray_triangle_intersect(p, dir, v0, v1, v2) + hits += 1 + end + end + else + ray_bbox_intersect(p, invdir, node.bbox) + end + end + + return isodd(hits) +end + +function point_in_mesh( + geometry::STLGeometry, + tree, idxs, dists, p::SVector{3, T}; - radius=T(Inf) + # k_nearest_neighbors = 3, + radius = T(Inf) ) where T <: Number # all allocations are coming from this guy below - idxs = inrange(tree, p, radius) + # idxs = inrange(tree, p, radius) + # inrange!(idxs, tree, p, radius) + knn!(idxs, dists, tree, p, length(idxs)) + # @show size(idxs) face_pts = face_points(geometry) hits = 0 for i in idxs @@ -109,6 +157,20 @@ function point_in_mesh( isodd(hits) end +@inline function ray_bbox_intersect( + origin::Vec3{Float32}, + invdir::Vec3{Float32}, + bbox::ImplicitBVH.BBox{Float32}, +)::Bool + t1 = (bbox.min .- origin) .* invdir + t2 = (bbox.max .- origin) .* invdir + + tmin = maximum(min.(t1, t2)) + tmax = minimum(max.(t1, t2)) + + return tmax ≥ max(tmin, 0f0) +end + @inline function ray_intersects_triangle( orig::SVector{3, T}, dir::SVector{3, T}, diff --git a/src/Settings.jl b/src/Settings.jl new file mode 100644 index 0000000..9056231 --- /dev/null +++ b/src/Settings.jl @@ -0,0 +1,49 @@ +module Settings + +using JSON + +struct Units + length::String + speed::String + time::String +end + +struct BuildStage + length::Float32 + width::Float32 +end + +struct Geometry + default_type::String + default_parallelism_method::String + default_search_method::String +end + +struct Projector + resolution::String +end + +struct DLPPrinterSettings + build_stage::BuildStage + projector::Projector +end + +const PrinterSettings = Union{ + DLPPrinterSettings +} + +struct SlicerSettings + geometry::Geometry + printer::PrinterSettings + units::Units +end + +function parse_settings(settings_file::String) + return open(settings_file) do f + str = read(f, String) + settings = JSON.parse(str, SlicerSettings) + return settings + end +end + +end # module PrinterSettings diff --git a/src/Slicey.jl b/src/Slicey.jl index f6ee62a..14076e4 100644 --- a/src/Slicey.jl +++ b/src/Slicey.jl @@ -4,6 +4,7 @@ module Slicey # read_geometry(...) instead of # Slicey.read_geometry(...) # you have to be careful of conflicting names though +export BuildStage export read_geometry using Adapt @@ -12,12 +13,13 @@ using GeometryBasics using ImageCore using ImageIO using ImageMagick +using ImplicitBVH using JSON using LinearAlgebra using MeshIO using NearestNeighbors using StaticArrays -using StructTypes +using TimerOutputs # will this fully evaluate the file? will the init_mssg show? include("Geometry.jl") @@ -25,7 +27,7 @@ include("InitializeMessage.jl") # re-incroporate TODO # include("Slicing_Prep.jl") # include("Slicing.jl") -include("parsers/PrinterSettings.jl") +include("Settings.jl") include("slicers/Slicers.jl") # Need to include a settings file, something with all the diff --git a/src/parsers/PrinterSettings.jl b/src/parsers/PrinterSettings.jl deleted file mode 100644 index 95b3de7..0000000 --- a/src/parsers/PrinterSettings.jl +++ /dev/null @@ -1,47 +0,0 @@ -struct Units - length::String - speed::String - time::String -end - -struct BuildStage - length::Float32 - width::Float32 -end - -struct Projector - resolution::String -end - -struct DLPPrinterSettings - build_stage::BuildStage - projector::Projector -end - -const AllPrinterSettings = Union{ - DLPPrinterSettings -} - -struct PrinterSettings - settings::DLPPrinterSettings - units::Units -end - -struct MyTestStruct - units::Units -end - -StructTypes.StructType(::Type{Units}) = StructTypes.Struct() -StructTypes.StructType(::Type{BuildStage}) = StructTypes.Struct() -StructTypes.StructType(::Type{Projector}) = StructTypes.Struct() -StructTypes.StructType(::Type{DLPPrinterSettings}) = StructTypes.Struct() -StructTypes.StructType(::Type{PrinterSettings}) = StructTypes.Struct() - -function parse_printer_settings(settings_file::String) - return open(settings_file) do f - str = read(f, String) - @show str - settings = JSON.parse(str, PrinterSettings) - return settings - end -end diff --git a/src/slicers/DLPSlicer.jl b/src/slicers/DLPSlicer.jl index ac8087a..5133a6b 100644 --- a/src/slicers/DLPSlicer.jl +++ b/src/slicers/DLPSlicer.jl @@ -1,48 +1,51 @@ -# TODO we should really break out the projector -# stuff from the rest of the build stage... -struct DLPBuildStage{I <: Integer, T} <: AbstractBuildStage{T} - bb::GeometryBasics.HyperRectangle{3, T} +struct DLPProjector{I <: Integer, T <: Number} pixel_count::NTuple{2, I} resolution::NTuple{2, T} end -function DLPBuildStage( - x_bounds::NTuple{2, T}, - y_bounds::NTuple{2, T}, +function DLPProjector( + build_stage::BuildStage{T}, pixel_count::NTuple{2, I} -) where {I <: Integer, T <: Number} - xmin, xmax = x_bounds - ymin, ymax = y_bounds - xlen, ylen = Float32(xmax - xmin), Float32(ymax - ymin) - xres = xlen / pixel_count[1] - yres = ylen / pixel_count[2] +) where {T <: Number, I <: Integer} + xres = build_stage.bb.widths[1] / pixel_count[1] + yres = build_stage.bb.widths[2] / pixel_count[2] + return DLPProjector(pixel_count, (xres, yres)) +end + + + +struct DLPPrinter{I <: Integer, T <: Number} <: AbstractPrinter{T} + build_stage::BuildStage{T} + projector::DLPProjector{I, T} +end - bb = GeometryBasics.HyperRectangle{3, Float32}((xmin, ymin, 0), (xlen, ylen, 0)) - return DLPBuildStage(bb, pixel_count, (xres, yres)) +function build_pixel_centroids(printer::DLPPrinter, z) + projector = printer.projector + centroids = zeros(Float32, 3, projector.pixel_count[2] * projector.pixel_count[1]) + n = 1 + for j in 1:projector.pixel_count[1] + for i in projector.pixel_count[2] + centroids[:, n] .= pixel_centroid(printer, i, j, z) + n = n + 1 + end + end + return centroids end # this z height might not be right -function pixel_centroid(build_stage::DLPBuildStage, i::Int, j::Int, z::T) where T <: Number +function pixel_centroid(printer::DLPPrinter, i::Int, j::Int, z::T) where T <: Number # x = build_stage.bb.origin[1] + (j + 0.5f0) * build_stage.resolution[1] # y = build_stage.bb.origin[2] + (i + 0.5f0) * build_stage.resolution[2] - x = build_stage.bb.origin[1] + (j - 0.5f0) * build_stage.resolution[1] - y = build_stage.bb.origin[2] + (i - 0.5f0) * build_stage.resolution[2] + bs, pro = printer.build_stage, printer.projector + x = bs.bb.origin[1] + (j - 0.5f0) * pro.resolution[1] + y = bs.bb.origin[2] + (i - 0.5f0) * pro.resolution[2] return Vec3{T}(x, y, z) end -struct DLPPrinter{B <: DLPBuildStage} <: AbstractPrinter{B} - build_stage::B -end - struct DLPSlicer{G <: AbstractPartGeometry, P <: DLPPrinter} <: AbstractSlicer{G, P} geometry::G printer::P -end - -function DLPSlicer(build_stage::DLPBuildStage, geometry_file::String) - geometry = STLGeometry(geometry_file) - printer = DLPPrinter(build_stage) - return DLPSlicer(geometry, printer) + timer::TimerOutput end function slice( @@ -50,47 +53,117 @@ function slice( sample_method = :brute_force, parallelism = :serial ) - z_heights = layer_heights(slicer, z_height) - slice( - slicer, z_heights, Val{parallelism}(); - sample_method = sample_method - ) + @timeit slicer.timer "Z height calculate" begin + z_heights = layer_heights(slicer, z_height) + end + @timeit slicer.timer "Slice" begin + slice( + slicer, z_heights, Val{parallelism}(); + sample_method = sample_method + ) + end end function slice( slicer::DLPSlicer, z_heights, ::Val{:serial}; + k_nearest_neighbors = 12, sample_method = :brute_force ) img_file_base = splitext(slicer.geometry.file_name)[1] - image = zeros(Bool, slicer.printer.build_stage.pixel_count[2], slicer.printer.build_stage.pixel_count[1]) + projector = slicer.printer.projector + image = zeros(Bool, projector.pixel_count[2], projector.pixel_count[1]) + img = Gray.(image) sample_method = Val{sample_method}() - for (n, z) in enumerate(z_heights) - @info "Slicing layer $n" - fill!(image, 0) - _process_layer!(image, slicer, z, sample_method) - img_file = img_file_base * "_$(lpad(n, 6, "0")).bmp" - img = Gray.(image) - save(img_file, img) + + if sample_method == Val{:kdtree}() + @timeit slicer.timer "Build kdtree" begin + # idxs = Int32[] + # dists = Float32[] + + # creating one for each thread so we don't overwrite + idxs = [zeros(Int32, k_nearest_neighbors) for n in 1:Threads.nthreads()] + dists = [zeros(Float32, k_nearest_neighbors) for n in 1:Threads.nthreads()] + tree = build_face_centroid_kdtree(slicer.geometry) + end + elseif sample_method == Val{:bvh}() + bvh = build_triangle_bvh(slicer.geometry) + end + Threads.@threads for (n, z) in collect(enumerate(z_heights)) + @timeit slicer.timer "Layer slice" begin + @info "Slicing layer $n" + fill!(image, 0) + @timeit slicer.timer "Process layer" begin + if sample_method == Val{:kdtree}() + tid = Threads.threadid() + _process_layer!(image, slicer, z, sample_method, tree, idxs[tid], dists[tid]) + elseif sample_method == Val{:bvh}() + _process_layer!(image, slicer, z, sample_method, bvh) + else + _process_layer!(image, slicer, z, sample_method) + end + end + @timeit slicer.timer "Create image" begin + map!(Gray, img, image) + end + @timeit slicer.timer "Save image" begin + img_file = img_file_base * "_$(lpad(n, 6, "0")).bmp" + save(img_file, img) + end + end end end # some internals function _process_layer!(image, slicer::DLPSlicer, z::T, ::Val{:brute_force}) where T <: Number + projector = slicer.printer.projector + # create image thing to sample over + for i in 1:projector.pixel_count[2] + for j in 1:projector.pixel_count[1] + # @timeit slicer.timer "Pixel geometry search" begin + cent = pixel_centroid(slicer.printer, i, j, z) + cent = SVector{3, T}(cent...) + # now loop over STL triangles + hits = 0 + for tri in axes(slicer.geometry.face_points, 2) + v0 = SVector{3, T}(slicer.geometry.face_points[1, tri]...) + v1 = SVector{3, T}(slicer.geometry.face_points[2, tri]...) + v2 = SVector{3, T}(slicer.geometry.face_points[3, tri]...) + ray_intersects_triangle(cent, RAY_DIR, v0, v1, v2) && (hits += 1) + end + + image[i, j] = isodd(hits) + # end + end + end +end + +function _process_layer!(image, slicer::DLPSlicer, z::T, ::Val{:bvh}, bvh) where T <: Number + projector = slicer.printer.projector # create image thing to sample over - for i in 1:slicer.printer.build_stage.pixel_count[2] - for j in 1:slicer.printer.build_stage.pixel_count[1] - cent = pixel_centroid(slicer.printer.build_stage, i, j, z) + for i in 1:projector.pixel_count[2] + for j in 1:projector.pixel_count[1] + cent = pixel_centroid(slicer.printer, i, j, z) cent = SVector{3, T}(cent...) - # now loop over STL triangles - hits = 0 - for tri in axes(slicer.geometry.face_points, 2) - v0 = SVector{3, T}(slicer.geometry.face_points[1, tri]...) - v1 = SVector{3, T}(slicer.geometry.face_points[2, tri]...) - v2 = SVector{3, T}(slicer.geometry.face_points[3, tri]...) - ray_intersects_triangle(cent, RAY_DIR, v0, v1, v2) && (hits += 1) - end + # empty!(idxs) + # empty!(dists) + image[i, j] = point_in_mesh(slicer.geometry, bvh, cent) + end + end +end - image[i, j] = isodd(hits) +function _process_layer!(image, slicer::DLPSlicer, z::T, ::Val{:kdtree}, tree, idxs, dists) where T <: Number + projector = slicer.printer.projector + # create image thing to sample over + for i in 1:projector.pixel_count[2] + # @show i + for j in 1:projector.pixel_count[1] + # @timeit slicer.timer "Pixel geometry search" begin + cent = pixel_centroid(slicer.printer, i, j, z) + cent = SVector{3, T}(cent...) + # empty!(idxs) + # empty!(dists) + image[i, j] = point_in_mesh(slicer.geometry, tree, idxs, dists, cent) + # end end end end diff --git a/src/slicers/Slicers.jl b/src/slicers/Slicers.jl index a82448a..c3fec3f 100644 --- a/src/slicers/Slicers.jl +++ b/src/slicers/Slicers.jl @@ -1,16 +1,32 @@ -abstract type AbstractBuildStage{T <: Number} end +# abstract type AbstractBuildStage{T <: Number} end -function bounding_box(build_stage::AbstractBuildStage) +struct BuildStage{T <: Number} + bb::GeometryBasics.HyperRectangle{3, T} +end + +function BuildStage( + x_bounds::NTuple{2, T}, + y_bounds::NTuple{2, T} +) where T <: Number + xmin, xmax = x_bounds + ymin, ymax = y_bounds + xlen, ylen = Float32(xmax - xmin), Float32(ymax - ymin) + bb = GeometryBasics.HyperRectangle{3, Float32}((xmin, ymin, 0), (xlen, ylen, 0)) + return BuildStage{Float32}(bb) +end + +function bounding_box(build_stage::BuildStage) return build_stage.bb end -function centroid(build_stage::AbstractBuildStage) +function centroid(build_stage::BuildStage) bb = bounding_box(build_stage) return (bb.origin + (bb.origin + bb.widths)) / 2 end abstract type AbstractPrinter{ - B <: AbstractBuildStage + # B <: AbstractBuildStage + T <: Number } end # TODO also we can add some diff --git a/test/TestPrinterSettings.jl b/test/TestPrinterSettings.jl new file mode 100644 index 0000000..39bdf2a --- /dev/null +++ b/test/TestPrinterSettings.jl @@ -0,0 +1,14 @@ +function test_printer_settings_units_block() + open("printer_settings/units_for_test.json") do f + str = read(f, String) + block = JSON.parse(str, Slicey.PrinterSettings.Units) + + @test block.length == "mm" + @test block.speed == "mm/s" + @test block.time == "s" + end +end + +@testset "PrinterSettings" begin + test_printer_settings_units_block() +end diff --git a/test/printer_settings/dlp_printter.json b/test/printer_settings/dlp_printter.json index 29720f3..1beef17 100644 --- a/test/printer_settings/dlp_printter.json +++ b/test/printer_settings/dlp_printter.json @@ -4,7 +4,12 @@ "speed": "mm/s", "time": "s" }, - "settings": { + "geometry": { + "default_type": "STLGeometry", + "default_parallelism_method": "serial", + "default_search_method": "brute_force" + }, + "printer": { "build_stage": { "length": 355.5555, "width": 200.0 diff --git a/test/printer_settings/units_for_test.json b/test/printer_settings/units_for_test.json index 9fc7c6c..74d781b 100644 --- a/test/printer_settings/units_for_test.json +++ b/test/printer_settings/units_for_test.json @@ -1,7 +1,5 @@ { - "units": { - "length": "mm", - "speed": "mm/s", - "time": "s" - } + "length": "mm", + "speed": "mm/s", + "time": "s" } diff --git a/test/runtests.jl b/test/runtests.jl index 09166c8..426fd13 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,10 @@ using Aqua +using JSON using Slicey using Test +include("TestPrinterSettings.jl") + @testset "Aqua Testing" begin Aqua.test_all(Slicey) end diff --git a/test/stl_files/Stanford_Bunny.stl b/test/stl_files/Stanford_Bunny.stl new file mode 100644 index 0000000..bec02fc Binary files /dev/null and b/test/stl_files/Stanford_Bunny.stl differ