Skip to content

Instantly share code, notes, and snippets.

@perusio

perusio/no_geom_MBR.jl

Forked from caseykneale/no_geom_MBR.jl
Created July 13, 2022 13:02
Show Gist options
  • Save perusio/d0737e94e7d5d41c94d0fde9c224538b to your computer and use it in GitHub Desktop.
Save perusio/d0737e94e7d5d41c94d0fde9c224538b to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. @caseykneale caseykneale created this gist Aug 1, 2020.
    69 changes: 69 additions & 0 deletions no_geom_MBR.jl
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,69 @@
    using GeometryBasics

    x(a::AbstractVector) = first(a)
    y(a::AbstractVector) = last(a)
    x(a::AbstractMatrix) = a[ :, 1 ]
    y(a::AbstractMatrix) = a[ :, 2 ]

    function orientation(p::AbstractVector, q::AbstractVector, r::AbstractVector)::Int
    val = ( y(q) - y(p) ) * ( x(r) - x(q) ) - ( x(q) - x(p) ) * ( y(r) - y(q) )
    return (val ≈ 0) ? 0 : ( (val > 0) ? 1 : 2 )
    end

    """
    Classic Jarvis/Gift Wrapping Convex Hull
    """
    function convexhull( points::Matrix )::Vector
    n = size(points, 1)
    @assert (n > 2) "Convex Hull requires at least 3 points."
    hull = []
    p, q = argmin( x( points ) ), 0
    init = p
    while ( p != init ) || ( length( hull ) == 0 )
    push!( hull, p )
    q = (( p + 1) % n ) + 1
    for i in 1:n
    if orientation(points[p,:], points[i,:], points[q,:]) == 2
    q = i
    end
    end
    p = q;
    end
    return hull
    end

    function minimum_bounding_rectangle( points, hull_idxs )
    n = length(hull_idxs)
    # calculate edge angles
    edges = points[ hull_idxs[ 2:end ] ] - points[ hull_idxs[ 1:(end-1) ] ]
    angles = unique( abs.( atan.( y( edges ), x( edges ) ) .% ( pi/2 ) ) )
    # find rotation matrices
    rot_matrices = zeros( 2, 2, length( angles ) )
    rot_matrices[1,1,:] = cos.(angles); rot_matrices[2,1,:] = sin.(angles )
    rot_matrices[1,2,:] = -rot_matrices[2,1,:]; rot_matrices[2,2,:] = rot_matrices[1,1,:]
    const_view = points[ hull_idxs,: ]
    rot_points = [ rot_matrices[:,:,z] * const_view' for z in 1:size( rot_matrices, 3 ) ]
    # find the bounding points
    min_x, max_x, min_y, max_y = Vector{eltype(points)}(undef, length(rot_points)),
    Vector{eltype(points)}(undef, length(rot_points)),
    Vector{eltype(points)}(undef, length(rot_points)),
    Vector{eltype(points)}(undef, length(rot_points))
    for (n, rp) in enumerate( rot_points )
    min_x[n], max_x[n] = extrema( rp[1,:] )
    min_y[n], max_y[n] = extrema( rp[2,:] )
    end
    # find the box with the best area
    smallest_box = argmin( (max_x .- min_x) .* (max_y .- min_y) )
    x1,x2 = min_x[smallest_box], max_x[smallest_box]
    y1,y2 = min_y[smallest_box], max_y[smallest_box]
    return ( rot_matrices[:,:,smallest_box]' * [ x1 y1; x1 y2; x2 y2; x2 y1 ]')'
    end

    p1 = randn(1000,2)
    hull_inds = convexhull( p1 )
    mbr = minimum_bounding_rectangle( p1, hull_inds )

    using Plots
    scatter( x(p1), y(p1), color = "black", legend = false);
    scatter!( x(p1[hull_inds,:]), y(p1[hull_inds,:]), color = "purple");
    scatter!( x(mbr), y(mbr), color = "pink")