Isovists should be computed even with fully-occluded rays

NEWS.md

@@ -9,6 +9,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added Zenodo and CRAN badges #15
 
+## Fixed
+
+- Isovist can now be computed also if some rays are fully occluded #17
+
 # Version 0.1.0 - 2025-04-03
 
 ## Added

R/isovist.R

@@ -159,19 +159,29 @@ occlude_rays <- function(rays, occluders = NULL) {
   # Filter the only features intersecting each others
   rays_intersect <- ray_geoms[ray_intersects_occluders]
   occluders_intersect <- occluder_geoms[occluder_intersects_rays]
+  occluders_intersect <- sf::st_union(occluders_intersect)
+
+  # Determine which rays are fully within the occluders and will thus be dropped
+  rays_within_occluders <- sf::st_within(rays_intersect,
+                                         occluders_intersect,
+                                         sparse = FALSE)
 
   # Determine the ray segments that do not overlap with the occluders. For each
-  # ray the segments form (multi)-linestrings
-  diffs <- sf::st_difference(rays_intersect, sf::st_union(occluders_intersect))
+  # ray the segments form (multi)-linestrings. Fully-occluded rays are dropped
+  diffs <- sf::st_difference(rays_intersect, occluders_intersect)
 
   # By casting `diffs` to linestrings, we keep only the first segment of each
   # ray, which, by construction, runs from the original viewpoint to the first
   # occluder (if any)
   rays_occluded <- sf::st_cast(diffs, "LINESTRING", group_or_split = FALSE) |>
     suppressWarnings()
 
+  # Determine the indices of the partially-occluded rays, i.e. the ones that
+  # intersect the occluders and are not fully within the occluders
+  idx <- which(ray_intersects_occluders)[!rays_within_occluders]
+
   # Update the occluded ray geometries and return the sf(c) object
-  ray_geoms[ray_intersects_occluders] <- rays_occluded
+  ray_geoms[idx] <- rays_occluded
   sf::st_set_geometry(rays, ray_geoms)
 }
 

tests/testthat/test-isovist.R

@@ -186,3 +186,14 @@ test_that("merge_isovists returns the merged polygon, with or without holes", {
   expected <- sf::st_union(isovists)
   expect_equal(actual, expected)
 })
+
+test_that("visor works with fully occluded rays", {
+  viewpoints <- sf::st_sfc(sf::st_point(c(0, 0)), sf::st_point(c(2, 0)))
+  occluders <- sf::st_sfc(
+    create_occluder(1, 0, 2, 1),
+    create_occluder(3, 1, 1, 1)
+  )
+  expect_no_error(
+    get_isovist(viewpoints, occluders = occluders, ray_length = 1)
+  )
+})