Commit 1cf1e763 authored by Ultima's avatar Ultima
Browse files

implemented nifty_los.py and nifty_los_implementation.pyx

parent bf83cec9
# -*- coding: utf-8 -*-
import pyximport; pyximport.install()
from nifty_los_implementation import los_integrator
# -*- coding: utf-8 -*-
import numpy as np
cimport numpy as np
cimport cython
class los_integrator(object):
def __init__(self, shape, start, end):
self.dtype = np.dtype('float')
self.shape = tuple(shape)
self.start = np.array(start, dtype=self.dtype)
self.end = np.array(end, dtype=self.dtype)
assert(np.all(np.array(self.shape) != 0))
assert(len(self.shape) == len(self.start) == len(self.end))
assert(len(self.start.shape) == len(self.end.shape) == 1)
def integrate(self):
if np.all(self.start == self.end):
return self._empty_results()
try:
projected_start = self._project_to_cuboid(mode='start')
projected_end = self._project_to_cuboid(mode='end')
except ValueError:
return self._empty_results()
(indices, weights) = self._integrate_through_cuboid(projected_start,
projected_end)
return (indices, weights)
def _empty_results(self):
return ([np.array([], dtype=np.dtype('int'))] * len(self.shape),
np.array([], dtype=np.dtype('float')))
def _project_to_cuboid(self, mode):
if mode == 'start':
a = self.start
b = self.end
elif mode == 'end':
a = self.end
b = self.start
else:
raise ValueError
if np.all(np.zeros_like(a) <= a) and np.all(a <= np.array(self.shape)):
return a
c = b - a
surface_list = []
for i in xrange(len(self.shape)):
surface_list += [[i, 0]]
surface_list += [[i, self.shape[i]]]
translator_list = map(lambda z:
self._get_translator_to_surface(a, c, *z),
surface_list)
# sort the translators according to their norm, save the sorted indices
translator_index_list = np.argsort(map(np.linalg.norm,
translator_list))
# iterate through the indices -from short to long translators- and
# take the first translator which brings a to the actual surface of
# the cuboid and not just to one of the parallel planes
found = False
for i in translator_index_list:
p = a + translator_list[i]
if np.all(np.zeros_like(p) <= p) and \
np.all(p <= np.array(self.shape)):
found = True
break
if not found:
raise ValueError(
"ERROR: Line-of-sight does not go through cuboid.")
return p
def _get_translator_to_surface(self, point, full_direction,
dimension_index, surface):
"""
translates 'point' along the vector 'direction' such that the
dimension with index 'dimension_index' has the value 'surface'
"""
direction_scaler = np.divide((surface - point[dimension_index]),
full_direction[dimension_index])
if direction_scaler < 0 or direction_scaler > 1:
return point * np.nan
scaled_direction = full_direction * direction_scaler
return scaled_direction
def _integrate_through_cuboid(self, start, end):
# estimate the maximum number of cells that could be hit
# the current estimator is: norm of the vector times number of dims
num_estimate = np.ceil(np.linalg.norm(end - start))*len(start)
index_list = np.empty((num_estimate, len(start)),
dtype=np.dtype('int'))
weight_list = np.empty((num_estimate), self.dtype)
current_position = start
i = 0
while True:
next_position, weight = self._get_next_position(current_position,
end)
floor_current_position = np.floor(current_position)
index_list[i] = floor_current_position
weight_list[i] = weight
if np.all(np.floor(current_position) == np.floor(end)):
break
current_position = next_position
i += 1
return list(index_list[:i].T), weight_list[:i]
def _get_next_position(self, position, end_position):
full_direction = end_position - position
surface_list = []
for i in xrange(len(position)):
surface_list += [[i, strong_floor(position[i])]]
surface_list += [[i, strong_ceil(position[i])]]
translator_list = map(lambda z: self._get_translator_to_surface(
position,
full_direction,
*z),
surface_list)
# index_of_best_translator = np.argsort(np.linalg.norm(translator_list,
# axis=1))[0]
translator_list = np.array(translator_list)
index_of_best_translator = np.linalg.norm(translator_list, axis=1)
index_of_best_translator = np.argsort(index_of_best_translator)[0]
best_translator = translator_list[index_of_best_translator]
# if the surounding surfaces are not reachable, it must be the case
# that the current position is in the same cell as the endpoint
if np.isnan(np.linalg.norm(best_translator)):
floor_position = np.floor(position)
# check if position is in the same cell as the endpoint
assert(np.all(floor_position == np.floor(end_position)))
weight = np.linalg.norm(end_position - position)
next_position = None
else:
next_position = position + best_translator
weight = np.linalg.norm(best_translator)
return (next_position, weight)
def strong_floor(x):
floor_x = np.floor(x)
if floor_x == x:
return floor_x-1
else:
return floor_x
def strong_ceil(x):
ceil_x = np.ceil(x)
if ceil_x == x:
return ceil_x+1
else:
return ceil_x
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment