Attempt to handle AbstractRange types for coordinate vectors

[Kevin-Mattheus-Moerman]

@t-bltg thanks for developing MarchingCubes.jl! I am working on using it in my Comodo.jl package.

I wanted to extend MarchingCubes to allow for range type coordinate input, i.e. such that when using this:

mc = MarchingCubes.MC(A,Int; x=xr, y=yr, z=zr)

the entries xr etc can be of the type AbstractRange.

Currently I think I implemented a very basic fix, but perhaps you can point me to improving it if what I propose is not efficient.

In MarchingCubes.jl I changed the types from AbstractVector{F} = F[], to Union{AbstractVector{F}, AbstractRange{F}} = F[],
Next in the definition of MC I added stuff like the below to collect the range into a vector before it is passed to Ref(x) :

       if isa(x,AbstractRange)
            x = collect(x)
        end

I did the simple collect mainly because I could not quickly figure out how to use/fix that Ref and RefValue business, so any advise on that is welcome.

Below is the code I used for testing (which also contains Makie based visualisation):

using LinearAlgebra
using GLMakie
using GeometryBasics
using MarchingCubes

# Define image matrix
nSteps = 50 
xr,yr,zr = ntuple(_->range(-1.0,1.0,nSteps),3) # Coordinate ranges 
A = [norm((x,y,z)) for x in xr, y in yr, z in zr] # Image 

# Using MarchingCubes to get the isosurface 
level = 0.5 # Level for isosurface 
mc = MarchingCubes.MC(A,Int; x=xr, y=yr, z=zr)
MarchingCubes.march(mc,level)

# Access faces and vertices and poor in GeometryBasics format to enable Makie based visualisation
F = [TriangleFace{Int64}(f) for f in mc.triangles]
V = [Point{3,Float64}(p) for p in mc.vertices]

# Visualization
fig = Figure(size=(800,800))
ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", limits=(minimum(xr),maximum(xr),minimum(yr),maximum(yr),minimum(zr),maximum(zr)))
hp1 = poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=1, strokecolor=:black, color=:white,shading=FastShading,transparency=false)
fig

Let me know if this is going in a direction you like or if the feature is not at all of interest (or what changes would be needed to have it fit).

Thanks!

[t-bltg]

Hi, @Kevin-Mattheus-Moerman.
Great idea !

An AbstractRange is a subtype of an AbstractVector:

julia> supertypes(AbstractRange)
(AbstractRange, AbstractVector, Any)

, so I think x::AbstractVector{F} = F[] should already work if given an AbstractRange, no ? If no, maybe you can omit the {F} part which is probably too restrictive.

Regarding the test failures, I think we must adapt to JuliaGeometry/GeometryBasics.jl#226: normals -> normal, can you have a look into that too ?

[Kevin-Mattheus-Moerman]

An AbstractRange is a subtype of an AbstractVector , so I think x::AbstractVector{F} = F[] should already work if given an AbstractRange, no ? If no, maybe you can omit the {F} part which is probably too restrictive.

Yes, sorry, I missed that. I'll revert that so.

@t-bltg so the only "issue" in terms of handling ranges then seems to be relating to that Ref and RefValue stuff (which I am completely new to and currently confused about). Is there an easy way for that to handle ranges? Or is the collect(x) approach I took fine?

Regarding the test failures, I think we must adapt to JuliaGeometry/GeometryBasics.jl#226: normals -> normal, can you have a look into that too ?

I can have a look, but I take it this issue does not relate to my PR.

[t-bltg]

This seems to work:

    x::RefValue{<:AbstractVector{F}}
    y::RefValue{<:AbstractVector{F}}
    z::RefValue{<:AbstractVector{F}}

The RefValue thing was added so that the struct attributes x,y,z can be mutated: mc.x[] = <some new vector>, without making the struct mutable.

Can you add a small test too ?

[t-bltg]

I can have a look, but I take it this issue does not relate to my PR.

I don't mind you fixing the tests in this PR, since the MarchingCubes.makemesh helper is not guaranteed to be stable.

[Kevin-Mattheus-Moerman]

Cool, thanks. I implemented the x::RefValue{<:AbstractVector{F}} and removed the need for convert with that.

[Kevin-Mattheus-Moerman]

@t-bltg there seems to be a new testing error now where bytes = @allocated march(mc) is not zero. However, when I test this on my machine I cannot reproduce the error and this bytes tests work fine.

edit: I now sometimes see this error locally too. Not sure what is causing it though.

[Kevin-Mattheus-Moerman]

@t-bltg I think the latest commit fixes the GeometryBasics normals -> normal bit (and also Mesh(meta(kwargs...), faces) -> Mesh(faces; kwargs...)). I think I have to call it a day now. Also I am unsure about the above errors. What are your thoughts on these?

[t-bltg]

@t-bltg there seems to be a new testing error now where bytes = @allocated march(mc) is not zero. However, when I test this on my machine I cannot reproduce the error and this bytes tests work fine.

Heck, there seems to be a StaticArrays regression here:

julia> using MarchingCubes, AllocCheck
julia> mc = MarchingCubes.scenario(case = :cushin)
julia> @check_allocs foo(mc) = march(mc)
julia> try
    foo(mc)
catch err
    foreach(display, err.errors)
end
Allocation of Memory{StaticArraysCore.SVector{3, Float32}} in ./boot.jl:516
  | ccall(:jl_alloc_genericmemory, Ref{GenericMemory{kind,T,addrspace}}, (Any, Int), self, m)

Stacktrace:
 [1] GenericMemory
   @ ./boot.jl:516 [inlined]
 [2] array_new_memory
   @ ./array.jl:1058 [inlined]
 [3] (::Base.var"#133#134"{Vector{StaticArraysCore.SVector{3, Float32}}, Int64, Int64, Int64, Int64, Int64, Memory{StaticArraysCore.SVector{3, Float32}}, MemoryRef{StaticArraysCore.SVector{3, Float32}}})()
   @ Base ./array.jl:1129
 [4] multiple call sites
   @ unknown:0

Allocation of Memory{StaticArraysCore.SVector{3, Int64}} in ./boot.jl:516
  | ccall(:jl_alloc_genericmemory, Ref{GenericMemory{kind,T,addrspace}}, (Any, Int), self, m)

Stacktrace:
 [1] GenericMemory
   @ ./boot.jl:516 [inlined]
 [2] array_new_memory
   @ ./array.jl:1058 [inlined]
 [3] (::Base.var"#133#134"{Vector{StaticArraysCore.SVector{3, Int64}}, Int64, Int64, Int64, Int64, Int64, Memory{StaticArraysCore.SVector{3, Int64}}, MemoryRef{StaticArraysCore.SVector{3, Int64}}})()
   @ Base ./array.jl:1129
 [4] multiple call sites
   @ unknown:0

# other errors related to `@warn` or `@error` calls omitted ...

[t-bltg]

Heck, there seems to be a StaticArrays regression here:

Hum, no I was wrong, I think that RefValue{...} has to take a concrete type for the compiler to use static arrays and avoid allocations.

So I think we should restore RefValue{Vector{F}}, and use collect on AbstractRanges in the constructor as you did before.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 98.75%. Comparing base (b3e2ca1) to head (68b78d8).
Report is 11 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #29      +/-   ##
==========================================
- Coverage   99.58%   98.75%   -0.84%     
==========================================
  Files           2        2              
  Lines         483      482       -1     
==========================================
- Hits          481      476       -5     
- Misses          2        6       +4     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[t-bltg]

The tests are now passing.
Adding some new tests for AbstractRanges would be appreciated !

[Kevin-Mattheus-Moerman]

@t-bltg okay I added testing for ranged coordinate input just now. The tests seem to pass mostly. Note my added tests might be a bit verbose with comments, as I noticed the other tests did not have comments. Feel free to purge it of the comments if you feel that is cleaner. Let me know if you feel the tests I added are appropriate or if more or some changes are needed. The tests I added check if for vector or range based coordinate input the following are the same:

• The output vertices and output triangles (if these are the same the normals and other metrics follow too). I used a simple == here, I hope that is okay.
• The coordinate input is used in the same way, i.e. the output geometry is as expected (presents with a zero mean and a known mean radius for the spherical isosurface). These tests are perhaps not needed if the standard coordinate output is already validated and the ranged version is tested to be the same. However, I did not see a specific test to see if the coordinate input is treated properly yet so I added it. Feel free to amend/edit/delete.

I tried addressing the format issue there but not sure what it relates to. Also not sure why the Julia nightly tests relating to bytes = @allocated march(mc) fail.

Side note, I removed using BenchmarkTools from runtests.jl as the package appeared unused.

[t-bltg]

Thanks for the added tests !
I will now format the code, and use my convention for the comments, and merge when the tests pass.

For the nightly failures, me don't mind about them, we will when an rc or a beta version will be out.