In [212]:
using PyPlot
In [211]:
include("JuliaDe3DFitting.jl")
Out[211]:
In [221]:
function makeSphere(a, b, c, rin, n)
xs = Float64[]
ys = Float64[]
zs = Float64[]
for i in 0:n
z = rin * sin(pi*(i-(n/2))/n) + c
r = rin * cos(pi*(i-(n/2))/n)
for j in 0:2*n
x = r * sin(2*pi*j/(2*n)) + a
y = r * cos(2*pi*j/(2*n)) + b
push!(xs, x)
push!(ys, y)
push!(zs, z)
end
end
return hcat(xs, ys, zs)
end
function makeNoisedData(X)
divider = 5
sampleSize = size(X)[1]
noised_xs = X[:, 1] - ((rand(sampleSize) - 0.5) / divider)
noised_ys = X[:, 2] - ((rand(sampleSize) - 0.5) / divider)
noised_zs = X[:, 3] - ((rand(sampleSize) - 0.5) / divider)
return hcat(noised_xs, noised_ys, noised_zs)
end
function plots(X)
for i in 1:3
for j in 1:2
elev = 20 * i
azim = 20 * j
fig = PyPlot.figure(figsize=(5, 5))
ax = PyPlot.Axes3D(fig, elev=elev, azim=azim)
PyPlot.scatter3D(X[:, 1], X[:, 2], X[:, 3], s=1, color="b")
end
end
end
function plots(X1, X2)
for i in 1:3
for j in 1:2
elev = 20 * i
azim = 20 * j
fig = PyPlot.figure(figsize=(5, 5))
ax = PyPlot.Axes3D(fig, elev=elev, azim=azim)
PyPlot.scatter3D(X1[:, 1], X1[:, 2], X1[:, 3], s=1, color="y")
PyPlot.scatter3D(X2[:, 1], X2[:, 2], X2[:, 3], s=1, color="b")
end
end
end
function makePlane(a, b, c)
function plane(x, y, a, b, c)
return a * x + b * y + c
end
data = zeros(100, 3)
for i in 1:10
for j in 1:10
z = plane(i/10.0, j/10.0, a, b, c)
data[10*(i-1)+j, :] = [i/10.0, j/10.0, z]
end
end
return data
end
Out[221]:
In [225]:
plots(makeNoisedData(makePlane(3, 2, 1)))
In [226]:
noisedPlane = makeNoisedData(makePlane(3, 2, 1))
result = planeFit(noisedPlane)
Out[226]:
In [227]:
plots(noisedPlane, makePlane(result[1][1], result[1][2], result[1][3]))
In [214]:
X = makeSphere(1, 1, 1, 2, 30)
Out[214]:
In [215]:
noisedX = makeNoisedData(X)
plots(noisedX)
In [216]:
result = sphereFit(noisedX, 100)
Out[216]:
In [217]:
fittedX = makeSphere(result[1], result[2], result[3], result[4], 60)
plots(noisedX, fittedX)
In [218]:
function makeCylinder(W, C, r, theta0, theta1)
sin0 = sin(theta0)
cos0 = cos(theta0)
sin1 = sin(theta1)
cos1 = cos(theta1)
orthVec1 = vec([-sin0, cos0, 0.0])
orthVec2 = vec([-cos0*cos1, -sin0*cos1, sin1])
orthVec3 = -orthVec1
orthVec4 = -orthVec2
function pointsOnCylinder(V1, V2, W, C, r)
X = [0 0 0]
for i in 1:10
VecOnCircle = i * V1 + (10 - i) * V2
VecOnCircle /= sqrt(dot(VecOnCircle, VecOnCircle))
VecOnCircle *= r
VecOnCircle += C
for j in -10:10
if j == 0
X = vcat(X, [VecOnCircle[1] VecOnCircle[2] VecOnCircle[3]])
else
VecOnCylinder = VecOnCircle + (j * W)/10.0
X = vcat(X, [VecOnCylinder[1] VecOnCylinder[2] VecOnCylinder[3]])
end
end
end
return X
end
X1 = pointsOnCylinder(orthVec1, orthVec2, W, C, r)
X2 = pointsOnCylinder(orthVec2, orthVec3, W, C, r)
X3 = pointsOnCylinder(orthVec3, orthVec4, W, C, r)
X4 = pointsOnCylinder(orthVec4, orthVec1, W, C, r)
X = vcat(X1, X2, X3, X4)
return X[1:end, :]
end
Out[218]:
In [219]:
Wcor = [cos(0.2)*sin(0.2), sin(0.2)*sin(0.2), cos(0.2)]
Xcyl = makeCylinder(Wcor, [0, 0, 0], 1, 0.2, 0.2)
Out[219]:
In [222]:
noisedXcyl = makeNoisedData(Xcyl)
plots(noisedXcyl)
In [223]:
W, C, r, theta0, theta1 = cylinderFit(noisedXcyl)
Out[223]:
In [199]:
Wcor
Out[199]:
In [224]:
plots(noisedXcyl, makeCylinder(W, C, r, theta0, theta1))
In [173]:
plots(makeNoisedData(makePlane(3, 2, 1)))
In [184]:
noisedPlane = makeNoisedData(makePlane(3, 2, 1))
result = planeFit(noisedPlane)
Out[184]:
In [190]:
plots(noisedPlane, makePlane(result[1][1], result[1][2], result[1][3]))
In [231]:
X
Out[231]:
In [ ]: