Merge pull request #5 from zygmuntszpak/demo_fundamental

Adds demonstration of fundamental matrix estimation

demo/example_fundamental_matrix.jl

@@ -0,0 +1,65 @@
+using MultipleViewGeometry, Images
+using BenchmarkTools, Compat
+using StaticArrays, MAT, Plots
+
+# Load MATLAB matrices that represent a pair of images and that contain
+# a set of manually matched corresponding points.
+file = matopen("./data/sene-adelaideRMF.mat")
+img1 = colorview(RGB, normedview(permutedims(read(file,"img1"),[3,1,2])))
+img2 = colorview(RGB, normedview(permutedims(read(file,"img2"),[3,1,2])))
+inlier_pts1 = read(file,"inlierPts1")
+inlier_pts2 = read(file,"inlierPts2")
+close(file)
+
+# Conver the images to Grayscale. 
+img1g = Gray.(img1)
+img2g = Gray.(img2)
+
+plotlyjs()
+
+_, npts1 = size(inlier_pts1)
+_, npts2 = size(inlier_pts2)
+
+# Construct a set of corresponding points expressed in homogeneous coordinates.
+ℳ = [ Point2DH(vcat(inlier_pts1[:,n],1)) for n in 1:npts1]
+ℳʹ =  [ Point2DH(vcat(inlier_pts2[:,n],1)) for n in 1:npts2]
+
+# Estimate the Fundamental matrix from the corresponding points using the
+# Direct Linear Transform algorithm.
+𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(),  (ℳ, ℳʹ))
+
+# Plot Keypoints, epipole and concomitant epipolar lines in the first view.
+p1 = Plots.plot(img1g,grid = false, box = :none, legend = false, size = (455,341))
+for n = 1:25:npts1
+    m = ℳ[n]
+    # Epipolar line in the first image.
+    l = 𝑛(𝐅₀'*ℳʹ[n])
+    draw!(EpipolarLineGraphic(), l, size(img1), p1)
+    # Keypoint.
+    Plots.plot!([m[1]],[m[2]], grid = false, box = :none, legend = false,
+                    seriestype = :scatter, w = 5, aspect_ratio = :equal)
+end
+e = epipole(𝐅₀)
+Plots.plot!([e[1]],[e[2]], grid = false, box = :none, legend = false,
+            seriestype = :scatter, markershape = :diamond, markerstrokecolor = :red,
+            markersize = 5, aspect_ratio = :equal)
+
+# Plot Keypoints and concomitant epipolar lines in the second view.
+p2 = Plots.plot(img2g,grid = false, box = :none, legend = false)
+for n = 1:25:npts2
+    mʹ = ℳʹ[n]
+    # Epipolar line in the second image.
+    l = 𝑛(𝐅₀*ℳ[n])
+    draw!(EpipolarLineGraphic(), l, size(img2), p2)
+    Plots.plot!([mʹ[1]],[mʹ[2]], grid = false, box = :none, legend = false,
+                    seriestype = :scatter, w = 5, aspect_ratio = :equal)
+end
+eʹ = epipole(𝐅₀')
+Plots.plot!([eʹ[1]],[eʹ[2]], grid = false, box = :none, legend = false,
+            seriestype = :scatter, markershape = :diamond, markerstrokecolor = :red,
+            markersize = 5, aspect_ratio = :equal)
+
+
+# Display both views simultaneously.
+p3 = Plots.plot(p1,p2,layout=(1,2), legend = false)
+display(p3)

src/MultipleViewGeometry.jl

@@ -42,14 +42,17 @@ export construct
 export construct
 
 # Functions exported from `project.jl`
-export project
+export project, epipole
 
 # Functions exported from `rotations.jl`
 export rotx, roty, rotz, rotxyz, rodrigues2matrix
 
 # Functions exported from `cost_functions.jl`
 export cost, X, covariance_matrix, covariance_matrix_debug
 
+# Functions exported from `draw.jl`
+export draw!, EpipolarLineGraphic
+
 include("math_aliases/ModuleMathAliases.jl")
 include("types/ModuleTypes.jl")
 include("operators/ModuleOperators.jl")
@@ -62,6 +65,7 @@ include("moments/ModuleMoments.jl")
 include("cost_function/ModuleCostFunction.jl")
 include("estimate/ModuleEstimation.jl")
 include("construct/ModuleConstruct.jl")
+include("draw/ModuleDraw.jl")
 
 
 using .ModuleMathAliases
@@ -76,6 +80,7 @@ using .ModuleEstimation
 using .ModuleMoments
 using .ModuleCostFunction
 using .ModuleConstruct
+using .ModuleDraw
 
 
 # package code goes here

src/draw/ModuleDraw.jl

@@ -0,0 +1,9 @@
+module ModuleDraw
+using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleOperators, MultipleViewGeometry.ModuleMathAliases
+using StaticArrays
+using Plots, Plotly
+using Juno
+
+export draw!, EpipolarLineGraphic
+include("draw.jl")
+end

src/draw/draw.jl

@@ -0,0 +1,48 @@
+abstract type GraphicEntity end
+
+type EpipolarLineGraphic <: GraphicEntity
+end
+
+
+function draw!(g::EpipolarLineGraphic, l::AbstractVector, dim::Tuple{<:Number,<:Number}, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend})
+
+    top = intersection(l, [0 ; -1  ; 0])
+    bottom = intersection(l, [0 ; -1 ; dim[1]])
+    left = intersection(l, [-1 ; 0 ; 0])
+    right = intersection(l, [-1 ; 0 ; dim[2]])
+
+    x = Float64[]
+    y = Float64[]
+
+    if is_inbounds(top,dim)
+        push!(x, top[1])
+        push!(y, top[2])
+    end
+
+    if is_inbounds(bottom,dim)
+        push!(x, bottom[1])
+        push!(y, bottom[2])
+    end
+
+    if is_inbounds(left,dim)
+        push!(x, left[1])
+        push!(y, left[2])
+    end
+
+    if is_inbounds(right,dim)
+        push!(x, right[1])
+        push!(y, right[2])
+    end
+
+    Plots.plot!(x,y,w=3)
+end
+
+function intersection(l1::AbstractArray, l2::AbstractArray)
+    l = 𝑛(cross(l1,l2))
+    l[1:2]
+end
+
+function is_inbounds(pt::AbstractVector, dim::Tuple{<:Number,<:Number})
+    nrow, ncol = dim
+    pt[1] >= -1.5 && pt[1] < ncol+1.5 && pt[2] >= -1.5 && pt[2] <= nrow + 1.5
+end

src/projection/epipole.jl

@@ -1,5 +1,5 @@
 function epipole(𝐅::AbstractArray)
     𝐔, 𝐒, 𝐕 = svd(𝐅)
-    𝐞 = 𝑛(MVector(𝐕[:,end]))
+    𝐞 = 𝑛(MVector{3}(𝐕[:,end]))
     Point2DH(𝐞)
 end