Merge pull request #32 from rmcaixeta/dev

Improve unfold methods and correct neighborhood bug for few data

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"

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

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

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)

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
@@ -27,60 +26,59 @@ matrices (unfolded domain and unfolded samples points).
 
 ### `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
   (for search=:radius) to local optimization.
 * `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)
-	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)
 
 	# 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=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
+  (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=50, 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

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
@@ -36,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