strangeseg.jl

using Arrows
using NamedTuples
import JLD: load
# import Arrows.BenchmarkArrows: STD_ROTATION_MATRIX, render
import Images: colorview, Gray

using NamedTuples
const STD_ROTATION_MATRIX = [0.94071758  -0.33430171 -0.05738258
                             -0.33835238 -0.91297877 -0.2280076
                             0.02383425  0.2339063 -0.97196698]


"(width, height, 2) array, res[i,j] = [i,j] - ray dir based on an increment"
function gen_fragcoords(width::Integer, height::Integer)
  raster_space = zeros(width, height, 2)
  for i = 1:width
    for j = 1:height
      raster_space[i, j, :] = [i, j] - 0.5
    end
  end
  return raster_space
end

"Render rays starting with raster_space according to geometry"
function make_ro(r, raster_space, width, height)
  nmatrices = 1
  resolution = [width, height]
  # Normalise it to be bound between 0 1
  norm_raster_space = raster_space ./ reshape(resolution, 1, 1, 2)
  # Put it in NDC space, -1, 1
  screen_space = -1.0 + 2.0 * norm_raster_space
  # Make pixels square by mul by aspect ratio
  aspect_ratio = resolution[1] / resolution[2]
  ndc_space = screen_space .* reshape([aspect_ratio, 1.0], (1, 1, 2))
  # Ray Direction

  # Position on z-plane
  scalars = ones(width, height, 1) * 1.0
  ndc_xyz = cat(3, ndc_space, scalars) * 0.5  # Change focal length

  # Put the origin farther along z-axis
  ro = [0, 0, 1.5]

  # Rotate both by same rotation matrix
  ro_t = reshape(ro, (1, 3)) * r
  ndc_t = Array{Float64}(width, height, nmatrices, 3)
  for w = 1:width, h = 1:height
    ndc_t[w, h, 1, :] = ndc_xyz[w, h, :]' * r
  end
  ndc_t = reshape(ndc_t, (width, height, nmatrices, 3))
  ndc_t = permutedims(ndc_t, (3, 1, 2, 4))

  # Increment by 0.5 since voxels are in [0, 1]
  ro_t = ro_t + 0.5
  ndc_t = ndc_t + 0.5
  # Find normalise ray dirs from origin to image plane
  unnorm_rd = ndc_t .- reshape(ro_t, (nmatrices, 1, 1, 3))
  norms = [norm(unnorm_rd[:,w,h,:]) for w = 1:size(unnorm_rd, 2), h = 1:size(unnorm_rd, 3)]

  rd = unnorm_rd ./ reshape(norms, (nmatrices, width, height, 1))
  return rd, ro_t
end

opt = @NT(width = 32, height = 32, nsteps = 3, res = 32, batch_size = 1,
          phong = false, density = 2)
width, height, res = opt.width, opt.height, opt.res
@show raster_space = gen_fragcoords(width, height)
@show rd, ro = make_ro(STD_ROTATION_MATRIX, raster_space, width, height)