Copyright (C) 2017-2018, Colin P Stark and Gavin J Stark. All rights reserved. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * * @file preprocess.ml * @brief Preprocessing *
type
t
=
{
props : Streamlines.Properties.t_props_preprocess; |
mutable
roi_gradx_array : Streamlines.Globals.t_ba_floats; |
mutable
roi_grady_array : Streamlines.Globals.t_ba_floats; |
mutable
where_looped : (int * int)
list; |
pad_width : int; |
}
t
Structure for the Preprocess workflow
val
pv_noisy : t ‑> (unit ‑>
unit) ‑> unit
pv_noisy t
Shortcut to use Preprocess.t verbosity for Properties.pv_noisy
val
pv_debug : t ‑> (unit ‑>
unit) ‑> unit
pv_debug t
Shortcut to use Preprocess.t verbosity for Properties.pv_debug
val pv_info
: t ‑> (unit ‑>
unit) ‑> unit
pv_info t
Shortcut to use Preprocess.t verbosity for Properties.pv_info
val
pv_verbose : t ‑> (unit ‑>
unit) ‑> unit
pv_verbose t
Shortcut to use Preprocess.t verbosity for Properties.pv_verbose
val
vector_length : float ‑> float
‑> float
vector_length a b
Calculate the length of the vector (a,b)
val
compute_topo_gradient_field : (float, 'a) ODM.t ‑> (float, BA.float32_elt) ODM.t *
(float, BA.float32_elt) ODM.t
compute_topo_gradient_field
roi_array
Differentiate ROI topo in x and y directions to create a gradient vector field
val
normalize_arrays : (float, 'a) ODM.t ‑> (float, 'a) ODM.t ‑> (float, 'a) ODM.t *
(float, 'a) ODM.t
normalize_arrays u_array v_array
Normalize a pair of arrays, by scaling them by 1/mod(u,v)
val
compute_gradient_velocity_field : (float, 'a) ODM.t ‑> (float, 'a) ODM.t ‑> (float, 'a) ODM.t *
(float, 'a) ODM.t
compute_gradient_velocity_field
roi_grady_array roi_grady_array
Compute normalized gradient velocity vector field from ROI topo grid.
This is effectively an array of unit (dh/dx, dh/dy) vectors, but held as two scalar arrays dh/dx = u and dh/dy = v
val
get_flow_vector : char ‑> int
* int
get_flow_vector n
Convert a char direction indicator to an int vector pair
val
get_unit_flow_vector : char ‑>
float * float
get_unit_flow_vector nn
Get a float vector pair of length 1 from a char direction indicator
val
has_just_diagonal_outflows : (int ‑> int ‑>
'a) ‑> int ‑> int
‑> char
has_just_diagonal_outflows get x y
A map for a 2D bigarray, hence provided with a get x y -> value function.
Returns mask of w,x,y,z otherwise for SW, SE, NW, NE outflows but 0 if there are any N, S, E or W outflows.
Consider a gradient field of the topology; it is composed using dh/dx and dh/dy. Hence diagonal outflows are missed at the vector flow, which is generated from that gradient field
Hence if a point only has diagonal outflows it is a 'blockage'.
Consider in particular:
roi -> roi_gradx_array -> roi_grady_array -> normalized velocity
6 6 6 6 6
6 7 9 9 6 +1.5 +1.0 -1.5 -1.5 -1.0 -1.5 -,+ -,+ +,+
6 9 8 9 6 -> +1.0 0 -1.0 ->-1.0 0 0 -> -,+ 0,0 +,0
6 9 9 9 6 +1.5 0 -1.5 +1.5 +1.0 +1.5 -,- 0,- +,-
6 6 6 6 6
This is mapped to blockages and blocked_neighbors
6 6 6 6 6
6 7 9 9 6 0 0 0 0 0 0 NW NW NE
6 9 8 9 6 -> 0 64 0 -> 0 0 1 -> NW NW W
6 9 9 9 6 0 0 0 0 8 64 SW N NW
6 6 6 6 6
The DTM will have been preconditioned by a GIS to drain to the edge using eight different directions
val
upstream_of_diagonal_outflow : (int ‑> int ‑> char)
‑> int ‑> int ‑>
char
upstream_of_diagonal_outflow get x
y
A map for a 2D bigarray, hence provided with a get x y -> value function.
The upstream pixels of a blockage need to be found and fixed; why is not quite clear.
This helps to remove loops introduced by tweaking the vector field at the blockages
val
calc_speed_div_curl : int ‑> int
‑> (float, 'a, 'b) Bigarray.Genarray.t ‑> (float, 'c, 'd) Bigarray.Genarray.t ‑> float * float
* float
calc_speed_div_curl x y u v
Get UV of x,y and its +1 neighbors, the speed=mod(avg uv), divergence and curl
Geograpically:
uv01 uv11
uv00 uv10
speed = sqrt((+u00+u01+u10+u11)^2 + (+v00+v01+v10+v11)^2 )/4
div = -u00-u01+u10+u11 -v00+v01-v10+v11a
curl = +u00-u01+u10-u11 -v00-v01+v10+v11
Another way to look at uv is uvNM = (dh/dx (x+N,y+M), dh/dy (x+N,y+M))
Hence:
uv00 . (-1,-1) = -dh/dx(x,y) - dh/dy(x,y)
etc
Hence div is
uv01.NW + uv11.NE
+ uv00.SW + uv10.SE
And curl is
uv01.SW + uv11.NW
+ uv00.SE + uv10.NE
As the vector field is a gradient vector field it ought to have a curl of zero (it would in a perfect world) and divergence should be greater than zero (as there is always a downhill);
val
find_blockages : t ‑> Streamlines.Core.t_core_data
‑>
(char, Bigarray.int8_unsigned_elt) ODN.t *
(char, Bigarray.int8_unsigned_elt) ODN.t
find_blockages t data
The DTM is preconditioned for draining to the edge in 8 directions; the vector gradient field generated by preprocess is 4 direction though. Hence there may be 'blockages' and 'loops.
Find blockages (pixels where there are only diagonal outflows, so the vector gradients of neighbors all lead in) and blocked neighbors
val
get_blockages_lists : ((char, 'a) Streamlines__.Globals.ODM.t
* (char, 'b) Streamlines__.Globals.ODM.t)
‑> (int * int) list * (int
* int) list
get_blockages_lists (blockages_array,
blocked_neighbors_array)
Convert blockage big arrays to lists of (x,y) locations where there are blockages or blocked neighbors
val
show_blockages : ((char, 'a) Streamlines__.Globals.ODM.t
* (char, 'b) Streamlines__.Globals.ODM.t)
‑> unit
show_blockages (blockages_array,
blocked_neighbors_array)
Show blockages from the big arrays
val
fix_blockages : t ‑> ((char, 'a) Streamlines__.Globals.ODM.t
* (char, 'b) Streamlines__.Globals.ODM.t)
‑> Streamlines.Core.t_core_data
‑> unit
fix_blockages t (blockages_array,
blocked_neighbors_array) data
Fix blockages and blocked neighbors by making them use their diagonal neighbor as specified by the big arrays
val
break_out_of_loop : Streamlines.Core.t_core_data
‑> (int * int) ‑>
unit
break_out_of_loop data (x,y)
Break out of a loop at a given (x,y), by forcing the pixel to flow to the lowest of the 8 neighbors
val
find_and_fix_loops : t ‑> Streamlines.Core.t_core_data
‑> (int * int) list
find_and_fix_loops t data
Find loops (sets of 4 elements where curl is much above zero, divergence is low - i.e. are likely to be very flat and poorly conditioned) by making them flow towards their lowest neighbors.
val
conditioned_gradient_vector_field : t ‑>
Streamlines.Core.t_core_data
‑> unit
conditioned_gradient_vector_field t
data
Compute topographic gradient vector field on a preconditioned DTM.
The preconditioning steps are:
1. Find blockages in gradient vector field 2. Calculate surface derivatives (gradient & 'curvature') 3. Set gradient vector field magnitudes ('speeds') 4. Find and fix loops in gradient vector field 5. Fix blockages in gradient vector field 6. Set initial streamline points ('seeds')
val
mask_nan_uv : Streamlines.Core.t_core_data
‑> unit
mask_nan_uv data
Mask out (in the mask big arrays) any elements that have a U or V value of NaN
raw_gradient_vector_field ?? data
Not implemented yet
Compute topographic gradient vector field without preconditioning the DTM. let raw_gradient_vector_field data = self.print('Compute raw gradient vector field') (self.roi_gradx_array,self.roi_grady_array) = derivatives(self.geodata.roi_array) self.speed_array = set_speed_field(self.u_array, self.v_array)
val create
: Streamlines.Properties.t_props
‑> t
create props
Create a Preprocess.t from its properties