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)
```