Code source de Seath

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This file is part of PHYMOBAT 2.0.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# 
# PHYMOBAT 2.0 is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# PHYMOBAT 2.0 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 General Public License
# along with PHYMOBAT 2.0.  If not, see <http://www.gnu.org/licenses/>.

import numpy, math, sys
import numpy as np # sort, ...
from collections import * # defaultdict

[docs]class Seath(): """ Get the optimal threshold and Bhattacharyya distance for a separability between two classes Source article : SEaTH–A new tool for automated feature extraction in the context of object-based image analysis S. Nussbaum et al. Source Info : Kenji Ose (IRSTEA) et Nathalie St-Geours (IRSTEA) :param value_1: List of index mean by polygons (sample 1) :type value_1: list :param value_2: List of index mean by polygons (sample 2) :type value_2: list :param threshold: Optimal threshold under this form *>0.56* :type threshold: str :param J: Jeffries-Matusita distance (measure separability between 2 classes on a 0 to 2 scale) :type J: list :Example: >>> import Seath >>> a = Seath() >>> a.value_1 = b[0].stats_dict >>> a.value_2 = b[1].stats_dict >>> a.separability_and_threshold() >>> a.threshold[0] '>0.56' >>> a.J[0] 1.86523428 """ def __init__(self): """Create a new 'Seath' instance """ self.value_1 = [] self.value_2 = [] self.threshold = [] self.J = []
[docs] def separability_and_threshold(self, **kwargs): """ Function to extract the optimal threshold for a separability between two classes :kwargs: **index** (str) - The processing will prints the string """ ind = kwargs['index'] if kwargs.get('index') else 'the index' field_class = [self.value_1, self.value_2] indict = defaultdict(list) for i in range(len(field_class)): for f in field_class[i].items(): if math.isnan(f[1][0]) == False: # if math.isnan(Terrain[i][f][0]) == False:# indict[i].append(f[1][0]) # Indict[Nommark[i]].append(Terrain[i][f][0])# ### Compute Bhattacharyya distance ### ############################### # Compute mean and variance m = defaultdict(list) # Average v = defaultdict(list) # Variance p = defaultdict(list) # Likelihood C = defaultdict(list) # Transpose classes B = [] # Bhattacharyya distance # Optimal threshold seuil1 = [] seuil2 = [] seuil = [] for mark in range(len(field_class)): C[mark] = np.array(indict[mark]).transpose() m[mark].append(np.mean(C[mark])) v[mark].append(np.var(C[mark])) p[mark].append(1 / float(len(C[mark]))) print "m : ", m print "v : ", v # Mean, standard deviation and likelihood initialisation phase for 2 classes m1 = m[0] m2 = m[1] v1 = v[0] v2 = v[1] p1 = p[0] p2 = p[1] for i in range(len(m[0])): B.append(( (1/float(8)) * ( (m1[i] - m2[i])**2 ) * (2 / ( v1[i] + v2[i] )) ) + ( (1/float(2)) * np.log( ( v1[i] + v2[i] ) / ( 2 * np.sqrt(v1[i] * v2[i]) ) ) )) self.J.append(2 * ( 1 - np.exp( -B[i] ) )) ### Optimal threshold calculation ### ###################### # Bayes theorem solution A = np.log( np.sqrt( v1[i] / v2[i] ) * ( p2[i] / p1[i] )) D = np.sqrt( v1[i] * v2[i] ) * np.sqrt( ( m1[i] - m2[i] )**2 + 2 * A * ( v1[i] - v2[i] ) ) seuil1.append(( 1 / ( v1[i] - v2[i] ) ) * ( m2[i] * v1[i] - m1[i] * v2[i] + D )) seuil2.append(( 1 / ( v1[i] - v2[i] ) ) * ( m2[i] * v1[i] - m1[i] * v2[i] - D )) # Optimal threshold # Logical condition depending on article figure 2 if ( seuil1[i] > m2[i] and seuil1[i] < m1[i] ) or ( seuil1[i] > m1[i] and seuil1[i] < m2[i] ) : print "Valid threshold !" else: seuil1[i] = "" if ( seuil2[i] > m2[i] and seuil2[i] < m1[i] ) or ( seuil2[i] > m1[i] and seuil2[i] < m2[i] ) : print "Valid threshold !" else: seuil2[i] = "" # Final condition if ( seuil1[i] == "" and seuil2[i] == "" ) or ( seuil1[i] != "" and seuil2[i] != "" ): seuil.append("") elif ( seuil1[i] != "" and seuil2[i] == "" ): seuil.append(seuil1[i]) elif ( seuil1[i] == "" and seuil2[i] != "" ): seuil.append(seuil2[i]) print("Bhattacharyya distance ", B) print("J : ", self.J) print("Threshold 1 : ", seuil1) print("Threshold 2 : ", seuil2) print("Optimal threshold :", seuil) for i in range(len(seuil)): if seuil[i] != "" and m1[i] > m2[i]: print('For ' + ind + ', the class 1 > ' + str(seuil[i])) self.threshold.append('>' + str(seuil[i])) elif seuil[i] != "" and m1[i] < m2[i]: print('For ' + ind + ', the class 1 < ' + str(seuil[i])) self.threshold.append('<' + str(seuil[i])) else: print('For ' + ind + ', not discrimination !') sys.exit(1)
# self.threshold.append('')