2021-08-17-1st

PBR_filter.py

@@ -0,0 +1,164 @@
+# Written by Dr Daniel Buscombe, Marda Science LLC
+# for the USGS Coastal Change Hazards Program
+#
+# MIT License
+#
+# Copyright (c) 2021, Marda Science LLC
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from tkinter import filedialog
+from tkinter import *
+import os, sys, getopt
+import numpy as np
+from imageio import imread, imwrite
+import matplotlib.pyplot as plt
+import numpy as np
+from skimage.color import rgb2hsv, hsv2rgb
+# from glob import glob
+from skimage.restoration import denoise_wavelet, estimate_sigma, rolling_ball
+
+from functools import partial
+# rescale_sigma=True required to silence deprecation warnings
+_denoise_wavelet = partial(denoise_wavelet, rescale_sigma=True)
+from skimage import util
+
+# import warnings filter
+from warnings import simplefilter
+# ignore all future warnings
+simplefilter(action='ignore', category=RuntimeWarning)
+
+# =========================================================
+def rescale(dat,mn,mx):
+    """
+    rescales an input dat between mn and mx
+    """
+    m = np.nanmin(dat.flatten())
+    M = np.nanmax(dat.flatten())
+    return (mx-mn)*(dat-m)/(M-m)+mn
+
+# =========================================================
+def sharpen(Z, radius, do_plot):
+
+    sigma_est = estimate_sigma(Z, multichannel=True, average_sigmas=False)
+    region = denoise_wavelet(Z, multichannel=True, rescale_sigma=True, wavelet_levels=6, convert2ycbcr=True,
+                              method='BayesShrink', mode='soft', sigma=np.max(sigma_est)*5)
+    original = rescale(region,0,255)
+
+    Zo = np.ma.filled(original, fill_value=np.nan).copy()
+    hsv = rgb2hsv(Zo)
+    im = (0.299 * Zo[:,:,0] + 0.5870*Zo[:,:,1] + 0.114*Zo[:,:,2])
+    im[Z[:,:,0]==0]=0
+
+    ##https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_rolling_ball.html#sphx-glr-auto-examples-segmentation-plot-rolling-ball-py
+    ##background = rolling_ball(im, radius=100)
+    image_inverted = util.invert(im)
+    background_inverted = rolling_ball(image_inverted, radius=radius)
+    filtered_image_inverted = image_inverted - background_inverted
+    filtered_image = util.invert(filtered_image_inverted)
+    background = util.invert(background_inverted)
+
+    background[np.isnan(background)] = 0
+    background[np.isinf(background)] = 0
+    intensity = (im/background)
+    intensity[np.isnan(intensity)] = 0
+    intensity[np.isinf(intensity)] = 0
+    intensity = (255*intensity).astype('uint8')
+
+    sharpened = hsv2rgb(np.dstack([hsv[:,:,0], hsv[:,:,1], intensity]))
+
+    sharpened[:,:,0] = rescale(sharpened[:,:,0],Z[:,:,0].min(), Z[:,:,0].max())
+    sharpened[:,:,1] = rescale(sharpened[:,:,1],Z[:,:,1].min(), Z[:,:,2].max())
+    sharpened[:,:,2] = rescale(sharpened[:,:,2],Z[:,:,2].min(), Z[:,:,1].max())
+    sharpened = (sharpened).astype('uint8')
+
+    if do_plot:
+        plt.figure(figsize=(12,12))
+        plt.subplot(231); plt.imshow(Z); plt.axis('off'); plt.title('a)', loc='left')
+        plt.subplot(232); plt.imshow(original.astype('uint8')); plt.axis('off'); plt.title('b)', loc='left')
+
+        plt.subplot(233); plt.imshow(background, cmap='gray'); plt.axis('off'); plt.title('c)', loc='left')
+
+        plt.subplot(234); plt.imshow(im, cmap='gray'); plt.axis('off'); plt.title('d)', loc='left')
+        plt.subplot(235); plt.imshow(intensity, cmap='gray'); plt.axis('off'); plt.title('e)', loc='left')
+
+        plt.subplot(236); plt.imshow(sharpened, cmap='gray'); plt.axis('off'); plt.title('f)', loc='left')
+
+        # plt.show()
+        plt.savefig(f.replace('.jpg','_filt_fig_breakdown.png'), dpi=300, bbox_inches='tight')
+        plt.close()
+
+    return sharpened.astype('uint8')
+
+
+#============================================================
+# =========================================================
+
+def do_filter(f, radius, do_plot):
+    Z = imread(f)
+    #radius = 3
+    sharpened = sharpen(Z, radius, do_plot)
+    imwrite(f.replace('.jpg','_filt.png'),sharpened)
+
+    if do_plot:
+        plt.subplot(221); plt.imshow(Z); plt.axis('off')
+        plt.subplot(222); plt.imshow(sharpened); plt.axis('off')
+        plt.subplot(223); plt.imshow(Z[:250,:250,:]); plt.axis('off')
+        plt.subplot(224); plt.imshow(sharpened[:250,:250,:]); plt.axis('off')
+        # plt.show()
+        plt.savefig(f.replace('.jpg','_filt_fig.png'), dpi=300, bbox_inches='tight')
+        plt.close()
+
+###==================================================================
+#===============================================================
+if __name__ == '__main__':
+
+    argv = sys.argv[1:]
+    try:
+        opts, args = getopt.getopt(argv,"h:r:p:")
+    except getopt.GetoptError:
+        print('python PBR_filter.py -r radius (px) -doplot 0 (1)')
+        sys.exit(2)
+    for opt, arg in opts:
+        if opt == '-h':
+            print('Example usage python PBR_filter.py -r 5')
+            print('Example usage python PBR_filter.py -r 6 -p 1')
+            sys.exit()
+        elif opt in ("-r"):
+            radius = arg
+            radius = int(radius)
+        elif opt in ("-p"):
+            do_plot = arg
+            do_plot = int(do_plot)
+
+    if 'do_plot' not in locals():
+        do_plot = 0
+    if 'radius' not in locals():
+        radius = 3
+    print("PBR: pan-sharpen with background subtraction and radius %i" % (radius))
+
+    #files = glob('*.jpg')
+    root = Tk()
+    files = filedialog.askopenfilenames(initialdir = "./",title = "Select image file",filetypes = (("image file","*.jpg"),("all files","*.*")))
+    root.withdraw()
+    print("%s files selected" % (len(files)))
+
+    for f in files:
+        print(f)
+        do_filter(f,radius, do_plot)

README.md

@@ -1,2 +1,72 @@
-# PBR_filter
+# "PBR" filter
 {P}ansharpening by {B}ackground {R}emoval algorithm for sharpening RGB images
+
+
+1. Read image, wavelet denoise, and convert to HSV
+2. Do 'inverted background subtraction' on the V (intensity) channel
+3. combine with HS, convert to RGB
+
+The effect is to sharpen details of object boundaries/transitions, brighten, and recolour
+
+But without changing the overall distribution of values within the image (i.e. the range), and recolouring in an internally consistent (deterministic) way
+
+Sample 'Madeira' imagery comes from https://coastal.er.usgs.gov/data-release/doi-P9L474WC/
+
+![](PBR.jpg)
+
+Sample 'OBX' imagery comes from https://www.sciencebase.gov/catalog/item/6037cca0d34eb12031175133
+
+![](PBR2.jpg)
+
+0. Create a conda environment
+
+A. Conda housekeeping
+
+`conda clean --all`
+`conda update -n base -c defaults conda`
+
+B. Create new `pbr` conda environment
+
+We'll create a new conda environment and install packages into it from conda-forge
+
+`conda env create -f install/pbr.yml`
+`conda activate pbr`
+
+1. Run the program
+
+`python PBR_filter.py`
+
+To get thicker dark lines use
+`python PBR_filter.py  -r 5`
+
+To get even thicker dark lines use something like
+`python PBR_filter.py  -r 7`
+
+If you want to print out two more plots showing the process step-by-step for each image, use :
+
+`python PBR_filter.py -p 1`
+
+navigate to a folder of images you wish to filter, and it will step through them one by one
+
+2. Visualise your results
+see the png images the program makes inside the same directory as the input images
+
+## How does this work?
+
+Let's take this image, which is a tile cropped out of a larger orthomosaic - see [here](https://coastal.er.usgs.gov/data-release/doi-P9L474WC/):
+![][example/20180619_MadeiraBeachFL_ortho_5cm_10_12.jpg]
+
+The filter with default radius of 3 pixels creates this PBR image
+![][example/20180619_MadeiraBeachFL_ortho_5cm_10_12_filt.png]
+
+In the figure below, the process is broken into stages
+![][example/20180619_MadeiraBeachFL_ortho_5cm_10_12_filt_fig_breakdown.png]
+
+a) original image
+b) wavelet denoised image, where noise over a range of scales is removed and mostly affects very small scale (pixel level) noise. This step isnt crucial but I always like to denoise imagery if I can as a precaution
+c) the greyscale background image that has been created with a rolling ball filter with ball of radius [whatever]
+d) greyscale version of the denoised image
+e) the intensity image that is the greyscale divided by the greyscale background image
+f) the filtered RGB image that is the result of swapping the greyscale with the intensity image in the HSV stack of the original RGB image, then converting that into RGB colorspace
+
+*disclaimer*: I do not know if I have reinvented the wheel - I have not searched for similar implementations. Please tell me by opening an Issue if this technique has previously been proposed

install/pbr.yml

@@ -0,0 +1,12 @@
+name: pbr
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+ - python
+ - scipy
+ - numpy
+ - scikit-image>=0.18.2
+ - ipython
+ - joblib
+ - tqdm