relatedcircumscribingcircle

Project.toml

@@ -5,12 +5,14 @@ version = "0.0.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 
 [compat]
 DelimitedFiles = "1.9"
+Random = "1.11.0"
 StatsBase = "0.34"
 Test = "1"
 TestItems = "1"

src/LandscapeMetrics.jl

@@ -4,6 +4,7 @@ using TestItems
 using Test
 using StatsBase
 import DelimitedFiles
+using Random
 
 # Patches, generic utilies, and overloads
 include("types.jl")
@@ -26,6 +27,22 @@ export patches!, patches
 include("utilities/cellcenters.jl")
 export cellcenters
 
+#Function to get the contiguity value of the cells
+include("utilities/contiguityvalue.jl")
+export contiguityvalue
+# Function to identify the minimum enclosing circle of a set of points
+include("utilities/welzlalgorithm.jl")
+export welzl_algorithm
+export minimum_enclosing_circle
+export Circle
+export is_point_in_circle
+export distance
+
+
+
+
+
+
 # Function to get the core areas
 include("utilities/core.jl")
 
@@ -47,13 +64,6 @@ export largestpatchindex
 include("area_and_edge/radiusofgyration.jl")
 export radiusofgyration
 
-include("area_and_edge/edgedensity.jl")
-export edgedensity
-
-include("area_and_edge/totaledge.jl")
-export totaledge
-
-
 # Shape
 include("shape/paratio.jl")
 export paratio, perimeterarearatio
@@ -62,5 +72,13 @@ export shapeindex
 include("shape/fractal.jl")
 export fractaldimensionindex
 
+include("shape/contiguityindex.jl")
+export contig_index
+include("shape/relatedcircumscribingcircle.jl")
+export related_circumscribing_circle
+export patch_cell_corners
+
+
+
 end # module LandscapeMetrics
 

src/shape/contiguityindex.jl

@@ -0,0 +1,40 @@
+function contig_index(l::Landscape, patch, template)
+
+    # Create a mask for the given patch
+    patch_mask = patches(l) .== patch
+
+    # Keep all values in the patch, set others to zero
+    patch_grid = zeros(Int, size(l))
+
+    # Fill in the values for the patch
+    patch_grid[patch_mask] .= l.grid[patch_mask]
+
+    # Get the contiguity values for the given template
+    contig_vals = contiguityvalue(patch_grid, template)
+
+    # Get the cells that belong to the patch
+    patch_cells = findall(patches(l) .== patch)
+
+    # Calculate the contiguity index
+    sum_vals = sum(contig_vals[patch_cells])
+    n_cells = length(patch_cells)
+    avg_val = sum_vals / n_cells
+    template_sum = sum(template)
+
+    return (avg_val - 1) / (template_sum)
+end
+
+@testitem "We can compute the contiguity index of a patch" begin
+    A = [
+        1 1;
+        1 1;
+        1 3
+    ]
+
+    B = [1 2 1;
+         2 1 2;
+         1 2 1]
+    L = Landscape(A)
+    @test contig_index(L, 1, B) == 0.40
+end
+

src/shape/relatedcircumscribingcircle.jl

@@ -0,0 +1,59 @@
+"""
+    patch_cell_corners(l::Landscape, patch)
+
+Returns the coordinates of the corners of all cells in the specified patch.
+
+"""
+function patch_cell_corners(l::Landscape, patch)
+
+    # We get the length of the side of one cell
+    s = side(l)
+
+    # We get the cell centers
+    centers = cellcenters(l)
+
+    # We find all cells in the patch
+    patch_cells = findall(patches(l) .== patch)
+
+    # We compute the corners of each cell in the patch
+    corners = Tuple{Float64,Float64}[]
+
+    # Loop through each cell in the patch
+    for idx in patch_cells
+        c = centers[idx]
+        push!(corners, (c[1] - s/2, c[2] - s/2))
+        push!(corners, (c[1] - s/2, c[2] + s/2))
+        push!(corners, (c[1] + s/2, c[2] - s/2))
+        push!(corners, (c[1] + s/2, c[2] + s/2))
+    end
+
+    return unique(corners)
+end
+
+"""
+    related_circumscribing_circle(l::Landscape, patch)
+
+Returns the related circumscribing circle metric for the specified patch.
+"""
+function related_circumscribing_circle(l::Landscape, patch)
+
+    # We get the area of the patch
+    patch_area = area(l, patch)
+
+    # We get the corners of the cells in the patch
+    patch_points = patch_cell_corners(l, patch)
+
+    # We compute the minimum enclosing circle
+    circle = minimum_enclosing_circle(copy(patch_points))
+
+    # We compute the area of the circle
+    circle_area = π * circle.radius^2
+
+    # We return the related circumscribing circle metric
+    return 1 - (patch_area / circle_area)
+end
+
+
+
+
+

src/utilities/contiguityvalue.jl

@@ -0,0 +1,40 @@
+B = [
+    1 2 1;
+    2 1 2;
+    1 2 1
+]
+
+function contiguityvalue(A::AbstractMatrix, B::AbstractMatrix)
+    nrows, ncols = size(A)
+    T = eltype(A)
+    result = zeros(T, nrows, ncols)
+    # Pad the matrix with zeros on all sides
+    padded = zeros(T, nrows+2, ncols+2)
+    padded[2:end-1, 2:end-1] .= A
+    for i in 1:nrows
+        for j in 1:ncols
+            neighborhood = padded[i:i+2, j:j+2]
+            result[i, j] = sum(neighborhood .* B)
+        end
+    end
+    return result
+end
+
+@testitem "We can calculate the contiguity value of a landscape" begin
+    A = [
+        1 1;
+        1 1;
+        1 3
+    ]
+    B = [
+    1 2 1;
+    2 1 2;
+    1 2 1
+]
+    expected = [
+        6 6;
+        11 13;
+        10 8
+    ]
+    @test contiguityvalue(A, B) == expected
+end

src/utilities/welzlalgorithm.jl

@@ -0,0 +1,125 @@
+using Random
+
+# Structure to represent a circle
+struct Circle
+    center::Tuple{Float64, Float64}
+    radius::Float64
+end
+
+
+"""
+    distance(p1::Tuple{Float64, Float64}, p2::Tuple{Float64, Float64})
+
+    Computes the Euclidian distance between two points.
+
+"""
+distance(p1::Tuple{Float64, Float64}, p2::Tuple{Float64, Float64}) =
+    sqrt((p1[1] - p2[1])^2 + (p1[2] - p2[2])^2)
+
+"""
+    is_point_in_circle(p::Tuple{Float64, Float64}, c::Circle)
+
+    Checks if point `p` is inside or on the boundary of circle `c`.
+
+"""
+is_point_in_circle(p::Tuple{Float64, Float64}, c::Circle) =
+    distance(p, c.center) <= c.radius + 1e-9
+
+"""
+
+    circumcircle(p1, p2, p3)
+
+    Computes the circumcircle of the triangle formed by points `p1`, `p2`, and `p3`.
+
+"""
+function circumcircle(p1, p2, p3)
+    A = 2 * (p2[1] - p1[1])
+    B = 2 * (p2[2] - p1[2])
+    C = 2 * (p3[1] - p1[1])
+    D = 2 * (p3[2] - p1[2])
+    E = p2[1]^2 + p2[2]^2 - p1[1]^2 - p1[2]^2
+    F = p3[1]^2 + p3[2]^2 - p1[1]^2 - p1[2]^2
+    denom = A * D - B * C
+    if abs(denom) < 1e-12
+        # Points are colinear; circle encloses all three points
+        center = ((p1[1] + p3[1]) / 2, (p1[2] + p3[2]) / 2)
+        radius = maximum([distance(center, p) for p in (p1, p2, p3)])
+        return Circle(center, radius)
+    end
+    cx = (E * D - B * F) / denom
+    cy = (A * F - E * C) / denom
+    center = (cx, cy)
+    radius = distance(center, p1)
+    return Circle(center, radius)
+end
+
+"""
+    trivial_circle(P)
+
+    Computes the minimum enclosing circle for 0, 1, 2, or 3 points in `P`.
+"""
+function trivial_circle(P)
+    if length(P) == 0
+        return Circle((0.0, 0.0), 0.0)
+    elseif length(P) == 1
+        return Circle(P[1], 0.0)
+    elseif length(P) == 2
+        center = ((P[1][1] + P[2][1]) / 2, (P[1][2] + P[2][2]) / 2)
+        radius = distance(P[1], P[2]) / 2
+        return Circle(center, radius)
+    elseif length(P) == 3
+        return circumcircle(P[1], P[2], P[3])
+    end
+end
+
+"""
+    welzl_helper(P, R, n)
+
+    Recursive helper function for Welzl's algorithm.
+    `P` is the list of points, `R` is the list of points on the boundary,
+    and `n` is the number of points in `P` to consider.
+"""
+function welzl_helper(P::Vector{Tuple{Float64, Float64}}, R::Vector{Tuple{Float64, Float64}}, n::Int)
+    if n == 0 || length(R) == 3
+        return trivial_circle(R)
+    end
+
+    p = P[n]
+    D = welzl_helper(P, R, n - 1)
+
+    if is_point_in_circle(p, D)
+        return D
+    else
+        return welzl_helper(P, vcat(R, [p]), n - 1)
+    end
+end
+
+"""
+    minimum_enclosing_circle(points)
+
+    Computes the minimum enclosing circle for a set of points using Welzl's algorithm.
+"""
+function minimum_enclosing_circle(points::Vector{Tuple{Float64, Float64}})
+
+    # Shuffle the points to ensure random order
+    Pshuffled = shuffle(copy(points))
+
+    # Calls the recursive helper function
+    return welzl_helper(Pshuffled, Tuple{Float64, Float64}[], length(Pshuffled))
+end
+
+
+@testitem "MEC of one point is a circle of radius 0 at that point" begin
+    p = (1.0, 2.0)
+    mec = minimum_enclosing_circle([p])
+    @test mec.radius == 0.0
+    @test mec.center == p
+end
+
+@testitem "MEC of two points is a circle with diameter the segment between them" begin
+    p1 = (0.0, 0.0)
+    p2 = (2.0, 0.0)
+    mec = minimum_enclosing_circle([p1, p2])
+    @test mec.radius == 1.0
+    @test mec.center == (1.0, 0.0)
+end