Hough Transform_.idea_.name

Code

Hough Transform_.idea_Code.iml

<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
  <component name="NewModuleRootManager">
    <content url="file://$MODULE_DIR$" />
    <orderEntry type="inheritedJdk" />
    <orderEntry type="sourceFolder" forTests="false" />
  </component>
</module>

Hough Transform_.idea_encodings.xml

<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
  <component name="Encoding" useUTFGuessing="true" native2AsciiForPropertiesFiles="false" />
</project>

Hough Transform_.idea_misc.xml

<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
  <component name="ProjectRootManager" version="2" project-jdk-name="Python 2.7.3 (C:\Python27\python.exe)" project-jdk-type="Python SDK" />
</project>

Hough Transform_.idea_modules.xml

<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
  <component name="ProjectModuleManager">
    <modules>
      <module fileurl="file://$PROJECT_DIR$/.idea/Code.iml" filepath="$PROJECT_DIR$/.idea/Code.iml" />
    </modules>
  </component>
</project>

Hough Transform_.idea_scopes_scope_settings.xml

<component name="DependencyValidationManager">
  <state>
    <option name="SKIP_IMPORT_STATEMENTS" value="false" />
  </state>
</component>

Hough Transform_.idea_vcs.xml

<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
  <component name="VcsDirectoryMappings">
    <mapping directory="" vcs="" />
  </component>
</project>

Hough Transform_Edge detection.py

import numpy as np
import scipy.ndimage as sci
import numpy.linalg as la
import cv2
import math
import functions

def detectBall(img, gray):
    """
    :param img: colored image of opencv image type
    :param gray: grayscale image of opencv image type
    :return: image with circles marked. Can also return (radius, x, y) coordinates of circles by returning circle_indices
    """
    # =============================
    # Read ball image
    # =============================
    RESIZE = 1                                          # Resize width and height to width/RESIZE and
                                                        # height/RESIZE. Value of 1 means no resizing.

    ROWS, COLS = gray.shape                             # Get size of image
    img = cv2.resize(img, (COLS/RESIZE,ROWS/RESIZE))    # Resize images
    gray = cv2.resize(gray,(COLS/RESIZE,ROWS/RESIZE))
    b,g,r = cv2.split(img)                              # split image into b, g, r channels.
    ROWS, COLS = gray.shape                             # Get new size of image

    # =============================
    # Initialize variables
    # =============================
    BLURSIZE = 5
    THRESH = 10
    blur = functions.blur(img, BLURSIZE)
    # =============================
    # Detect Edges in Different channels and then recombine channel edges
    # =============================
    isEdgeGray, angleGradGray = functions.detect_edge(gray,THRESH, BLURSIZE)
    isEdgeBlue, angleGradBlue = functions.detect_edge(b,THRESH, BLURSIZE)
    isEdgeGreen, angleGradGreen = functions.detect_edge(g, THRESH, BLURSIZE)
    isEdgeRed, angleGradRed = functions.detect_edge(r, THRESH, BLURSIZE)

    isEdge = np.logical_or(isEdgeBlue, np.logical_or(isEdgeGreen, isEdgeRed))

    # Convert edgeImg from boolean to int. Every True pixel becomes 255.
    edgeImg = np.array(isEdge, dtype = np.uint8)
    edgeImg = np.multiply(edgeImg, 255)

    # ===============================
    # Set parameters for hough transform
    radMin = ROWS/9                 # minimum radius of circles detected
    radMax = (ROWS + COLS)/2        # maximum radius of circles detected
    distMin = ROWS/8                # minimum distance apart for circles detected
    votesMin = 12                    # threshold for hough many votes a (radius, x, y) triplet needs to have before it's
                                    # considered a circle

    img, circle_indices = functions.hough_circle(img, gray, isEdge, angleGradGray, radMin, radMax, distMin, votesMin)
    return img

gray = cv2.imread("overlapping ball.jpg", 0)
img = cv2.imread("overlapping ball.jpg", -1)
detectedImg = detectBall(img, gray)
print type(detectedImg)
print "hi"
cv2.imwrite("overlapping ball pre reduce.jpg", detectedImg)
cv2.imshow("wtf", detectedImg)

Hough Transform_functions.py

import numpy as np
import scipy.ndimage.filters as filter
import scipy.signal as signal
import numpy.linalg as la
import cv2
import math

np.set_printoptions(edgeitems = 4)

def grad_to_edge_angle(angleGrad):
    """
    :param angleGrad: angle of gradient
    :return: approximate angle of edge direction. Directions can be 0, pi/4, pi/2, 3pi/4.

    Angles ~0 or pi become pi/2
    Angles ~pi/4 or -3pi/4 become 3pi/4
    Angles ~ pi/2 or -pi/2 become 0
    Angles ~ 3pi/4 or -pi/4 become pi/4
    """
    angleGrad = angleGrad * 4 / math.pi + 2
    angleGrad = round(angleGrad)
    return (angleGrad % 4)* math.pi / 4

def grad_to_approx_angle(angleGrad):
    """
    :param angleGrad: angle of gradient
    :return: approximate angle of gradient direction. Directions can be 0, pi/4, pi/2, 3pi/4.

    Angles ~0 or pi become 0
    Angles ~pi/4 or -3pi/4 become pi/4
    Angles ~ pi/2 or -pi/2 become pi/2
    Angles ~ 3pi/4 or -pi/4 become 3pi/4
    """
    angleGrad = angleGrad * 4 / math.pi
    angleGrad = round(angleGrad)
    return (angleGrad % 4)* math.pi / 4

def canny_edge(img, HIGHTHRESH, LOWTHRESH, BLURSIZE):
    """
    Detects edges with the Canny edge detector. Unused in this project because run time was too high.
    :param img:
    :param HIGHTHRESH: gradient magnitude threshold value needed to start an edge
    :param LOWTHRESH: gradient magnitude threshold value needed to continue an edge
    :param BLURSIZE: size of Gaussian kernel
    :return: boolean matrix of edge pixels, and matrix of gradient angle each pixel
    """
    ROWS, COLS = img.shape

    # =============================
    # "Generate Sobel and Gaussian kernels"
    # =============================

    sobelx = np.array([[-1.0, 0.0, 1.0],
              [-2.0, 0.0, 2.0],
              [-1.0, 0.0, 1.0]])

    sobely = np.array([[-1.0, -2.0, -1.0],
              [0.0, 0.0, 0.0],
              [1.0, 2.0, 1.0]])

    gauss = cv2.getGaussianKernel(BLURSIZE, BLURSIZE/3)
    gauss = np.outer(gauss, gauss)

    # =============================
    # "Apply Gaussian mask"
    # =============================
    img = cv2.filter2D(img, -1, gauss)
    gradImgx = np.copy(img)
    gradImgy = np.copy(img)
    gradImgy = signal.convolve2d(gradImgy, sobely, mode = "same")/4.0
    gradImgx = signal.convolve2d(gradImgx, sobelx, mode = "same")/4.0

    # =============================
    # "Calculate gradient angle and magnitude at each pixel"
    # =============================
    angleGrad = np.arctan2(gradImgy, gradImgx)
    magGrad = np.sqrt(np.square(gradImgy) + np.square(gradImgx))
    # =====================
    # "Trace along edges"
    # =====================
    # Go to every pixel. If  >= HIGHTHRESH, set it as an edge. Trace along edge direction to find more edges.

    isEdge = np.greater_equal(magGrad, HIGHTHRESH)
    isEdge = np.array(isEdge, dtype = np.uint8)
    # temp = np.multiply(isEdge, 255)
    # cv2.imshow("Edges with thresholding", temp)
    for r in range(0,ROWS):
        for c in range(0,COLS):
            # Hysteresis
            if isEdge[r][c]:
                trace_edge(r, c, LOWTHRESH, isEdge, magGrad, angleGrad)
                trace_edge(r, c, LOWTHRESH, isEdge, magGrad, np.multiply(angleGrad, -1))

    # temp = np.multiply(isEdge, 255)
    # cv2.imshow('Edge w/ Hysteresis', temp)

    # =====================
    # Nonmaximal suppression
    # =====================
    for r in range(0,ROWS):
        for c in range(0,COLS):
            if isEdge[r][c]:
                suppress_edge(r, c, magGrad, angleGrad, isEdge)

    # isEdge = np.multiply(isEdge, 255)

    isEdge = np.greater_equal(magGrad, LOWTHRESH)
    isEdge = np.array(isEdge, dtype = np.uint8)
    return isEdge, angleGrad
    # cv2.imshow('Edge with hysteresis and nonmaximal suppression', isEdge)

def trace_edge(r, c, LOWTHRESH, isEdge, magGrad, angleGrad):
    """
    Helper function for canny_edge. Takes pixel at r,c and traces along its gradient for edges
    :param r: row of edge
    :param c: column of edge
    :param LOWTHRESH: gradient magnitude threshold value needed to continue an edge
    :param isEdge: boolean matrix of whether a pixel is an edge or not
    :param magGrad: matrix of gradient magnitude at each pixel
    :param angleGrad: matrix of gradient angle at each pixel
    :return: None. (All operations are performed on the given isEdge matrix)
    """
    edgeAngle = grad_to_edge_angle(angleGrad[r][c])
    deltaX = round(math.cos(edgeAngle))
    deltaY = round(math.sin(edgeAngle))

    r = r+deltaY
    c = c+deltaX

    if r >= 0 and r < np.shape(isEdge)[0] and c > 0 and c < np.shape(isEdge)[1] and magGrad[r][c] >= LOWTHRESH:
        isEdge[r][c] = 1
        # trace_edge(r, c, LOWTHRESH, isEdge, magGrad, angleGrad)

def suppress_edge(r, c, magGrad, angleGrad, isEdge):
    """
    Helper function for canny_edge. Suppresses edges that are not maximal edges.
    :param r: row of edge
    :param c: column of edge
    :param isEdge: boolean matrix of whether a pixel is an edge or not
    :param magGrad: matrix of gradient magnitude at each pixel
    :param angleGrad: matrix of gradient angle at each pixel
    :return: None. (All operations are performed on the given isEdge matrix)
    """
    currentGrad= grad_to_approx_angle(angleGrad[r][c])
    deltaX = round(math.cos(currentGrad))
    deltaY = round(math.sin(currentGrad))

    rstar = r + deltaY
    cstar = c + deltaX
    
    if rstar >= 0 and rstar < np.shape(isEdge)[0] and cstar > 0 and cstar < np.shape(isEdge)[1] and isEdge[rstar][cstar]:
        nextGrad = angleGrad[r][c]
        if nextGrad > currentGrad:
            isEdge[r][c] = 0
        else:
            isEdge[rstar][cstar] = 0

        suppress_edge(rstar, cstar, magGrad, angleGrad,isEdge)

def blur(img, blursize):
    """
    :param img: image to be blurred
    :param blursize: size of gaussian kernel to blur with
    :return: blurred image
    """
    # =============================
    # Generate Gaussian kernel
    # =============================
    # first parameter is kernel size, second parameter is standard deviation
    gauss = cv2.getGaussianKernel(blursize, blursize/3)
    gauss = np.outer(gauss, gauss)
    # =============================
    # Filter with Gaussian kernel
    # =============================
    img = cv2.filter2D(img, -1, gauss)
    return img

def detect_edge(img, thresh, blursize):
    """
    Detects edges in an image
    :param img: image to detect edges in
    :param thresh: minimum magnitude of gradient before it's considered an edge
    :param blursize: size of Gaussian blur filter
    :return: boolean matrix of edges in image
    """
    # =============================
    # "Generate Sobel kernel"
    # =============================

    img = blur(img, blursize)
    sobelx = np.array([[-1, 0, 1],
              [-2, 0, 2],
              [-1, 0, 1]])

    sobely = np.array([[-1, -2, -1],
              [0, 0, 0],
              [1, 2, 1]])

    # find partial derivatives of image by convolving with Sobel kernels
    gradImgx = np.copy(img)
    gradImgy = np.copy(img)
    gradImgy = signal.convolve2d(gradImgy, sobely, mode = "same")/4.0
    gradImgx = signal.convolve2d(gradImgx, sobelx, mode = "same")/4.0

    # =============================
    # "Calculate gradient angle and magnitude at each pixel"
    # =============================
    angleGrad = np.arctan2(gradImgy, gradImgx)
    magGrad = np.sqrt(np.square(gradImgy) + np.square(gradImgx))

    isEdge = np.greater_equal(magGrad, thresh)
    return isEdge, angleGrad


def hough_circle(colorImg, img, isEdge, angleGrad, radMin, radMax, distMin, votesMin):
    """
    Uses the Hough circle transform to detect circles in image

    :param img: Image to be operated on. size rows x cols. cv2 image format
    :param isEdge: rows x col size boolean matrix. True if pixel is edge pixel.
    :param angleGrad: rows x col size float matrix. gradient direction at each pixel.
    :param radMin: int, minimum radius of circles detected
    :param radMax: int, maximum radius of circles detected
    :param distMin: int, minimum distance between circles detected
    :param votesMin: int, minimum threshold for accumulator value. The higher this value is,
                    the less false circles detected

    :return: numpy array of circles detected. Stored as a = [[radii], [rows], [cols]]
            where
            a[0][i] is the radius of the ith circle
            a[1][i] is the row where the center of the ith circle is located
            a[2][i] is the col where the center of the ith circle is located
    """
    rows, cols = img.shape
    # Initialize accumulator matrix
    param_votes = np.zeros((radMax, rows, cols), dtype = int)
    # Initialize accumulator matrix of local maximums
    param_maxes = np.zeros((radMax, rows, cols), dtype = int)

    # Perform Hough transform on each edge point with radii ranging from radMin to radMax
    for r in range(0, rows):
        for c in range(0, cols):
            # Get gradient angle at pixel (r,c)
            # (r, c) is the proposed center for the circle
            theta = angleGrad[r][c]
            if isEdge[r][c]:
                for rad in range(radMin, radMax):
                    # For each proposed, center,
                    # Test each possible radius and give a vote to the resulting (rad, center) coordinates
                    xstar = round(c - rad * math.cos(theta))
                    ystar = round(r - rad * math.sin(theta))
                    if xstar >=0 and xstar < cols and ystar >=0 and ystar < rows:
                        param_votes[rad][ystar][xstar] += 1

                    # Add pi to the gradient angle so the angle also points in the other direction because
                    # who knows is the gradient is pointing towards or away from the center of the circle.
                    theta = theta + math.pi
                    xstar = round(c - rad * math.cos(theta))
                    ystar = round(r - rad * math.sin(theta))

                    # Add a vote to (rad, center) if the center is within the bounds of the image
                    if xstar >=0 and xstar < cols and ystar >=0 and ystar < rows:
                        param_votes[rad][ystar][xstar] += 1

    # Find local maximums in the accumulator matrix.
    param_votes[param_votes < votesMin] = 0
    param_max = filter.maximum_filter(param_votes, size = (distMin, distMin, distMin))
    param_max = np.logical_and(param_votes, param_max)

    # Find indices of maximums
    circle_indices = np.where(param_max)

    # Will most likely find too many circles, so use helper function to average clusters of circles
    # Also removes false positives

    circle_indices = reduce_circles(circle_indices, distMin, param_votes)

    total_circles = len(circle_indices[0])

    # draw detected circles onto img
    for i in range(0, total_circles):
        rad = circle_indices[0][i]
        r = circle_indices[1][i]
        c = circle_indices[2][i]
        cv2.circle(colorImg, (c, r), rad, (0,255,0), thickness = 3)
        cv2.circle(colorImg, (c, r), 1, (0, 255, 0), thickness = 2)

    return colorImg, circle_indices

def reduce_circles(circle_indices, distMin, param_votes):
    """
    Helper function for hough_circle
    Averages clusters of circles

    :param circle_indices: numpy array of circles detected. Stored as a = [[radii], [rows], [cols]]
            where
            a[0][i] is the radius of the ith circle
            a[1][i] is the row where the center of the ith circle is located
            a[2][i] is the col where the center of the ith circle is located
    :param distMin: int, minimum distance between circles detected. (The function averages circles clustered within
                    a distMin x distMin square)
    :param votesMin: int, minimum threshold for accumulator value. The higher this value is,
                    the less false circles detected
    :return: numpy array of circles
    """

    rad_indices = circle_indices[0]
    r_indices = circle_indices[1]
    c_indices = circle_indices[2]
    print rad_indices
    print r_indices
    print c_indices

    total_circles = len(rad_indices)
    print "total circles ", total_circles
    new_circle_indices = [[],[],[]]


    radtemp = 0                 # accumulator value for radii
    rtemp = 0                   # accumulator value for rows
    ctemp = 0                   # accumulator value for columns
    weightcount = 0             # accumulator value for weight of each circle.
                                # A proposed circle with n votes contributes n^2 weight
    count = 0                   # total number of circles detected in a cluster.

    minThresh = total_circles / 5   # Each cluster has a count that records number of circles in that cluster.
                                    # If count is below minThresh, the cluster is disregarded as a false positive.

    # Remove an (radius, center) once it's been determined part of a cluster. Iterate until no more proposed circles to
    # iterate through. Iterate backwards so deleting things doesn't mess up indexing.
    while np.size(rad_indices)!= 0:
        # get center
        r0 = r_indices[0]
        c0 = c_indices[0]
        for i in reversed(range(0, total_circles)):
            # If iterate through all other proposed circles, and find ones with center within a
            # distMin x distMin square.
            if abs(r0 - r_indices[i]) < distMin and abs(c0 - c_indices[i] < distMin):
                # (radius, center) values of proposed circle wwe've iterated to
                rad = rad_indices[i]
                r = r_indices[i]
                c = c_indices[i]

                # increase accumulator values, weighted by weight
                weight = param_votes[rad][r][c]*param_votes[rad][r][c]
                radtemp += rad * weight
                rtemp += r * weight
                ctemp += c * weight

                # delete proposed circle from rad_indices
                rad_indices = np.delete(rad_indices, i)
                r_indices = np.delete(r_indices, i)
                c_indices = np.delete(c_indices, i)
                weightcount += weight
                count += 1
        total_circles = len(rad_indices)

        # if more circles in cluster than minThresh, average the radii and center values to combine the cluster into
        # one circle
        if count >= minThresh:
            new_circle_indices[0].append(radtemp / weightcount)
            new_circle_indices[1].append(rtemp / weightcount)
            new_circle_indices[2].append(ctemp /weightcount)

        # Optional code. If you know you won't have circles within circles, uncomment this code.
        # This removes circles with centers inside of other circles.

        # for i in reversed(range(0, total_circles)):
        #     if abs(rtemp / weightcount - r_indices[i]) < radtemp / weightcount and abs(ctemp /weightcount - c_indices[i]) < radtemp / weightcount:
        #         rad_indices = np.delete(rad_indices, i)
        #         r_indices = np.delete(r_indices, i)
        #         c_indices = np.delete(c_indices, i)
        #     total_circles = len(rad_indices)

        # reset accumulators to 0 for the next cluster
        radtemp = 0
        rtemp = 0
        ctemp = 0
        count = 0
        weightcount = 0

    return np.array(new_circle_indices)