@stevengj we had a conversation (at meepcon) about some of the potential pitfalls of the current implementation of GeometryPrimitives.jl compared to the C version of libctl.

You mentioned the primary issue had to do with dispatch. libctl has a custom dispatch routine that is rather optimized for many function calls, whereas the current GeometryPrimitives.jl implementation tries to generalize for multiple dispatch (but I may be misremembering...).

Could you talk a bit more about that here and maybe discuss what it would take to overcome some of these issues?

I read a recent article by @ChrisRackauckas that describes how they cut down on compile times by limiting dispatch flexibility. The issue here isn't with compile time, but I'm wondering if there are similar lessons to be learned.

In order to be more consistent with the definition of Base.findfirst(), wouldn't it be better to return the index of the shape, rather than the shape itself?

Hi, we are leveraging this package for one of our packages, MESTI.jl, to perform subpixel smoothing. Thanks for your great package!

Since our package utilizes your package and has the dependency. When we recently registered our updated version of MESTI.jl. We found an error message from GeometryPrimitives.jl during tests in the registration.

Error: LoadError: UndefVarError: SurfaceLike not defined

The error seems to come from recent updates in Makie.jl. This error also appears when we try to install the GeometryPrimitives.jl package on a new machine. By default, Julia utilizes the latest Makie.jl (v0.20.3).

This error can be fixed by adding the restriction on the old version of Makie (v0.19.12) when we register our package or install the GeometryPrimitives.jl package.

I am wondering if you have any plans to fix this possible error by updating your GeometryPrimitives.jl to utilize the latest Makie.jl? Thanks.

I was leveraging this useful little package for one of my codes which I'm now trying to update to support Julia 0.7/1.0.

Do you guys have any plans to update this to support Julia 0.7/1.0?

Currently surfpt_nearby is type-stable for every concrete subtype of Shape{N} and returns NTuple{2,SVector{N,Float64}}. However, when surfpt_nearby is called with an abstract Shape{N} object, it is type-unstable. For example, consider the following code:

julia> VERSION
v"0.6.1-pre.0"

julia> using GeometryPrimitives, StaticArrays

julia> svec = Vector{Shape{3}}(2);

julia> function test(p, svec::Vector{Shape{3}})
           rsum = @SVector zeros(3)
           nsum = @SVector zeros(3)
           for s = svec
               r, n = surfpt_nearby(p, s)
               rsum += r
               nsum += n
           end
           return rsum, nsum
       end
test (generic function with 1 method)

The above test function is type-unstable, because the second argument s of surfpt_nearby is an element of Vector{Shape{3}} and thus typed Shape{3}:

julia> @code_warntype test(@SVector(zeros(3)), svec)
Variables:
  ...
Body:
  begin
      ...
      SSAValue(2) = (Main.surfpt_nearby)(p::SVector{3,Float64}, s::GeometryPrimitives.Shape{3,D} where D)::Any
      ...
  end::Tuple{Any,Any}

This makes sense. When surfpt_nearby is called with an abstract Shape object, the following method defined in GeometryPrimitives.jl is called:

surfpt_nearby(x::AbstractVector, o::Shape{N}) where {N} = surfpt_nearby(SVector{N}(x), o)

Because the specific type of o is not known, the compiler cannot further infer the return type in this case.

To resolve this type instability, I tried to change the above method definition using return type annotation:

(surfpt_nearby(x::AbstractVector, o::Shape{N})::NTuple{2,SVector{N,Float64}}) where {N} = surfpt_nearby(SVector{N}(x), o)

I thought this would remove the type instability, because now the compiler would know the return type of the above method. Strangely, the type instability remains:

julia> @code_warntype test(@SVector(zeros(3)), svec)
Variables:
  ...
Body:
  begin
      ...
      SSAValue(2) = (Main.surfpt_nearby)(p::SVector{3,Float64}, s::GeometryPrimitives.Shape{3,D} where D)::Any
      ...
  end::Tuple{Any,Any}

I'm not sure which method of surfpt_nearby is being called here. To me, the only method with a matching argument types is the one defined in GeometryPrimitives.jl, which is now annotated with a return type. Still, the compiler fails to pick up this annotated return type.

In fact, if we call surfpt_nearby with a regular Vector instead of an SVector in the first argument p, the return type annotation works as intended and test becomes type-stable:

julia> @code_warntype test(zeros(3), svec)
Variables:
  ...
Body:
  begin
      ...
      SSAValue(2) = (Main.surfpt_nearby)(p::Array{Float64,1}, s::GeometryPrimitives.Shape{3,D} where D)::Tuple{SVector{3,Float64},SVector{3,Float64}}
      ...
  end::Tuple{SVector{3,Float64},SVector{3,Float64}}

I don't understand why the return type annotation does not kick in for SVector p. How can we resolve this problem? Any suggestions?

Currently all the concrete shapes inherit from Shape{N}. I hope them to inherit from Shape{N,D}, where D specifies the type of the data field of concrete shape types. (This will require changing the definition of Shape to abstract type Shape{N,D} end.)

The reason for this change is as follows. In my another module that uses GeometryPrimitives, I need to force all shapes to take an additional SVector, which can be conveniently stored in data. Then, when I create an array of these shapes, I would want to type it ::Vector{Shape{N,SVector{N}}}, which is only doable when all the concrete shapes inherit from Shape{N,D} rather than Shape{N}.

Because Shape{N,D} is a subtype of Shape{N}, all the concrete types will remain subtypes of Shape{N}. Therefore, this change will not do much harm on existing code (e.g., functions taking Shape{N} can remain taking Shape{N}).

KDTree constructs the left and right child trees recursively. There are cases where this recursion occurs endlessly, leading to StackOverflowError. Such as a case occurs when attempting to construct a KDTree with the array of 2D Spheres:

s = [Sphere([0,0], 1), Sphere([0,0], 1), Sphere([0,0], 1), Sphere([0,0], 1),  # all centered at origin
     Sphere([2,0], 1), Sphere([3,0], 1), Sphere([4,0], 1)]

Here, we have first four spheres centered at the origin and the next three spheres on the positive side of the x-axis. Because we have more than a half of the entire shapes at the origin, x = 0 is taken as the median separating the left and right trees. Then, the left tree is constructed with the first 4 spheres because their left boundaries are below x = 0, but the right tree is constructed with the entires spheres because all of them have right boundaries above x = 0. Because the right tree is constructed with the entires spheres, it faces the same situation as the parent tree, and the recursive construction repeats forever.

Such infinite recursion can be prevented by giving up branching when either one of the left and right child trees has a similar number of shapes as the parent tree.

@stevengj do you have any examples that use GeometryPrimitives (e.g. maybe a notebook from one of your classes)?

The tests are pretty comprehensive, but it would be nice to see the entire "pipeline" in action (i.e. different primitives that correspond to different materials and visualizing the net geometry).

Something with subpixel smoothing (i.e. using volfrac, surfpt_nearby , and, normal ) would be really cool too. It seems somewhat straightforward, but I figured I'd ask first so as to not reinvent the wheel. Thanks!

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Currently, the Ellipsoid constructor takes the diameters d, whereas the Sphere constructor takes the radius r.

Because sphere is a special case of ellipsoid, to me it seems more consistent to change the Ellipsoid constructor to take the radii as well (or to change the Sphere constructor to take the diameter).

Can I create a PR on this, or it might be better to leave them as it is?