Highlight a section of a line between two points (Leaflet/R)

Question

Given a line and two points along that line, how can I highlight the section of the line between the two points?

The code below provides an example of a line (tr) and two points (sg) that fall on it, with the output (R/Leaflet) looking like this.

I have found examples of connecting points (e.g. here), but isn't really what I'm looking for: I need to draw/highlight the segment of the line between the two points, and not a straight line between the two points. The closest thing I could find is highlighting a section of a street (e.g. here), but I can't really apply that here (working with railway lines + I won't be able to use Google Maps).

library(tidyverse)
library(sf) 
library(leaflet)

tr <- structure(
  list(
    SHAPE_Length = 22210.9047689, 
    SHAPE = structure(list(
      structure(c(-3.447431, -3.426467, -3.3717497, -3.3504529, 
                  -3.34384, -3.3063404, -3.2875194, -3.2388496, -3.1727098, 
                  -3.136474, -3.1031447, 58.51414, 58.5140392, 58.5066418, 
                  58.5013098, 58.4950488, 58.4901998, 58.4837109, 58.4745371, 
                  58.4552422, 58.4512823, 58.4421872), .Dim = c(11L, 2L),
                class = c("XY", "LINESTRING", "sfg"))), 
      n_empty = 0L, crs = structure(list(
        input = "+init=epsg:4326", 
        wkt = "GEOGCRS[\"WGS 84\",\n    
        DATUM[\"World Geodetic System 1984\",\n        
        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            
        LENGTHUNIT[\"metre\",1]],\n        
        ID[\"EPSG\",6326]],\n    
        PRIMEM[\"Greenwich\",0,\n        
        ANGLEUNIT[\"degree\",0.0174532925199433],\n        
        ID[\"EPSG\",8901]],\n    
        CS[ellipsoidal,2],\n        
        AXIS[\"longitude\",east,\n            
        ORDER[1],\n            
        ANGLEUNIT[\"degree\",0.0174532925199433,\n                
        ID[\"EPSG\",9122]]],\n        
        AXIS[\"latitude\",north,\n            
        ORDER[2],\n            
        ANGLEUNIT[\"degree\",0.0174532925199433,\n                
        ID[\"EPSG\",9122]]],\n    
        USAGE[\n        
        SCOPE[\"unknown\"],\n        
        AREA[\"World.\"],\n        
        BBOX[-90,-180,90,180]]]"), 
        class = "crs"), 
      class = c("sfc_LINESTRING", "sfc"), 
      precision = 0, 
      bbox = structure(c(xmin = -3.447431, ymin = 58.4421872, 
                         xmax = -3.1031447, ymax = 58.51414), class = "bbox"))), 
  row.names = 1L, class = c("sf", "data.frame"), 
  sf_column = "SHAPE", 
  agr = structure(c(SHAPE_Length = NA_integer_), 
                  .Label = c("constant", "aggregate", "identity"), class = "factor"))


sg <- structure(list(
  geometry = structure(list(structure(c(-3.12213008366773, 58.4467460240336), 
                                      class = c("XY", "POINT", "sfg")), 
                            structure(c(-3.28773441769833, 58.483750111291), 
                                      class = c("XY", "POINT", "sfg"))), 
                       class = c("sfc_POINT", "sfc"), 
                       precision = 0, 
                       bbox = structure(c(xmin = -3.28773441769833, ymin = 58.4467460240336, 
                                          xmax = -3.12213008366773, ymax = 58.483750111291), 
                                        class = "bbox"), 
                       crs = structure(list(input = "+init=epsg:4326",
                                            wkt = "GEOGCRS[\"WGS 84\",\n    
                                            DATUM[\"World Geodetic System 1984\",\n        
                                            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            
                                            LENGTHUNIT[\"metre\",1]],\n        
                                            ID[\"EPSG\",6326]],\n    
                                            PRIMEM[\"Greenwich\",0,\n        
                                            ANGLEUNIT[\"degree\",0.0174532925199433],\n        
                                            ID[\"EPSG\",8901]],\n    
                                            CS[ellipsoidal,2],\n        
                                            AXIS[\"longitude\",east,\n            
                                            ORDER[1],\n            
                                            ANGLEUNIT[\"degree\",0.0174532925199433,\n                
                                            ID[\"EPSG\",9122]]],\n        
                                            AXIS[\"latitude\",north,\n            
                                            ORDER[2],\n            
                                            ANGLEUNIT[\"degree\",0.0174532925199433,\n                
                                            ID[\"EPSG\",9122]]],\n    
                                            USAGE[\n        
                                            SCOPE[\"unknown\"],\n        
                                            AREA[\"World.\"],\n        
                                            BBOX[-90,-180,90,180]]]"), 
                                       class = "crs"), n_empty = 0L)), 
  row.names = 1:2, class = c("sf","data.frame"), 
  sf_column = "geometry", 
  agr = structure(integer(0), .Label = c("constant","aggregate", "identity"), 
                  class = "factor", .Names = character(0)))


# Plot the stuff
leaflet() %>%
  addTiles() %>%
  addPolylines(data = tr) %>%
  addMarkers(data = sg)

sessionInfo():

R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
tidyverse_1.3.1
sf_1.0-3
leaflet_2.0.4.1