Convert image coordinate system to correct coordinate system in R

I'm doing object detection of trees on a raster image in TIFF format. The output of my model is a series of bounding boxes that provide xmin, xmax, ymin, ymax for each bounding box for each tree in my image. These values are provided in the IMAGE coordinate system, which assumes that the top left corner is (0,0). However, my raster is in NAD83(CSRS) / UTM zone 18N (EPSG:2959).

How do I convert these values into the correct projection so that I can overlay them on my raster?

Here is the raster.

Here is the data that I get from my model. If you're curious, I'm using the deepforest model which you can find here: https://deepforest.readthedocs.io/en/latest/getting_started.html#training

Here is the code to get the csv data:

library(deepforestr)
library(terra)
library(sf)

model = df_model()

model$use_release()

raster_path = get_data("sample2.tif") # Gets a path to an example raster tile

bounding_boxes = model$predict_tile(raster_path, return_plot=TRUE)

The CSV data looks like this:

The raster looks like this:

First I create some example data. A SpatRaster x and data.frame with cell numbers d (instead of files and screenshots, please use code- generated data like this when asking questions).

library(terra)
x <- rast(ncols=1999, nrows=1249, nlyrs=3, xmin=331567.84, xmax=331887.68, 
     ymin=5078301.12, ymax=5078500.96, crs='EPSG:2959')

d <- data.frame(xmin=c(1521, 799), ymin=c(910, 1121), 
                xmax=c(1550, 832), ymax=c(876, 1087))

You can apply some simple math:

fun <- function(r, e) {
    cbind(
        xmin(r) + e[, c("xmin", "xmax")] * xres(r) + 0.5 * xres(r),
        ymax(r) - e[, c("ymin", "ymax")] * yres(r) - 0.5 * yres(r)
    )
}

crds <- fun(x, d)
crds
#      xmin     xmax    ymin    ymax
#1 331811.3 331815.9 5078355 5078361
#2 331695.8 331701.0 5078322 5078327

To create polygons, you could now do

p <- apply(crds, 1, \(x) as.polygons(ext(x))) |> vect(crs=crs(x))
p
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 2, 0  (geometries, attributes)
# extent      : 331695.8, 331815.9, 5078322, 5078361  (xmin, xmax, ymin, ymax)
# coord. ref. : NAD83(CSRS) / UTM zone 18N (EPSG:2959) 

And with that goal in mind, here is an alternative shorter route that first computes the cell numbers

cmn <- cellFromRowCol(x, d$ymax+1, d$xmin+1) 
cmx <- cellFromRowCol(x, d$ymin+1, d$xmax+1) 
v <- apply(cbind(cmn,cmx), 1, \(i) as.polygons(ext(x, i))) |> vect(crs=crs(x))
v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 2, 0  (geometries, attributes)
# extent      : 331695.7, 331816, 5078321, 5078361  (xmin, xmax, ymin, ymax)
# coord. ref. : NAD83(CSRS) / UTM zone 18N (EPSG:2959)