Thresholding

import glob
import scipy
import scikitplot
import matplotlib
import numpy as np
import imageio
import rasterio
import skimage
import imagecodecs
from matplotlib import pyplot as plt
from skimage import data
from skimage import io as io
from skimage.filters import threshold_multiotsu, threshold_otsu, threshold_yen
from skimage import io
from skimage import viewer
from skimage import color
from skimage import img_as_ubyte
from skimage import img_as_float64
from skimage import img_as_uint
from skimage import img_as_int


# input_dir = "/mnt/DD/Bonds_Project/Aurea_Project/Input/S1A/"
images_list = []
s1list = []


# OOB algorithm doing multi otsu thresholding on multiple image objects Sentinel 1


class Image:
    def __init__(self, path):
        # initialize mandatory variables
        # path is /path/to/file
        # array is the numpy array
        # timestamp is for managing chronogically files
        self.path = path
        # self.array = io.imread(self.path)
        self.timestamp = path.split('_')[6:7]
        intermediatevalue = ''.join(self.timestamp)
        self.timestamp = intermediatevalue.split('.')[:1]

    def givePath(self, path):
        # give path to object Image eg : a.givePath(path)
        self.path = path
        return self.path

    def giveTimestamp(self, timestamp):
        # give timestamp to object Image eg : im1.giveTimestamp(timestamp)
        self.timestamp = timestamp
        return self.timestamp

    def createTimestamp(self, path):
        self.timestamp = path.split('_')[6:7]
        return self.timestamp

    def readImage(self, path):
        # transform image from a string list (path) to a numpy array readable (array) eg : im1.readImage(path)
        if not self.path:
            # If Path was set beforehand , just read it
            self.array = io.imread(path)
        else:
            # otherwise set it
            self.array = io.imread(self.path)
        return self.array

    def determineTimeStamp(self):
        if not self.timestamp:
            intermediate_timestamp = self.path.split('dB')[-1:]
            self.timestamp = ''.join(intermediate_timestamp).split('.')[:1]
            print(self.timestamp)
        elif not (len(self.timestamp) == 15):
            print('correcting values for time stamp')
            intermediate_timestamp = self.path.split('dB')[-1:]
            self.timestamp = ''.join(intermediate_timestamp).split('.')[:1]
            print(self.timestamp)
        else:
            print('time stamp already set')
            pass

    def getPath(self):
        return self.path  # method to return a string like variable ('/mnt/DD/etc)

    def getArray(self):
        return self.array  # method to return a numpy array

    def getTimestamp(self):
        # get timestamp from object Image
        return self.timestamp

    def do_otsu(self):
        # Single otsu binary threshold
        self.getTimestamp()
        threshold = threshold_otsu(io.imread(self.path, as_gray=True), 256)
        print('threshold is otsu ' + str(threshold))
        regions = np.digitize(io.imread(self.path), bins=256)
        output = img_as_float64(regions)
        imageio.imwrite('/mnt/DD/Bonds_Project/Aurea_Project/Input/' + ''.join(self.getTimestamp()) + "binaryotsu.tiff",
                   output)
        print(''.join(self.getTimestamp()) + ' is done')

    # def do_multi_otsu(self,classe=3):
    #     # First look at images
    #     # inputIm = rgb2gray(io.imread(self.path, plugin='gdal'))
    #     print('image has been read, determining thresholds')
    #     threshold = threshold_multiotsu(io.imread(self.path, as_gray=True), classe, 256)
    #     print('threshold is multi otsu ' + str(threshold))
    #     regions = np.digitize(io.imread(self.path), bins=threshold)
    #     print('image has been digitized')
    #     outputname = '/mnt/DD/Bonds_Project/Aurea_Project/Input/S1B/' + ''.join(self.timestamp) + '_otsu_thresh.tiff'
    #     print('outputname: ' + outputname)
    #     output = img_as_uint(regions)
    #     imageio.imwrite(outputname, output)
    #     # io.imsave(outputname, output)
    #     print(''.join(self.getTimestamp()) + ' is done')

    def do_multi_otsu(self, classe=3):
        # First look at images
        # inputIm = rgb2gray(io.imread(self.path, plugin='gdal'))
        # Use rasterio to keep georeferencing
        # First read image
        with rasterio.open(self.path) as src:
            img = src.read()
            profile = src.profile
        # Determining thresholds using scikit and multi otsu
        print('image has been read, determining thresholds')
        threshold = threshold_multiotsu(img, classe, 256)
        print('threshold is multi otsu ' + str(threshold))
        # Digitizing using numpy
        regions = np.digitize(img, bins=threshold)
        print('image has been digitized')
        # determining output name
        outputname = '/mnt/DD/Bonds_Project/Aurea_Project/Input/S1B/' + ''.join(self.timestamp) + \
                     '_'+str(classe)+'_otsu_thresh.tiff'
        print('outputname: ' + outputname)
        output = img_as_ubyte(regions)
        # Writing output with rasterio with same CRS as input
        with rasterio.open(outputname, 'w', **profile) as dst:
            dst.write(output)
        # io.imsave(outputname, output)
        print(''.join(self.getTimestamp()) + ' is done')

    def do_binarize(self):
        input_array = io.imread(self.path)


def multiotsu(image, classe=3, nbins=256):
    print('reading image')
    img = io.imread(image,as_gray=True)
    print('image is ready')
    print('doing thresholding')
    # img[img==0]=[np.nan]
    # img = img[np.logical_not(np.isnan(img))]
    threshold = threshold_multiotsu(img, classe, nbins)
    print('number of class ' + str(classe) + '  number of bins ' + str(nbins) + ' ' + str(threshold))
    regions = np.digitize(img, bins=threshold)
    output = img_as_int(regions)
    #plt.imsave('/mnt/DD/Bonds_Project/Aurea_Project/Input/multiotsu6.tiff', output)
    fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(10, 3.5))
    ax[0].imshow(img, cmap='gray')
    ax[0].set_title('Original')
    ax[0].axis('off')

    # Plotting the histogram and the two thresholds obtained from
    # multi-Otsu.
    print('ploting the histogram')
    ax[1].hist(img.ravel(), bins=256)
    ax[1].set_title('Histogram')
    for thresh in threshold:
        ax[1].axvline(thresh, color='r')

    # Plotting the Multi Otsu result.i
    print('ploting the multi otsu result')
    ax[2].imshow(regions, cmap='Accent')
    ax[2].set_title('Multi-Otsu result')
    ax[2].axis('off')
    # viewer = skimage.viewer.ImageViewer(regions)
    # viewer.show()
    plt.subplots_adjust()
    plt.show()
    print('saving image')
    #outname= image.split('_2020')
    #plt.imsave('/home/marielouise/Images/'+image+'_MapOtsu.tiff', output)
    imageio.imwrite('/home/marielouise/Images/MapOtsu.tiff', output)

def imageSearchtif(folder):
    for images in glob.glob(folder + "*.tif"):
        images_list.append(images)
    print("listes d'images: " + str(images_list))
    return images_list


def run(input_dir=''):
    # We want to have 2 lists in order to create Image objects using path values ,
    arraylist = []
    res = []
    if input_dir:
        folder = input_dir
    else:
        folder = "/mnt/DD/Bonds_Project/Aurea_Project/Input/2017A/"
    imageSearch(folder, extension='tif')  # The goal here is to fetch a list of path and give it to our variable object
    for path in images_list:
        # a list  named arraylist will be created and will get line of path for each image into the path attribute
        arraylist.append(Image(path))
    # Setting attributes for the multiple object we are going to otsu thresh
    for objet in arraylist:
        # objet.readImage(objet.getPath())
        objet.determineTimeStamp()
        print(objet.getPath())
        res.append(objet)
    # doing the multi otsu thresholding
    for array in res:
        print('Doing this one : ' + ''.join(array.getTimestamp()))
        array.do_multi_otsu()


def imageSearch(folder, extension):
    # images_list = []
    for images in glob.glob(folder + '/*.' + extension):
        images_list.append(images)
    print("listes d'images: " + str(images_list))
    return images_list


def createS1Imglist(inputfolder='', extension='zip'):
    s1list.clear()
    if inputfolder:
        folder = inputfolder
    else:
        folder = "/mnt/DD/Bonds_Project/Aurea_Project/Input/S1A/"
    imageSearch(folder, extension)  # The goal here is to fetch a list of path and give it to our variable object
    for path in images_list:
        # a list  named s1list will be created and will get line of path for each image into the path attribute
        s1list.append(Image(path))
    return s1list, images_list


def tableau(folder=''):
    if not folder:
        folder = "/mnt/DD/Bonds_Project/Aurea_Project/Input/S1A/"
    print("folder is " + str(folder))
    print("looking for images ")
    imageslist = []
    timestamp = []
    for images in glob.glob(folder + '*.zip'):
        imageslist.append(images)
    for images in imageslist:
        timestamp.append(images.split('_')[6:7])