In [1]:

    

using RegionTrees
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using Revise


    

In [2]:

    

import RegionTrees: needs_refinement, refine_data


    

In [3]:

    

struct AdaptiveDistanceField{C <: Cell} <: Function
    root::C
end

function AdaptiveDistanceField(signed_distance::Function, origin::AbstractArray,
              widths::AbstractArray,
              rtol=1e-2,
              atol=1e-2)
    refinery = SignedDistanceRefinery(signed_distance, atol, rtol)
    boundary = HyperRectangle(origin, widths)
    root = Cell(boundary, refine_data(refinery, boundary))
    adaptivesampling!(root, refinery)
    AdaptiveDistanceField(root)
end


    

    
Out[3]:





AdaptiveDistanceField

In [27]:

    

struct SignedDistanceRefinery{F <: Function} <: AbstractRefinery
    signed_distance_func::F
    atol::Float64
    rtol::Float64
end

function _get_cubeindex(iso_vals, iso)
    cubeindex = iso_vals[1] < iso ? 0x01 : 0x00
    iso_vals[2] < iso && (cubeindex |= 0x02)
    iso_vals[3] < iso && (cubeindex |= 0x04)
    iso_vals[4] < iso && (cubeindex |= 0x08)
    cubeindex
end

function needs_refinement(refinery::SignedDistanceRefinery, cell::Cell)
    minimum(cell.boundary.widths) > refinery.atol || return false

    val = _get_cubeindex(cell.data, 0)
    #@show val
    bounds = cell.boundary
    mid = (bounds.origin.+bounds.widths)./2
    
    value_interp = sum(cell.data)
    value_true = refinery.signed_distance_func(mid)
    if (val == 0x0f && value_true < 0) || (iszero(val) && value_true > 0)
        return false
    elseif !isapprox(value_interp, value_true, rtol=refinery.rtol, atol=refinery.atol)
        return true
    else
        return false
    end
end

function refine_data(refinery::SignedDistanceRefinery, cell::Cell, indices)
    #@show indices
    refine_data(refinery, child_boundary(cell,indices))
end

function refine_data(refinery::SignedDistanceRefinery, bounds::HyperRectangle)
    pts = (bounds.origin, bounds.origin.+bounds.widths.*SVector(1,0),
           bounds.origin.+bounds.widths.*SVector(1,1),bounds.origin.+bounds.widths)
    map(refinery.signed_distance_func, pts)
end


    

    
Out[27]:





refine_data (generic function with 3 methods)

In [49]:

    

true_signed_distance(x) = sqrt(sum(x.^2)) - 1


    

    
Out[49]:





true_signed_distance (generic function with 1 method)

In [50]:

    

origin = SVector(-1, -1)
widths = SVector(2, 2)
@benchmark adf = AdaptiveDistanceField(true_signed_distance, origin, widths, 0.001,0.001)


    

    
Out[50]:





BenchmarkTools.Trial: 
  memory estimate:  183.92 MiB
  allocs estimate:  3013414
  --------------
  minimum time:     99.183 ms (0.00% GC)
  median time:      144.120 ms (32.14% GC)
  mean time:        158.957 ms (37.31% GC)
  maximum time:     340.307 ms (71.36% GC)
  --------------
  samples:          32
  evals/sample:     1

In [51]:

    

@show length(collect(allleaves(adf.root)))


    

    




length(collect(allleaves(adf.root))) = 904021






    
Out[51]:





904021

In [ ]:

    

# compare to uniform sampling