How to identify intersections of hollow geometries while ignoring others?

I'm working with spatial data in R using the sf package, and I need to identify geometries that are completely within a larger geometry (not just intersecting or touching boundaries).

In other words, I need the red square (which has a hollow center) to merge with the purple, while leaving the blue unaffected.

Does anyone know how to do this?

Here is an data example:

# Geometry 1 (geom1) is a polygon with a hole (an outer boundary and an inner "empty" ring).
# Geometry 2 (geom2) is a simple polygon that aligns perfectly with the inner boundary of geom1 (it is completely inside the "hole").
# Geometry 3 (geom3) is a separate polygon that intersects the boundary of geom1 but is not fully contained by it.

library(sf)

# Create geometry 1 as a polygon with a hole
geom1 <- st_polygon(list(
  rbind(c(0, 0), c(0, 5), c(5, 5), c(5, 0), c(0, 0)),           # Outer boundary
  rbind(c(1, 1), c(1, 4), c(4, 4), c(4, 1), c(1, 1))            # Inner "hole"
))

# Create geometry 2 as a simple polygon
geom2 <- st_polygon(list(
  rbind(c(1, 1), c(1, 4), c(4, 4), c(4, 1), c(1, 1))            # Polygon inside the hole
))

# Create geometry 3 as a separate polygon
geom3 <- st_polygon(list(
  rbind(c(5, 5), c(5, 9), c(9, 9), c(9, 5), c(5, 5))            # Polygon outside geom1
))

# Combine geometries into a simple feature collection
sf_data <- st_sf(geometry = st_sfc(geom1, geom2, geom3))

Note that st_intersects returns both geometry 2 (internal) and geometry 3. In this case, it would be necessary to return only geometry 2. st_contains also not work since geometry 1 has a blank space.

I was expecting to find a method to identify geometries that are fully within the outer boundary of another geometry.

Included example suggests that you could drop holes and operate with exterior rings to check what gets covered (by what).
With sfheaders::sf_remove_holes():

library(sf)

(exterior_rings <- sfheaders::sf_remove_holes(sf_data))
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 9 ymax: 9
#> CRS:           NA
#>                         geometry
#> 1 POLYGON ((0 0, 0 5, 5 5, 5 ...
#> 2 POLYGON ((1 1, 1 4, 4 4, 4 ...
#> 3 POLYGON ((5 5, 5 9, 9 9, 9 ...

st_covered_by(exterior_rings)
#> Sparse geometry binary predicate list of length 3, where the predicate
#> was `covered_by'
#>  1: 1
#>  2: 1, 2
#>  3: 3

which(lengths(st_covered_by(exterior_rings)) > 1)
#> [1] 2