[GRASS-SVN] r60876 - sandbox/annakrat/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 19 21:46:57 PDT 2014
Author: annakrat
Date: 2014-06-19 21:46:57 -0700 (Thu, 19 Jun 2014)
New Revision: 60876
Modified:
sandbox/annakrat/r3.flow/r3.flow.py
sandbox/annakrat/r3.flow/rast3d_functions.py
sandbox/annakrat/r3.flow/voxel_traversal.py
Log:
r3.flow: fixed voxel traversal algorithm and added test for it
Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py 2014-06-20 04:30:33 UTC (rev 60875)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-06-20 04:46:57 UTC (rev 60876)
@@ -99,7 +99,8 @@
from grass.pygrass.vector.geometry import Line, Point
from integrate import integrate_next, get_time_step
-from rast3d_functions import get_velocity, norm, compute_gradient, rast3d_location2coord
+from rast3d_functions import get_velocity, norm, compute_gradient, rast3d_location2coord, rast3d_is_valid_location
+from voxel_traversal import traverse
EPSILON = 1e-8
@@ -231,16 +232,23 @@
break # zero velocity means end of propagation
# convert to time
delta_T = get_time_step(unit, step, velocity_norm, cell_size)
- point = integrate_next(map_info, velocity, point, delta_T * direction)
- if point is None:
+ new_point = integrate_next(map_info, velocity, point, delta_T * direction)
+ if new_point is None or \
+ not rast3d_is_valid_location(map_info, new_point[1], new_point[0], new_point[2]):
break
if seed.flowline:
- line.append(Point(point[0], point[1], point[2]))
+ line.append(Point(new_point[0], new_point[1], new_point[2]))
if seed.flowaccum:
- col, row, depth = rast3d_location2coord(map_info, point[1], point[0], point[2])
+ col, row, depth = rast3d_location2coord(map_info, new_point[1], new_point[0], new_point[2])
if last_coords != (col, row, depth):
flowacc[depth][row][col] += 1
+ if last_coords[0] is not None and \
+ sum([abs(coord[0] - coord[1]) for coord in zip(last_coords, (col, row, depth))]) > 1:
+ additional = traverse(point, new_point, map_info)
+ for p in additional:
+ flowacc[p[2]][p[1]][p[0]] += 1
last_coords = (col, row, depth)
+ point = new_point
count += 1
if seed.flowline and len(line) > 1:
Modified: sandbox/annakrat/r3.flow/rast3d_functions.py
===================================================================
--- sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-20 04:30:33 UTC (rev 60875)
+++ sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-20 04:46:57 UTC (rev 60876)
@@ -31,6 +31,12 @@
# except for trilinear interpolation which could be added there.
+def rast3d_is_valid_location(rast3d_region, north, east, top):
+ return (north >= rast3d_region['south'] and north <= rast3d_region['north']) and \
+ (east >= rast3d_region['west'] and east <= rast3d_region['east']) and \
+ (top >= rast3d_region['bottom'] and top <= rast3d_region['top'])
+
+
def rast3d_location2coord(rast3d_region, north, east, top):
"""!Returns column, row and depth of a cell where the input
coordinates (north, east, top) lie."""
Modified: sandbox/annakrat/r3.flow/voxel_traversal.py
===================================================================
--- sandbox/annakrat/r3.flow/voxel_traversal.py 2014-06-20 04:30:33 UTC (rev 60875)
+++ sandbox/annakrat/r3.flow/voxel_traversal.py 2014-06-20 04:46:57 UTC (rev 60876)
@@ -6,8 +6,12 @@
import numpy as np
import math
+import unittest
+from random import uniform
+# based on Amanatides, John, and Andrew Woo. "A fast voxel traversal algorithm for ray tracing."
+# Proceedings of EUROGRAPHICS. Vol. 87. 1987.
def traverse(start, end, region_info):
"""!Computes which voxel cells are traversed by a line.
To each transected cell we add 1.
@@ -15,13 +19,17 @@
:param start: start point (3 coordinates)
:param end: end point (3 coordinates)
:param dict region_info: information about region
- :return: coordinates of transected voxel cells
+ :return: coordinates of transected voxel cells from north west bottom
"""
# initialize
dx = end[0] - start[0]
dy = end[1] - start[1]
dz = end[2] - start[2]
+ stepX = 1 if start[0] < end[0] else -1
+ stepY = 1 if start[1] < end[1] else -1
+ stepZ = 1 if start[2] < end[2] else -1
+
X = int(math.floor((start[0] - region_info['west']) / region_info['ewres']))
Y = int(math.floor((start[1] - region_info['south']) / region_info['nsres']))
Z = int(math.floor((start[2] - region_info['bottom']) / region_info['tbres']))
@@ -29,78 +37,90 @@
Y_end = int(math.floor((end[1] - region_info['south']) / region_info['nsres']))
Z_end = int(math.floor((end[2] - region_info['bottom']) / region_info['tbres']))
- tDeltaX = region_info['ewres'] / float(dx) if dx != 0 else 1e6
- tDeltaY = region_info['nsres'] / float(dy) if dy != 0 else 1e6
- tDeltaZ = region_info['tbres'] / float(dz) if dz != 0 else 1e6
+ tDeltaX = abs(region_info['ewres'] / float(dx) if dx != 0 else 1e6)
+ tDeltaY = abs(region_info['nsres'] / float(dy) if dy != 0 else 1e6)
+ tDeltaZ = abs(region_info['tbres'] / float(dz) if dz != 0 else 1e6)
- xtmp = start[0] / region_info['ewres']
- ytmp = start[1] / region_info['nsres']
- ztmp = start[2] / region_info['tbres']
- xtmp = xtmp if xtmp > 0 else -xtmp
- ytmp = ytmp if ytmp > 0 else -ytmp
- ztmp = ztmp if ztmp > 0 else -ztmp
- tMaxX = tDeltaX * (1.0 - (xtmp - int(xtmp)))
- tMaxY = tDeltaY * (1.0 - (ytmp - int(ytmp)))
- tMaxZ = tDeltaZ * (1.0 - (ztmp - int(ztmp)))
+ xtmp = (start[0] - region_info['west']) / region_info['ewres']
+ ytmp = (start[1] - region_info['south']) / region_info['nsres']
+ ztmp = (start[2] - region_info['bottom']) / region_info['tbres']
- stepX = 1 if start[0] < end[0] else -1
- stepY = 1 if start[1] < end[1] else -1
- stepZ = 1 if start[2] < end[2] else -1
+ if stepX > 0:
+ tMaxX = tDeltaX * (1.0 - (xtmp - math.floor(xtmp)))
+ else:
+ tMaxX = tDeltaX * (xtmp - math.floor(xtmp))
+ if stepY > 0:
+ tMaxY = tDeltaY * (1.0 - (ytmp - math.floor(ytmp)))
+ else:
+ tMaxY = tDeltaY * (ytmp - math.floor(ytmp))
+ if stepZ > 0:
+ tMaxZ = tDeltaZ * (1.0 - (ztmp - math.floor(ztmp)))
+ else:
+ tMaxZ = tDeltaZ * (ztmp - math.floor(ztmp))
coordinates = []
while True:
- if abs(tMaxX) < abs(tMaxY):
- if abs(tMaxX) < abs(tMaxZ):
+ if tMaxX < tMaxY:
+ if tMaxX < tMaxZ:
tMaxX = tMaxX + tDeltaX
X = X + stepX
- if (X - X_end) * stepX > 0:
- break
+
else:
tMaxZ = tMaxZ + tDeltaZ
Z = Z + stepZ
- if (Z - Z_end) * stepZ > 0:
- break
+
else:
- if abs(tMaxY) < abs(tMaxZ):
+ if tMaxY < tMaxZ:
tMaxY = tMaxY + tDeltaY
Y = Y + stepY
- if (Y - Y_end) * stepY > 0:
- break
+
else:
tMaxZ = tMaxZ + tDeltaZ
Z = Z + stepZ
- if (Z - Z_end) * stepZ > 0:
- break
- coordinates.append((X, Y, Z))
+ if (X, Y, Z) == (X_end, Y_end, Z_end):
+ break
+ coordinates.append((X, region_info['rows'] - Y - 1, Z))
return coordinates
+
+class TestTraverseVoxel(unittest.TestCase):
+
+ def test_traverse(self):
+ for i in range(100):
+ region_info = {'rows': int(uniform(1, 50)), 'cols': int(uniform(1, 50)), 'depths': int(uniform(1, 50)),
+ 'nsres': uniform(0.5, 10), 'ewres': uniform(0.5, 10), 'tbres': uniform(0.5, 10),
+ 'south': uniform(-100, 100), 'west': uniform(-100, 100), 'bottom': uniform(-100, 100)}
+
+ start = [(uniform(region_info['west'], region_info['west'] + region_info['cols'] * region_info['ewres'])),
+ (uniform(region_info['south'], region_info['south'] + region_info['rows'] * region_info['nsres'])),
+ (uniform(region_info['bottom'], region_info['bottom'] + region_info['depths'] * region_info['tbres']))]
+ end = [(uniform(region_info['west'], region_info['west'] + region_info['cols'] * region_info['ewres'])),
+ (uniform(region_info['south'], region_info['south'] + region_info['rows'] * region_info['nsres'])),
+ (uniform(region_info['bottom'], region_info['bottom'] + region_info['depths'] * region_info['tbres']))]
+
+ coords = traverse(start, end, region_info)
+
+ differ = np.array(end) - np.array(start)
+ t = np.linspace(0, 1, 10000)
+ x = start[0] + differ[0] * t
+ y = start[1] + differ[1] * t
+ z = start[2] + differ[2] * t
+
+ X = np.floor((x - region_info['west']) / region_info['ewres']).astype(int)
+ Y = np.floor((y - region_info['south']) / region_info['nsres']).astype(int)
+ Y = region_info['rows'] - Y - 1
+ Z = np.floor((z - region_info['bottom']) / region_info['tbres']).astype(int)
+ coords2 = [tuple(row) for row in np.column_stack((X, Y, Z))]
+
+ seen = set()
+ # we don't use unique because we want to keep the order
+ unique = [k for k in coords2 if k not in seen and not seen.add(k)][1:-1]
+# message = "start: {start}\nend: {end}\nregion info: {info}".format(start=start,
+# end=end, info=region_info)
+ # we accept 'First list contains 1 additional elements.' because comparing is not as precise
+ self.assertListEqual(coords, unique)
+
+
if __name__ == '__main__':
- info = {}
- info['cols'] = 10
- info['rows'] = 10
- info['depths'] = 10
- info['nsres'] = 4
- info['ewres'] = 8
- info['tbres'] = 1
- info['west'] = -80
- info['south'] = -40
- info['bottom'] = -10
- start = [-66.2,-38.9,-9.5]
- end = [-25,-24.5,-9.5]
- X = int(math.floor((start[0] - info['west']) / info['ewres']))
- Y = int(math.floor((start[1] - info['south']) / info['nsres']))
- Z = int(math.floor((start[2] - info['bottom']) / info['tbres']))
- X_end = int(math.floor((end[0] - info['west']) / info['ewres']))
- Y_end = int(math.floor((end[1] - info['south']) / info['nsres']))
- Z_end = int(math.floor((end[2] - info['bottom']) / info['tbres']))
-# start = [1, 1, 0.5]
-# end = [74, 42, 0.5]
-
- array = np.zeros((info['depths'], info['rows'], info['cols']), dtype=np.int8)
- coords = traverse(start, end, info)
- array[Z, array.shape[1] - Y - 1, X] = 1
- array[Z_end, array.shape[1] - Y_end - 1, X_end] = 3
- for coord in coords:
- array[coord[2], array.shape[1] - coord[1] - 1, coord[0]] = 2
- q = array[0]
+ unittest.main()
More information about the grass-commit
mailing list