diff --git a/src/ImageFeatures.jl b/src/ImageFeatures.jl
index f1d2173..82c2861 100644
--- a/src/ImageFeatures.jl
+++ b/src/ImageFeatures.jl
@@ -66,6 +66,7 @@ export
 
     #Lines
     hough_transform_standard,
+    hough_line_probabilistic,
 
     #Circles
     hough_circle_gradient
diff --git a/src/houghtransform.jl b/src/houghtransform.jl
index 85f9aaf..891f946 100644
--- a/src/houghtransform.jl
+++ b/src/houghtransform.jl
@@ -103,6 +103,292 @@ threshold::Integer, linesMax::Integer) where T<:Union{Bool,Gray{Bool}}
 
 end
 
+"""
+```
+lines = hough_line_probabilistic(image, ρ, θ, threshold, lineLength, lineGap, linesMax)
+```
+
+Returns lines :
+      vector of lines identified, lines in format ((r0, c0), (r1, c1))
+      indicating line start and end.
+
+The lines are generated by applying hough transform on the image.
+
+Parameters:
+-    `image`       = Image to be transformed (eltype should be `Bool`)
+-    `ρ`           = Discrete step size for perpendicular length of line
+-    `θ`           = List of angles for which the transform is computed
+-    `threshold`   = Accumulator threshold for line detection
+-    'lineLength'  = minimum length of a good_line
+-    'lineGap'     = minimum gap between two different lines.
+-    `linesMax`    = Maximum no of lines to return
+
+# Example
+```julia
+julia> img = load("line.jpg");
+
+julia> img_edges = canny(img, (Percentile(0.99), Percentile(0.97)), 1);
+
+julia> lines = hough_line_probabilistic(img_edges, 1, linspace(0,π,180),30,30,10,10)
+10-element Array{NTuple{4,Int64},1}:
+ (186, 283, 20, 283) 
+ (186, 20, 20, 20)   
+ (200, 218, 200, 291)
+ (20, 68, 20, 180)   
+ (186, 85, 186, 197) 
+ (48, 59, 69, 199)   
+ (50, 58, 65, 160)   
+ (200, 35, 200, 147) 
+ (20, 186, 20, 282)  
+ (155, 138, 75, 198)  
+```
+May use LineSegment of ImageDraw to draw lines.
+
+References
+    ----------
+    .. [1] C. Galamhos, J. Matas and J. Kittler, "Progressive probabilistic
+           Hough transform for line detection", in IEEE Computer Society
+           Conference on Computer Vision and Pattern Recognition, 1999.
+"""
+type Param
+       numangle::Int64
+       constadd::Float64
+       accumulator_matrix::AbstractArray{Int64, 2}
+       mask::AbstractArray{Bool, 2}
+       nzloc::Vector{Tuple{Int64,Int64}}
+       lines::Vector{Tuple{Int64, Int64, Int64, Int64}}
+       shift::Int64
+       threshold::Int64
+       lineLength::Int64
+       lineGap::Int64
+end
+
+type Sample
+    x0::Int64
+    y0::Int64
+    dx0::Int64
+    dy0::Int64
+end    
+
+#function to mark and collect all non zero points
+function collect_points(params, img)
+    for pix in CartesianRange(size(img))
+        pix1 = (pix[1], pix[2])
+        if(img[pix])
+            push!(params.nzloc, pix1)
+            params.mask[pix] = true
+        else
+            params.mask[pix] = false
+        end
+    end
+end    
+
+#function to update the accumulator matrix for every point selected
+function update_accumulator(params, point, sinθ, cosθ)
+    max_n = 1
+    max_val = params.threshold-1
+    for n in 0:params.numangle-1
+            dist = point[2]*cosθ[n+1] + point[1]*sinθ[n+1]
+            dist += params.constadd
+            dist = Int64(floor(dist))
+            params.accumulator_matrix[n+1 , dist + 1] += 1
+            val = params.accumulator_matrix[n+1 , dist + 1]
+            if(max_val < val)
+                max_val = val
+                max_n = n+1
+            end    
+    end
+    return max_n, max_val
+end    
+
+#function to detect the line segment after merging lines within lineGap
+function pass_1(img, params, sample, xflag)
+    line_end = [[0,0],[0,0]]
+    h, w = size(img)
+    for k = 1:2
+        gap = 0
+        x = sample.x0
+        y = sample.y0
+        dx = sample.dx0
+        dy = sample.dy0
+
+        if k>1
+            dx = -dx
+            dy = -dy
+        end
+        
+        while(true)
+            i1 = 0
+            j1 = 0
+            if(xflag==1)
+                j1 = x
+                i1 = y>>params.shift
+            else
+                j1 = x>>params.shift
+                i1 = y
+            end
+
+            # check when line exits image boundary
+            if( j1 < 0 || j1 >= w || i1 < 0 || i1 >= h )
+                break;
+            end
+            gap+=1
+
+            # if non-zero point found, continue the line
+            if(params.mask[i1+1, j1+1])
+                gap = 0
+                line_end[k][1] = i1+1
+                line_end[k][2] = j1+1
+             # if gap to this point was too large, end the line    
+            elseif(gap > params.lineGap)
+                break
+            end
+            x = Int64(x+dx)
+            y = Int64(y+dy)
+        end
+    end
+    return line_end
+end
+
+#function to reset the mask and accumulator_matrix 
+function pass_2(params, sample, xflag, good_line, line_end, sinθ, cosθ)
+    for k = 1:2
+        x = sample.x0
+        y = sample.y0
+        dx = sample.dx0
+        dy = sample.dy0
+
+        if k>1
+            dx = -dx
+            dy = -dy
+        end
+
+        # walk along the line using fixed-point arithmetics,
+        while(true)
+            i1, j1 = 0,0
+
+            if (xflag==1)
+                j1 = x
+                i1 = y >> params.shift
+            else
+                j1 = x >> params.shift
+                i1 = y
+            end
+
+            # if non-zero point found, continue the line
+            if(params.mask[i1+1, j1+1])
+                if(good_line)
+                    for n = 0:params.numangle-1
+                        r = ((j1+1)*cosθ[n+1] + (i1+1)*sinθ[n+1])
+                        r = Int64(floor(r+params.constadd))
+                        params.accumulator_matrix[n+1, r+1]-=1
+                        params.mask[i1+1, j1+1] = false
+                    end
+                end
+            end
+            # exit when the point is the line end                
+            if((i1+1) == line_end[k][1] && (j1+1) == line_end[k][2])
+                break
+            end
+            x = Int64(x+dx)
+            y = Int64(y+dy)              
+        end
+    end
+end    
+
+function hough_line_probabilistic(
+img::AbstractArray{T,2},
+ρ::Real, θ::Range,
+threshold::Integer, lineLength::Integer, lineGap::Integer, linesMax::Integer) where T<:Union{Bool,Gray{Bool}}
+  
+    ρ > 0 || throw(ArgumentError("Discrete step size must be positive"))
+    indsy, indsx = indices(img)
+    ρinv = 1 / ρ
+    numangle = length(θ)
+    numrho = round(Int,(2(length(indsx) + length(indsy)) + 1)*ρinv) 
+    constadd = (numrho-1)/2
+    accumulator_matrix = zeros(Int, numangle + 2, numrho + 2)
+    h, w = size(img)
+    mask = zeros(Bool, h, w)
+    #Pre-Computed sines and cosines in tables
+    sinθ, cosθ = sin.(θ).*ρinv, cos.(θ).*ρinv
+    nzloc = Vector{Tuple{Int64,Int64}}(0)
+    lines = Vector{Tuple{Int64, Int64, Int64, Int64}}(0)
+    params = Param(numangle, constadd, accumulator_matrix, mask, nzloc, lines, 16, threshold, lineLength, lineGap)
+    sample = Sample(0,0,0,0)
+
+    #collect non-zero image points
+    collect_points(params, img)
+
+    count_ = size(nzloc)[1]+1
+
+    # stage 2. process all the points in random order
+    while(count_>1)
+        count_-=1
+        good_line = false
+        # choose random point out of the remaining ones
+        idx = rand(1:count_)
+        max_n = 1
+        point = nzloc[idx]
+        i, j = point[1]-1, point[2]-1    
+        sample.x0, sample.y0, sample.dx0, sample.dy0, xflag = 0, 0, 0, 0, 0
+        max_n = 1
+
+        # "remove" it by overriding it with the last element
+        params.nzloc[idx] = params.nzloc[count_]
+
+        if(!(params.mask[point[1], point[2]]))
+            continue
+        end
+        
+        # update accumulator, find the most probable line
+        max_n, max_val = update_accumulator(params, point, sinθ, cosθ)
+
+        # if it is too "weak" candidate, continue with another point
+        if(max_val < params.threshold)
+            continue
+        end
+        
+        # from the current point walk in each direction along the found line
+        a = -sinθ[max_n]
+        b = cosθ[max_n]
+        sample.x0 = j
+        sample.y0 = i
+        good_line = false
+
+        if(abs(a) > abs(b))
+            xflag = 1
+            sample.dx0 = a > 0 ? 1 : -1
+            sample.dy0 = round(b*(1 << params.shift)/abs(a))
+            sample.y0 = (sample.y0 << params.shift) + (1 << (params.shift-1))
+        else
+            xflag = 0
+            sample.dy0 = b > 0 ? 1 : -1
+            sample.dx0 = round( a*(1 << params.shift)/abs(b) );
+            sample.x0 = (sample.x0 << params.shift) + (1 << (params.shift-1));    
+        end   
+
+        # pass 1: walk the line, merging lines less than specified gap length
+        line_end = pass_1(img, params, sample, xflag)
+
+        # confirm line length is sufficient
+        good_line = abs(line_end[2][1] - line_end[1][1]) >= lineLength || abs(line_end[2][2] - line_end[1][2]) >= lineLength              
+
+        # pass 2: walk the line again and reset accumulator and mask
+        pass_2(params, sample, xflag, good_line, line_end, sinθ, cosθ) 
+
+        # add line to the result
+        if(good_line)
+            push!(lines, (line_end[1][1], line_end[1][2], line_end[2][1], line_end[2][2]))
+
+            if(size(lines)[1] >= linesMax)
+                return lines
+            end
+        end         
+    end
+    return lines        
+end
+
 """
 ```
 circle_centers, circle_radius = hough_circle_gradient(img_edges, img_phase, scale, min_dist, vote_thres, min_radius:max_radius)  
diff --git a/test/houghtransform.jl b/test/houghtransform.jl
index a665471..8a4b2cb 100644
--- a/test/houghtransform.jl
+++ b/test/houghtransform.jl
@@ -37,6 +37,46 @@ using ImageFeatures
         @test er <= 0.1
     end
 
+    @testset "Hough Line Probabilistic" begin
+        #one horizontal line segment with few scattered points 
+        img = zeros(Bool,20,20)
+        for j in 5:15
+            img[15,j] = true
+        end
+
+        srand(1234)
+        for k in 1:5
+            img[rand(1:13),rand(6:14)] = true
+        end
+        lines = hough_line_probabilistic(img, 1, linspace(0,π,180),7,5,10,4)
+        @test length(lines) == 1
+        @test lines[1] == (15,5,15,15)
+
+        #for a square image
+        img = zeros(Bool, 20, 20)
+
+        for i in (3,14)
+            for j in 3:14
+                img[i,j] = true
+            end
+        end     
+
+        for i in (3,14)
+            for j in 3:14
+                img[j,i] = true
+            end
+        end
+        srand(1234)
+        for k in 1:4
+            img[rand(6:10),rand(6:10)] = true
+        end
+        lines = hough_line_probabilistic(img, 1, linspace(0,π,180),9,7,15,4)
+        lines_ = [line for line in lines]
+        @test length(lines) == 4
+        @test all(lines_ == [(3, 3, 3, 14), (14, 3, 14, 14), (13, 14, 4, 14), (13, 3, 4, 3)])
+        
+    end
+
     @testset "Hough Circle Gradient" begin
 
     dist(a, b) = sqrt(sum(abs2, (a-b).I))