# -*- coding: utf-8 -*-
"""
Spyder Editor


"""
import getopt
import sys
import ogr
import gdalconst
import numpy as np
import rtree
import os
import xlsxwriter
import six
import packaging
import packaging.version
import packaging.specifiers
import packaging.requirements



def ROW(path,field_name):
    source = ogr.Open(path, 1) # 1 for read an write
    layer = source.GetLayer()
    
    
    layer_defn = layer.GetLayerDefn()
    field_names = [layer_defn.GetFieldDefn(i).GetName() for i in range(layer_defn.GetFieldCount())]
    
    l = len(field_names)
    
    field =np.zeros(l)
    
 
    for i in range(l):
        if field_names[i] == field_name:
            field[i] = 1
        else:
            field[i]=0
            
    
    
    if sum(field) == 0:
        
        new_field = ogr.FieldDefn(field_name, ogr.OFTInteger)
    
        layer.CreateField(new_field)
       
    
        i = 0
    
        for feature in layer:
            i = i +1
            feature.SetField(field_name, i)
            layer.SetFeature(feature)
         

    
    else:
       
    
        i = 0
    
        for feature in layer:
            i = i +1
            feature.SetField(field_name, i)
            layer.SetFeature(feature)
         

        
    source = None
    return




def intersect(index,path_1,path_2,field_name,p):
    
    source_1 = ogr.Open(path_1, gdalconst.GA_ReadOnly) 

    layer_1 = source_1.GetLayer()
    
   
    source_2 = ogr.Open(path_2, gdalconst.GA_ReadOnly)
    

    layer_2 = source_2.GetLayer()
    
    p = float(p)/100
    
  
    
    Tab=[[],[],[],[],[],[],[]]
    
    for j in range(0, layer_2.GetFeatureCount()):
        feature2 = layer_2.GetFeature(j)
        geometry2 = feature2.GetGeometryRef()
        xmin, xmax, ymin, ymax = geometry2.GetEnvelope()
        
        for i in list(index.intersection((xmin, xmax, ymin, ymax))):
            
            feature1 = layer_1.GetFeature(i)
            geometry1 = feature1.GetGeometryRef()
            
            
            
            
            inte = geometry2.Intersects(geometry1)
            
            if inte== True:
                
                inter= geometry2.Intersection(geometry1)
                                
                if inter.GetArea() != 0:
                    
                    area1 = geometry1.GetArea()                    
                    ch_area1 = area1*p
                    area2 = geometry2.GetArea()  
                    ch_area2 = area2*p
                    area = inter.GetArea()
                    
                    if ch_area1 < area:
                                           
                        row_gt = feature1.GetField('ROW_gt')
                        row_seg = feature2.GetField(field_name)
                        
                        C=[row_gt, row_seg, area, ch_area1, ch_area2, area1, area2]
                        Tab= np.c_[Tab,C]
                        
                    elif ch_area2 < area:
                        
                        row_gt = feature1.GetField('ROW_gt')
                        row_seg = feature2.GetField(field_name)
                        
                        C=[row_gt, row_seg, area, ch_area1, ch_area2, area1, area2]
                        Tab= np.c_[Tab,C]
                        
     
    return Tab



def ED2(tab, number_gt):

    ED2_value = np.zeros(10)
    
    row_gt = tab[0,:]
    row_gt_or= np.unique(row_gt)
    n_gt= len(row_gt_or)
    
    row_seg = tab[1,:]
    row_seg_or =np.unique(row_seg)
    
    n_seg = len(row_seg_or)
    
   
    if n_gt != number_gt:
        frequency =np.bincount(row_gt.astype(int))
        
        penality_nsr = np.max(frequency)
        n_seg = n_seg + penality_nsr*np.absolute(number_gt-n_gt)
        nsr = np.absolute(float(n_gt)-float(n_seg))/float(n_gt)
        
    else:
        nsr = np.absolute(float(n_gt)-float(n_seg))/float(n_gt)
    
    
    area_gt = 0
    
    under_area = 0
    
    
    if n_gt == number_gt:
        for i in range(0,n_gt):
            
            gt = row_gt_or[i]
            
            A = row_gt == gt
            
            col = np.column_stack(np.where(A))
            
            l = len(col)
            
            area_gt = area_gt + tab[5,col[0]]
            
            for j in range(0,l):
                
                under_area = under_area + np.absolute(tab[6,col[j]] - tab[2,col[j]])
                
        
        pse = under_area/area_gt
        
    else:
        
        under_area = []
        
        
        for i in range(0,n_gt):
            
            under_area_gt=0
            
            gt = row_gt_or[i]
            
            A = row_gt == gt
            
            col = np.column_stack(np.where(A))
            
            l = len(col)
            
            area_gt = area_gt + tab[5,col[0]]
            
            for j in range(0,l):
                
                under_area_gt = under_area_gt + np.absolute(tab[6,col[j]] - tab[2,col[j]])
             
            under_area =np.append(under_area,under_area_gt)
        
    
        
        pse_penality = np.max(under_area)
        
        under_area = np.sum(under_area)

        pse = (under_area + pse_penality*np.absolute(number_gt-n_gt))/area_gt
        
    
    ED2_index = np.sqrt(np.power(nsr,2) + np.power(pse,2))
    
    ED2_value[3] = n_gt
    ED2_value[4] = n_seg
    ED2_value[5] = area_gt
    ED2_value[6] = under_area
    ED2_value[7] = nsr
    ED2_value[8] = pse
    ED2_value[9] = ED2_index
    
    
    return ED2_value




def index_gt(path):
    
    source_1 = ogr.Open(path, gdalconst.GA_ReadOnly) 

    layer_1 = source_1.GetLayer()
    
    number_gt = layer_1.GetFeatureCount()
    
    index = rtree.index.Index(interleaved=False)

    for i in range(0,layer_1.GetFeatureCount()):
        feature1 = layer_1.GetFeature(i)
        geometry1 = feature1.GetGeometryRef()
        xmin, xmax, ymin, ymax = geometry1.GetEnvelope()
        index.insert(i, (xmin, xmax, ymin, ymax))
    
    return (index, number_gt)    



def path_segmentation(path_2):
    
    input_seg = os.listdir(path_2)
    
    seg=[]
    seg_fold = []

    path =[]
    
    for i in input_seg:
        

    
        path.append(os.path.join(path_2,i))
   
        
        
    for i in path:
        
        files = os.listdir(i)
        seg_fold =[]       
        for j in files:
            if j.endswith('shp'):
                seg_fold.append(j)
        seg.append(seg_fold)
            

            
    return (seg,path,input_seg)


###############################################################################
###############################################################################
######################## main #################################################
###############################################################################
###############################################################################


def main(argv):
   working_folder = ''
   outputfile = ''
   p = 0
   try:
      opts, args = getopt.getopt(argv,"hw:o:p:",["w_path=","ofile=","percetage"])
      
     
      if len(opts) != 3:
          print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
          sys.exit(2)          
      
   except getopt.GetoptError:
      print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
      sys.exit(2)
   for opt, arg in opts:
      if opt == '-h':
         print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
         sys.exit()
      elif opt in ("-w", "--w_path"):
         working_folder = arg
         if working_folder == '':
             print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
             sys.exit()
              
      elif opt in ("-o", "--ofile"):
         outputfile = arg
         if outputfile == '':
             print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
             sys.exit()         
         
      elif opt in ("-p", "--percentage"):
         p = arg
         if p == '':
             print( 'AssesSeg -w <working folder> -o <outputfile> -p <overlapping percentage>')
             sys.exit()           
          
         
   return(working_folder, p, outputfile)     
   






(working_folder, p, outputfile) = main(sys.argv[1:])  



if outputfile.endswith('.xlsx'):
	outputfile = outputfile
	
elif outputfile.endswith('.xls'):
	outputfile = outputfile + 'x'
	
else: 
	outputfile = outputfile +'.xlsx'

os.chdir(working_folder)    



os.chdir(working_folder)

	
path_2 = './seg_test'
name_2= 'ROW_seg'
path_1= './gt'
name_1 = 'ROW_gt'


input_gt = os.listdir(path_1)
gt = []


for i in input_gt:
	if i.endswith('shp'):
		gt.append(i)


path_gt = os.path.join(path_1,gt[0])

		
ROW(path_gt,name_1)



(index, number_gt) =  index_gt(path_gt)



 
seg, path , input_seg= path_segmentation(path_2)       

l_path = len(path)




path_gt = os.path.join(path_1,gt[0])

   
title = ('name','scale','shape','compactness','n. gt geometry',
		 'n. seg geometry','area gt','under seg area',
		  'nsr','pse','ED2_index')



col =('B2','C2','D2','E2','F2','G2','H2','I2','J2','K2')



workbook = xlsxwriter.Workbook(outputfile)



for i in range(0,l_path):
	
	shp_names = []
	excel= [[],[],[],[],[],[],[],[],[],[]]
	worksheet = workbook.add_worksheet(input_seg[i])
	worksheet.write_row('A1', title)

		
	for segments in seg[i]:
		
	   try: 
		
			if len(segments)==24:
				
			
				scale = float(segments[3:5])
				shape = float(segments[9:12])
				compactness = float(segments[17:20])
				
				path_seg = os.path.join(path[i],segments)
				
				print path_seg, ' processing' 
				
				ROW(path_seg,name_2)
				
				
				
				tab = intersect(index,path_gt,path_seg,name_2,p)
				
				ED2_value =ED2(tab,number_gt)
				
				ED2_value[0] = scale
				ED2_value[1] = shape
				ED2_value[2] = compactness
				
				excel= np.c_[excel,ED2_value]
				
				shp_names.append(segments)
					
		   
				
			elif len(segments)==25:
				
				
				scale = float(segments[3:6])
				shape = float(segments[10:13])
				compactness = float(segments[18:21])
				
				path_seg = os.path.join(path[i],segments)
				
				print path_seg, ' processing' 
				
				ROW(path_seg,name_2)
				
				
				
				tab= intersect(index,path_gt,path_seg,name_2,p)
				
				ED2_value =ED2(tab, number_gt)
				
				ED2_value[0] = scale
				ED2_value[1] = shape
				ED2_value[2] = compactness
				
				excel= np.c_[excel,ED2_value]
				shp_names.append(segments)
				
		
			else:        
							   
				scale = 0
				shape = 0
				compactness = 0
				
				path_seg = os.path.join(path[i],segments)
				
				print path_seg, ' processing'
				
				ROW(path_seg,name_2)
				   
				
				tab= intersect(index,path_gt,path_seg,name_2,p)
				
				ED2_value =ED2(tab, number_gt)
				
				ED2_value[0] = scale
				ED2_value[1] = shape
				ED2_value[2] = compactness
				
				excel= np.c_[excel,ED2_value]
				shp_names.append(segments)
			
	   except ValueError:
			
			scale = 0
			shape = 0
			compactness = 0
			
			path_seg = os.path.join(path[i],segments)
			
			print path_seg, ' processing'
			
			ROW(path_seg,name_2)
			   
			
			tab= intersect(index,path_gt,path_seg,name_2,p)
			
			ED2_value =ED2(tab, number_gt)
			
			ED2_value[0] = scale
			ED2_value[1] = shape
			ED2_value[2] = compactness
			
			excel= np.c_[excel,ED2_value]
			shp_names.append(segments)           
		   
		   
		   

		   
	worksheet.write_column('A2',shp_names)    
			
	for j in range(0,10):
		
		worksheet.write_column(col[j], excel[j,:])
		
workbook.close()       
print('Done!!!')
