Weighted Voronoi polygons in R

Eventually I managed to solve my problem, building from the excellent answer https://gis.stackexchange.com/a/17290/103010 from user "whuber" in the already cited Create weighted Thiessen polygons?.

This is meant to be a simple summary of the code, not of the science behind the weighted Voronoi polygons:

1. Calculate the standard Euclidean distance of every pixel in a raster to each point.
2. Apply a weighting function that change the distance value: in the example provided the value increases exponentially with the distance but decreases with higher size of the points (this works for my specific application). Bigger trees will have larger areas of influence but will not override completely smaller ones.
3. Identify which tree claims influence on a pixel: the one with the lower weighted distance.
4. With this info, build a raster, polygon or simply a data.table for your need

    # necessary packages
    library(data.table)
    library(raster)
    library(rgeos)

    # create base raster
    rez= 1
    r <- raster(xmn=0, xmx=40, ymn=0, ymx=40, 
                crs = '+proj=utm +zone=01 +datum=WGS84', resolution=rez)
    
    # create example data
    treez <- data.table(treeID=c("a","b","c","d"),
                        x=c(5,35,25,10),
                        y=c(35,35,15,10),
                        size=c(5,20,10,1))
    
    # create an empty datatable with the length equal to the raster size
    dt <- data.table(V1=seq(1,length.out = 40^2*(1/rez)^2))

    # calculate distance grid for each tree, apply weigth and store the data
    for (i in treez[,unique(treeID)]){
      # 1) standard Euclidean distance
      d <- distanceFromPoints(r, treez[treeID==i,.(x,y)])
      # 2) apply your own custom weight
      d <- ((d@data@values)^2)/treez[treeID==i,size]
      # 3) store the results with the correct name
      dt <- cbind(dt,d)
      setnames(dt, "d", i)
    }
    
    # carry out comparison
    # find the column with the lowest value (+1 to account to column V1 at the beginning)
    dt[,minval:=which.min(.SD)+1,by=V1,.SDcols=2:length(names(dt))]
    # find the name of that column, that is the id of the point
    dt[,minlab:=names(dt)[minval]]
    
    # create resulting raster
    out <- raster(xmn=0, xmx=40, ymn=0, ymx=40, 
                crs = '+proj=utm +zone=01 +datum=WGS84', resolution=rez,
                vals=dt[,as.factor(minlab)])
    
    # Transform into polygons
    outpoly <- rasterToPolygons(out, dissolve=T)
    # Plot results with starting points
    plot(outpoly);points(treez[,.(x,y)],col="red",pch=16, cex=sqrt(treez[,size]))