From 93d8049811a921c4600d644b2203c75a8345682e Mon Sep 17 00:00:00 2001
From: tejus-gupta <tejusgupta14@gmail.com>
Date: Tue, 13 Nov 2018 21:14:37 +0530
Subject: [PATCH] Added method for constructing RAG from image (#31)

---
 src/core.jl  | 59 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 test/core.jl | 23 ++++++++++++++++++++
 2 files changed, 82 insertions(+)

diff --git a/src/core.jl b/src/core.jl
index ecf7158..53b54d6 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -264,3 +264,62 @@ end
 # Once Base has colon defined here we can replace this
 _colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
     map((i,j) -> i:j, Tuple(I), Tuple(J))
+
+
+"""
+    G, vertex2cartesian = region_adjacency_graph(img, weight_fn, R)
+
+Constructs a region adjacency graph (RAG) from a N-D image `img`. It returns the RAG along
+with a mapping from vertex index in RAG to cartesian index in the image.
+
+`weight_fn` is used to assign weights to the edges using pixel similarity and spatial proximity,
+where higher weight means greater similarity and thus stronger association. Zero weight is assigned
+to edges between any pair of nodes that are more than `R` pixels apart. `R` can be specified
+as a N-dimensional `CartesianIndex`. Alternatively, `R` can be an integer, in which a 
+N-dimensional `CartesianIndex` with value `R` along each dimension is used. `weight_fn` should have
+signature - 
+
+    edge_weight = weight_fn(p1::Pair{CartesianIndex{N},T}, p2::Pair{CartesianIndex{N},T}) where {N,T}
+
+Any graph clustering technique can be used with the constructed RAG to segment the image.
+
+# Example
+
+```julia
+    julia> using ImageSegmentation, SimpleWeightedGraphs
+    julia> img = fill(1.0, (10,10))
+    julia> img[4:6, 4:6] .= 0
+
+    julia> weight_fn(I, J) = 1-abs(I.second - J.second)
+    julia> G, vertex2cartesian = region_adjacency_graph(img, weight_fn, 1)
+```
+
+"""
+
+function region_adjacency_graph(img::AbstractArray{CT,N}, weight_fn::Function, R::CartesianIndex{N}) where {CT<:Union{Colorant,Real}, N}
+    cartesian2vertex = LinearIndices(img)
+    vertex2cartesian = CartesianIndices(img)
+
+    sources = Vector{Int}()
+    destinations = Vector{Int}()
+    weights = Vector{Float64}()
+
+    indices = CartesianIndices(axes(img));
+    Istart, Iend = first(indices), last(indices)
+    for I in indices
+        for J in CartesianIndices(map((i,j)->i:j, Tuple(max(Istart, I-R)), Tuple(min(Iend, I+R))))
+            if I <= J 
+                continue
+            end
+            push!(sources, cartesian2vertex[I])
+            push!(destinations, cartesian2vertex[J])
+            push!(weights, weight_fn(I=>img[I], J=>img[J]))
+        end
+    end
+
+    G = SimpleWeightedGraph(sources, destinations, weights)
+
+    return G, vertex2cartesian
+end
+
+region_adjacency_graph(img::AbstractArray{CT,N}, weight_fn::Function, R::Int) where {CT<:Union{Colorant,Real}, N} = region_adjacency_graph(img, weight_fn, R * CartesianIndex{N}())
diff --git a/test/core.jl b/test/core.jl
index 549a81f..c0d62b7 100644
--- a/test/core.jl
+++ b/test/core.jl
@@ -117,4 +117,27 @@
   @test all(label->(expected_means[label] ≈ new_seg.segment_means[label]), new_seg.segment_labels)
   @test new_seg.image_indexmap == expected
 
+  # region_adjacency_graph
+  img = fill(1.0, (10,10))
+  img[4:6, 4:6] .= 0
+
+  weight_fn(I, J) = exp(-abs(I.second - J.second)) / ((I.first[1]-J.first[1])^2 + (I.first[2]-J.first[2])^2)
+  G, vertex2cartesian = region_adjacency_graph(img, weight_fn, 2)
+
+  @test nv(G) == 100
+  @test G.weights[45,46] == 1
+  @test G.weights[43,44] == exp(-1)
+
+  #test that all edges weights equal weights computed using weight_fn
+  @test all(map(edge-> edge.weight == weight_fn(vertex2cartesian[edge.src]=>img[vertex2cartesian[edge.src]], vertex2cartesian[edge.dst]=>img[vertex2cartesian[edge.dst]]), edges(G)))
+
+  #test that all edges are between nodes less than r pixels away
+  @test all(map(edge-> max(abs(vertex2cartesian[edge.src][1] - vertex2cartesian[edge.dst][1]), abs(vertex2cartesian[edge.src][2] - vertex2cartesian[edge.dst][2])) <= 2 , edges(G)))
+
+  @test region_adjacency_graph(img, weight_fn, 2) == region_adjacency_graph(img, weight_fn, CartesianIndex(2,2))
+
+  G, vertex2cartesian = region_adjacency_graph(img, weight_fn, CartesianIndex(1,2))
+  @test G.weights[45,47] == 0
+  @test G.weights[45,65] != 0
+
 end