I’m afraid I have no idea how to reduce this to a snippet. I’ll post a script which reproduces the error, for anyone kind enough to have a look. Hmm, no code files as attachments. Inline it is then. Apologies for the huge dump. The function findlenscalibration and the closure generaterationallenscostfn contain the things I’m trying to understand.
using JuMP, Parameters
import Base: dot
@with_kw immutable LensCalibSolveParams
a_init = 1.0
b_init = 2*10^7
w_init = 0.5
imagewidth = 1920
imageheight = 1080
maxedgelsperline = 20
end
@with_kw immutable Conic
a = 0.0
b = 0.0
c = 0.0
d = 0.0
e = 0.0
f = 0.0
end
function dot{T<:AbstractArray}(C::Conic, b::T)
return C.a*b[1]+C.b*b[2]+C.c*b[3]+C.d*b[4]+C.e*b[5]+C.f*b[6]
end
function conicdesignmatrix!(matrix, xs, ys)
for r in 1:size(matrix,1)
x = xs[r]
y = ys[r]
matrix[r,1] = x^2
matrix[r,2] = x*y
matrix[r,3] = y^2
matrix[r,4] = x
matrix[r,5] = y
matrix[r,6] = 1.0
end
matrix
end
function fitconic(xs, ys, design_matrix)
conicdesignmatrix!(design_matrix, xs, ys)
U, S, Vs = svd(design_matrix, thin=false)
# The null space is the column corresponding to the smallest singular value.
# The SVD returns the singular values sorted, so this is just the final column of V:
Conic(Vs'[:,6]...)
end
"A closure to create the cost function"
function generaterationallenscostfn(params, linesxs, linesys)
@unpack_LensCalibSolveParams params
# allocate memory
design_matrix = zeros(5,6)
# Fit a conic to each line, then find the dot product of the conic
# with each corresponding lifted point. Store the dot products in an array.
# This needs a 2D array with dimensions NoOfLinesxMaxNoOfEdgelsPerLine
conics = Array(Conic, length(linesxs))
for i in 1:length(linesxs)
conics[i] = fitconic(linesxs[i], linesys[i], design_matrix)
end
theta_dot_chis = zeros(length(linesxs), maxedgelsperline)
for j in 1:length(linesxs), i in 1:min(maxedgelsperline, length(linesxs[j]))
x = linesxs[j][i]
y = linesys[j][i]
chi = [x^2, x*y, y^2, x, y, 1]
theta_dot_chis[j,i] = dot(conics[j], chi)
end
nlines = length(linesxs)
A = zeros(3, 6)
Axx = 0.0
Axy = 0.0
Ax = 0.0
Ayy = 0.0
Ay = 0.0
return function(v) # v = [a, b, w, lines...]
result = 0.0
asquared = v[1]^2
wsquared = v[3]^2
A31 = 4*asquared*wsquared
# Each iteration:
# For each line:
#for (lidx, line) in enumerate(lines)
for linen in 1:nlines
line = v[linen*3+1:linen*3+3] # 4, 7, 10, 13...
Axx = sum(line[1]*1.0, line[2]*1.0, line[3]*A31) #dot(A[:,1], line)
Axy = 0.0#dot(A[:,2], line)
Ayy = sum(line[1]/asquared, line[2]/asquared, line[3]*4*wsquared)#dot(A[:,3], line)
Ax = sum((-2*v[2]-imagewidth)*line[1], -imagewidth*line[2], -A31*imagewidth*line[3])#dot(A[:,4], line)
Ay = sum(line[2]*(-2*v[2]/asquared-imageheight/asquared), line[3]*(-4*wsquared*imageheight))#dot(A[:,5], line)
# Find A_xx, A_xy... for the line
# For each point on the line:
for pidx in 1:nlines
# Take the conic-point product for this point, and divide it by the denom made
# up using this line and A_xx A_xy...
# Accumulate these & return.
x = linesxs[lidx][pidx]
y = linesys[lidx][pidx]
denom = sqrt((2Axx*x+Axy*y+Ax)^2+(2Ayy*y+Axy*x+Ay)^2)
result += (theta_dot_chis[lidx, pidx]/denom)^2
end
end
return result
end
end
function initline(xs, ys)
# Fit a straight line to the points (x,y) in xs, ys to initialise.
A = zeros(2,2)
b = zeros(2)
c = zeros(2)
A[1,1] = length(xs)
# BLASification
A[1,2] = sum(xs)
A[2,2] = BLAS.dot(length(xs), xs, 1, xs, 1)
A[2,1] = A[1,2]
b[1] = sum(ys)
b[2] = BLAS.dot(length(xs), xs, 1, ys, 1)
try
c = inv(A)*b
catch
c = pinv(A)*b
end
[c... 1]
end
# Computes the initial estimates of straight lines for each of the curves.
function init_lines(linesxs, linesys, params)
@unpack_LensCalibSolveParams params
lines = zeros(3,length(linesxs))
for i in 1:length(linesxs)
lines[:,i] = initline(linesxs[i], linesys[i])
end
lines
end
function findlenscalibration(linesxs, linesys, params)
@unpack_LensCalibSolveParams params
nlines = length(linesxs)
m = Model()
# estimate initial values:
lines_init = init_lines(linesxs, linesys, params)
lines_init = reshape(lines_init, 3*nlines)
initvals = [a_init b_init w_init lines_init...]
@variable(m, v[i=1:length(initvals)], start=initvals[i])
# Get a function to compute the cost, given the constant parameters.
rationallenscost = generaterationallenscostfn(params, linesxs, linesys)
JuMP.register(m, :rationallenscost, 1, rationallenscost, autodiff=true)
@NLobjective(m, Min, rationallenscost(v))
status = solve(m)
computeA(getvalue(a), getvalue(b), getvalue(w))
end
line1x = [403.5, 641.5, 926.5, 1220.0, 1476.5]
line1y = [201.0, 169.5, 152.0, 156.5, 179.5]
line2x = [391.0, 633.0, 928.0, 1232.5, 1495.0]
line2y = [465.0, 456.0, 450.0, 449.5, 454.0]
line3x = [405.0, 643.5, 931.0, 1228.5, 1486.0]
line3y = [733.5, 749.5, 756.5, 750.5, 734.5]
# "Vertical" lines:
line4x = [515.0, 507.0, 504.0, 508.0, 516.5, 531.0]
line4y = [184.0, 319.0, 460.0, 603.0, 742.0, 873.0]
line5x =[779.5, 776.0, 776.0, 778.5, 783.0, 790.0]
line5y = [158.0, 301.0, 452.5, 607.0, 754.5, 893.0]
line6x = [1074.0, 1079.0, 1081.5, 1083.0, 1081.5, 1078.5]
line6y = [151.0, 295.0, 449.0, 605.0, 755.5, 895.5]
line7x = [1355.0, 1365.5, 1371.5, 1370.5, 1364.5, 1352.5]
line7y = [166.0, 305.0, 451.0, 600.0, 743.5, 879.5]
# Collect these together into arrays of arrays:
linesxs = [line1x, line2x, line3x, line4x, line5x, line6x, line7x]
linesys = [line1y, line2y, line3y, line4y, line5y, line6y, line7y]
params = LensCalibSolveParams()
findlenscalibration(linesxs, linesys, params)