push!(LOAD_PATH, "./../src/")

using Images, TestImages
using ApApproximation

img = testimage("lena_gray_256") |> x->convert(Array{Float64,2}, data(x))
convert(Image, img)

circle(i,j) = (i-size(img,1)/2)^2 + (j-size(img,1)/2)^2 <= (min(size(img)...)/2)^2*.95
global C = eltype(img)[circle(i,j) for i in 1:size(img,1), j in 1:size(img,2)]
img = C.*img + (1-C)

convert(Image, img)

Q, N = 20, 15

def = Q/2

v = [  (1, 1:1.:def/2),
    (3, def/2+1:1.:def)
    ];

Qv = sum([v[i][1]*length(v[i][2]) for i in 1:length(v)])

println("Does v agree with Q? $(Qv == Q)")

E = BispectralSet(N, v);

println("Size of the set E: (Q,N) = $(size(E))")

# plot(E)

# Since the image is fixed we conveniently define the interpolation functions 
# to take into account the size of `img`.
ApApproximation.bispectral2cart(f; clamped = true) = bispectral2cart(f, size(img)...; clamped = clamped)

# approximation of the set E
f1 = cartesian2bispectral(img, E); 

# interpolation on a cartesian grid
img2 = pol2cart(bispectral2pol(f1), 256, 256)

convert(Image, bispectral2cart(f1))

Q, N = 340, 64

max_ϱ = Q/4
step_ϱ = 1.

v = [ (4, 1:step_ϱ:max_ϱ) ];

Qv = sum([v[i][1]*length(v[i][2]) for i in 1:length(v)])

println("Does v agree with Q? $(Qv == Q)")

E = BispectralSet(N, v);

# approximation of the set E
f = cartesian2bispectral(img, E); 

convert(Image, bispectral2cart(f))

BesselMatrix(E) ## Dry run for compilation
@time J = BesselMatrix(E)

println("Dimension of J is (N,Q,Q): $(size(J.matrix)==(N,Q,Q))")

ap(f, J) ## Dry run for compilation
@time af = ap(f, J);

println("The size of af is (Q,N): $(size(af)==(Q,N))")

println()
println("Norm: $(norm(vec(af.f),2) |> x->round(x,2))")
println()

diff = norm(vec(af.f))/norm(vec(f.f)) |> x->round(x,2) #|> log10 |> x->round(x,1)
# println("L^2 norm of the resulting AP approximation $(round(Int,norm(vec(af.f),2)/norm(f.f,2)*100))% of the norm of the image")
println("The norm of the AP interpolation is $diff time bigger than that of the image")

# Function to visualize the magnitude on the bispectral set
function ApApproximation.bispectral2cart{N,T<:Complex}(af::BispInterpolation{N,T}) 
    normalize(x) = C.*x/maximum(C.*x) + (1-C)
    bispectral2cart(abs(af)) |> normalize
end

convert(Image,bispectral2cart(af))

@time f2 = iap(af, J);

println()
println("Norm: $(norm(vec(f2.f),2) |> x->round(x,2))")
println()

extr = extrema(bispectral2cart(f2, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

convert(Image, bispectral2cart(f2))

Base.exp(x::Frequency) =  exp(im*x.λ*cos(angle(x)))

function Base.exp(x::Frequency, y::Frequency)
    ApApproximation.composition(x,y) |> exp
end

ev(af::BispInterpolation) = ev(af, af.E)
function ev{N,T}(af::BispInterpolation, F::BispectralSet{N,T})
    f = Array{Complex128, 2}(size(F)...)
    for n in 1:size(f,2), j in 1:size(f,1)
        f[j,n] = sum([ af.f[k,m]*exp(af.E[k,m],F[j,n]) for k in 1:size(af.E,1), m in 1:size(af.E,2) ])
    end
    BispInterpolation{N,Float64}(real(f), F)
end

af1 = ap(f1)
convert(Image, bispectral2cart(af1))

ev(f1) # Dry run for compilation
iap(f1) # Dry run for compilation

tic() 
f1_ap = iap(f1)
t_ap = toq()

tic()
f1_ev = ev(f1)
t_ev = toq()

println("AP interpolation timing: $t_ap")
println("Direct computation timing: $t_ev")
println("Direct computation is $(round(Int,t_ev/t_ap)) times slower")
println("Relative L^2 error between the two: $(norm(vec(f1_ev.f - f1_ap.f) ,2)/norm(vec(f1_ap.f),2))")

# Rotation angle
rot = 2pi/N * 10

println("Rotation of $(round(2/N*10,2)) π radiants")
@time f2 = rotate(af, rot) |> x -> iap(x, J);

extr = extrema(bispectral2cart(f2, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

println()
println("Norm: $(norm(vec(f2.f),2) |> x->round(x,2))")
println()

convert(Image, bispectral2cart(f2))

# Translation (in polar coordinates)
ρ = 30.
θ = pi/3
println("Translation of $(map(x->round(x,1),(ρ*cos(θ), ρ*sin(θ))))")
@time f2 = translate(af, ρ, θ) |> x -> iap(x, J);

extr = extrema(bispectral2cart(f2, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

println()
println("Norm: $(norm(vec(f2.f),2) |> x->round(x,2))")
println()

convert(Image, bispectral2cart(f2))

P, M = Q, N

max_λ = 2. 
step_λ = max_λ / (Q/4)

v = [ (4, step_λ:step_λ:max_λ) ];

Pv = sum([v[i][1]*length(v[i][2]) for i in 1:length(v)])
println("Does v agree with P? $(Pv == P)")

F = BispectralSet(N, v);

quarter = round(Int,size(F,1)/4)
α = 100
d = [   ones(2quarter)*α/10; 
        ones(quarter)*α; 
        ones(size(F,1) - 3*quarter)*1e2α; ]

println("Generating Bessel matrices: ")
@time J_int = BesselMatrix(E,F, weights = d);
println("AP approximation: ")
@time af_int = ap(f, F, J_int)

println()
println("Norm: $(norm(vec(af_int.f),2) |> x->round(x,2))")
println()

diff = norm(vec(af_int.f))/norm(vec(f.f)) |> x->round(x,4) 
println("The norm of the AP interpolation is $diff time bigger than that of the image")

diff = norm(vec(af.f))/norm(vec(af_int.f)) |> log10 |> x->round(x,1)
println("The norm of the AP approximation is 10^($diff) time bigger than that of the AP interpolation")

convert(Image, bispectral2cart(af_int))

println("Inverse AP approximation: ")
@time f_int = iap(af_int, E, J_int)

println()
println("Norm: $(norm(vec(f_int.f),2) |> x->round(x,2))")
println()

extr = extrema(bispectral2cart(f_int, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

convert(Image, bispectral2cart(f_int))

# Rotation angle
rot = 2pi/N * 10

println("Rotation of $(round(2/N*10,2)) π radiants")
@time f_rot = rotate(af_int, rot) |> x -> iap(x, E, J_int);

extr = extrema(bispectral2cart(f_rot, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

println()
println("Norm: $(norm(vec(f_rot.f),2) |> x->round(x,2))")
println()

convert(Image, bispectral2cart(f_rot))

@time f_trans = translate(af_int, ρ, θ) |> x -> iap(x, E, J_int)

println()
println("Norm: $(norm(vec(f_trans.f),2) |> x->round(x,2))")
println()

extr = extrema(bispectral2cart(f_trans, clamped = false))
println("Maximal and minimum values of the resulting image: $extr")

convert(Image, bispectral2cart(f_trans))