From 77c260a5fb423a433ab4cf738f53060ecb68c44a Mon Sep 17 00:00:00 2001 From: Rafael Caixeta <8386288+rmcaixeta@users.noreply.github.com> Date: Sat, 24 Apr 2021 09:27:57 -0300 Subject: [PATCH 1/3] Additional calls to unfold --- src/coords_optimization.jl | 4 +- src/unfold.jl | 143 +++++++++++++++++++++++++++++-------- test/runtests.jl | 5 +- 3 files changed, 117 insertions(+), 35 deletions(-) diff --git a/src/coords_optimization.jl b/src/coords_optimization.jl index a44b81a..6ccd01d 100644 --- a/src/coords_optimization.jl +++ b/src/coords_optimization.jl @@ -1,13 +1,13 @@ # Optimization coordinates -function opt(known_pts, known_unf, to_unf, search, neigh, guess=nothing) +function opt(known_pts, known_unf, to_unf, search, neigh, guess) idxs, dists = get_neighbors(known_pts, to_unf, search, neigh, true) out_coords = zeros(Float64, size(to_unf)) Threads.@threads for i in 1:length(idxs) locs = view(known_unf, :, idxs[i]) - initguess = isnothing(guess) ? known_unf[:,idxs[i][1]] : guess[:,i] + initguess = guess[:,i] opt = optimize(x->mse_coords(x, locs, dists[i]), initguess) res = Optim.minimizer(opt) out_coords[:,i] .= res diff --git a/src/unfold.jl b/src/unfold.jl index 66ed947..88151db 100644 --- a/src/unfold.jl +++ b/src/unfold.jl @@ -1,15 +1,14 @@ """ - unfold(ref_pts, domain, samps=nothing; isomap=:default, optim=:default) + unfold(ref, to_unf, samps=nothing; isomap=:default, optim=:default) -Unfold the input points based on the reference points informed. Returns -a coordinate matrix with the unfolded domain points. Or a tuple of two -matrices (unfolded domain and unfolded samples points). +Unfold the input points `to_unf` based on the reference points informed in `ref`. +Returns a coordinate matrix with the unfolded points. Or a tuple of matrices if +`to_unf` is a list of points. ## Parameters: -* `ref_pts` - coordinate matrix with the reference points for unfolding -* `domain - coordinate matrix with domain points for unfolding -* `samps` - coordinate matrix of the sample points (optional) +* `to_unf` - coordinate matrix with the points that will be unfolded +* `ref` - coordinate matrix with the reference points for unfolding * `isomap` - additional parameters for isomap step; either :default or a NamedTuple with the keys shown below * `optim` - additional parameters for optimization step; either :default or a @@ -35,8 +34,8 @@ matrices (unfolded domain and unfolded samples points). * `maxerr` - the maximum accepted distortion for unfolding optimization. * `nchunks` - number of rounds of optimization. """ -function unfold(ref_pts::AbstractMatrix, domain::AbstractMatrix, samps=nothing; - isomap=:default, optim=:default) +function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix; + isomap=:default, optim=:default) # read isomap and optim parameters or assign the defaults below ipars = (search=:knn, neigh=16, anchors=1500, reftol=0.01) @@ -46,41 +45,40 @@ function unfold(ref_pts::AbstractMatrix, domain::AbstractMatrix, samps=nothing; opars = updatepars(opars, optim) # conversions if necessary - !(ref_pts[1] isa Float64) && (ref_pts = Float64.(ref_pts)) - !(domain[1] isa Float64) && (domain = Float64.(domain)) - !isnothing(samps) && !(samps[1] isa Float64) && (samps = Float64.(samps)) + !(ref[1] isa Float64) && (ref = Float64.(ref)) + !(to_unf[1] isa Float64) && (to_unf = Float64.(to_unf)) # random seed Random.seed!(1234567890) # pre-process reference points - ref_pts = remove_duplicates(ref_pts, tol=ipars.reftol) - resol = get_resolution(ref_pts) + ref = remove_duplicates(ref, tol=ipars.reftol) + resol = get_resolution(ref) # do landmark isomap at reference points - unf_ref = landmark_isomap(ref_pts, ipars.search, ipars.neigh, ipars.anchors) + unf_ref = landmark_isomap(ref, ipars.search, ipars.neigh, ipars.anchors) # get initial guess - normals, good = getnormals(ref_pts, ipars.search, ipars.neigh) - unf_dom = firstguess(domain, view(ref_pts,:,good), view(unf_ref,:,good), view(normals,:,good)) + normals, good = getnormals(ref, ipars.search, ipars.neigh) + unf = firstguess(to_unf, view(ref,:,good), view(unf_ref,:,good), view(normals,:,good)) # unfold points in random chunks - ndom = size(domain,2) + ndom = size(to_unf,2) shuffled_ids = shuffle(1:ndom) nchunks = ceil(Int, ndom/opars.nchunks) chunks = collect(Iterators.partition(shuffled_ids, nchunks)) for (i, ids) in enumerate(chunks) # only reference surface as conditioner - known_org = ref_pts[:,good] + known_org = ref[:,good] known_unf = unf_ref[:,good] ids_to_unf = ids # add unfolded from previous chunks if i > 1 extra_ids = vcat(chunks[1:(i-1)]...) - extra_org = view(domain,:,extra_ids) - extra_unf = view(unf_dom,:,extra_ids) + extra_org = view(to_unf,:,extra_ids) + extra_unf = view(unf,:,extra_ids) good_, bad_ = error_ids(extra_org, extra_unf, opars.search, opars.neigh, opars.maxerr) @@ -90,22 +88,105 @@ function unfold(ref_pts::AbstractMatrix, domain::AbstractMatrix, samps=nothing; end # unfold points - initguess = view(unf_dom, :, ids_to_unf) - to_unf = view(domain, :, ids_to_unf) - unf_dom[:,ids_to_unf] .= opt(known_org, known_unf, to_unf, + initguess = view(unf, :, ids_to_unf) + chunk_to_unf = view(to_unf, :, ids_to_unf) + unf[:,ids_to_unf] .= opt(known_org, known_unf, chunk_to_unf, opars.search, opars.neigh, initguess) end - # unfold samples if informed; otherwise, return just the domain unfolded - if !isnothing(samps) - initguess = firstguess(samps, domain, unf_dom) - unf_samps = opt(domain, unf_dom, samps, opars.search, opars.neigh, initguess) - unf_dom, unf_samps - else - unf_dom + unf +end + +function unfold(to_unf::AbstractVector, ref::AbstractMatrix; + isomap=:default, optim=:default) + + merged_to_unf = reduce(hcat, to_unf) + unf = unfold(merged_to_unf, ref, isomap=isomap, optim=optim) + j_ = cumsum(size.(to_unf, 2)) + i_ = [i == 1 ? 1 : j_[i-1]+1 for i in 1:length(j_)] + [unf[:,i:j] for (i,j) in zip(i_,j_)] +end + +""" + unfold(to_unf, ref, unf_ref; optim=:default) + +Unfold the input points `to_unf` based on points already unfolded `unf_ref` and +their correspondent in the original space `ref`. Useful to unfold samples after +the domain is already unfolded, for example. Returns a coordinate matrix with +the unfolded points. Or a tuple of matrices if `to_unf` is a list of points. + +## Parameters: + +* `to_unf` - coordinate matrix with the points that will be unfolded +* `ref` - coordinate matrix with the reference points for unfolding +* `unf_ref` - coordinate matrix with the unfolded reference points +* `optim` - additional parameters for optimization step; either :default or a + NamedTuple with the keys shown below + +### `optim` parameters: + +*Default:* optim = (search=:knn, neigh=16, maxerr=5, nchunks=1) +* `search` - search type to optimize locally (:knn for k-nearest neighbor or + :radius for radius search). +* `neigh` - number of neighbors (for search=:knn) or radius distance + (for search=:radius) to local optimization. +* `maxerr` - the maximum accepted distortion for unfolding optimization. +* `nchunks` - number of rounds of optimization. +""" +function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix, unf_ref::AbstractMatrix; + optim=:default) + + # read optim parameters or assign the defaults below + opars = (search=:knn, neigh=16, maxerr=5, nchunks=2) + opars = updatepars(opars, optim) + + # conversions if necessary + !(to_unf[1] isa Float64) && (to_unf = Float64.(to_unf)) + + # unfold points in random chunks + ndom = size(to_unf,2) + shuffled_ids = shuffle(1:ndom) + nchunks = ceil(Int, ndom/opars.nchunks) + chunks = collect(Iterators.partition(shuffled_ids, nchunks)) + unf = firstguess(to_unf, ref, unf_ref) + + for (i, ids) in enumerate(chunks) + # only reference surface as conditioner + known_org = ref + known_unf = unf_ref + ids_to_unf = ids + + # add unfolded from previous chunks + if i > 1 + extra_ids = vcat(chunks[1:(i-1)]...) + extra_org = view(to_unf,:,extra_ids) + extra_unf = view(unf,:,extra_ids) + good_, bad_ = error_ids(extra_org, extra_unf, opars.search, + opars.neigh, opars.maxerr) + + known_org = hcat(known_org, view(extra_org, :, good_)) + known_unf = hcat(known_unf, view(extra_unf, :, good_)) + ids_to_unf = Int.(union(ids, view(extra_ids, bad_))) + end + + # unfold points + initguess = view(unf, :, ids_to_unf) + chunk_to_unf = view(to_unf, :, ids_to_unf) + unf[:,ids_to_unf] .= opt(known_org, known_unf, chunk_to_unf, + opars.search, opars.neigh, initguess) end + unf end +function unfold(to_unf::AbstractVector, ref::AbstractMatrix, + unf_ref::AbstractMatrix; optim=:default) + + merged_to_unf = reduce(hcat, to_unf) + unf = unfold(merged_to_unf, ref, unf_ref, optim=optim) + j_ = cumsum(size.(to_unf, 2)) + i_ = [i == 1 ? 1 : j_[i-1]+1 for i in 1:length(j_)] + [unf[:,i:j] for (i,j) in zip(i_,j_)] +end function updatepars(default, newpars) if newpars != :default diff --git a/test/runtests.jl b/test/runtests.jl index caf5649..b40fc68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,8 +18,9 @@ using Test println("- Reference surface extracted") # Get transformed coordinates of blocks and samples after unfolding - optim = (neigh=15,) - unf_block, unf_samps = unfold(ref_surface, input_block, input_samps, optim=optim) + optim = (neigh=70,) + unf_block = unfold(input_block, ref_surface, optim=optim) + unf_samps = unfold(input_samps, input_block, unf_block, optim=optim) println("- Unfolding finished") # Write new XT, YT and ZT columns with the transformed coordinates From 391727160a66d7e715e338ed4446f35217c7ebff Mon Sep 17 00:00:00 2001 From: Rafael Caixeta <8386288+rmcaixeta@users.noreply.github.com> Date: Sun, 25 Apr 2021 13:12:25 -0300 Subject: [PATCH 2/3] Adjust README and tests --- Project.toml | 2 +- README.md | 16 ++++++++-------- src/unfold.jl | 8 ++++---- test/runtests.jl | 4 ++++ 4 files changed, 17 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index a784a68..79743b2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Unfolding" uuid = "c5d2ca1b-4127-462e-b5fa-7f6e5ca16f83" authors = ["Rafael Moniz Caixeta"] -version = "0.2.0" +version = "0.2.1" [deps] Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" diff --git a/README.md b/README.md index 880ba31..14b378e 100644 --- a/README.md +++ b/README.md @@ -52,7 +52,7 @@ input_samp = coords( df_samp, columns=["X","Y","Z"] ) # Get reference surface points for unfolding ref_surface = getreference(input_block) # Get transformed coordinates of blocks and samples after unfolding -unf_block, unf_samp = unfold(ref_surface,input_block,input_samp) +unf_block, unf_samp = unfold([input_block,input_samp], ref_surface) # Write new XT, YT and ZT columns with the transformed coordinates for (i,c) in enumerate([:XT,:YT,:ZT]) @@ -61,12 +61,12 @@ for (i,c) in enumerate([:XT,:YT,:ZT]) end # Write output to CSV -CSV.write( "out_dh.csv", df_samp ) -CSV.write( "out_blks.csv", df_block ) +CSV.write("out_dh.csv", df_samp ) +CSV.write("out_blks.csv", df_block ) # Write output to VTK format -to_vtk(unf_block,"out_blks") -to_vtk(unf_samp,"out_dh") +to_vtk("out_blks", unf_block) +to_vtk("out_dh", unf_samp) ``` The code can be saved in a textfile with `.jl` extension and be called in a terminal: `julia file.jl` or `julia -t 4 file.jl` to run faster using 4 threads (or any number of threads you want). Or you can organize it in notebooks. @@ -104,7 +104,7 @@ input_samp = df_samp.to_numpy().T # Get reference surface points for unfolding ref_surface = unf.getreference(input_block) # Get transformed coordinates of blocks and samples after unfolding -unf_block, unf_samp = unf.unfold(ref_surface,input_block,input_samp) +unf_block, unf_samp = unf.unfold([input_block,input_samp], ref_surface) # Write new XT, YT and ZT columns with the transformed coordinates for i,c in enumerate(["XT","YT","ZT"]): @@ -116,8 +116,8 @@ df_block.to_csv("out_blks.csv",index=False) df_samp.to_csv("out_samp.csv",index=False) # Write output to VTK format -unf.to_vtk(unf_block,"out_blks") -unf.to_vtk(unf_samp,"out_dh") +unf.to_vtk("out_blks", unf_block) +unf.to_vtk("out_dh", unf_samp) ``` ### Results diff --git a/src/unfold.jl b/src/unfold.jl index 88151db..4b60c0f 100644 --- a/src/unfold.jl +++ b/src/unfold.jl @@ -26,7 +26,7 @@ Returns a coordinate matrix with the unfolded points. Or a tuple of matrices if ### `optim` parameters: -*Default:* optim = (search=:knn, neigh=16, maxerr=5, nchunks=4) +*Default:* optim = (search=:knn, neigh=50, maxerr=5, nchunks=4) * `search` - search type to optimize locally (:knn for k-nearest neighbor or :radius for radius search). * `neigh` - number of neighbors (for search=:knn) or radius distance @@ -39,7 +39,7 @@ function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix; # read isomap and optim parameters or assign the defaults below ipars = (search=:knn, neigh=16, anchors=1500, reftol=0.01) - opars = (search=:knn, neigh=16, maxerr=5, nchunks=4) + opars = (search=:knn, neigh=50, maxerr=5, nchunks=4) ipars = updatepars(ipars, isomap) opars = updatepars(opars, optim) @@ -125,7 +125,7 @@ the unfolded points. Or a tuple of matrices if `to_unf` is a list of points. ### `optim` parameters: -*Default:* optim = (search=:knn, neigh=16, maxerr=5, nchunks=1) +*Default:* optim = (search=:knn, neigh=50, maxerr=5, nchunks=1) * `search` - search type to optimize locally (:knn for k-nearest neighbor or :radius for radius search). * `neigh` - number of neighbors (for search=:knn) or radius distance @@ -137,7 +137,7 @@ function unfold(to_unf::AbstractMatrix, ref::AbstractMatrix, unf_ref::AbstractMa optim=:default) # read optim parameters or assign the defaults below - opars = (search=:knn, neigh=16, maxerr=5, nchunks=2) + opars = (search=:knn, neigh=50, maxerr=5, nchunks=2) opars = updatepars(opars, optim) # conversions if necessary diff --git a/test/runtests.jl b/test/runtests.jl index b40fc68..a3b4543 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,4 +37,8 @@ using Test good, bad = errors(:ids, all_input, all_unf) @test length(good) > (30 * length(bad)) + + println("- Quick test of unfold with multiple inputs") + s1, s2 = unfold([input_samps[:,1:80],input_samps[:,80:end]], ref_surface) + s1, s2 = unfold([input_samps[:,1:80],input_samps[:,80:end]], input_block, unf_block) end From d64aede32a426ca297be7adb91612488453450ff Mon Sep 17 00:00:00 2001 From: Rafael Caixeta <8386288+rmcaixeta@users.noreply.github.com> Date: Sun, 25 Apr 2021 13:37:58 -0300 Subject: [PATCH 3/3] Clip max neighbors if n < k --- src/data_handling.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/data_handling.jl b/src/data_handling.jl index 1e903a7..d66a810 100644 --- a/src/data_handling.jl +++ b/src/data_handling.jl @@ -45,7 +45,8 @@ function get_neighbors(origin::AbstractMatrix, target::AbstractMatrix, @assert nhood in [:knn,:radius] "invalid neighborhood type" tree = KDTree(origin) if nhood==:knn - idxs, dists = knn(tree, target, neigh_val, true) + nneigh = minimum([neigh_val, size(origin,2)]) + idxs, dists = knn(tree, target, nneigh, true) idxs, dists elseif nhood==:radius idxs = inrange(tree, target, neigh_val)