<s> import ee <EOL> import math <EOL> import numpy <EOL> import scipy <EOL> import scipy . special <EOL> import scipy . optimize <EOL> import matplotlib <EOL> import matplotlib . pyplot as plt <EOL> '''<STR_LIT>''' <EOL> class RadarHistogram ( object ) : <EOL> '''<STR_LIT>''' <EOL> BACKSCATTER_MODEL_GAMMA = <NUM_LIT:1> <EOL> BACKSCATTER_MODEL_GAUSSIAN = <NUM_LIT:2> <EOL> BACKSCATTER_MODEL_DIP = <NUM_LIT:3> <EOL> BACKSCATTER_MODEL_PEAK = <NUM_LIT:4> <EOL> def __init__ ( self , domain , sensor , backscatter_model = None ) : <EOL> self . domain = domain <EOL> self . sensor = sensor <EOL> self . backscatter_model = [ ] <EOL> for i in range ( len ( sensor . band_names ) ) : <EOL> try : <EOL> if backscatter_model != None : <EOL> self . backscatter_model . append ( backscatter_model ) <EOL> else : <EOL> if sensor . water_distributions [ sensor . band_names [ i ] ] [ '<STR_LIT>' ] == '<STR_LIT>' : <EOL> self . backscatter_model . append ( RadarHistogram . BACKSCATTER_MODEL_PEAK ) <EOL> elif sensor . water_distributions [ sensor . band_names [ i ] ] [ '<STR_LIT>' ] == '<STR_LIT>' : <EOL> self . backscatter_model . append ( RadarHistogram . BACKSCATTER_MODEL_DIP ) <EOL> elif sensor . water_distributions [ sensor . band_names [ i ] ] [ '<STR_LIT>' ] == '<STR_LIT>' : <EOL> self . backscatter_model . append ( RadarHistogram . BACKSCATTER_MODEL_GAMMA ) <EOL> else : <EOL> self . backscatter_model . append ( RadarHistogram . BACKSCATTER_MODEL_GAUSSIAN ) <EOL> except KeyError : <EOL> pass <EOL> self . hist_image = self . __preprocess_image ( sensor ) <EOL> self . histograms = self . __compute_histogram ( self . hist_image ) <EOL> self . __find_thresholds ( ) <EOL> def __preprocess_image ( self , sensor ) : <EOL> image = sensor . image <EOL> if self . sensor . log_scale : <EOL> image = image . log10 ( ) <EOL> for b in sensor . band_names : <EOL> if '<STR_LIT>' in sensor . water_distributions [ b ] : <EOL> other_bands = list ( sensor . band_names ) <EOL> other_bands . remove ( b ) <EOL> r = sensor . water_distributions [ b ] [ '<STR_LIT>' ] <EOL> image = image . select ( [ b ] , [ b ] ) . clamp ( r [ <NUM_LIT:0> ] , r [ <NUM_LIT:1> ] ) . addBands ( image . select ( other_bands , other_bands ) ) <EOL> return image <EOL> def __compute_histogram ( self , image ) : <EOL> buckets = <NUM_LIT> if '<STR_LIT>' not in self . sensor . water_distributions [ self . sensor . band_names [ <NUM_LIT:0> ] ] else self . sensor . water_distributions [ self . sensor . band_names [ <NUM_LIT:0> ] ] [ '<STR_LIT>' ] <EOL> histogram = image . reduceRegion ( ee . Reducer . histogram ( buckets , None , None ) , self . domain . bounds , <NUM_LIT:30> , None , None , True ) . getInfo ( ) <EOL> h = [ ] <EOL> for c in range ( len ( self . sensor . band_names ) ) : <EOL> ch = self . sensor . band_names [ c ] <EOL> total = sum ( histogram [ ch ] [ '<STR_LIT>' ] ) <EOL> histogram [ ch ] [ '<STR_LIT>' ] = map ( lambda x : x / total , histogram [ ch ] [ '<STR_LIT>' ] ) <EOL> h . append ( ( histogram [ ch ] [ '<STR_LIT>' ] , histogram [ ch ] [ '<STR_LIT>' ] , histogram [ ch ] [ '<STR_LIT>' ] ) ) <EOL> return h <EOL> def __cdf ( self , params , x , backscatter_model ) : <EOL> mode = params [ <NUM_LIT:0> ] <EOL> k = params [ <NUM_LIT:1> ] <EOL> offset = params [ <NUM_LIT:2> ] <EOL> if backscatter_model == RadarHistogram . BACKSCATTER_MODEL_GAUSSIAN : <EOL> return <NUM_LIT:0.5> * ( <NUM_LIT:1> + scipy . special . erf ( ( x - mode ) / ( k * math . sqrt ( <NUM_LIT:2> ) ) ) ) <EOL> theta = ( mode - offset ) / ( k - <NUM_LIT:1> ) <EOL> return scipy . special . gammainc ( k , ( x - offset ) / theta ) <EOL> def __cdf_percentile ( self , params , percentile , backscatter_model ) : <EOL> mode = params [ <NUM_LIT:0> ] <EOL> k = params [ <NUM_LIT:1> ] <EOL> offset = params [ <NUM_LIT:2> ] <EOL> if backscatter_model == RadarHistogram . BACKSCATTER_MODEL_GAUSSIAN : <EOL> return scipy . special . erfinv ( percentile / <NUM_LIT:0.5> - <NUM_LIT:1> ) * k * math . sqrt ( <NUM_LIT:2> ) + mode <EOL> theta = ( mode - offset ) / ( k - <NUM_LIT:1> ) <EOL> v = scipy . special . gammaincinv ( k , percentile ) * theta + offset <EOL> return v <EOL> def __show_histogram ( self , channel ) : <EOL> start = self . histograms [ channel ] [ <NUM_LIT:0> ] <EOL> width = self . histograms [ channel ] [ <NUM_LIT:1> ] <EOL> values = self . histograms [ channel ] [ <NUM_LIT:2> ] <EOL> ind = numpy . arange ( start = start , stop = start + width * len ( values ) , step = width ) [ : - <NUM_LIT:1> ] <EOL> plt . bar ( ind , height = values [ : - <NUM_LIT:1> ] , width = width , color = '<STR_LIT:b>' ) <EOL> plt . ylabel ( '<STR_LIT>' ) <EOL> plt . xlabel ( self . sensor . band_names [ channel ] + '<STR_LIT>' ) <EOL> bm = self . backscatter_model [ channel ] <EOL> ( threshold , params ) = self . __find_threshold_histogram ( channel ) <EOL> if params != None : <EOL> mid = int ( ( params [ <NUM_LIT:0> ] - start ) / width ) <EOL> cumulative = sum ( values [ : mid ] ) + values [ mid ] / <NUM_LIT:2> <EOL> scale = cumulative / self . __cdf ( params , params [ <NUM_LIT:0> ] , bm ) <EOL> plt . bar ( ind , map ( lambda x : scale * ( self . __cdf ( params , x + width , bm ) - self . __cdf ( params , x , bm ) ) , ind ) , width = width , color = '<STR_LIT:r>' , alpha = <NUM_LIT:0.5> ) <EOL> plt . plot ( ( threshold , threshold ) , ( <NUM_LIT:0> , <NUM_LIT> ) , '<STR_LIT>' ) <EOL> def __gamma_function_errors ( self , p , mode , fit_end , offset , channel ) : <EOL> start = self . histograms [ channel ] [ <NUM_LIT:0> ] <EOL> width = self . histograms [ channel ] [ <NUM_LIT:1> ] <EOL> values = self . histograms [ channel ] [ <NUM_LIT:2> ] <EOL> k = p [ <NUM_LIT:0> ] <EOL> bm = self . backscatter_model [ channel ] <EOL> if ( bm == RadarHistogram . BACKSCATTER_MODEL_GAMMA and k <= <NUM_LIT:1.0> ) : <EOL> return [ float ( '<STR_LIT>' ) ] * len ( values ) <EOL> if ( bm == RadarHistogram . BACKSCATTER_MODEL_GAUSSIAN and k <= <NUM_LIT:0.0> ) : <EOL> return [ float ( '<STR_LIT>' ) ] * len ( values ) <EOL> error = <NUM_LIT:0.0> <EOL> last_cdf = <NUM_LIT:0.0> <EOL> errors = numpy . zeros ( len ( values ) ) <EOL> mid = int ( ( mode - start ) / width ) <EOL> cumulative = sum ( values [ : mid ] ) + values [ mid ] / <NUM_LIT:2> <EOL> scale = cumulative / self . __cdf ( ( mode , k , offset ) , mode , bm ) <EOL> for i in range ( len ( values ) ) : <EOL> if start + i * width - offset > fit_end : <EOL> break <EOL> cdf = scale * self . __cdf ( ( mode , k , offset ) , start + i * width , bm ) <EOL> errors [ i ] = ( cdf - last_cdf ) - values [ i ] <EOL> last_cdf = cdf <EOL> return errors <EOL> def __find_threshold_histogram ( self , channel ) : <EOL> start = self . histograms [ channel ] [ <NUM_LIT:0> ] <EOL> width = self . histograms [ channel ] [ <NUM_LIT:1> ] <EOL> values = self . histograms [ channel ] [ <NUM_LIT:2> ] <EOL> temp = self . sensor . water_distributions [ self . sensor . band_names [ channel ] ] [ '<STR_LIT>' ] <EOL> ( minv , maxv ) = ( temp [ '<STR_LIT>' ] , temp [ '<STR_LIT>' ] ) <EOL> i = int ( ( minv - start ) / width ) <EOL> biggest_bin = i <EOL> while i < len ( values ) and start + i * width <= maxv : <EOL> if values [ i ] > values [ biggest_bin ] : <EOL> biggest_bin = i <EOL> i += <NUM_LIT:1> <EOL> mode = start + biggest_bin * width <EOL> if self . backscatter_model [ channel ] == RadarHistogram . BACKSCATTER_MODEL_PEAK : <EOL> return ( mode , None ) <EOL> i = biggest_bin + <NUM_LIT:1> <EOL> while i < len ( values ) - <NUM_LIT:3> and start + i * width <= maxv : <EOL> if values [ i ] < values [ i + <NUM_LIT:1> ] and values [ i ] < values [ i + <NUM_LIT:2> ] and values [ i ] < values [ i + <NUM_LIT:3> ] : <EOL> break <EOL> i += <NUM_LIT:1> <EOL> local_min = start + i * width <EOL> if self . backscatter_model [ channel ] == RadarHistogram . BACKSCATTER_MODEL_DIP : <EOL> return ( local_min , None ) <EOL> m = self . sensor . minimum_value <EOL> if self . sensor . log_scale : <EOL> m = math . log10 ( m ) <EOL> ( value , result ) = scipy . optimize . leastsq ( self . __gamma_function_errors , [ <NUM_LIT:10> ] , factor = <NUM_LIT:1.0> , args = ( mode , mode , m , channel ) ) <EOL> params = ( mode , value [ <NUM_LIT:0> ] , m ) <EOL> mid = int ( ( mode - start ) / width ) <EOL> cumulative = sum ( values [ : mid ] ) + values [ mid ] / <NUM_LIT:2> <EOL> scale = cumulative / self . __cdf ( params , params [ <NUM_LIT:0> ] , self . backscatter_model [ channel ] ) <EOL> i = mid <EOL> while i < len ( values ) : <EOL> cumulative += values [ i ] <EOL> diff = cumulative - scale * self . __cdf ( params , start + i * width , self . backscatter_model [ channel ] ) <EOL> if diff > <NUM_LIT> : <EOL> break <EOL> i += <NUM_LIT:1> <EOL> threshold = start + i * width <EOL> return ( threshold , params ) <EOL> def __find_thresholds ( self ) : <EOL> self . thresholds = [ ] <EOL> self . distributions = [ ] <EOL> for c in range ( len ( self . sensor . band_names ) ) : <EOL> ( threshold , params ) = self . __find_threshold_histogram ( c ) <EOL> if self . sensor . log_scale : <EOL> threshold = <NUM_LIT:10> ** threshold <EOL> self . thresholds . append ( threshold ) <EOL> self . distributions . append ( params ) <EOL> def get_thresholds ( self ) : <EOL> return self . thresholds <EOL> def find_loose_thresholds ( self , percentile = <NUM_LIT> ) : <EOL> results = [ ] <EOL> for c in range ( len ( self . sensor . band_names ) ) : <EOL> ( threshold , params ) = ( self . thresholds [ c ] , self . distributions [ c ] ) <EOL> if params != None : <EOL> t = self . __cdf_percentile ( params , percentile , self . backscatter_model [ c ] ) <EOL> if self . sensor . log_scale : <EOL> t = <NUM_LIT:10> ** t <EOL> results . append ( t ) <EOL> else : <EOL> if self . sensor . log_scale : <EOL> threshold = math . log10 ( threshold ) <EOL> start = self . histograms [ c ] [ <NUM_LIT:0> ] <EOL> width = self . histograms [ c ] [ <NUM_LIT:1> ] <EOL> values = self . histograms [ c ] [ <NUM_LIT:2> ] <EOL> i = int ( ( threshold - start ) / width ) + <NUM_LIT:1> <EOL> if self . backscatter_model [ c ] == RadarHistogram . BACKSCATTER_MODEL_DIP : <EOL> while i < len ( values ) - <NUM_LIT:3> : <EOL> if values [ i ] > values [ i + <NUM_LIT:1> ] and values [ i ] > values [ i + <NUM_LIT:2> ] and values [ i ] > values [ i + <NUM_LIT:3> ] : <EOL> break <EOL> i += <NUM_LIT:1> <EOL> else : <EOL> while i < len ( values ) - <NUM_LIT:3> : <EOL> if values [ i ] < values [ i + <NUM_LIT:1> ] and values [ i ] < values [ i + <NUM_LIT:2> ] and values [ i ] < values [ i + <NUM_LIT:3> ] : <EOL> break <EOL> i += <NUM_LIT:1> <EOL> t = start + i * width <EOL> if self . sensor . log_scale : <EOL> t = <NUM_LIT:10> ** t <EOL> results . append ( t ) <EOL> return results <EOL> def show_histogram ( self ) : <EOL> plt . figure ( <NUM_LIT:1> ) <EOL> for c in range ( len ( self . sensor . band_names ) ) : <EOL> ch = self . sensor . band_names [ c ] <EOL> plt . subplot ( <NUM_LIT:100> * len ( self . sensor . band_names ) + <NUM_LIT:10> + c + <NUM_LIT:1> ) <EOL> self . __show_histogram ( c ) <EOL> plt . show ( ) <EOL> def __computeJT ( relHist , binVals , T ) : <EOL> '''<STR_LIT>''' <EOL> FAIL_VAL = <NUM_LIT> <EOL> histogram1 = relHist [ <NUM_LIT:0> : T ] <EOL> histogram2 = relHist [ T + <NUM_LIT:1> : ] <EOL> P1 = sum ( histogram1 ) <EOL> P2 = sum ( histogram2 ) <EOL> if ( P1 <= <NUM_LIT:0> ) or ( P2 <= <NUM_LIT:0> ) : <EOL> return FAIL_VAL <EOL> weightedBins1 = numpy . multiply ( histogram1 , binVals [ <NUM_LIT:0> : T ] ) <EOL> weightedBins2 = numpy . multiply ( histogram2 , binVals [ T + <NUM_LIT:1> : ] ) <EOL> mean1 = sum ( weightedBins1 ) / P1 ; <EOL> mean2 = sum ( weightedBins2 ) / P2 ; <EOL> diffs1 = numpy . subtract ( binVals [ <NUM_LIT:0> : T ] , mean1 ) <EOL> diffs2 = numpy . subtract ( binVals [ T + <NUM_LIT:1> : ] , mean2 ) <EOL> diffsSq1 = [ d * d for d in diffs1 ] <EOL> diffsSq2 = [ d * d for d in diffs2 ] <EOL> weightedBins1 = numpy . multiply ( histogram1 , diffsSq1 ) <EOL> weightedBins2 = numpy . multiply ( histogram2 , diffsSq2 ) <EOL> sigma1 = math . sqrt ( sum ( weightedBins1 ) / P1 ) <EOL> sigma2 = math . sqrt ( sum ( weightedBins2 ) / P2 ) <EOL> if ( sigma1 <= <NUM_LIT:0> ) or ( sigma2 <= <NUM_LIT:0> ) : <EOL> return FAIL_VAL <EOL> J = <NUM_LIT:1> + <NUM_LIT:2> * ( P1 * math . log ( sigma1 ) + P2 * math . log ( sigma2 ) ) - <NUM_LIT:2> * ( P1 * math . log ( P1 ) + P2 * math . log ( P2 ) ) <EOL> return J <EOL> def splitHistogramKittlerIllingworth ( histogram , binVals ) : <EOL> '''<STR_LIT>''' <EOL> histSum = float ( sum ( histogram ) ) <EOL> relativeHist = numpy . divide ( histogram , histSum ) <EOL> numBins = len ( binVals ) <EOL> J = [ ] <EOL> for i in range ( <NUM_LIT:0> , numBins ) : <EOL> J . append ( __computeJT ( relativeHist , binVals , i ) ) <EOL> minJ = J [ <NUM_LIT:1> ] <EOL> threshold = binVals [ <NUM_LIT:0> ] <EOL> for i in range ( <NUM_LIT:1> , numBins ) : <EOL> if J [ i ] < minJ : <EOL> minJ = J [ i ] <EOL> threshold = ( binVals [ i ] + binVals [ i - <NUM_LIT:1> ] ) / <NUM_LIT:2> <EOL> return threshold </s>
