Update to v1.0.0 (#13)

* Breaking change to NLLSOptions
* Generalize lens distortion conversion function
* Several bug fixes and stability improvements

Project.toml

@@ -1,7 +1,7 @@
 name = "NLLSsolver"
 uuid = "50b5cf9e-bf7a-4e29-8981-4280cbba70cb"
 authors = ["Oliver Woodford"]
-version = "0.1.1"
+version = "1.0.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

src/iterators.jl

@@ -1,10 +1,14 @@
 using LDLFactorizations
 export iterate!, NewtonData, DoglegData, LevMarData
 
-function symmetricsolve(A::AbstractMatrix, b::AbstractVector, options)
+function symmetricsolve(A::SparseMatrixCSC, b::AbstractVector, options)
     return ldl(A) \ b
 end
 
+function symmetricsolve(A::AbstractMatrix, b::AbstractVector, options)
+    return Symmetric(A) \ b
+end
+
 function linearsolve(A::AbstractMatrix, b::AbstractVector, options)
     return solve(LinearProblem(A, b)).u
 end
@@ -141,7 +145,7 @@ function iterate!(levmardata::LevMarData, data::NLLSInternal, problem::NLLSProbl
             # Success (or convergence) - update lambda
             uniformscaling!(hessian, -lastlambda)
             step_quality = (cost_ - data.bestcost) / (((data.step' * hessian) * 0.5 + gradient') * data.step)
-            levmardata.lambda *= max(0.1, 1 - (step_quality - 1) ^ 3)
+            levmardata.lambda *= step_quality < 1.966 ? 1 - (step_quality - 1) ^ 3 : 0.1
             # Return the cost
             return cost_
         end

src/optimize.jl

@@ -58,7 +58,7 @@ function optimize!(problem::NLLSProblem{VarTypes}, options::NLLSOptions, data::N
                 push!(trajectory, copy(data.step))
             end
             # Check for convergence
-            if options.callback(problem, data, cost) || !(dcost >= data.bestcost * options.dcost) || (maximum(abs, data.step) < options.dstep) || (fails > options.maxfails) || data.iternum >= options.maxiters
+            if options.callback(problem, data, cost) || !(dcost >= data.bestcost * options.reldcost) || !(dcost >= options.absdcost) || (maximum(abs, data.step) < options.dstep) || (fails > options.maxfails) || data.iternum >= options.maxiters
                 break
             end
             # Update the variables

src/structs.jl

@@ -17,7 +17,8 @@ function Base.String(iterator::NLLSIterator)
 end
 
 struct NLLSOptions
-    dcost::Float64              # Minimum relative reduction in cost required to avoid termination
+    reldcost::Float64           # Minimum relative reduction in cost required to avoid termination
+    absdcost::Float64           # Minimum relative reduction in cost required to avoid termination
     dstep::Float64              # Minimum L-infinity norm of the update vector required to avoid termination
     maxfails::Int               # Maximum number of consecutive iterations that have a higher cost than the current best before termination
     maxiters::Int               # Maximum number of outer iterations
@@ -26,8 +27,8 @@ struct NLLSOptions
     storecosts::Bool            # Indicates whether the cost per outer iteration should be stored
     storetrajectory::Bool       # Indicates whether the step per outer iteration should be stored
 end
-function NLLSOptions(; maxiters=100, dcost=1.e-6, dstep=1.e-6, maxfails=3, iterator=gaussnewton, callback=(args...)->false, storecosts=false, storetrajectory=false)
-    NLLSOptions(dcost, dstep, maxfails, maxiters, iterator, callback, storecosts, storetrajectory)
+function NLLSOptions(; maxiters=100, reldcost=1.e-6, absdcost=1.e-15, dstep=1.e-6, maxfails=3, iterator=levenbergmarquardt, callback=(args...)->false, storecosts=false, storetrajectory=false)
+    NLLSOptions(reldcost, absdcost, dstep, maxfails, maxiters, iterator, callback, storecosts, storetrajectory)
 end
 
 struct NLLSResult
@@ -100,7 +101,7 @@ mutable struct NLLSInternal{VarTypes}
             end
             varlen = UInt(nvars(problem.variables[unfixed]))
             linsystem = UniVariateLS(unfixed, varlen, computehessian ? varlen : UInt(lengthresiduals(problem.residuals)))
-            return new(copy(problem.variables), copy(problem.variables), linsystem, Vector{Float64}(undef, varlen), 0., 0., 0., 0., 0, 0, 0, starttimens)
+            return new(copy(problem.variables), copy(problem.variables), linsystem, Vector{Float64}(undef, varlen), 0., 0., 0., 0., 0, 0, 0, 0, starttimens)
         end
 
         # Multiple variables. Use a block sparse matrix

src/visualgeometry/camera.jl

@@ -1,5 +1,5 @@
-export SimpleCamera, NoDistortionCamera, ExtendedUnifiedCamera, BarrelDistortion
-export ideal2image, image2ideal, pixel2image, image2pixel, update, nvars, nres, varindices, getvars, computeresidual, barrel2eulens
+export SimpleCamera, NoDistortionCamera, ExtendedUnifiedCamera, BarrelDistortion, EULensDistortion
+export ideal2image, image2ideal, pixel2image, image2pixel, ideal2distorted, distorted2ideal, update, nvars, nres, varindices, getvars, computeresidual, convertlens
 using StaticArrays, LinearAlgebra
 
 
@@ -135,24 +135,23 @@ struct LensDistortResidual{T} <: AbstractResidual
 end
 nvars(@objtype(LensDistortResidual)) = 1 # The residual depends on one variable
 nres(@objtype(LensDistortResidual)) = 1 # The residual vector has length one
-varindices(::LensDistortResidual) = [1]
+varindices(::LensDistortResidual) = SVector(1)
 function computeresidual(residual::LensDistortResidual, lens::EULensDistortion)
     return SVector(residual.rdistort - ideal2distorted(lens, residual.rlinear))
 end
 function getvars(::LensDistortResidual{T}, vars::Vector) where T
     return (vars[1]::EULensDistortion{T},)
 end
 
-function barrel2eulens(k1, k2, halfimsz)
+function convertlens(tolens, fromlens, halfimsz)
     # Create an optimization problem to convert the lens distortion
-    problem = NLLSsolver.NLLSProblem{EULensDistortion}()
-    NLLSsolver.addvariable!(problem, EULensDistortion(0.001, 1.))
-    for x in range(0, halfimsz, 100)
-        x2 = x ^ 2
-        NLLSsolver.addresidual!(problem, LensDistortResidual(x, x * (1 + x2 * (k1 + x2 * k2))))
+    problem = NLLSsolver.NLLSProblem{typeof(tolens)}()
+    NLLSsolver.addvariable!(problem, tolens)
+    for x in LinRange(0., convert(Float64, halfimsz), 100)
+        NLLSsolver.addresidual!(problem, LensDistortResidual(x, ideal2distorted(fromlens, x)))
     end
     # Optimize
-    NLLSsolver.optimize!(problem, NLLSsolver.NLLSOptions(dcost=1.e-6, iterator=NLLSsolver.levenbergmarquardt))
+    NLLSsolver.optimize!(problem, NLLSsolver.NLLSOptions(iterator=NLLSsolver.dogleg))
     # Return the lens
     return problem.variables[1]
-end
+end

test/camera.jl

@@ -36,4 +36,14 @@ end
     testcamera(NoDistortionCamera(f, c), x)
     testcamera(ExtendedUnifiedCamera(f, c, 0.1, 0.1), x)
     testcamera(ExtendedUnifiedCamera(f, c, 0.5, 0.2), x)
+
+    # Test lens distortion model conversion
+    halfimsz = 1.
+    fromlens = BarrelDistortion(-1.e-5, -1.e-7)
+    tolens = convertlens(EULensDistortion(0.1, 0.1), fromlens, halfimsz)
+    maxerror = 0.
+    for x in LinRange(0., halfimsz, 100)
+        maxerror = max(maxerror, abs(ideal2distorted(tolens, x) - ideal2distorted(fromlens, x)))
+    end
+    @test maxerror < 1.e-6
 end