In [19]:
import numpy as np
import imageio
import math
import cv2
import matplotlib.pyplot as plt
import tkinter as tk
from tkinter import simpledialog
import random
import itertools
#Question 1
def FindLines1(image_path_name,alpha,threshold):
    img_orig = cv2.imread(image_path_name);
    img_pts_display = cv2.imread(image_path_name);
    img = cv2.imread(image_path_name);
    #Threshold justification: Images considered for this experiment range from 640x480 to 1600x1200 and threshold represents
    #the minimum length of line segment and I feel 150 pixel line would be a good choice since anything less than that could be a
    #part of curve most likely
    width, height = img.shape[:2];
    #Rho value could range upto the length of the diagonal of the image
    diag_len = int(round(math.sqrt(width * width + height * height)))  
    cos_t = math.cos(alpha)
    sin_t = math.sin(alpha)
    # 1D Hough accumulator array of rho values
    accumulator = np.zeros((diag_len), dtype=np.uint8)
    #Convert image to grayscale for edge detection
    gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    #Canny edge detection
    are_edges = cv2.Canny(gray, 50, 150);
    #Choosing only non-zero edges which means all the significant lines found by the algorithm are found by hand 
    y_idxs, x_idxs = np.nonzero(are_edges)
    # Vote in the hough accumulator
    for i in range(len(x_idxs)):
        x = x_idxs[i]
        y = y_idxs[i]
        rho = int(round(x * cos_t + y * sin_t))
        accumulator[rho] += 1

    #Finding peaks in 1d accumulator array
    d_peaks_array = np.where( accumulator >= threshold);
    d_peaks_array = d_peaks_array[0];
    for rho in d_peaks_array:
        a = cos_t;
        b = sin_t;
        #Points that contributed to selected lines
        x0, y0 = cos_t*rho, sin_t*rho;
        pt1 = ( int(x0+1000*(-b)), int(y0+1000*(a)) )
        pt2 = ( int(x0-1000*(-b)), int(y0-1000*(a)) )
        cv2.line(img, pt1, pt2, (255, 0, 0), 2, cv2.LINE_AA);
    plt.subplot(121),plt.imshow(img_orig)
    plt.title('Original Image'), plt.xticks([]), plt.yticks([]);
    plt.subplot(122),plt.imshow(img)
    plt.title('Hough Image'), plt.xticks([]), plt.yticks([]);
    plt.show()
    plt.subplot(121),plt.imshow(gray,cmap = 'gray');
    plt.title('Significant lines w/o algo'), plt.xticks([]), plt.yticks([]);
    plt.show()
    
    
#Question 2    
def compute_image_gradient(Edge_Gradient,Gradient_angle,sobelx,sobely,array_length):
    for i in range(0,array_length):
        Gx = sobelx[i];
        Gy = sobely[i];
        G = np.square(Gx) + np.square(Gy)
        G = np.sqrt(G);
        Edge_Gradient[i] = G;
        Theta = np.zeros((len(Gy)),dtype=np.float)
        for j in range(0,len(Gy)):
            if(Gx[j] != 0):
                temp_deg = math.degrees(math.atan2(Gy[j],Gx[j]));
                #Putting angles that differ by 180 degree or pi into one bin
                if(not((temp_deg >= 0) and (temp_deg <= 180))):
                    if(temp_deg < 0):
                        temp_deg = temp_deg + 180;
                Theta[j] = temp_deg;
        Gradient_angle[i] = Theta;
    return Edge_Gradient,Gradient_angle

#Question 3
def CompareHoughP(image_path_name):
    img = cv2.imread(image_path_name);
    #Lists that store the point values((x1,y1),(x2,y2)) of the lines found through Hough and HoughP so that it can be used for
    #vanishing point detection
    hough_lines = [];
    houghP_lines = [];
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    dst = cv2.Canny(gray, 100, 200)
    lines= cv2.HoughLines(dst, 1, math.pi/180.0, 250, np.array([]), 0, 0)
    a,b,c = lines.shape
    for i in range(a):
        rho = lines[i][0][0]
        theta = lines[i][0][1]
        a = math.cos(theta)
        b = math.sin(theta)
        x0, y0 = a*rho, b*rho
        pt1 = ( int(x0+1000*(-b)), int(y0+1000*(a)) )
        pt2 = ( int(x0-1000*(-b)), int(y0-1000*(a)) );
        x_1 = int(x0 + 1000 * (-b))
        y_1 = int(y0 + 1000 * (a))
        x_2 = int(x0 - 1000 * (-b))
        y_2 = int(y0 - 1000 * (a))
        hough_lines.append(((x_1, y_1), (x_2, y_2)));
        cv2.line(img, pt1, pt2, (255, 0, 0), 2, cv2.LINE_AA)
    plt.subplot(121),plt.imshow(img)
    plt.title('Hough Image'), plt.xticks([]), plt.yticks([]);
    minLineLength = 50
    maxLineGap = 10
    img_1 = cv2.imread(image_path_name)
    gray = cv2.cvtColor(img_1,cv2.COLOR_BGR2GRAY)
    dst = cv2.Canny(gray, 100, 200)
    lines1 = cv2.HoughLinesP(dst,1,np.pi/180,50,minLineLength,maxLineGap)
    for line in lines1:    
        for x1,y1,x2,y2 in line:
            houghP_lines.append(((x1,y1),(x2,y2)));
            cv2.line(img_1, (x1,y1), (x2,y2), (0, 255, 0), 2, cv2.LINE_AA);
    plt.subplot(122),plt.imshow(img_1)
    plt.title('Probabilstic Hough Image'), plt.xticks([]), plt.yticks([]);
    plt.show();
    return hough_lines,houghP_lines;
 
#Question 2(finding lines from peak alpha values of edge histograms)
def FindLines2(image_path_name,angle_list,threshold):
    #List that stores the point values((x1,y1),(x2,y2)) of the lines found through edge histogram peaks so that it can be used for
    #vanishing point detection
    part2_lines = []
    img_orig = cv2.imread(image_path_name);
    img = cv2.imread(image_path_name);
    width, height = img.shape[:2]
    diag_len = int(round(math.sqrt(width * width + height * height)))
    cos_t = np.cos(angle_list)
    sin_t = np.sin(angle_list)
    num_thetas = len(angle_list)
    # 2D Hough accumulator array of theta vs rho
    accumulator = np.zeros((diag_len, num_thetas), dtype=np.uint8)
    #Convert image to grayscale for edge detection step
    gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    #Canny edge detection
    are_edges = cv2.Canny(gray, 50, 150);
    y_idxs, x_idxs = np.nonzero(are_edges)

    # Vote in the hough accumulator
    for i in range(len(x_idxs)):
        x = x_idxs[i]
        y = y_idxs[i]
        for t_idx in range(num_thetas):
            # Calculate rho for each theta value(alpha peaks) and increment the counter by 1 for that rho,theta pair
            rho = int(round(x * cos_t[t_idx] + y * sin_t[t_idx]));
            if((rho >= 0 and rho < diag_len)):
                accumulator[rho][t_idx] += 1 
    accumulator = np.argmax(accumulator,axis=0);
    #print(accumulator);
    for i in range(0,num_thetas):
        rho = accumulator[i]; 
        a = cos_t[i];
        b = sin_t[i];
        #Points that contributed to selected lines
        x0, y0 = a*rho, b*rho;
        pt1 = ( int(x0+1000*(-b)), int(y0+1000*(a)) );
        pt2 = ( int(x0-1000*(-b)), int(y0-1000*(a)) );
        x_1 = int(x0 + 1000 * (-b))
        y_1 = int(y0 + 1000 * (a))
        x_2 = int(x0 - 1000 * (-b))
        y_2 = int(y0 - 1000 * (a))
        part2_lines.append(((x_1, y_1), (x_2, y_2)));
        cv2.line(img, pt1, pt2, (255, 0, 0), 2, cv2.LINE_AA);
    plt.subplot(121),plt.imshow(img_orig)
    plt.title('Original Image'), plt.xticks([]), plt.yticks([]);
    plt.subplot(122),plt.imshow(img)
    plt.title('Edge histogram lines'), plt.xticks([]), plt.yticks([]);
    plt.show()
    return part2_lines;

#Question 2(Edge histogram and peak alpha value list computation)
def EdgeHistogramFindLines(image_path_name):
    part2_lines = [];
    img = cv2.imread(image_path_name);
    #Convert the image to grayscale for gray edge detection
    gray_image = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY);
    height, width = gray_image.shape[:2];
    #Apply Sobel operator, includes Guassian smoothing
    sobelx = cv2.Sobel(gray_image,cv2.CV_64F,1,0,ksize=5)
    sobely = cv2.Sobel(gray_image,cv2.CV_64F,0,1,ksize=5)
    #After sobel operator , find the image gradients Edge_Gradient(G)=sqrt(Gx^2+Gy^2), Angle(θ)=tan−1(Gy/Gx);
    Edge_Gradient_Gray = np.zeros((len(sobelx),len(sobelx[0])),dtype=np.float)
    Gradient_angle_gray = np.zeros((len(sobelx),len(sobelx[0])),dtype=np.float);
    array_length = len(sobelx);
    Edge_Gradient_Gray,Gradient_angle_gray = compute_image_gradient(Edge_Gradient_Gray,Gradient_angle_gray,sobelx,sobely,array_length);
    #Create gray edge histograms(36 bin)
    Gradient_angle_gray = Gradient_angle_gray.astype(np.uint8);
    histr = cv2.calcHist([Gradient_angle_gray],[0],None,[36],[0,180]);
    #Finding peaks of the histogram, returns indices of peak values
    peak_alpha_list = np.argmax(histr, axis=0)
    #Since each bin is 5 degrees, value in peak_alpha_list refers to index of peak angles , for ex: if peak_alpha_list has value 1
    #it means that most of the edges have orientations in the range of 5 to 10 degrees. Hence we consider alpha = 7(median value of 5-10), and range = [7-3,7+3] incrementing/decrementing 1 at a time i.e
    #Our significant angles are 7, 6, 5, 4, 8, 9, 10. Negative values(if any) are converted to absolute.
    #Finding all the significant edge orientations()
    angle_list = np.zeros((len(peak_alpha_list) * 7), dtype=np.uint8);
    ctr = 0;
    for val in peak_alpha_list:
        temp = ((val+1)*5) - 3;
        angle_list[ctr:ctr+7] = [temp,abs(temp-1),abs(temp-2),abs(temp-3),temp+1,temp+2,temp+3];
        ctr+=7;
    
    #Angles converted to radians 
    angle_list = np.radians(angle_list);
    #Function to find lines using 2D accumualtor with the angle list determined from edge histogram peaks
    part2_lines = FindLines2(image_path_name,angle_list,100);
    return part2_lines;

#Question 4(Helper function to detect vanishing points)
#Function to find determinant given 2 columns of a matrix
def det(a, b):
    return a[0] * b[1] - a[1] * b[0]

#Question 4
# Find intersection point of two lines
def line_intersection(line1, line2):
    x_diff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
    y_diff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])

    div = det(x_diff, y_diff)
    if div == 0:
        return None  # Lines don't cross

    d = (det(*line1), det(*line2))
    x = det(d, x_diff) / div
    y = det(d, y_diff) / div

    return x, y

#Question 4
# Find intersections between multiple lines
def find_intersections(lines):
    intersections = []
    for i, line_1 in enumerate(lines):
        for line_2 in lines[i + 1:]:
            if not line_1 == line_2:
                intersection = line_intersection(line_1, line_2)
                if intersection:  # If lines cross, then add
                    intersections.append(intersection)

    return intersections

#Question 4(Detect vanishing point, using the intersections and the grid size)
# Given intersections, find the grid where most intersections occur and treat as vanishing point
def find_vanishing_point(img, grid_size, intersections):
    # Image dimensions
    image_height = img.shape[0]
    image_width = img.shape[1]

    # Grid dimensions
    grid_rows = (image_height // grid_size) + 1
    grid_columns = (image_width // grid_size) + 1

    # Current cell with most intersection points
    max_intersections = 0
    best_cell = (0.0, 0.0)

    for i, j in itertools.product(range(grid_rows), range(grid_columns)):
        cell_left = i * grid_size
        cell_right = (i + 1) * grid_size
        cell_bottom = j * grid_size
        cell_top = (j + 1) * grid_size
        
        current_intersections = 0  # Number of intersections in the current cell
        for x, y in intersections:
            if cell_left < x < cell_right and cell_bottom < y < cell_top:
                current_intersections += 1

        # Current cell has more intersections that previous cell (better)
        if current_intersections > max_intersections:
            max_intersections = current_intersections
            best_cell = ((cell_left + cell_right) / 2, (cell_bottom + cell_top) / 2)
            
    if best_cell[0] != None and best_cell[1] != None:
        rx1 = int(best_cell[0] - grid_size / 2)
        ry1 = int(best_cell[1] - grid_size / 2)
        rx2 = int(best_cell[0] + grid_size / 2)
        ry2 = int(best_cell[1] + grid_size / 2)
        #Marking the vanishing point as a small rectangle on the image
        cv2.rectangle(img, (rx1, ry1), (rx2, ry2), (0, 255, 0), 10)
        plt.subplot(121),plt.imshow(img)
        plt.title('Vanishing point'), plt.xticks([]), plt.yticks([]);
        plt.show();
    return best_cell


if __name__ == '__main__':
    root = tk.Tk() # Create an instance of tkinter
    root.withdraw();# Close the parent window
    Question1_response = simpledialog.askstring(title = "Question 1 ",
    prompt = "Question 1: Enter N/n for seeing demo of results on pre-stored images in this file or Y/y for entering the image path of your choice\n For input other than Y/y demo results shown");
    #If you choose demo, then place all the images that are submitted in the same folder as the python source file,
    #Else enter the correct image path name for which you wish to see the Radon transform/Question 1 results
    image_path_name = "";
    demo_flag = True;
    if(Question1_response == "Y" or Question1_response == "y"):
        demo_flag = False;
    if(demo_flag == False):
        root = tk.Tk() # Create an instance of tkinter
        root.withdraw();# Close the parent window
        image_path_name = simpledialog.askstring(title = "Question 1 ",
                                    prompt = "Enter the image path for which you wish to see Radon Transform");
        alpha_str = simpledialog.askstring(title = "Question 1 ",
                                    prompt = "Enter the alpha value in degrees(range preferred 0-90) For ex: if your alpha is 90 degrees enter 90");
        alpha_f = float(alpha_str);
        #Mapping the angle to the correct value between 0-90 for out of range values
        if(alpha_f < 0):
            alpha_f = alpha_f + 360;
        if(not((alpha_f >= 0) and (alpha_f <= 90.0))):
            if(alpha_f <= 180):
                alpha_f = alpha_f - 90;
            elif(alpha_f <= 270):
                alpha_f = alpha_f - 180;
            elif(alpha_f <= 360):
                alpha_f = alpha_f - 270;
        alpha_r = (alpha_f) * (math.pi)/180;
        #Entering threshold as 150
        FindLines1(image_path_name,alpha_r,150);
    else: 
        #Demo results for 2 images for 2 angles
        #Image 1 : Alpha = 0 degrees
        image_path_name = "ST2MainHall4001_HW2_1.jpg";
        alpha_r = 0 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
        #Image 1 : Alpha = 90 degrees
        alpha_r = 90 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
        #Image 2 : Alpha = 0 degrees
        image_path_name = "ST2MainHall4058_HW2_2.jpg";
        alpha_r = 0 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
        #Image 2 : Alpha = 90 degrees
        alpha_r = 90 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
        #Generating random angle for alpha between 0 to 90 degrees, using Image 2
        rand_alpha1 = random.randint(0,91);
        alpha_r = rand_alpha1 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
        rand_alpha2 = random.randint(0,91);
        while(rand_alpha2 == rand_alpha1):
            rand_alpha2 = random.randint(0,361);
        alpha_r = rand_alpha2 * (math.pi)/180;
        FindLines1(image_path_name,alpha_r,150);
    #Question 2 : Building edge histograms from edge orientations
    #Demonstrating on 3 images
    #Lists to store lines to detect vanishing points later
    part2_lines_img1 = [];
    part2_lines_img2 = [];
    part2_lines_img3 = [];
    hough_line_img1 = [];
    hough_line_img2 = [];
    hough_line_img3 = [];
    houghP_lines_img1 = [];
    houghP_lines_img2 = [];
    houghP_lines_img3 = [];
    #Image 1
    image_path_name = "ST2MainHall4001_HW2_1.jpg";
    part2_lines_img1 = EdgeHistogramFindLines(image_path_name);
    #Image 2
    image_path_name = "ST2MainHall4058_HW2_2.jpg";
    part2_lines_img2 = EdgeHistogramFindLines(image_path_name);
    #Image 3
    image_path_name = "ST2MainHall4079_HW2_3.jpg";
    part2_lines_img3 = EdgeHistogramFindLines(image_path_name);
    #Question 3 : Comparing Hough and Probabilistic Hough transform for above 3 images
    #Image 1
    image_path_name = "ST2MainHall4001_HW2_1.jpg";
    hough_line_img1,houghP_lines_img1 = CompareHoughP(image_path_name);
    #Image 2
    image_path_name = "ST2MainHall4058_HW2_2.jpg";
    hough_line_img2,houghP_lines_img2 = CompareHoughP(image_path_name);
    #Image 3
    image_path_name = "ST2MainHall4079_HW2_3.jpg";
    hough_line_img3,houghP_lines_img3 = CompareHoughP(image_path_name);
   
    #Question 4: Finding vanishing points on the lines found in part 2 and part 3
    #Part 2 Image 1 
    print("Vanishing points for lines found in part 2");
    print("Image 1");
    image_path_name = "ST2MainHall4001_HW2_1.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(part2_lines_img1);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Part 2 Image 2
    print("Image 2");
    image_path_name = "ST2MainHall4058_HW2_2.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(part2_lines_img2);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Part 2 Image 3
    print("Image 3");
    image_path_name = "ST2MainHall4079_HW2_3.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(part2_lines_img3);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Part 3, Lines found from Hough Transform
    #Image 1
    print("Vanishing points for lines found in part 3 hough transform");
    print("Image 1");
    image_path_name = "ST2MainHall4001_HW2_1.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(hough_line_img1);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Image 2
    print("Image 2");
    image_path_name = "ST2MainHall4058_HW2_2.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(hough_line_img2);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Image 3
    print("Image 3");
    image_path_name = "ST2MainHall4079_HW2_3.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(hough_line_img3);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Part 3, Lines found from Probabilistic Hough Transform
    #Image 1
    print("Vanishing points for lines found in part 3 probabilstic hough transform");
    print("Image 1");
    image_path_name = "ST2MainHall4001_HW2_1.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(houghP_lines_img1);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Image 2
    print("Image 2");
    image_path_name = "ST2MainHall4058_HW2_2.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(houghP_lines_img2);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
    #Image 3
    print("Image 3");
    image_path_name = "ST2MainHall4079_HW2_3.jpg";
    img = cv2.imread(image_path_name);
    intersections = find_intersections(houghP_lines_img3);
    if intersections:
        grid_size = min(img.shape[0], img.shape[1]) // 27
        vanishing_point = find_vanishing_point(img, grid_size, intersections)
Vanishing points for lines found in part 2
Image 1
Image 2
Image 3
Vanishing points for lines found in part 3 hough transform
Image 1
Image 2
Image 3
Vanishing points for lines found in part 3 probabilstic hough transform
Image 1
Image 2
Image 3
In [ ]:
 
In [ ]: