add files

Canadian-Severe-Storms-Laboratory · Jun 10, 2024 · a7fe85b · a7fe85b
1 parent 433da2f
commit a7fe85b
Show file tree

Hide file tree

Showing 7 changed files with 265 additions and 52 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,52 +1,4 @@
-# Prerequisites
-*.d
-
-# Object files
-*.o
-*.ko
-*.obj
-*.elf
-
-# Linker output
-*.ilk
-*.map
-*.exp
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Libraries
-*.lib
-*.a
-*.la
-*.lo
-
-# Shared objects (inc. Windows DLLs)
-*.dll
-*.so
-*.so.*
-*.dylib
-
-# Executables
-*.exe
-*.out
-*.app
-*.i*86
-*.x86_64
-*.hex
-
-# Debug files
-*.dSYM/
-*.su
-*.idb
-*.pdb
-
-# Kernel Module Compile Results
-*.mod*
-*.cmd
-.tmp_versions/
-modules.order
-Module.symvers
-Mkfile.old
-dkms.conf
+test_images/
+simple_images/
+python/__pycache__
+python/.idea/
diff --git a/python/colour_swaths.py b/python/colour_swaths.py
@@ -0,0 +1,29 @@
+import os
+import numpy as np
+import cv2
+from PIL import Image
+from glob import glob
+from tqdm import tqdm
+
+colour_table = np.array([[255, 255, 255], [255, 204, 153], [255, 153, 0], [102, 255, 0], [102, 255, 0], [0, 204, 0], [0, 204, 0], [0, 153, 0], [0, 153, 0], [0, 102, 0], [0, 102, 0], [51, 255, 255], [51, 255, 255], [51, 255, 255], [51, 255, 255], [51, 255, 255], [0, 204, 255], [0, 204, 255], [0, 204, 255], [0, 204, 255], [0, 204, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 153, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 102, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [0, 0, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [153, 2, 255], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [204, 51, 153], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [153, 0, 102], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100], [100, 100, 100]]
+                        , np.uint8)
+
+
+def colour_swath(img):
+    if not isinstance(img, np.ndarray) or img.dtype != np.uint8 or len(img.shape) != 2:
+        raise TypeError('image must be single chanel np.uint8')
+
+    return colour_table[img]
+
+
+if __name__ == '__main__':
+
+    files = glob('simple_images/*.tif')
+
+    for f in tqdm(files):
+        name = os.path.basename(f).split(".")[0]
+        img = np.array(Image.open(f), np.uint8)
+
+        #img = colour_swath(img)
+
+        cv2.imwrite("simple_images/" + name + ".png", img)
diff --git a/python/denoise.py b/python/denoise.py
@@ -0,0 +1,46 @@
+import sys
+
+import numpy as np
+import cv2
+from main import imread_tif32f
+from colour_swaths import colour_swath
+
+
+def remove_artifacts(im):
+
+    artifact_polygons = [[[2084, 1653], [2284, 1634], [2281, 1614]],
+                         [[4608, 2175], [4654, 2284], [4623, 2187]],
+                         [[4752, 2576], [4769, 2651], [4778, 2648]]]
+
+    for polygon in artifact_polygons:
+        cv2.fillPoly(im, [np.array(polygon)], 0)
+
+
+def denoise_mesh(im):
+    img = im.copy()
+    remove_artifacts(img)
+
+    mask = cv2.blur(img, (3, 3))
+
+    mask = cv2.inRange(mask, 6, 255)
+
+    mask2 = np.zeros_like(mask)
+
+    num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(mask, connectivity=8)
+
+    for i in range(1, num_labels):  # Start from 1 to skip the background component
+        if stats[i, cv2.CC_STAT_AREA] >= 16:
+            mask2[labels == i] = 255
+
+    return cv2.bitwise_and(img, img, mask=mask2)
+
+
+if __name__ == '__main__':
+    img = imread_tif32f("simple_images/20230715.tif")
+
+    img = denoise_mesh(img)
+
+    cimg = colour_swath(img)
+
+    cv2.imwrite("denoised.png", cimg)
+
diff --git a/python/denoised.png b/python/denoised.png
diff --git a/python/found.png b/python/found.png
diff --git a/python/img_name.py b/python/img_name.py
@@ -0,0 +1,10 @@
+import os
+from glob import glob
+
+if __name__ == '__main__':
+    files = glob('simple_images/*.tif')
+
+    for f in files:
+        #print("simple_images/" + f.split('\\')[1].split('_')[0] + ".tif")
+        os.rename(f, "simple_images/" + f.split('\\')[1].split('.')[0] + ".tif")
+
diff --git a/python/main.py b/python/main.py
@@ -0,0 +1,176 @@
+# MESH-MAP, (Maximum Estimated Size of Hail - Monitoring and Analysis Program)
+
+import math
+import os
+
+import numpy as np
+import cv2
+import networkx as nx
+from colour_swaths import colour_swath
+from denoise import denoise_mesh
+from PIL import Image
+from glob import glob
+from tqdm import tqdm
+
+
+def imread_tif32f(path):
+    return np.array(Image.open(path), dtype=np.uint8)
+
+
+def extract_cc_images(img, bimg, connectivity, area_threshold):
+    (num_labels, labels_im, stats, centroids) = cv2.connectedComponentsWithStats(bimg, connectivity=connectivity,
+                                                                                 ltype=cv2.CV_32S)
+    # List to store cropped component images
+    cropped_images = []
+    cropped_stats = []
+    cropped_centroids = []
+
+    # Iterate through the components
+    for i in range(1, num_labels):  # Start from 1 to skip the background
+        area = stats[i, cv2.CC_STAT_AREA]
+
+        if area >= area_threshold:
+            x = stats[i, cv2.CC_STAT_LEFT]
+            y = stats[i, cv2.CC_STAT_TOP]
+            w = stats[i, cv2.CC_STAT_WIDTH]
+            h = stats[i, cv2.CC_STAT_HEIGHT]
+
+            # Crop the bounding box of the component
+            component_image = img[y:y + h, x:x + w]
+
+            # Create a mask for the current component
+            component_mask = (labels_im[y:y + h, x:x + w] == i).astype(np.uint8) * 255
+
+            # Apply the mask to the cropped image
+            filtered_component_image = cv2.bitwise_and(component_image, component_image, mask=component_mask)
+
+            cropped_images.append(filtered_component_image)
+            cropped_stats.append(stats[i])
+            cropped_centroids.append(centroids[i])
+
+    return cropped_images, cropped_stats, cropped_centroids
+
+
+def create_cc_graph(img):
+    graph = nx.Graph()
+
+    for i in range(img.shape[0]):
+        for j in range(img.shape[1]):
+
+            if img[i, j] != 0:
+
+                for y in range(max(0, i - 1), min(img.shape[0], i + 2)):
+                    for x in range(max(0, j - 1), min(img.shape[1], j + 2)):
+
+                        if y == i and x == j:
+                            continue
+
+                        if img[y, x] != 0:
+                            graph.add_edge((i, j), (y, x), weight=math.hypot(x - j, y - i))
+
+    return graph
+
+
+def path_length(graph, path):
+    return sum(graph[u][v]['weight'] for u, v in zip(path[:-1], path[1:]))
+
+
+def weighted_diameter(G):
+    shortest_paths = dict(nx.all_pairs_dijkstra_path(G))
+
+    # Find the longest of the shortest paths
+    longest_shortest_path_length = 0
+    longest_shortest_path = None
+
+    for source, paths in shortest_paths.items():
+        for target, path in paths.items():
+            length = path_length(G, path)
+            if length > longest_shortest_path_length:
+                longest_shortest_path_length = length
+                longest_shortest_path = path
+
+    return longest_shortest_path_length, longest_shortest_path
+
+
+def long_enough_subsection(image, lb, ub, length):
+    bimage = cv2.inRange(image, lb, ub)
+    subcomp_images, _, _ = extract_cc_images(image, bimage, 8, length)
+
+    for sub_image in subcomp_images:
+        g = create_cc_graph(sub_image)
+
+        wd, path = weighted_diameter(g)
+
+        if wd >= length:
+            return path
+
+    return None
+
+
+def draw_referenced_contour(im, comp_im, stats):
+    contour = cv2.findContours(comp_im, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[0][0]
+    x = stats[cv2.CC_STAT_LEFT]
+    y = stats[cv2.CC_STAT_TOP]
+
+    offset_contours = [np.array(contour, int) + np.array([x, y])]
+
+    cv2.drawContours(im, offset_contours, -1, (0, 0, 0), 1)
+
+
+def draw_referenced_box(im, stats):
+    x = stats[cv2.CC_STAT_LEFT]
+    y = stats[cv2.CC_STAT_TOP]
+    w = stats[cv2.CC_STAT_WIDTH]
+    h = stats[cv2.CC_STAT_HEIGHT]
+
+    cv2.rectangle(im, (x - 1, y - 1), (x + w + 1, y + h + 1), (0, 0, 0), 1)
+
+
+def draw_referenced_path(im, path, stats):
+    x = stats[cv2.CC_STAT_LEFT]
+    y = stats[cv2.CC_STAT_TOP]
+
+    # [y, x] -> [x, y]
+    path = [[p[1], p[0]] for p in path]
+
+    path = np.array(path) + np.array([x, y])
+
+    cv2.polylines(im, [path], False, (0, 0, 0), 1)
+
+
+# TODO: Overall length rounded to .5
+# TODO: Bounding box coords
+# TODO: Max MESH value in swath
+# TODO: Max swath width
+# TODO: Max swath width path coords
+if __name__ == '__main__':
+
+    for file in tqdm(glob("simple_images/*.tif")[21:]):
+        name = os.path.basename(file)
+        im = imread_tif32f(file)
+        im = denoise_mesh(im)
+        cimg = colour_swath(im)
+
+        binary_im = np.zeros_like(im, dtype=np.uint8)
+        cv2.threshold(im, 10, 255, cv2.THRESH_BINARY, binary_im)
+
+        comp_images, stats, centroids = extract_cc_images(im, binary_im, 8, 40)
+
+        swath_idx = []
+        swath_paths = []
+
+        for i, image in enumerate(comp_images):
+
+            req1 = long_enough_subsection(image, 10, 255, 40)
+            req2 = long_enough_subsection(image, 30, 255, 2)
+
+            if req1 is not None and req2 is not None:
+                swath_idx.append(i)
+                swath_paths.append([req1, req2])
+
+        for i, s in enumerate(swath_idx):
+            #draw_referenced_path(cimg, swath_paths[i][0], stats[s])
+            draw_referenced_contour(cimg, comp_images[s], stats[s])
+            draw_referenced_box(cimg, stats[s])
+
+        cv2.imwrite("simple_images/results1/" + name, cimg)