Commit 76d7be01 authored by Jayant Khatkar's avatar Jayant Khatkar
Browse files

calculate contour segments in a voxel (tested for most cases)

parent e698be2c
......@@ -3,6 +3,7 @@ using CSV
using JSON
using LightGraphs
using NearestNeighbors
using Statistics
struct contour
......@@ -85,9 +86,9 @@ end
struct voxmap
seglen::Array{Number}
segoffset::Array{Number}
segcontours::Array{Int}
seglen::Vector{Float64}
segoffset::Vector{Float64}
segcontours::Vector{Int}
end
......@@ -122,19 +123,31 @@ function seg_intersect(p1,q1,p2,q2)
return false
end
dist(p1, p2) = √sum((p1 -p2).^2)
interpolate(p1, p2, xi, axis) = p1 + (p2-p1)*(xi-p1[axis])/(p2[axis]-p1[axis])
interpolate(p1, p2, x1, xi, x2) = p1 + (p2-p1)*(xi-x1)/(x2-x1)
function voxmap(vox, vox_d, problem)
function voxmap(vox, vox_d, prob::problem)
# for one vox, get all contours which pass through it
# only need to search contours in its layer
l = Int((vox[3] + problem.layer_height/2)/problem.layer_height)
l = Int((vox[3] + prob.layer_height/2)/prob.layer_height)
voxx1 = vox[1] + vox_d/2
voxx2 = vox[1] - vox_d/2
voxy1 = vox[2] + vox_d/2
voxy2 = vox[2] - vox_d/2
seg_now = false
seglen = Vector{Number}()
segoffset = Vector{Number}()
segcontours = Vector{Int}()
seglen_sofar = 0
t_start = 0
for cid in prob.layers[l]
c = prob.contours[cid]
for c in problem.layers[l]:
# check if contour passes thorough this vox
for i in 2:size(c.pos)[1]
......@@ -142,41 +155,154 @@ function voxmap(vox, vox_d, problem)
if c.pos[i, 1] > voxx1 && c.pos[i-1, 1] > voxx1 ||
c.pos[i,1] < voxx2 && c.pos[i-1, 1] < voxx2 ||
c.pos[i,2] < voxy2 && c.pos[i-1, 2] < voxy2 ||
c.pos[i,2] < voxy1 && c.pos[i-1, 2] > voxy1
c.pos[i,2] > voxy1 && c.pos[i-1, 2] > voxy1
# segment outside loop entirely
# segment outside vox entirely
if seg_now
print("Something's gone wrong: segment entirely outside voxel, but last segment inside")
end
continue
end
if c.pos[i, 1] > voxx1 && c.pos[i-1, 1] > voxx1 ||
c.pos[i,1] < voxx2 && c.pos[i-1, 1] < voxx2 ||
c.pos[i,2] < voxy2 && c.pos[i-1, 2] < voxy2 ||
c.pos[i,2] < voxy1 && c.pos[i-1, 2] > voxy1
p1inside = c.pos[i-1, 1] < voxx1 && c.pos[i-1, 1] > voxx2 && c.pos[i-1, 2] > voxy2 && c.pos[i-1, 2] < voxy1
p2inside = c.pos[i, 1] < voxx1 && c.pos[i,1] > voxx2 && c.pos[i,2] > voxy2 && c.pos[i,2] < voxy1
# is this line segment completely inside vox?
if p1inside && p2inside
# segment outside loop entirely
if seg_now
# TODO add this to segment
seglen_sofar += dist(c.pos[i], c.pos[i-1]) # append to existing contour
if !seg_now # start new seg
t_start = 0 # 0 bc contour must be starting for this case
seg_now = true
if i!=2
print("Whole segment inside but something wrong")
end
continue
end
end
cross_side1 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx1, voxy1], [voxx1, voxy2])
cross_side2 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx1, voxy1], [voxx2, voxy1])
cross_side3 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx2, voxy1], [voxx2, voxy2])
cross_side4 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx2, voxy2], [voxx1, voxy2])
println(cid)
println([cross_side1, cross_side2, cross_side3, cross_side4])
# does this line segment intersect with vox only once
if p1inside p2inside
# find intersection point
if cross_side1 || cross_side3
# intersection with x
xi = [voxx1, voxx2][[cross_side1, cross_side3]][1]
p_i = interpolate(c.pos[i-1,:], c.pos[i,:], xi, 1)
t_i = interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], xi, c.pos[i,1])
elseif cross_side3 || cross_side4
# intersection with y
yi = [voxy1, voxy2][[cross_side2, cross_side4]][1]
p_i = interpolate(c.pos[i-1,:], c.pos[i,:], yi, 2)
t_i = interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], yi, c.pos[i,2])
end
if p1inside
# end existing segment
if !seg_now
print("Something's gone wrong")
end
seglen_sofar += dist(c.pos[i-1, :], p_i)
push!(segcontours, cid)
push!(seglen, seglen_sofar)
push!(segoffset, (t_i + t_start)/2)
seglen_sofar = 0
seg_now = false
else
# TODO START NEW SEGMENT
# start new contour
t_start = t_i
seglen_sofar = dist(p_i, c.pos[i, :])
seg_now = true
end
continue
elseif sum([cross_side1, cross_side2, cross_side3, cross_side4]) >= 2
# intersects twice
p_is = []
t_is = []
if cross_side1
push!(p_is, interpolate(c.pos[i-1,:], c.pos[i,:], voxx1, 1))
push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], voxx1, c.pos[i,1]))
end
if cross_side2
push!(p_is, interpolate(c.pos[i-1,:], c.pos[i,:], voxy1, 2))
push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], voxy1, c.pos[i,2]))
end
if cross_side3
push!(p_is, interpolate(c.pos[i-1,:], c.pos[i,:], voxx2, 1))
push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], voxx2, c.pos[i,1]))
end
if cross_side4
push!(p_is, interpolate(c.pos[i-1,:], c.pos[i,:], voxy2, 2))
push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], voxy2, c.pos[i,2]))
end
if seg_now
print("Something's wrong")
end
# check line segment i and i-1
elseif seg_intersect(c.pos[i,:], c.pos[i,:], [voxx1, voxy1], [voxx1, voxy2]) ||
seg_intersect(c.pos[i,:], c.pos[i,:], [voxx1, voxy1], [voxx2, voxy1]) ||
seg_intersect(c.pos[i,:], c.pos[i,:], [voxx2, voxy1], [voxx2, voxy2]) ||
seg_intersect(c.pos[i,:], c.pos[i,:], [voxx2, voxy2], [voxx1, voxy2])
push!(segoffset, mean(t_is))
push!(segcontours, cid)
n_intersections = sum([cross_side1, cross_side2, cross_side3, cross_side4])
if n_intersections ==2
push!(seglen, dist(p_is[1], p_is[2]))
else
push!(seglen, dist(p_is[1], p_is[3]))
end
end
print("intersecting")
end
# if contour ends inside the voxel
if seg_now
# end segment
push!(segcontours, cid)
push!(seglen, seglen_sofar)
push!(segoffset, (t_start + last(c.time))/2)
seglen_sofar = 0
t_start = 0
seg_now = false
end
# for those contours find exact segments
end
return voxmap(seglen, segoffset, segcontours)
end
voxels = DataFrames.DataFrame(CSV.File("tensile-1-1.csv"))
contours = contour.(JSON.parse(open("tensilecontours.json")))
dt = problem(contours, 5, 5)
#voxels = DataFrames.DataFrame(CSV.File("tensile-1-1.csv"))
#contours = contour.(JSON.parse(open("tensilecontours.json")))
#dt = problem(contours, 5, 5)
function test_voxmap()
# create vox
vox = [0,0,0.5]
vox_d = 2
pos1 = [[0.5, -0.5] ones(2)*0.5 ones(2)]
time1 = [0,1]
pos2 = [[1.5, 0.5, -0.5, -1.5] ones(4)*0.5 ones(4)]
time2 = Vector(0:3)
pos3 = [[-0.5, -0.5] [2, -2] ones(2)]
time3 = [0, 2.5]
pos4 = [[0.5, 2.5] [-1.5, -0.5] ones(2)]
time4 = [0,1]
pos5 = [[-2,2] [-2,2] ones(2)]
time5 = [0,1]
contour1 = contour(pos1, time1)
contour2 = contour(pos2,time2)
contour3 = contour(pos3,time3)
contour4 = contour(pos4,time4)
contour5 = contour(pos5,time5)
contours = [contour1, contour2, contour3, contour4, contour5]
prob = problem(contours, 1, 1)
vm = voxmap(vox, vox_d, prob)
return vm
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment