lidar_spin.py

#!/usr/bin/env python

import rasterio as rio
from sys import argv, exit
from skimage import io
import numpy as np
from math import pi, radians

def px_to_rad(n):
  return n / 200.0 # 1 px = 5 mrad

smallest_inc = 0
largest_inc = 0

def rincaz_to_xyz(r, inc, az):
  r += 8.241897
  r *= 20 # arbitrary
  
  # going from the range
  # zenith..low in pixels to
  # low..zenith, horizon = 0 in rad
  inc = px_to_rad(454 - inc) # 454 because that's the image height
  inc += 0.8726646259971648 # 0.87... is the image height in rad
  
  az = px_to_rad(az)
  
  # standard spherical->xyz conversion:
  x = r * np.sin(inc) * np.cos(az)
  y = r * np.sin(inc) * np.sin(az)
  z = r * np.cos(inc)
  
  x, y, z = map(int, [x, y, z])
  return (x, y, z)


with rio.open(argv[1]) as src:
  razinc = src.read()

volume = np.zeros((1080, 1080, 1080), dtype=np.uint8)

# most of this is terrible ad-hoc nonsense:

for step in range(360):
  volume[...] = 0
  for az in range(razinc.shape[1]):
    for inc in range(razinc.shape[2]):
      r = razinc[:, az-step, inc][0]
      if r == 0:
        continue

      x, y, z = rincaz_to_xyz(r, inc, az)

      try:
        volume[x+540, y+540, z+540] += 1
      except IndexError:
        pass

  flatter = np.sum(volume.astype(np.float64), axis=0)
  flatter = ((flatter / 16) * (2**16-1)).astype(np.uint16)
  io.imsave('degrees/%s.tif' % (step), flatter)