#!/usr/bin/python3
# tutorial: https://medium.com/@jamesthesken/crop-row-detection-using-python-and-opencv-93e456ddd974


import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('/home/caroline/Pictures/trees1.png')
#cv2.imshow("img",img)
b,g,r = cv2.split(img) # split image in blue, green and red channels
#get grayscale image by 'amplifying' the color green and getting rid of red & blue
gscale = 2*g-r-b

#Canny edge detection: 
#1. noise reduction with 5x5 Gaussian filter 
#2. Finding Intensity Gradient: filtering with sobel kernel in horizintal 
#   & vertical direction to get 1st derivative in horizontal & vertical direction
#   calculation of edge gradient & direction for each pixel
#   => getting gradient, magnitude & direction
#3. Non-maximum Suppression: scan to remove unwanted pixels (not constituting edge)
#   checking if pixel is local maximum in its neighbourhood in direction of gradient
#   => binary image with hin edges
#4. Hysteresis Thresholding: threshold values maxVal & minVal to determine edges & non-edges (discarded)
#   edges between values are only edges if connected to to "sure-edges"
#   also removes small pixels as edges are assumed to be long lines 
gscale = cv2.Canny(gscale,280,290,apertureSize = 3) #(input image, minVal, maxVal, aperture_size=size of Sobel kernel to find image gradients: default=3)
# default L2gradient=False => Edge_Gradient(G)=|G_x|+|G_y|; True => Edge_Gradient(G)=sqrt(G_x²+G_y²)
gscale2 = cv2.Canny(gscale,100,290,apertureSize = 3)
gscale3 = cv2.Canny(gscale,250,290,apertureSize = 3)
gscale4 = cv2.Canny(gscale,280,350,apertureSize = 3)

""" plt.figure()
plt.plot(), plt.imshow(gscale)
plt.title('Canny Edge-Detection Results')
plt.xticks([]), plt.yticks([])
plt.show()  """

##
plt.subplot(221)
plt.plot(),plt.imshow(gscale) 
plt.title('Canny Edge-Detection 1'), plt.xticks([]), plt.yticks([])
plt.subplot(222)
plt.plot(),plt.imshow(gscale2)
plt.title('Canny Edge-Detection 2'), plt.xticks([]), plt.yticks([])
plt.subplot(223)
plt.plot(),plt.imshow(gscale3)
plt.title('Canny Edge-Detection 3'), plt.xticks([]), plt.yticks([])
plt.subplot(224)
plt.plot(),plt.imshow(gscale4)
plt.title('Canny Edge-Detection 4'), plt.xticks([]), plt.yticks([])
plt.show()
##


size = np.size(gscale) #returns the product of the array dimensions
skel = np.zeros(gscale.shape,np.uint8) #array of zeros
ret,gscale = cv2.threshold(gscale,128,255,0) #thresholding the image
#(grayscale source image, threshod value to classify pixel values, max value, type of thresholding)
#output: threshold that was used, thresholded image

element = cv2.getStructuringElement(cv2.MORPH_CROSS,(3,3)) # get kernel: (shape,size,anchor=new cv.Point(-1.-1))
# anchor position [-1,-1] at center of image, only cross-shaped element depends on anchor position
done = False
while not done:
    eroded = cv2.erode(gscale,element) # erodes away boundaries of foreground object & removes small noise
    temp = cv2.dilate(eroded,element) # increase object area 
    temp = cv2.subtract(gscale,temp) # subtract images of same type & depth
    skel = cv2.bitwise_or(skel,temp)
    gscale = eroded.copy()
    zeros = size - cv2.countNonZero(gscale)
    if zeros==size:
        done = True


lines = cv2.HoughLines(skel,1,np.pi/180,130) #lines=vector storing lines (r,theta)
#(output of edge detector, resolution of r in pixels, resolution of theta in rad, threshold: min # of intersections, srn & stn)
a,b,c = lines.shape

#draw the lines 
for i in range(a):
    rho = lines[i][0][0]
    theta = lines[i][0][1]    
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))    
    
    cv2.line(img,(x1,y1),(x2,y2),(0,0,255),2, cv2.LINE_AA)
    #display lines: image, point1, point2, color, thickness int, line type, shift)

plt.subplot(121)#OpenCV reads images as BGR, this corrects so it is displayed as RGB
plt.plot(),plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB)) 
plt.title('Row Detection'), plt.xticks([]), plt.yticks([])
plt.subplot(122)
plt.plot(),plt.imshow(skel,cmap='gray')
plt.title('Skeletal Image'), plt.xticks([]), plt.yticks([])
plt.show()