Module Streamlines__Preprocess

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
 *
module ODM = Owl.Dense.Matrix.Generic
module ODN = Owl.Dense.Ndarray.Generic
module BA = Bigarray

Types

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

pv_verbosity functions

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

Vector field creation

val vector_length : float ‑> float ‑> float

vector_length a b

Calculate the length of the vector (a,b)

val compute_topo_gradient_field : (float, 'aODM.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, 'aODM.t ‑> (float, 'aODM.t ‑> (float, 'aODM.t * (float, 'aODM.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, 'aODM.t ‑> (float, 'aODM.t ‑> (float, 'aODM.t * (float, 'aODM.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

Blockages and loop finding and fixing

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.

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, 'aStreamlines__.Globals.ODM.t * (char, 'bStreamlines__.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, 'aStreamlines__.Globals.ODM.t * (char, 'bStreamlines__.Globals.ODM.t) ‑> unit

show_blockages (blockages_array, blocked_neighbors_array)

Show blockages from the big arrays

val fix_blockages : t ‑> ((char, 'aStreamlines__.Globals.ODM.t * (char, 'bStreamlines__.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.

Vector field creation

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)

Workflow functions

val create : Streamlines.Properties.t_props ‑> t

create props

Create a Preprocess.t from its properties

val process : t ‑> Streamlines.Core.t_core_data ‑> unit

process t data

Run the Preprocess workflow