Merge pull request #112 from GO-Eratosthenes/test-shadow-detection

implemented test for shadow classification

dhdt/generic/handler_dat.py

@@ -0,0 +1,21 @@
+import os
+
+def get_list_files(dat_path, ending):
+    """ generate a list with files, with a given ending/extension
+
+    Parameters
+    ----------
+    dat_path : string
+        path of interest.
+    ending : string
+        ending of a file.
+
+    Returns
+    -------
+    list
+        collection of file names.
+    """
+    if not os.path.isdir(dat_path):
+        dat_path = os.path.dirname(dat_path)
+    file_list = [x for x in os.listdir(dat_path) if x.endswith(ending)]
+    return file_list

dhdt/preprocessing/image_transforms.py

@@ -6,7 +6,7 @@
 
 
 def mat_to_gray(Z, notZ=None, vmin=None, vmax=None):
-    """ transform matix array  to float, omitting nodata values
+    """ transform matix array to float, omitting nodata values
 
     Parameters
     ----------

dhdt/preprocessing/shadow_geometry.py

@@ -221,7 +221,7 @@ def get_shadow_polygon(M, t_siz):
 
 
 # threshholding and classification functions
-def sturge(M):
+def sturge(M, neighbors=2):
     """ Transform intensity array to labelled image
 
     Parameters
@@ -235,7 +235,7 @@ def sturge(M):
         array with numbered labels
     """
     values, base = normalized_sampling_histogram(M)
-    dips = find_valley(values, base, 2)
+    dips = find_valley(values, base, neighbors=neighbors)[0]
     val = max(dips)
     imSeparation = M > val  # or <
 
@@ -280,11 +280,11 @@ def find_valley(values, base, neighbors=2):
     wall_plu = np.vstack((values, wall_plu))
     wall_min = np.vstack((values, wall_min))
 
+    concav_plu = np.sign(np.diff(wall_plu, n=1, axis=0)) == +1
+    concav_min = np.sign(np.diff(wall_min, n=1, axis=0)) == +1
     if neighbors > 1:
-        concav_plu = np.all(np.sign(np.diff(wall_plu, n=1, axis=0)) == +1,
-                            axis=0)
-        concav_min = np.all(np.sign(np.diff(wall_min, n=1, axis=0)) == +1,
-                            axis=0)
+        concav_plu = np.all(concav_plu, axis=0)
+        concav_min = np.all(concav_min, axis=0)
 
     # select the dips
     selec = np.all(np.vstack((concav_plu, concav_min)), axis=0)

dhdt/preprocessing/shadow_transforms.py

@@ -1,14 +1,38 @@
 import numpy as np
 
-from dhdt.generic.data_tools import s_curve
+from inspect import signature
 
+from ..generic.data_tools import s_curve
 from ..generic.unit_check import are_three_arrays_equal, are_two_arrays_equal
 from .color_transforms import (erdas2hsi, lab2lch, lms2lab, rgb2hcv, rgb2hsi,
                                rgb2lms, rgb2xyz, rgb2ycbcr, rgb2yiq, xyz2lab)
 from .multispec_transforms import (pca_rgb_preparation,
                                    principle_component_analysis)
 
 
+def _get_shadow_transforms():
+    methods_dict = {
+        'siz': 'shadow index first version',
+        'isi': 'improved shadow index',
+        'sei': 'shadow enhancement index',
+        'fcsdi': 'false color shadow difference index',
+        'csi': 'combinational shadow index',
+        'nsvdi': 'normalized saturation value difference index',
+        'sil': 'shadow index second',
+        'sr': 'specthem ratio',
+        'sdi': 'shadow detector index',
+        'sisdi': 'saturation intensity shadow detector index',
+        'mixed': 'mixed property based shadow index',
+        'msf': 'modified shadow fraction',
+        'c3': 'color invariant',
+        'entropy': 'entropy shade removal',
+        'shi': 'shade index',
+        'sp': 'shadow probabilities',
+        'nri': 'normalized range shadow index'
+    }
+    return methods_dict
+
+
 def apply_shadow_transform(method, Blue, Green, Red, RedEdge, Near, Shw,
                            **kwargs):
     """ Given a specific method, employ shadow transform
@@ -62,6 +86,9 @@ def apply_shadow_transform(method, Blue, Green, Red, RedEdge, Near, Shw,
         reads the specific bands given
     """
     method = method.lower()
+    methods = _get_shadow_transforms()
+    assert method in list(methods.keys()), \
+        'please provide correct method string'
 
     if method == 'siz':
         M = shadow_index_zhou(Blue, Green, Red)
@@ -86,20 +113,22 @@ def apply_shadow_transform(method, Blue, Green, Red, RedEdge, Near, Shw,
     elif method == 'mixed':
         M = mixed_property_based_shadow_index(Blue, Green, Red)
     elif method == 'msf':
-        p_s = kwargs['P_S'] if 'P_S' in kwargs else .95
+        p_s = kwargs['P_S']
+        if p_s is None:
+            p_s = signature(modified_shadow_fraction).parameters['p_s'].default
         M = modified_shadow_fraction(Blue, Green, Red, P_S=p_s)
     elif method == 'c3':
         M = color_invariant(Blue, Green, Red)
     elif method == 'entropy':
-        a = kwargs['a'] if 'a' in kwargs else 138.
+        a = kwargs['a']
+        if a is None:
+            a = signature(entropy_shade_removal).parameters['a'].default
         M, R = entropy_shade_removal(Green, Red, Near, a=a)
         return M, R
     elif method == 'shi':
         M = shade_index(Blue, Green, Red, Near)
     elif method == 'sp':
-        ae = kwargs['ae'] if 'ae' in kwargs else 1e+1
-        be = kwargs['be'] if 'be' in kwargs else 5e-1
-        M = shadow_probabilities(Blue, Green, Red, Near, ae=ae, be=be)
+        M = shadow_probabilities(Blue, Green, Red, Near, **kwargs)
     elif method in (
             'nri',
             'normalized_range_shadow_index',
@@ -220,6 +249,8 @@ def modified_shadow_fraction(Blue, Green, Red, P_S=.95):
         green band of satellite image
     Red : numpy.ndarray, size=(m,n)
         red band of satellite image
+    P_S: float, range=0...1
+        threshold value for classification, based on the quantile
 
     Returns
     -------
@@ -243,7 +274,7 @@ def modified_shadow_fraction(Blue, Green, Red, P_S=.95):
     T_S = np.quantile(r, P_S)
     sig = np.std(r)
 
-    SF = np.exp(-np.divide((r - T_S)**2, 4 * sig))
+    SF = np.exp(-np.divide((r - T_S) ** 2, 4 * sig))
     np.putmask(SF, SF >= T_S, 1.)
     return SF
 
@@ -520,7 +551,7 @@ def normalized_difference_water_index(Green, Near):
     return NDWI
 
 
-def modified_normalized_difference_water_index(Green, Short):
+def modified_normalized_difference_water_index(Green, Shortwave):
     """transform green and shortwave infrared arrays to modified NDW-index.
     See also [Xu06]_.
 
@@ -550,17 +581,17 @@ def modified_normalized_difference_water_index(Green, Short):
               International journal of remote sensing, vol.27(14) pp.3025–3033,
               2006.
     """
-    are_two_arrays_equal(Green, Short)
+    are_two_arrays_equal(Green, Shortwave)
 
-    denom = (Green + Short)
-    MNDWI = np.divide((Green - Short),
+    denom = (Green + Shortwave)
+    MNDWI = np.divide((Green - Shortwave),
                       denom,
                       out=np.zeros_like(denom),
                       where=denom != 0)
     return MNDWI
 
 
-def normalized_difference_moisture_index(Near, Short):
+def normalized_difference_moisture_index(Near, Shortwave):
     """transform near and shortwave infrared arrays to NDM-index.
     See also [Wi02]_.
 
@@ -588,10 +619,10 @@ def normalized_difference_moisture_index(Near, Short):
               Landsat TM imagery" Remote sensing of environment, vol.80
               pp.385–396, 2002.
     """
-    are_two_arrays_equal(Near, Short)
+    are_two_arrays_equal(Near, Shortwave)
 
-    denom = (Near + Short)
-    NDMI = np.divide((Near - Short),
+    denom = (Near + Shortwave)
+    NDMI = np.divide((Near - Shortwave),
                      denom,
                      out=np.zeros_like(denom),
                      where=denom != 0)
@@ -1052,8 +1083,8 @@ def entropy_shade_removal(Ia, Ib, Ic, a=None):
         lower frequency band of satellite image
     Ic : numpy.array, size=(m,n)
         lowest frequency band of satellite image
-    a : float, default=-1
-        angle to use [in degrees]
+    a : float, default=-1, unit=degrees
+        angle of the coordinate rotation
 
     Returns
     -------
@@ -1218,7 +1249,7 @@ def fractional_range_shadow_index(*args):
     return SI
 
 
-def shadow_probabilities(Blue, Green, Red, Near, ae=1e+1, be=5e-1):
+def shadow_probabilities(Blue, Green, Red, Near, **kwargs):
     """transform blue, green, red and near infrared to shade probabilities, see
     also [FS10] and [Rü13].
 
@@ -1250,6 +1281,13 @@ def shadow_probabilities(Blue, Green, Red, Near, ae=1e+1, be=5e-1):
     are_three_arrays_equal(Blue, Green, Red)
     are_two_arrays_equal(Blue, Near)
 
+    ae = kwargs.get('ae')
+    if ae is None:
+        ae = signature(s_curve).parameters['ae'].default
+    be = kwargs.get('ae')
+    if be is None:
+        be = signature(s_curve).parameters['be'].default
+
     Fk = np.amax(np.stack((Red, Green, Blue), axis=2), axis=2)
     F = np.divide(np.clip(Fk, 0, 2), 2)  # (10) in [FS10]
     L = np.divide(Red + Green + Blue, 3)  # (4) in [FS10]
@@ -1263,7 +1301,6 @@ def shadow_probabilities(Blue, Green, Red, Near, ae=1e+1, be=5e-1):
     M = np.multiply(D, (1 - F))
     return M
 
-
 # recovery - normalized color composite
 
 # todo: Makaru 2011, for reflectance values (L2 data)

dhdt/processing/matching_tools.py

@@ -72,7 +72,7 @@ def get_integer_peak_location(C, metric=None):
     C = np.atleast_2d(C)
     max_corr = np.argmax(C)
 
-    ij = np.unravel_index(max_corr, C.shape, order='F') # 'C'
+    ij = np.unravel_index(max_corr, C.shape, order='F')  # 'C'
     if C.ndim >= 3: ij = ij[:2]
     di, dj = ij[::-1]
     di -= C.shape[0] // 2
@@ -160,12 +160,12 @@ def get_template(Z, idx_1, idx_2, radius):
         array with Gaussian peak in the center
     """
     sub_idx = np.mgrid[idx_1 - radius:idx_1 + radius + 1,
-                       idx_2 - radius:idx_2 + radius + 1]
+              idx_2 - radius:idx_2 + radius + 1]
 
     sub_ids = np.ravel_multi_index(np.vstack(
         (sub_idx[0].flatten(), sub_idx[1].flatten())),
-                                   Z.shape,
-                                   mode='clip')
+        Z.shape,
+        mode='clip')
     Z_sub = np.take(Z, sub_ids, mode='clip')
     Z_sub = np.reshape(Z_sub, (2 * radius + 1, 2 * radius + 1))
     return Z_sub
@@ -275,7 +275,7 @@ def pad_radius(Z, radius, cval=0):
         data array
     radius : {positive integer, tuple}
         extra boundary to be added to the data array
-    cval = {integer, flaot}
+    cval : {integer, flaot}
         conastant value to be used in the padding region
 
     Returns
@@ -284,12 +284,12 @@ def pad_radius(Z, radius, cval=0):
         extended data array
     """
     assert type(Z) in (np.ma.core.MaskedArray, np.ndarray), \
-        ("please provide an array")
+        "please provide an array"
     if not type(radius) is tuple:
         radius = (radius, radius)
 
     if Z.ndim == 3:
-        if type(Z) in (np.ma.core.MaskedArray, ):
+        if type(Z) in (np.ma.core.MaskedArray,):
             Z_xtra = np.ma.array(np.pad(Z, ((radius[0], radius[1]),
                                             (radius[0], radius[1]), (0, 0)),
                                         'constant',
@@ -305,7 +305,7 @@ def pad_radius(Z, radius, cval=0):
                         'constant',
                         constant_values=cval)
         return Z_xtra
-    if type(Z) in (np.ma.core.MaskedArray, ):
+    if type(Z) in (np.ma.core.MaskedArray,):
         Z_xtra = np.ma.array(np.pad(Z, ((radius[0], radius[1]),
                                         (radius[0], radius[1])),
                                     'constant',
@@ -372,9 +372,9 @@ def prepare_grids(im_stack, ds, cval=0):
 
 def make_templates_same_size(I1, I2):
     assert type(I1) in (np.ma.core.MaskedArray, np.ndarray), \
-        ("please provide an array")
+        "please provide an array"
     assert type(I2) in (np.ma.core.MaskedArray, np.ndarray), \
-        ("please provide an array")
+        "please provide an array"
 
     mt, nt = I1.shape[0], I1.shape[1]  # dimenstion of the template
     ms, ns = I2.shape[0], I2.shape[1]  # dimension of the search space
@@ -423,7 +423,7 @@ def remove_posts_outside_image(Z, i, j):
     --------
     remove_posts_pairs_outside_image
     """
-    assert Z.ndim >= 2, ('please provide an 2D array')
+    assert Z.ndim >= 2, 'please provide an 2D array'
 
     IN = np.logical_and.reduce((i >= 0, i < (Z.shape[0] - 1), j >= 0, j
                                 < (Z.shape[1] - 1)))
@@ -454,8 +454,8 @@ def remove_posts_pairs_outside_image(I1, i1, j1, I2, i2, j2):
     --------
     remove_posts_outside_image
     """
-    assert I1.ndim >= 2, ('please provide I1 as an 2D array')
-    assert I2.ndim >= 2, ('please provide I2 as an 2D array')
+    assert I1.ndim >= 2, 'please provide I1 as an 2D array'
+    assert I2.ndim >= 2, 'please provide I2 as an 2D array'
 
     IN = np.logical_and.reduce(
         (i1 >= 0, i1 < (I1.shape[0] - 1), j1 >= 0, j1 < (I1.shape[1] - 1), i2
@@ -495,10 +495,10 @@ def reposition_templates_from_center(I1, I2, di, dj):
 
     if I1.ndim == 3:
         I2sub = I2[mc - (mt // 2) - di:mc + (mt // 2) - di,
-                   nc - (nt // 2) - dj:nc + (nt // 2) - dj, :]
+                nc - (nt // 2) - dj:nc + (nt // 2) - dj, :]
     else:
         I2sub = I2[mc - (mt // 2) - di:mc + (mt // 2) - di,
-                   nc - (nt // 2) - dj:nc + (nt // 2) - dj]
+                nc - (nt // 2) - dj:nc + (nt // 2) - dj]
     return I1, I2sub
 
 
@@ -557,7 +557,7 @@ def get_coordinates_of_template_centers(Grid, temp_size):
 
     radius = np.floor(temp_size / 2).astype('int')
     I_idx, J_idx = np.mgrid[radius:(m - radius):temp_size,
-                            radius:(n - radius):temp_size]
+                   radius:(n - radius):temp_size]
     return I_idx, J_idx
 
 
@@ -586,7 +586,7 @@ def get_value_at_template_centers(Grid, temp_size):
     assert isinstance(temp_size, int), "please provide an integer"
 
     Iidx, Jidx = get_coordinates_of_template_centers(Grid, temp_size)
-    if (Grid.ndim == 3):
+    if Grid.ndim == 3:
         m, n = Iidx.shape
         b = Grid.shape[2]
         Iidx = np.tile(np.atleast_3d(Iidx), (1, 1, b))

dhdt/testing/classification_tools.py

@@ -0,0 +1,7 @@
+import numpy as np
+
+
+def calculate_accuracy(Predicted, Truth):
+    acc = 1 - np.divide(np.sum(np.logical_xor(Predicted, Truth)),
+                        Truth.size)
+    return acc