confidence_map_exploitation.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Tool to generate reference cloud masks for validation of operational cloud masks.
The elaboration is performed using an active learning procedure.

The code was written by Louis Baetens during a training period at CESBIO, funded by CNES, under the direction of O.Hagolle

==================== Copyright
Software (confidence_map_exploitation.py)

Copyright© 2019 Centre National d’Etudes Spatiales

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License version 3
as published by the Free Software Foundation.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this program.  If not, see
https://www.gnu.org/licenses/gpl-3.0.fr.html
"""
import os
import sys
import os.path as op
import json
import glob
import otbApplication
import numpy as np
import tempfile
from collections import defaultdict
import matplotlib.pyplot as plt

import merge_shapefiles


def confidence_map_change(app_confidence_map, out_tif, median_radius=5):
    ''' 
    Change the confidence map to an enhanced one
    '''
    # If number is even, add 1 to be odd
    if median_radius % 2 == 0:
        median_radius = median_radius+1

    MedianFilter = otbApplication.Registry.CreateApplication("BandMathX")
    MedianFilter.AddImageToParameterInputImageList("il", app_confidence_map)
    MedianFilter.SetParameterString("out", str(out_tif))
    MedianFilter.SetParameterString(
        "exp", "(median(im1b1N{}x{}))".format(median_radius, median_radius))
    MedianFilter.UpdateParameters()
    MedianFilter.ExecuteAndWriteOutput()

    return


def shapefile_rasterization(raw_img_tif, in_shps, out_tif):
    '''
    Rasterize one or two shapefile. Use for the points classification 
    to a raster.
    If two shapefiles are entered, merge the two tif afterwards
    in_shps can also be a directory, containing the various class
    layers (water.shp, land.shp, etc)
    '''
    if isinstance(in_shps, list):
        if len(in_shps) == 1:
            Rasterization = otbApplication.Registry.CreateApplication("Rasterization")
            Rasterization.SetParameterString("in", str(in_shps[0]))
            Rasterization.SetParameterString("im", str(raw_img_tif))
            Rasterization.SetParameterString("mode", "attribute")
            Rasterization.UpdateParameters()
            Rasterization.SetParameterString("mode.attribute.field", "class")
            Rasterization.SetParameterString("out", str(out_tif))
            Rasterization.UpdateParameters()
            Rasterization.ExecuteAndWriteOutput()
        elif len(in_shps) == 2:
            Rasterization1 = otbApplication.Registry.CreateApplication("Rasterization")
            Rasterization1.SetParameterString("in", str(in_shps[0]))
            Rasterization1.SetParameterString("im", str(raw_img_tif))
            Rasterization1.SetParameterString("mode", "attribute")
            Rasterization1.UpdateParameters()
            Rasterization1.SetParameterString("mode.attribute.field", "class")
            Rasterization1.UpdateParameters()
            Rasterization1.Execute()

            Rasterization2 = otbApplication.Registry.CreateApplication("Rasterization")
            Rasterization2.SetParameterString("in", str(in_shps[1]))
            Rasterization2.SetParameterString("im", str(raw_img_tif))
            Rasterization2.SetParameterString("mode", "attribute")
            Rasterization2.UpdateParameters()
            Rasterization2.SetParameterString("mode.attribute.field", "class")
            Rasterization2.UpdateParameters()
            Rasterization2.Execute()

            Combination = otbApplication.Registry.CreateApplication("BandMathX")
            Combination.AddImageToParameterInputImageList(
                "il", Rasterization1.GetParameterOutputImage("out"))
            Combination.AddImageToParameterInputImageList(
                "il", Rasterization2.GetParameterOutputImage("out"))
            Combination.SetParameterString("out", str(out_tif))
            Combination.SetParameterString("exp", "im1b1 + im2b1")
            Combination.UpdateParameters()
            Combination.ExecuteAndWriteOutput()

        else:
            print('Please enter 1 or 2 shapefiles')

    else:
        tmp_shp = op.join('tmp', next(tempfile._get_candidate_names()) + '.shp')
        merge_shapefiles.merge_all_types(in_shps, tmp_shp)

        Rasterization = otbApplication.Registry.CreateApplication("Rasterization")
        Rasterization.SetParameterString("in", str(tmp_shp))
        Rasterization.SetParameterString("im", str(raw_img_tif))
        Rasterization.SetParameterString("mode", "attribute")
        Rasterization.UpdateParameters()
        Rasterization.SetParameterString("mode.attribute.field", "class")
        Rasterization.SetParameterString("out", str(out_tif))
        Rasterization.UpdateParameters()
        Rasterization.ExecuteAndWriteOutput()

        tmps_files = glob.glob(op.abspath(tmp_shp)[0:-4] + '*')
        for tmp_file in tmps_files:
            os.remove(tmp_file)


def confidence_map_mean(global_parameters, mode='all', samples_set='train', extended=True):
    '''
    Compute the mean of the confidence map based on some filtering parameters

    mode can be:
    - all : will select all the pixels of the image (in this case, samples_set has no importance)
    - all_classified_samples : select only the pixels corresponding to all the classified samples
    - well_classified_samples : select only the pixels corresponding to the well classified samples
    - misclassified_samples : select only the pixels corresponding to the misclassified samples

    samples_set can be:
    - both : both training and validation sample set
    - train : only training sample set
    - validation : only validation sample set    
    '''
    main_dir = global_parameters["user_choices"]["main_dir"]

    # If select only samples, have to create a raster before of the samples
    if mode != 'all':
        if extended:
            train_points_shp = global_parameters["general"]["training_shp_extended"]
            validation_points_shp = global_parameters["general"]["validation_shp_extended"]
        else:
            train_points_shp = global_parameters["general"]["training_shp"]
            validation_points_shp = global_parameters["general"]["validation_shp"]

        in_shps = []
        if samples_set == 'both':
            in_shps.append(train_points_shp)
            in_shps.append(validation_points_shp)
        elif samples_set == 'train':
            in_shps.append(train_points_shp)
        elif samples_set == 'validation':
            in_shps.append(validation_points_shp)

        in_shps = [str(op.join(main_dir, 'Intermediate', i)) for i in in_shps]
        raw_img_tif = op.join(main_dir, 'In_data', 'Image',
                              global_parameters["user_choices"]["raw_img"])
        # Following can and should be changed
        rasterized_selection_tif = op.join(main_dir, 'Intermediate', 'rasterized_samples.tif')
        print(in_shps)
        shapefile_rasterization(raw_img_tif, in_shps, rasterized_selection_tif)

    # Order of images: im1: confidence, im2: classification, im3: samples selection
    if mode == 'all':
        expression = "im1b1"  # all
    elif mode == 'all_classified_samples':
        expression = "im3b1 != 0 ? im1b1 : 0"  # all_samples
    elif mode == 'well_classified_samples':
        expression = "im3b1 != 0 and im3b1 == im2b1 ? im1b1 : 0"  # well_classified_samples
    elif mode == 'misclassified_samples':
        expression = "im3b1 != 0 and im3b1 != im2b1 ? im1b1 : 0"  # misclassified_samples
    else:
        print('Please enter a valid mode')

    confidence_map = op.join(main_dir, 'Out', 'confidence.tif')
    classification_map = op.join(
        main_dir, 'Out', global_parameters["general"]["img_labeled_regularized"])

    temp_out = op.join(main_dir, 'Intermediate', 'confidence_selected.tif')

    ConfMapSelection = otbApplication.Registry.CreateApplication("BandMathX")
    if mode == 'all':
        ConfMapSelection.SetParameterStringList("il",
                                                [str(confidence_map)])
    else:
        ConfMapSelection.SetParameterStringList("il",
                                                [str(confidence_map), str(classification_map), str(rasterized_selection_tif)])
    ConfMapSelection.SetParameterString("out", str(temp_out))
    ConfMapSelection.SetParameterString("exp", expression)
    ConfMapSelection.UpdateParameters()
    #~ ConfMapSelection.ExecuteAndWriteOutput()
    ConfMapSelection.Execute()
    extraction_output = ConfMapSelection.GetImageAsNumpyArray('out')

    confidence_pixels = np.array(extraction_output)

    # remove the zero values
    confidence_pixels = confidence_pixels[confidence_pixels != 0]
    confidence_pixels = confidence_pixels[np.isnan(confidence_pixels) == False]

    # compute the mean on the valid pixels
    mean_confidence = np.mean(confidence_pixels)
    std_confidence = np.std(confidence_pixels)

    print(mean_confidence)
    print(std_confidence)
    print(len(confidence_pixels))

    return mean_confidence, std_confidence, len(confidence_pixels)


def compute_all_confidence_stats(global_parameters, out_json=''):
    modes = ['all', 'all_classified_samples', 'well_classified_samples', 'misclassified_samples']
    samples_sets = ['both', 'train', 'validation']

    data = defaultdict(dict)
    innerdict = defaultdict(dict)
    for mode in modes:
        for samples_set in samples_sets:
            mean_confidence, std_confidence, nb_pixels = confidence_map_mean(
                global_parameters, mode, samples_set, extended=True)
            innerdict["mean"] = mean_confidence
            innerdict["std"] = std_confidence
            innerdict["nb_pixels"] = nb_pixels
            data[mode][samples_set] = defaultdict(dict)
            data[mode][samples_set] = dict(innerdict)

    # Save our changes to JSON file
    if out_json != '':
        json_path = out_json
    else:
        main_dir = global_parameters["user_choices"]["main_dir"]
        json_path = op.join(main_dir, 'Statistics', 'confidence_stats.json')

    jsonFile = open(json_path, "w+")
    jsonFile.write(json.dumps(data, indent=3, sort_keys=True))
    jsonFile.close()


def plot_confidence_evolution(global_parameters):
    '''
    Plot the evolution of the mean confidence for a scene
    '''
    main_dir = global_parameters["user_choices"]["main_dir"]
    previous_ite_dir = op.join(main_dir, 'Previous_iterations')
    save_dirs = glob.glob(op.join(previous_ite_dir, 'SAVE_*'))

    if len(save_dirs) > 0:
        all_pixels_mean = []
        all_samples_mean = []
        wellclassified_samples_mean = []
        misclassified_samples_mean = []

        for save_dir in save_dirs:
            json_path = op.join(save_dir, 'Statistics', 'confidence_stats.json')
            json_data = json.load(open(json_path))
            all_pixels_mean.append(json_data["all"]["both"]["mean"])
            all_samples_mean.append(json_data["all_classified_samples"]["both"]["mean"])
            misclassified_samples_mean.append(json_data["misclassified_samples"]["both"]["mean"])
            wellclassified_samples_mean.append(json_data["well_classified_samples"]["both"]["mean"])

        save_nb = range(len(save_dirs))
        plt.plot(save_nb, all_pixels_mean, color='g',
                 linestyle='-', marker='o', label='Of all pixels')
        plt.plot(save_nb, all_samples_mean, color='b',
                 linestyle='-', marker='o', label='Of the samples')

        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
        location = global_parameters["user_choices"]["location"]
        date = global_parameters["user_choices"]["current_date"]
        plt.title('Evolution of the mean confidence\n{}, {}'.format(location, date))
        plt.xlabel('Iteration')
        plt.ylabel('Mean confidence')

        out_fig = op.join(main_dir, 'Statistics', 'confidence_evolution.png')
        plt.savefig(out_fig, bbox_inches='tight')
        plt.close()


def plot_samples_evolution(global_parameters):
    '''
    Plot the evolution of the samples number for a scene
    '''
    main_dir = global_parameters["user_choices"]["main_dir"]
    previous_ite_dir = op.join(main_dir, 'Previous_iterations')
    save_dirs = glob.glob(op.join(previous_ite_dir, 'SAVE_*'))

    if len(save_dirs) > 0:
        all_samples_nb = []
        wellclassified_samples_nb = []
        misclassified_samples_nb = []

        for save_dir in save_dirs:
            json_path = op.join(save_dir, 'Statistics', 'confidence_stats.json')
            json_data = json.load(open(json_path))
            all_samples_nb.append(json_data["all_classified_samples"]["both"]["nb_pixels"])
            misclassified_samples_nb.append(json_data["misclassified_samples"]["both"]["nb_pixels"])
            wellclassified_samples_nb.append(
                json_data["well_classified_samples"]["both"]["nb_pixels"])

        nb_of_saves = len(save_dirs)
        wellclassified_samples_ratio = [
            float(wellclassified_samples_nb[k])/all_samples_nb[k] for k in range(nb_of_saves)]
        misclassified_samples_ratio = [
            float(misclassified_samples_nb[k])/all_samples_nb[k] for k in range(nb_of_saves)]
        all_samples_ratio = [float(all_samples_nb[k])/all_samples_nb[-1]
                             for k in range(nb_of_saves)]

        save_nb = range(nb_of_saves)
        plt.plot(save_nb, all_samples_ratio, color='b', linestyle='-', marker='o',
                 label='Proportion of samples compared\nto the last iteration')
        plt.plot(save_nb, wellclassified_samples_ratio, color='g', linestyle='-', marker='o',
                 label='Proportion of well-classified\nsamples at the i$^{th}$ iteration')
        #~ plt.plot(save_nb, misclassified_samples_ratio, color='r', linestyle='--', marker='o', label = 'Of the misclassified samples')

        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
        #~ plt.legend()
        location = global_parameters["user_choices"]["location"]
        date = global_parameters["user_choices"]["current_date"]
        plt.title('Evolution of samples\n{}, {}'.format(location, date))
        plt.xlabel('Iteration')
        plt.ylabel('Samples proportion')

        #~ plt.show()
        out_fig = op.join(main_dir, 'Statistics', 'samples_evolution.png')
        plt.savefig(out_fig, bbox_inches='tight')
        plt.close()


def main():

    global_parameters = json.load(open(op.join('parameters_files', 'global_parameters.json')))
    plot_confidence_evolution(global_parameters)
    plot_samples_evolution(global_parameters)

    return

    in_tif = '/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Arles_31TFJ_20170917/Out/confidence.tif'
    out_tif = '/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Arles_31TFJ_20170917/Out/confidence_modified_med.tif'
    confidence_map_change(in_tif, out_tif, median_radius=9)


if __name__ == '__main__':
    main()