applyTransformation to matrix

Question

applyTransformation to matrix

Silviculturalist opened this issue 3 years ago · 12 comments

I'm using the NiftyReg linear algorithm to find the scaling, translation and rotation (yaw) of a set of coordinate points (trees on inventory plots) between two times. I've been told to ignore shear. I first create an image representation using a gaussian function with the amplitude= a coordinate attribute (diameter), and the standard deviation = positional error (e.g. 0.5 m) with the pixel resolution 0.1x0.1 m :

2015 (source):

2019 (target):

Suggested transformation:

This looks good.

However, when I then try to apply the transformation to the (local) x, y coordinates of the trees, I get quite odd results - for example, including the suggested translation (suggestion*resolution for real unit?) makes for a significantly poorer result:

Target:

Result with only rotation (yaw) - quite good:

Result with (in order - scaling, yaw, translation of points):

Do you have any ideas of what can be causing this?

Answer 1 · 2022-03-02T14:19:37.000Z

Hi! It might be easier to work out what's going on if you could outline the code behind these operations, but are you essentially trying to reconstruct and apply the transform with the skew clamped to zero? If so, the problem may be that all the component transforms interrelate, so by discarding the skew/shear you leave the final transformation inconsistent with the estimated one. Ideally the estimation would be done under the same constraint, but unfortunately NiftyReg only allows for full affine and rigid (rotation and translation) transforms. If I'm misunderstanding perhaps you could clarify a bit further.

Answer 2 · 2022-03-02T15:11:36.000Z

Thanks Jon for your quick reply. That makes sense! I did have some questions when the list results from decomposeAffine were not readable directly from the affine matrix.

I have separate S3 classes and methods for handling inventory plots with local coordinates - but here's as minimal an example as I could give.

Named stand100
Data: (48 trees, measured 2015 and again 2019):
structure(list(TreeID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L), x = c(0.74, 0.4, -1.85, -6.1, -7.22, -2.85, -2.18, -3.15, -6.39, -6.21, -2.76, 0.26, 1.17, 3.61, 2.86, 3.51, 4.66, 2.15, 8, 8.98, -1.27, -3.12, -6.86, -7.51, -2.74, -1.95, -2.9, -10.11, -9.72, -5.63, -5.5, -1.43, -3.55, 1.06, 2.46, 2.33, 1.57, 1.28, 2.8, 4.03, 4.74, 5.03, 6.08, 8.54, 8.38, 2.86, 0.32, 2.28), y = c(1.41, 9.48, 8.05, 5.54, 2.88, -0.15, -1.1, -1.28, -4.47, -4.55, -8.55, -6.01, -7.43, -8.8, -7.78, -7.42, -8.57, -4.52, -4.46, 2.18, 9.42, 7.61, 4.46, 1.63, -0.56, -1.31, -1.59, -3.05, -3.83, -5.25, -5.19, -8.77, -10.92, -11.1, -9.67, -7.06, -6.24, -5.76, -4.02, -7.08, -6.57, -7.94, -7.53, -3.06, 3.71, 2.51, 1.43, 9.93), Species = structure(c(2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("10", "11", "12", "20", "21", "30", "40", "41", "70", "71", "81", "92", "94", "95"), class = "factor"), Diameter = c(12.6, 14.4, 14, 6.7, 12.8, 19.7, 22.4, 10.4, 4.4, 5.8, 8.1, 17.4, 8.2, 11.4, 5.2, 9.5, 11, 13.7, 9.2, 23.4, 15.55, 14.95, 7.1, 11.45, 21, 23.8, 9.75, 8.8, 16.65, 4.4, 4.8, 8.4, 18.1, 4.45, 7.2, 8, 3.6, 17.05, 14.2, 4.2, 9.3, 12.25, 10.95, 9.95, 24.05, 4.65, 13.95, 22.5), Year = c("2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019")), row.names = c(NA, -48L), class = "data.frame")

I at first attempted to apply the entire transformation at once using the RNiftyReg::applyTransform() function. But since it did not like me manually changing the values for the translation (I was probably doing it wrong), I had to resort to writing a function that would pick out the important parts from the decomposed affine and try to apply it to the data.frame instead - giving the above. But here's how I tried to use the built-in matrix method:

Transformation matrix Named stand100_2015_2019_affine
NiftyReg affine matrix: 0.9882789 -0.1661159 0.0000000 8.7488136 0.1656586 0.9933031 0.0000000 -8.2611609 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000

Here's how I at first applied the entire matrix using the built in function:

stand100[,c("x","y")] <- applyTransform(stand100_2015_2019_affine,as.matrix(stand100[,c("x","y")]))

It is thrown off centre (probably because of unit issues between translation ? ):

Answer 3 · 2022-03-02T15:14:21.000Z

(Essentially what I want to do is conduct the same transformation on the coordinates as the representative images - so that the trees which would be the same between the years are placed 'above' eachother). applyTransform does this to the coordinates, but the translation scale is off. Can I set some type of a scale to the image?

Answer 4 · 2022-03-04T17:11:47.000Z

A key issue here is that the coordinate systems aren't the same. I suspect there are more subtle effects too, because of the skew consideration, but first and foremost, RNiftyReg will assume point locations are like matrix indices – starting in the corner of the image and increasing from one – whereas your coordinates don't have a range of about 100 (like the image dimensions) and go negative and so have some kind of central origin.

If you can easily convert them back to index-type coordinates that's probably quickest, but otherwise you can tell RNiftyReg about real-world resolution and origin centres through attributes – I can help with that if necessary.

Answer 5 · 2022-03-15T12:38:42.000Z

EDIT: Rotation is working well after realising I mistakenly asked it to plot the transformed 2019 values rather than the source 2015. There is a translation issue between the transformed and the target. I apologise if this has created any issue.

Hi Jon, I've been working on transforming the local coordinates to some index-type, but have gotten stuck. It seems like the translation (?) is off? The images look like they work great, so something is wrong with my writing... are the cell coordinates as handled by R/RNiftyReg in cell centre..? Looks also as if rotation is in wrong direction?

I've tried to make it as easy for you to see what I've done as possible. I'm very thankful for any help - been stuck on this for quite some time now.

Best,
Carl

# Some data to use
stand100 <- structure(list(TreeID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L), x = c(0.74, 0.4, -1.85, -6.1, -7.22, -2.85, -2.18, -3.15, -6.39, -6.21, -2.76, 0.26, 1.17, 3.61, 2.86, 3.51, 4.66, 2.15, 8, 8.98, -1.27, -3.12, -6.86, -7.51, -2.74, -1.95, -2.9, -10.11, -9.72, -5.63, -5.5, -1.43, -3.55, 1.06, 2.46, 2.33, 1.57, 1.28, 2.8, 4.03, 4.74, 5.03, 6.08, 8.54, 8.38, 2.86, 0.32, 2.28), y = c(1.41, 9.48, 8.05, 5.54, 2.88, -0.15, -1.1, -1.28, -4.47, -4.55, -8.55, -6.01, -7.43, -8.8, -7.78, -7.42, -8.57, -4.52, -4.46, 2.18, 9.42, 7.61, 4.46, 1.63, -0.56, -1.31, -1.59, -3.05, -3.83, -5.25, -5.19, -8.77, -10.92, -11.1, -9.67, -7.06, -6.24, -5.76, -4.02, -7.08, -6.57, -7.94, -7.53, -3.06, 3.71, 2.51, 1.43, 9.93), Species = structure(c(2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("10", "11", "12", "20", "21", "30", "40", "41", "70", "71", "81", "92", "94", "95"), class = "factor"), Diameter = c(12.6, 14.4, 14, 6.7, 12.8, 19.7, 22.4, 10.4, 4.4, 5.8, 8.1, 17.4, 8.2, 11.4, 5.2, 9.5, 11, 13.7, 9.2, 23.4, 15.55, 14.95, 7.1, 11.45, 21, 23.8, 9.75, 8.8, 16.65, 4.4, 4.8, 8.4, 18.1, 4.45, 7.2, 8, 3.6, 17.05, 14.2, 4.2, 9.3, 12.25, 10.95, 9.95, 24.05, 4.65, 13.95, 22.5), Year = c("2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019")), row.names = c(NA, -48L), class = "data.frame")

# Transformation for images 
# Removed external pointers .nifti_image_ptr
# applyTransformation gets angry if it cannot use a transformation with image pointers.
# I realise this will not work for you, but might help.
transformation <- structure(c(0.988278865814209, 0.165658593177795, 0, 0, -0.16611585021019, 
                              0.993303120136261, 0, 0, 0, 0, 1, 0, 8.74881362915039, -8.2611608505249, 
                              0, 1), .Dim = c(4L, 4L), class = "affine", source = structure("NIfTI image", imagedim = c(100L, 
                                                                                                                        100L), pixdim = c(1, 1), pixunits = "Unknown", .nifti_image_ver = 2L, class = c("internalImage", 
                                                                                                                                                                                                        "niftiImage")), target = structure("NIfTI image", imagedim = c(100L, 
                                                                                                                                                                                                                                                                       100L), pixdim = c(1, 1), pixunits = "Unknown", .nifti_image_ver = 2L, class = c("internalImage", 
                                                                                                                                                                                                                                                                                                                                                       "niftiImage")))

coordinate_to_index_central_origin <- function(coordinate,axis="x",radius,res){
  if(axis=="x")return(findInterval(coordinate,seq(-radius,radius,res))) else return(((radius/res)*2)-findInterval(coordinate,seq(-radius,radius,res))+1)
}

coordinates_to_index_central_origin <- Vectorize(coordinates_to_index_central_origin, vectorize.args = 'index')


index_to_coordinates_central_origin <- function(index,axis="x",radius,res){
  if(((radius-(index*res))+(res/2))>radius) warning("Coordinate Point outside radius bounds")
  if(((radius-(index*res))+(res/2))<(-radius)) warning("Coordinate Point outside radius bounds")
  #add half resolution for px centre representation.
  
  if(axis=="y")
    return(
      (radius-(index*res))+(res/2)
    )
  
  if(axis=="x")
    return(
      (-radius+(index*res))-(res/2)
    )
  
}

index_to_coordinates_central_origin <- Vectorize(index_to_coordinates_central_origin,vectorize.args = 'index')


# Creating Transformation --
# How my function will load a transform that has been saved from batch registering images.
#  transformation <- RNiftyReg::loadTransform("matched_trees/100/100_year_2015_project_to_2019.rds")

# Local coord sys to index coordinates --
  data2 <- stand100 # As per above comment.
data2[,"x"] <- coordinate_to_index_central_origin(data2[,"x"],axis = "x",radius = 10, res=0.1)
data2[,"y"] <- coordinate_to_index_central_origin(data2[,"y"],axis = "y",radius = 10, res=0.1)

# Apply transform --
data2[,c("x","y")] <- applyTransform(transform=stand100_2015_2019_affine, x=as.matrix(data2[,c("x","y")]))


# Change from index coordinates to local coord system. --
data2[,"x"] <- index_to_coordinates_central_origin(index=data2[,"x"],axis="x",radius = 10, res=0.1)
data2[,"y"] <- index_to_coordinates_central_origin(index=data2[,"y"],axis="y",radius = 10, res=0.1)


# Start situation 2015: --
  stand100 %>% filter(Year==2015)  %>% ggplot2::ggplot(aes(x=x,y=y)) +
  ggplot2::geom_point(aes(size=Diameter,color=Species))+
  ggplot2::coord_fixed()+
  ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10)+ # add a circle
  ggtitle("Source")


# Plotting expected 2019: --
  
  stand100 %>% filter(Year==2019)  %>% ggplot2::ggplot(aes(x=x,y=y)) +
  ggplot2::geom_point(aes(size=Diameter,color=Species))+
  ggplot2::coord_fixed()+
  ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10)+ # add a circle
  ggtitle("Target")

# Plotting result: --
  data2 %>% filter(Year==2015)  %>% ggplot2::ggplot(aes(x=x,y=y)) + # Error was here plotting Year==2019. Resulting in off rotation. 
  ggplot2::geom_point(aes(size=Diameter,color=Species))+
  ggplot2::coord_fixed()+
  ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10)+ # add a circle.
  ggtitle("Result")

Answer 6 · 2022-03-15T12:40:20.000Z

EDIT: I realised I had plotted the wrong year, which is why the rotation was doubled. Plotting 2015 results in the correct rotation. However translation is off. Compare the result image with the radius filtered target for easy viewing.
Graphical presentation in case you run into issues.

Radius filtered target:

Original target

Answer 7 · 2022-03-16T22:18:31.000Z

Sorry, I'm still struggling to work out what's going on here, and the code doesn't run as-is so it's hard to debug. Can I ask for two further clarifications?

What is stand100_2015_2019_affine? Is it the same as transformation? Is transformation the forward transformation from the registration?
If I understand your conversion function correctly your space runs from -10 to 10 in steps of 0.1, but that would result in indices up to 200, which is twice as wide as the images. How should I reconcile this?

I would write the code along the following lines, which seems to do a reasonable job:

library(RNiftyReg)
library(loder)
library(ggplot2)

stand100 <- structure(list(TreeID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L), x = c(0.74, 0.4, -1.85, -6.1, -7.22, -2.85, -2.18, -3.15, -6.39, -6.21, -2.76, 0.26, 1.17, 3.61, 2.86, 3.51, 4.66, 2.15, 8, 8.98, -1.27, -3.12, -6.86, -7.51, -2.74, -1.95, -2.9, -10.11, -9.72, -5.63, -5.5, -1.43, -3.55, 1.06, 2.46, 2.33, 1.57, 1.28, 2.8, 4.03, 4.74, 5.03, 6.08, 8.54, 8.38, 2.86, 0.32, 2.28), y = c(1.41, 9.48, 8.05, 5.54, 2.88, -0.15, -1.1, -1.28, -4.47, -4.55, -8.55, -6.01, -7.43, -8.8, -7.78, -7.42, -8.57, -4.52, -4.46, 2.18, 9.42, 7.61, 4.46, 1.63, -0.56, -1.31, -1.59, -3.05, -3.83, -5.25, -5.19, -8.77, -10.92, -11.1, -9.67, -7.06, -6.24, -5.76, -4.02, -7.08, -6.57, -7.94, -7.53, -3.06, 3.71, 2.51, 1.43, 9.93), Species = structure(c(2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 6L, 2L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("10", "11", "12", "20", "21", "30", "40", "41", "70", "71", "81", "92", "94", "95"), class = "factor"), Diameter = c(12.6, 14.4, 14, 6.7, 12.8, 19.7, 22.4, 10.4, 4.4, 5.8, 8.1, 17.4, 8.2, 11.4, 5.2, 9.5, 11, 13.7, 9.2, 23.4, 15.55, 14.95, 7.1, 11.45, 21, 23.8, 9.75, 8.8, 16.65, 4.4, 4.8, 8.4, 18.1, 4.45, 7.2, 8, 3.6, 17.05, 14.2, 4.2, 9.3, 12.25, 10.95, 9.95, 24.05, 4.65, 13.95, 22.5), Year = c("2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2015", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019", "2019")), row.names = c(NA, -48L), class = "data.frame")

# This matrix contains information about the resolution (diagonal) and the offset/origin (last column)
xform <- matrix(c(0.1,    0, 0,  -5,
                    0, -0.1, 0,   5,
                    0,    0, 1,   0,
                    0,    0, 0,   1), 4, 4, byrow=TRUE)

# Read images and apply metadata
source <- readPng("source.png")
pixdim(source) <- c(0.1, 0.1, 1)
sform(source) <- structure(xform, code=2L)
target <- readPng("target.png")
pixdim(target) <- c(0.1, 0.1, 1)
sform(target) <- structure(xform, code=2L)

# Register
reg <- niftyreg(source, target)

# Transform the locations in three steps
# World-to-voxel and voxel-to-world convert between indices and "real space" according to the corresponding image's xform matrix
points <- RNifti::worldToVoxel(as.matrix(stand100[stand100$Year==2015,c("x","y")]), source)
points <- applyTransform(forward(reg), points)
points <- RNifti::voxelToWorld(points, target)
result <- data.frame(x=points[,1], y=points[,2], Diameter=stand100$Diameter[stand100$Year==2015], Species=stand100$Species[stand100$Year==2015])

print(ggplot(subset(stand100, Year==2015), aes(x=x, y=y)) + geom_point(aes(size=Diameter,color=Species)) + coord_fixed() + ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10) + ggtitle("Source"))

print(ggplot(subset(stand100, Year==2019), aes(x=x, y=y)) + geom_point(aes(size=Diameter,color=Species)) + coord_fixed() + ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10) + ggtitle("Target"))

print(ggplot(result, aes(x=x, y=y)) + geom_point(aes(size=Diameter,color=Species)) + coord_fixed() + ggforce::geom_circle(aes(x0=0,y0=0,r=10),inherit.aes=FALSE,size=1/10) + ggtitle("Result"))

Here I'm treating the origin offset as -5, not -10, following my comment above. I hope this helps.

Answer 8 · 2022-03-18T12:25:57.000Z

Many thanks for your kind help! This has certainly moved things forward. Is there any way I could at least pay you something for your time and help? I'm a phd student but would be happy to send you something for your troubles.

You are entirely correct, the images should be 200x200 px when there is a resolution of 0.1x0.1 for a 10m radius plot (10+10)/0.1=200, this was my fault.

I've tried running your code, but I get some kind of an error (same as at bottom).

I've tried to detail as clearly as possible now, hopefully it is readable.

Here are the output images, now corrected:

Source:

Target:

Here's exactly how I get the images from stand100:
NB This was much trial and error making the code until the images look right. Maybe this creates issues with conventions for image() or RNifti?
Function to produce matrix for the images (a gaussian noise /w amplitude equal to diameter or other <point_strength>).

#Gauss array Rotate Left, Rotate Left, flip vertical.
#Byrow ... y inverted?
#'@export
produce_gauss_array.default <- function(data,point_strength,position_error=0.5,radius=10,res=0.1){

# Grab the Diameters.
  point_strength_vector <- data %>% select({{point_strength}})

# Make an empty dataframe for the pixels 200x200, long format, columns x indices, y indices, fill with value 0.
  value_matrix <-  data.frame(x=seq((-radius+{res}/2),(radius-{res}/2),{res})) %>% merge(data.frame(y=seq((radius-{res}/2),(-radius+{res}/2),-{res}))) %>% mutate(value=0)

  value_matrix2 <- value_matrix

# Find the highest function value for each point by painting the gaussian noise from the diameters onto the image one by one.
  for(i in 1:nrow(data)){
    x0 <- data[i,"x"]
    y0 <- data[i,"y"]
    point_strength_list <- data %>% select(!!{point_strength})
    point_strength_list <- point_strength_list[i,]
    value_matrix2 <- value_matrix2 %>% mutate(value2=round(point_strength_list*exp(-((((x-x0)^2)/(2*position_error^2))+((y-y0)^2)/(2*position_error^2))),digits = 3))
    value_matrix2 <- value_matrix2 %>% rowwise() %>% mutate(value= max(value,value2,na.rm=TRUE)) %>% select(-value2) %>% ungroup()
  }

# reformat this into a matrix without col or row names. Value is the maximum value produced from the diameter noise process.
  value_matrix2 <- matrix(value_matrix2$value,ncol=2*radius/{res},nrow=2*radius/{res},dimnames=NULL,byrow=TRUE)

# This just used to make sure I can't process something weird.
  class(value_matrix2) <- c("gauss.coord.array","matrix")
  attr(value_matrix2,'res') <- {res}
  attr(value_matrix2,'radius') <- radius

  return(
    value_matrix2
  )

}




#'Convert a coord.data.frame to a matrix with points represented as a gaussian process.
#'@param data A coord.data.frame
#'@param point_strength An attribute to be used as the amplitude of the gaussian process. Typically 'Diameter'.
#'@param position_error Standard deviation of gaussian process.
#'@param radius Radius of plot.
#'@param resolution Resolution of output matrix. Do not recommend higher resolution (lower values) than 0.1. Significant computing time.
#'@export
#'@return A gauss.coord.array containing columns: x,y,value.
produce_gauss_array <- function(data,point_strength,position_error=0.5,radius=10,res=0.1){UseMethod("produce_gauss_array")}

Function to save matrix as png image (incorrect convention for RNiftyReg?) - much inspired by your mmand::display().

save_gauss.default <- function(x,filename){

# Setting the filename
  stopifnot(is.character(filename))
  if(!grepl(filename,pattern='.png$')) filename <- stringr::str_glue(filename,'.png')

  #Z (color range)
  x_range <- range(x[is.finite(x)])


  #Image has odd conventions... first flip matrix vertically and then transpose.
  x2 <- apply(x,2,rev)
  x2 <- t(x2)

  #Create graphics device
  png(filename,width = dim(x),height=dim(x),units = 'px',res=1)

  #No margins, grey background
  par(mai=c(0,0,0,0),bg='grey70')

  #Plot matrix as raster in greyscale.
  image(x2,asp=1,useRaster=TRUE,col=grey(0:255/255),zlim=sort(x_range))

  #Close graphics device.
  dev.off()

}

save_gauss <- function(x,filename){UseMethod("save_gauss")}

My wrapper for your suggested code:

#' Apply the transform between two images to a coord.data.frame
#' @param data Coord.data.frame with points to be transformed. Must contain vars: 'Diameter', 'Year', 'x','y'.
#' @param source_year Coordinate year to transform. From which to produce source image.
#' @param target_year Target coordinate year from which to produce target image.
#' @param radius Default (10 m radius plots) Plot radius, to calculate coordinate indices.
#' @param res Default (0.1) Plot resolution, to calculate coordinate indices.
#' @return A coord.data.frame
#' @export
coordinate_transform <- function(data,target_year,source_year,radius=10,res=0.1){
  #Ignore additional class 'coord.data.frame', just ensures that it has cols: 'id','Year','x','y' and attribute describing coord. type.
  #stopifnot("coord.data.frame" %in% class(data))
  #stopifnot(attr(data,"coordinate_type")=='world')

  # This matrix contains information about the resolution (diagonal) and the offset/origin (last column)
  xform <- matrix(c(res,    0, 0,  -radius,
                    0, -res, 0,   radius,
                    0,    0, 1,   0,
                    0,    0, 0,   1), 4, 4, byrow=TRUE)


  #create a temporary folder
  if(!dir.exists("temp_coordinate_transformation")) dir.create("temp_coordinate_transformation")

  #Target and source data
  target_data <- data %>% filter(Year==target_year)
  source_data <- data %>% filter(Year==source_year)

  #create images for target and source
  target_data %>% produce_gauss_array(radius={radius},point_strength='Diameter',res={res}) %>% save_gauss("temp_coordinate_transformation/target.png")
  source_data %>% produce_gauss_array(radius={radius},point_strength='Diameter',res={res}) %>% save_gauss("temp_coordinate_transformation/source.png")

  #read images
  source <- png::readPNG("temp_coordinate_transformation/source.png")
  target <- png::readPNG("temp_coordinate_transformation/target.png")
  RNiftyReg::pixdim(source) <- c({res},{res},1)
  RNiftyReg::sform(source) <- structure(xform,code=2L)
  RNiftyReg::pixdim(target) <- c({res},{res},1)
  RNiftyReg::sform(target) <- structure(xform,code=2L)

  #register
  # stop() Seems to work fine down to here.
  reg <- niftyreg(source,target)

  # Transform the locations in three steps
  # World-to-voxel and voxel-to-world convert between indices and "real space" according to the corresponding image's xform matrix
  points <- RNifti::worldToVoxel(as.matrix(source_data[,c("x","y")]), source)
  points <- applyTransform(forward(reg), points)
  points <- RNifti::voxelToWorld(points, target)
  result <- data.frame(cbind("x"=points[,1], "y"=points[,2], subset(source_data,select=-c(x,y))))

  out <- rbind(result,target_data) %>% as.coord.data.frame(id = id,x=x,y=y,Species,Diameter,Year)

  #delete temp folder.
  unlink("temp_coordinate_transformation",recursive=TRUE)

  return(
    out
  )

}

I am now getting an error message from running the following:

stand100 %>% coordinate_transform(target_year = 2019, source_year = 2015, radius = 10, res = 0.1)

> [NiftyReg ERROR] Function: initialise_block_matching_method()
> [NiftyReg ERROR] There are no active blocks
> Error during wrapup: [NiftyReg] Fatal error
> Error: no more error handlers available (recursive errors?); invoking 'abort' restart

Answer 9 · 2022-03-21T16:11:45.000Z

OK, this is quite a different problem. In effect NiftyReg is reporting that there isn't enough information in the images to perform the registration, and I have seen this in other contexts. This may be due to a bug in NiftyReg's downsampling routine – although if it is I can't easily see it at the moment – or it could just be because the image is very sparse.

Your original images have a blue cast, and adding a small baseline blue value seems to get the registration going again, viz.

  #read images
  source <- png::readPNG("temp_coordinate_transformation/source.png")
  source[,,3] <- pmax(source[,,3], 0.1)   # Set a minimum blue channel value of 0.1
  target <- png::readPNG("temp_coordinate_transformation/target.png")
  target[,,3] <- pmax(target[,,3], 0.1)
  RNiftyReg::pixdim(source) <- c({res},{res},1)
  RNiftyReg::sform(source) <- structure(xform,code=2L)
  RNiftyReg::pixdim(target) <- c({res},{res},1)
  RNiftyReg::sform(target) <- structure(xform,code=2L)

This is a total hack though, based only on the idea of trying to avoid a large fraction of the images having zero intensity, and I wouldn't really recommend it for use at scale.

Another route might be to only register at full resolution:

  #register
  reg <- niftyreg(source,target,nLevels=1)

It's possible that this might not produce as good an alignment as the usual multi-scale approach, but it is probably worth experimenting with.

By the way, it isn't necessary from RNiftyReg's point of view to save and reload the pixel data to/from file: you could use the result of produce_gauss_array() as the input niftyreg(), as long as you set the attributes as above.

Answer 10 · 2022-03-22T12:40:29.000Z

Registering at full resolution works for my demo! Also, thanks for the tip about not needing to save and reload the images.

Here are the results, absolutely beautiful to me:

Thank you so much for helping me get this far.

I will now try the solution at scale.

Answer 11 · 2022-03-22T12:42:36.000Z

By the way, for whatever reason I have to run the reverse transformation rather than the forward transformation. Unsure why that happens, but still works.

Answer 12 · 2022-03-22T12:56:23.000Z

Excellent! The forward transformation should have the source-to-target sense that you want, but if the reverse transformation works then don’t overthink it, I guess! I’ll close the issue as we seem to have a solution, but you can reopen if necessary.

By the way, your offer of compensation was very kind, but there’s no need. Best of luck with the rest of your PhD!