rd_target_compare.jl

function interp(x, ys, xs)
  ys[1] + (ys[2] - ys[1]) * ((x - xs[1]) / (xs[2] - xs[1]))
end

base_data = readdlm(ARGS[1])
comp_data = readdlm(ARGS[2])
rows = max(size(base_data, 1), size(comp_data, 1))

# Calculate bits per pixel
bpp = zeros(rows, 2)
bpp[1:size(base_data,1),1] = base_data[:,3] * 8 ./ base_data[:,2]
bpp[1:size(comp_data,1),2] = comp_data[:,3] * 8 ./ comp_data[:,2]

# Find the first indices that are less than 0.1 and 0.05 in the base data
idx = zeros(2)
idx[1] = find(bpp[:,1] .< 0.1)[1]
idx[2] = find(bpp[:,1] .< 0.05)[1]

# Interpolate to find the dB values for PSNR-HVS and FastSSIM
i = idx[1]-1:idx[1]
j = idx[2]-1:idx[2]
psnrhvs = [base_data[i,5]'; base_data[j,5]']
fastssim = [base_data[i,7]'; base_data[j,7]']
target_db = [interp(0.1, psnrhvs[1,:], bpp[i,1]) interp(0.1, fastssim[1,:], bpp[i,1]);
             interp(0.05, psnrhvs[2,:], bpp[j,1]) interp(0.05, fastssim[2,:], bpp[j,1])]

# Find first indices that are less than the target dBs in the comparison data
idx = zeros(2, 2)
idx[1,1] = find(comp_data[:,5] .< target_db[1,1])[1]
idx[1,2] = find(comp_data[:,7] .< target_db[1,2])[1]
idx[2,1] = find(comp_data[:,5] .< target_db[2,1])[1]
idx[2,2] = find(comp_data[:,7] .< target_db[2,2])[1]

# Interpolate to find the bpp used for target dB in comparison data
r = Array(FloatRange, (2,2))
for i in 1:2
  for j in 1:2
    r[i,j] = idx[i,j]-1:idx[i,j]
  end
end

println(comp_data[:,5:7])
println(bpp)
println(target_db)
println(r)

comp_bpp = zeros(2, 2)
for i in 1:2
  for j in 1:2
    metric = 5
    if j == 2
      metric = 7
    end
    comp_bpp[i,j] = interp(target_db[i,j], bpp[r[i,j],2], comp_data[r[i,j],metric])
  end
end

base_bpp = [0.1 0.1; 0.05 0.05]

delta_bpp = (comp_bpp ./ base_bpp - 1) * 100

println(base_bpp)
println(comp_bpp)
println(delta_bpp)

#println("PSNR-HVS: $delta_bpp[1,1] @ 0.1bpp    $delta_bpp[2,1] @ 0.05bpp")