#first we have a piece which imports all we need but only once
# so that when we run the whole file again it does not waste time
#reimporting. 
l=0
if 1==1:
 from scipy.ndimage import filters
 from scipy import ndimage as ndimage

 import matplotlib.pyplot as plt

 from PIL import Image
 from PIL import ImageFilter

 from scipy.ndimage import measurements,morphology

 import numpy as np 

 

#######################Next we make a list with our picture-dresses
L1=["00006.jpg","00026.jpg","00071.jpg","00117.jpg","00145.jpg","00183.jpg"]
L2=["00007.jpg","00027.jpg","00082.jpg","00126.jpg","00189.jpg"]
L3=["00008.jpg","00044.jpg","00085.jpg","00129.jpg","00195.jpg"]
L4=["00012.jpg","00054.jpg","00087.jpg","00132.jpg","00197.jpg"]
L5=["00015.jpg","00055.jpg","00095.jpg","00174.jpg","00198.jpg"]
List_cars=L1+L2+L3+L4+L5

#######################next we define a function which given as input 
#a binary, gives as output a binary picture
#with the biggest connected omponent

def biggest_connected_component(a_binary_image):

 labels,nbr_objects=measurements.label(a_binary_image) #image where each connected component has a different color
 list_of_labels=list(set(labels.flatten()))
 ONE=np.zeros(labels.shape)+1
 size_labels=ndimage.sum(ONE,labels, index=list_of_labels)
 size_labels[0]=0

 maximum_size_index=np.argmax(size_labels)
 max_size_object=(labels==(list_of_labels[maximum_size_index]))
 return(max_size_object)



#######################next we define a function which given as input 
#a binary, gives as output a binary picture
#with the biggest five connected omponent

def five_biggest_connected_components(a_binary_image,number=5):
 n,m=a_binary_image.shape
 labels,nbr_objects=measurements.label(a_binary_image) #image where each connected component has a different color
 list_of_labels=list(set(labels.flatten()))
 ONE=np.zeros(labels.shape)+1
 size_labels=ndimage.sum(ONE,labels, index=list_of_labels)
 size_labels=list(size_labels)
 size_labels=[int(round(a)) for a in size_labels]
 size_labels[0]=0
 
 order_size_index=np.argsort(size_labels,axis=0)
 order_size_index=np.array(order_size_index)
 list_of_labels=np.array(list_of_labels)
 
 order_size_index=np.flip(order_size_index)
 
 reduced_list_labels=list_of_labels[order_size_index[:number]]
 
 
 
 labels2=np.copy(labels)
 for i in range(n):

   for j in range(m):
       if not(labels[i,j] in reduced_list_labels):
           labels2[i,j]=0
 return(labels2)




#################################Next we define areas of little variation,
#that is areas where image is quite constant

def areas_of_low_variation(limage,size_window=10,pourcentage=30):
  n=limage.shape[0]
  m=limage.shape[1]
  y=np.copy(limage)
   #size of the window in which we want small variation of color pixel
  # epsilon is going to be the condition that we request the points in the window to have the 3 dimensional vector of color channels to not be further
 # from center or we replace center

  B=np.copy(y[:,:,0])
  counter=0
  for i in range(size_window,n-size_window):

   for j in range(size_window,m-size_window):
     
    A=np.copy(limage[(i-size_window+1):(i+size_window),(j-size_window+1):(j+size_window),:])
    
    
    v=np.array([0,0,0])
    v[0]=y[i,j,0]
    v[1]=y[i,j,1]
    v[2]=y[i,j,2]
     #A=[(a-y[i,j,:]) for a in A]
    A=A-v
    A=np.abs(A)
    value=np.mean(A.flatten())
    B[i,j]=value
    
  
  OUTPUT=np.zeros(B.shape)+np.max(B)
  for i in range(size_window,n-size_window):

   for j in range(size_window,m-size_window):
          OUTPUT[i,j]=B[i,j]
  quantile_small=np.percentile(B.flatten(),pourcentage)
  FINAL_OUTPUT=1*(OUTPUT<quantile_small)
  return(FINAL_OUTPUT)

#############################################################
 ###NExt we find connected components for a binary image
def connected_components(binary_image):
 labels,nbr_objects=measurements.label(binary_image)
 return(labels)


###################################################################
#the next algorithm gets the connected components and colors
#them in their average color from initial picture
def connected_components_colored(binary_image,image):
   component_image=np.copy(image)
   labels,nbr_objects=measurements.label(binary_image)
   list_labels=list(set(labels.flatten()))
   ONE=np.zeros(labels.shape)+1
   size_labels=ndimage.sum(ONE,labels, index=list_labels)
   for color in [0,1,2]:
    av_color_labels=ndimage.sum(image[:,:,color],labels, index=list_labels)
    av_color_labels=np.round(av_color_labels/size_labels)

   
    av_color_labels[0]=0
   
    for k in range(len(list_labels)):
         label=list_labels[k]
         
         component_image[:,:,color][binary_image==label]=int(av_color_labels[k])
   
   
   return(component_image)

################################################################
#the next algorithm finds the connected component which is closest to the 
#center in terms of center of gravity

def component_most_centered(binary_image):
 nb_of_components=5
 labels,nbr_objects=measurements.label(binary_image)
 list_of_labels=list(set(labels.flatten()))
 ONE=np.zeros(labels.shape)+1
 size_labels=ndimage.sum(ONE,labels, index=list_of_labels)
 size_labels[0]=1
 X_matrix=np.zeros(labels.shape)
 Y_matrix=np.zeros(labels.shape)
 for i in range(229):
  for j in range(229):
       X_matrix[i,j]=j
       Y_matrix[i,j]=i
 X_mean_labels=ndimage.sum(X_matrix,labels, index=list_of_labels)
 Y_mean_labels=ndimage.sum(Y_matrix,labels, index=list_of_labels)
 X_mean_labels=X_mean_labels/size_labels
 Y_mean_labels=Y_mean_labels/size_labels 
 OUTPUT=np.zeros(labels.shape)
 order_of_sizes=np.argsort(-np.array(size_labels))
 distance=[]
 for l in range(nb_of_components):
   
   k=order_of_sizes[l]  
   OUTPUT[int(round(Y_mean_labels[k])),int(round(X_mean_labels[k]))]=list_of_labels[k]
   OUTPUT[int(round(Y_mean_labels[k]))+1,int(round(X_mean_labels[k]))]=list_of_labels[k]
   OUTPUT[int(round(Y_mean_labels[k]))-1,int(round(X_mean_labels[k]))]=list_of_labels[k] 
   OUTPUT[int(round(Y_mean_labels[k])),1+int(round(X_mean_labels[k]))]=list_of_labels[k] 
   OUTPUT[int(round(Y_mean_labels[k])),int(round(X_mean_labels[k]))-1]=list_of_labels[k]
   DIST=np.power((X_mean_labels[k]-110),2)+np.power((Y_mean_labels[k]-110),2)
   distance.append(DIST)
 distance=np.array(distance)
 ARGMIN_dist=np.argmin(distance)
 label_of_closest_component=list_of_labels[order_of_sizes[ARGMIN_dist]] 
 Image_of_closest_component=(labels==label_of_closest_component)
 return(Image_of_closest_component)
###########################################################
###########################################################
##########################################################
###########THE next functions finds points which 
#are similar in coloring to a certain area given as input,
# in the sence that the color minus brightness
#rescaled is close to the areas average. 
def similar_coloring(binary_area,image):
  horizontal_level=100
  brightness=np.zeros((229,229))
  n,m=image[:,:,0].shape
  color_difference=np.zeros((229,229,3))
  for i in range(n):
   for j in range(m):
     brightness[i,j]=np.mean(image[i,j,:])
     color_difference[i,j,0]=image[i,j,0]-brightness[i,j]
     color_difference[i,j,1]=image[i,j,1]-brightness[i,j]
     color_difference[i,j,2]=image[i,j,2]-brightness[i,j]
   
  col_diff_intensity=np.sqrt(np.power(color_difference[:,:,0],2)+np.power(color_difference[:,:,1],2)+np.power(color_difference[:,:,2],2))
  IM_high_col_diff=(col_diff_intensity>10)
  col_diff_intensity=np.copy(col_diff_intensity)+1
  col_difference_rescaled=np.copy(color_difference)  
  col_difference_rescaled[:,:,0]=color_difference[:,:,0]/col_diff_intensity
  col_difference_rescaled[:,:,1]=color_difference[:,:,1]/col_diff_intensity
  col_difference_rescaled[:,:,2]=color_difference[:,:,2]/col_diff_intensity
  
  col_diff_mean=np.array([0.1,1.0,2.0])
  col_diff_mean[0]=np.mean((col_difference_rescaled[:,:,0][binary_area==True]).flatten())
  col_diff_mean[1]=np.mean(col_difference_rescaled[:,:,1][binary_area==True])
  col_diff_mean[2]=np.mean(col_difference_rescaled[:,:,2][binary_area==True])
  
  a=np.power(col_diff_mean[0],2)+np.power(col_diff_mean[1],2)+np.power(col_diff_mean[2],2)
  a=np.sqrt(a)
  col_diff_mean=col_diff_mean/a
  #print("col_diff_mean is:")
  #print(col_diff_mean)
  A=col_difference_rescaled-col_diff_mean
  distance_col_diff=np.sqrt(np.power(A[:,:,0],2)+np.power(A[:,:,1],2)+np.power(A[:,:,2],2))
  B=np.ndarray((229,229,1))
  B[:,:,0]=IM_high_col_diff
  #print(distance_col_diff[100,:])
  C=col_difference_rescaled*B
  #print("col_difference_mean: "+str(col_diff_mean))
  #print("size of binary area"+str(np.sum(binary_area==True)))
  return(np.dot(C,col_diff_mean)>0.85)


####Next filter detemines horizontal lines:
def horizontal_filter(image):
    
    vertical_radius=2
    horizontal_radius=2
    if len(image.shape)==3:
     n,m=image[:,:,0].shape
    else:
     n,m=image.shape   
    brightness=np.zeros((n,m))
    hor_filter=np.zeros((2*vertical_radius+1,2*horizontal_radius+1))
    #HORIZONTAL FILTER:
    #hor_filter[vertical_radius,:]=-1
    #hor_filter[:(vertical_radius),:]=1/(2*vertical_radius)
    #hor_filter[(vertical_radius+1):,:]=1/(2*vertical_radius)
    
    #MOVING AVERAGE FILTER:
    hor_filter=hor_filter+1/((2*vertical_radius+1)*(2*horizontal_radius+1))
    #print("hor_filter")
    #print(hor_filter)
    #print("type hor filter:")
    #print(type(hor_filter))
    #print("this is hor_filter")
    #
    #print(hor_filter)
    #print("type of hor_filter is"+str(type(hor_filter)))
    if len(image.shape)==3:
     for i in range(n):
       for j in range(m):
         brightness[i,j]=np.mean(image[i,j,:])
    else:
        brightness=np.copy(image)
    new_image=np.zeros_like(brightness)     
    for i in range(vertical_radius+1,n-vertical_radius-1):
       for j in range(horizontal_radius+1,m-horizontal_radius-1): 
           v1=i-vertical_radius
           v2=i+vertical_radius+1
           h1=j-horizontal_radius
           h2=j+horizontal_radius+1
           new_image[i,j]=np.sum(brightness[v1:v2,h1:h2]*hor_filter)
       # color_difference[i,j,0]=image[i,j,0]-brightness[i,j]
       #color_difference[i,j,1]=image[i,j,1]-brightness[i,j]
       #color_difference[i,j,2]=image[i,j,2]-brightness[i,j]
       epsilon=np.percentile(new_image.flatten(),90)
       
       new_image2=np.zeros_like(brightness)     
    for i in range(vertical_radius+1,n-vertical_radius-1):
       for j in range(horizontal_radius+1,m-horizontal_radius-1): 
           v1=i-vertical_radius
           v2=i+vertical_radius+1
           h1=j-horizontal_radius
           h2=j+horizontal_radius+1
           new_image2[i,j]=np.sum(brightness[v1:v2,h1:h2]*(-hor_filter))
       # color_difference[i,j,0]=image[i,j,0]-brightness[i,j]
       #color_difference[i,j,1]=image[i,j,1]-brightness[i,j]
       #color_difference[i,j,2]=image[i,j,2]-brightness[i,j]
       epsilon2=np.percentile(new_image.flatten(),90)
       #output_image=(((new_image2>=epsilon2)+(new_image>=epsilon))>0)
       
    return(new_image)

######################################################################
    
    
    

    
 ################################################################################
 ########### the next function given the "the body of the car" finds the object
 #inside like lights and so on. For this it finds all the connected components
 # which are "looked" inside the body of the car and sufficiently large. 
def inside_objects(body,image):   
    
    #body is binary image which is supposed to give the body of the car
    #image is the color picture of the car in numpy matrix form 
    d=5 #distance to the border
    size_of_object=30 # size of the objects we keep
    negative_body=1-body
    labels,nbr_objects=measurements.label(negative_body)
    list_of_labels=list(set(labels.flatten()))
    ONE=np.zeros(labels.shape)+1
    size_labels=ndimage.sum(ONE,labels, index=list_of_labels)
    size_labels[0]=0
    X_matrix=np.zeros(labels.shape)
    Y_matrix=np.zeros(labels.shape)
    for i in range(229):
      for j in range(229):
       X_matrix[i,j]=j
       Y_matrix[i,j]=i
    X_max_labels=ndimage.maximum(X_matrix,labels, index=list_of_labels)
    X_min_labels=ndimage.minimum(X_matrix,labels, index=list_of_labels)
    Y_max_labels=ndimage.maximum(Y_matrix,labels, index=list_of_labels)
    Y_min_labels=ndimage.minimum(Y_matrix,labels, index=list_of_labels)
    count=0
    for u in range(len(list_of_labels)):
        if size_labels[u]<size_of_object:
            labels[labels==list_of_labels[u]]=0
        elif ((X_max_labels[u]>(229-d)) or (X_min_labels[u]<d) or (Y_max_labels[u]>(229-d)) or (Y_min_labels[u]<d)):
            labels[labels==list_of_labels[u]]=0
        else:
            count=count+1
    print("number of sizeable objects inside body of car is:")
    print(count)        
    return(labels)
##############################################################################
    ########################## Next we find the lower border of a 0,1 object
    
def lower_border(image):
    
    n,m=image[:,:].shape
    y=np.zeros(image.shape)
    maximum_line=[0]*m
    for j in range(m):
        if np.sum(image[:,j])>0:
         MAX=np.max(np.arange(n)[image[:,j]==1])
         y[MAX,j]=1
         y[MAX-1,j]=1
         maximum_line[j]=MAX
        else:
         maximum_line[j]=0   
         
    maximum_line=np.array(maximum_line)    
    return(y,maximum_line)
 
###########################################################################
#next we analyse the lower border trying to find the points where weehls start
# and long segments the imput is the y cordinate of the lower border as a function
    # y(x)

def analyser_lower_border(line,BODY_OF_CAR):
    n,m=IM.shape[:2]
    dist_border=6
    Delta_y=np.copy(line)
    for j in range(len(Delta_y)-1):
     Delta_y[j]=line[j+1]-line[j]
    output_image=np.zeros(IM.shape)
    for i in range(n):
        for j in range(m):
            if BODY_OF_CAR[i,j]==1:
                output_image[i,j,:]=(30,30,30)
    list_of_decrease_points=np.argsort(Delta_y[dist_border:(m-dist_border)])[:4]
    
    for point in list_of_decrease_points:
     max_decrease_point=point+dist_border
     output_image[:,max_decrease_point]=(200,0,0)
    
    list_of_increase_points=np.flip(np.argsort(Delta_y[dist_border:(m-dist_border)]))[:4]
    for point2 in list_of_increase_points:
     max_increase_point=point2+dist_border
     output_image[:,max_increase_point]=(0,200,0) 
     
    for j in range(m):
        output_image[line[j],j,:]=(0,200,200)
        output_image[line[j],j,:]=(0,200,200)
        
    right_lower_point=np.max(list_of_decrease_points)+dist_border-5
    for i in range(5):
        for j in range(5):
         output_image[line[right_lower_point]-i,right_lower_point-j,:]=(200,0,0)
    left_lower_point=4+dist_border+np.max(list_of_increase_points[list_of_increase_points<right_lower_point])
    for i in range(5):
        for j in range(5):
         output_image[line[left_lower_point]-i,left_lower_point-j,:]=(200,0,0)
    end_right_tire=np.max(list_of_increase_points)+dist_border
    end_right_tire=1+max([end_right_tire,right_lower_point])
    max_right_tire=np.argmin(line[right_lower_point:end_right_tire])+right_lower_point
    for i in range(5):
        for j in range(5):
         output_image[line[max_right_tire]-i,max_right_tire-j,:]=(200,0,0)
    left_end_left_tire= np.min(list_of_increase_points[list_of_increase_points<left_lower_point])    
    max_left_tire=np.argmin(line[left_end_left_tire:left_lower_point])+left_end_left_tire
    for i in range(5):
        for j in range(5):
         output_image[line[max_left_tire]-i,max_left_tire-j,:]=(200,0,0)
    Delta=4*(line[left_lower_point]-line[max_left_tire])//10   
    predicted_point_y_coordinate=line[max_left_tire]-Delta
    
    for i in range(5):
        for j in range(5):
         output_image[predicted_point_y_coordinate-i,left_lower_point-j,:]=(200,0,0)
    
    return(output_image)
    

######################################################################
####################################################################
#Here is where we start MAKING DISPLAY:


#for i in range(len(List_cars)): 
for nb in range(8,20):
      IMAGE=Image.open("foto_cars/"+List_cars[nb])
      IM_black_white=np.array(IMAGE.convert("L"))
      #epsilon_dark=np.quantile(IM_black_white.flatten(),0.15)
      IM_dark_region=1*(IM_black_white/np.max(IM_black_white)<0.1)
      
      IM=np.array(IMAGE,dtype="float64")

      #next we create the figure with several images:
      fig, axs = plt.subplots(nrows=3, ncols=2)
      
      fig.set_size_inches(16, 16)
      axs = axs.flatten()
      fig.suptitle('INSIDE BODY PARTS FINDING',fontsize=50)
      fig.subplots_adjust(wspace=0.3)
      
      
      
      #next we go on with our calculation 
      C=areas_of_low_variation(np.array(IM))
      #axs[0].imshow(IM)
      #next we create all the pictures for our display
     
      IM_copy=np.copy(IMAGE)
      horizontal_level=150
      
      for j in range(229):
        IM_copy[horizontal_level,j,:]=(250,250,0)
        IM_copy[horizontal_level+1,j,:]=(250,250,0)   
      axs[0].imshow(IM_copy)
      axs[0].set_title("Picture with horizontal line")
      IM_copy=np.copy(IM)
      IM_colored_comp=connected_components_colored(connected_components(areas_of_low_variation(IM)),IM)
      #axs[1].imshow(np.array(IM_colored_comp,dtype="uint8"))
      #axs[1].set_title("Areas_of_low_variation")
      ##Next we take out intensity and then apply areas of low variation
      IM_intensity=np.zeros((229,229))
      IM_col_diff_0=np.zeros((229,229))
      IM_col_diff_1=np.zeros((229,229))
      IM_col_diff_2=np.zeros((229,229))
      for i in range(229):


        for j in range(229):
          IM_intensity[i,j]=(IM[i,j,0]+IM[i,j,1]+IM[i,j,2])/3
          
          IM_col_diff_0[i,j]=IM[i,j,0]-IM_intensity[i,j]
          IM_col_diff_1[i,j]=IM[i,j,1]-IM_intensity[i,j]
          IM_col_diff_2[i,j]=IM[i,j,2]-IM_intensity[i,j]
      
      # Next we show a horizontal line through picture and plot the picture ther
      IM_copy=np.copy(IM)

      for j in range(229):
        IM_copy[horizontal_level,j,:]=(1,1,0)
        IM_copy[horizontal_level+1,j,:]=(1,1,0)   
      
      #axs[2].imshow(biggest_connected_component(areas_of_low_variation(IM)))
      #axs[2].set_title("Biggest_components of_low_var")
      H_level_0=IM[horizontal_level,:,0]
      H_level_1=IM[horizontal_level,:,1]
      H_level_2=IM[horizontal_level,:,2]
      
      
      
      #axs[2].plot(range(229),IM_intensity[horizontal_level,:])
      #axs[2].set_title("horizontal_cut_brightness")

      CENTERED_COMP=component_most_centered(areas_of_low_variation(IM))
      SIMILAR_COLOR=similar_coloring(CENTERED_COMP,IM)
      BODY_OF_CAR=biggest_connected_component(SIMILAR_COLOR)
      
      axs[1].imshow(BODY_OF_CAR)
 
      axs[1].set_title("BODY OF CAR")

      #axs[4].plot(range(229),H_level_0,"r",range(229),H_level_1,"g",range(229),H_level_2,"b")
      #axs[4].set_title("Horizontal_cut_picture")
      
      IM_col_diff_1[horizontal_level,:]=IM[horizontal_level,:,1]-IM_intensity[horizontal_level,:]
      #axs[5].plot(range(229),IM_col_diff_0[horizontal_level,:],"r",range(229),IM_col_diff_1[horizontal_level,:],"g",range(229),IM_col_diff_2[horizontal_level,:],"b")
      #axs[5].set_title("horizontal_cut_image_minus_brightness")
      #axs[4].imshow(component_most_centered(areas_of_low_variation(IM)))
      #imageWithEdges = IMAGE.filter(ImageFilter.FIND_EDGES)
      
      IM_colored_comp_in_body=connected_components_colored(connected_components((1-BODY_OF_CAR)),IM)
      #picture_with_line,graph=black_line(IM)
      #axs[4].imshow(np.array(picture_with_line,dtype="uint8"))
      #axs[4].set_title("Color intensity through the cut")
      #axs[5].plot(range(len(graph)),graph)
      #axs[5].imshow(np.array(IM_colored_comp_in_body,dtype="uint8"))
      #
      
      axs[2].imshow(inside_objects(BODY_OF_CAR,IM))
      axs[2].set_title("Inside_objects")
      five_biggest_connected_components
      LOWER_BORDER_IMAGE,lower_border_string=lower_border(BODY_OF_CAR)
      axs[3].imshow(LOWER_BORDER_IMAGE)
      
      axs[3].set_title("lower_border") 
      contour_image=(1-areas_of_low_variation(IM,size_window=5,pourcentage=75))
      contour_image=five_biggest_connected_components(horizontal_filter(contour_image))
      BODY_CONV_FILTER=horizontal_filter((1.0*BODY_OF_CAR))
      CONTOUR=(1*(BODY_CONV_FILTER<1))*(1*(BODY_CONV_FILTER>0))
      axs[4].imshow(analyser_lower_border(lower_border_string,BODY_OF_CAR))
      
      axs[4].set_title("lower_border_analysed") 
      
      axs[5].imshow(CONTOUR)
      axs[5].set_title("CONTOUR") 
      
      file_name_to_save_picture="./OUTPUT/car"+str(nb)
      #plt.savefig(file_name_to_save_picture)