R: Finding the nearest raster cell within a thresh

I have two perfectly overlapping rasters (same extents and cell size). For every cell in one raster (i.e. for every XY), I would like to determine the Euclidean geographical distance to the closest cell within a given threshold difference between the rasters.

Put another way: raster1 and raster2 measure some variable Z. I have a threshold difference for Z values (t) which constitutes a "matching" value (or "close enough") between raster1 and raster2. For each reference cell in raster1, I need to 1) find all the cells in raster2 with a Z value of abs(Z2-Z1)

Each raster has ~26 million cells, ~10 million of which have non-NA values. I have come up with a non-raster-based work-around for this problem, but only by converting the rasters to XYZ tables/vectors and performing a looping function for each reference cell. This is much too computationally intensive for the data size that I'm dealing with (takes ~10 days to process!). To assist comprehension of my question, however, that code is as follows:

library(SDMTools)
c.in <- asc2dataframe("reference.asc"); names(c.in) <- c("X","Y","Z")
f.in <- asc2dataframe("destination.asc"); names(f.in) <- c("X","Y","Z")

x=c.in$X
y=c.in$Y
c=c.in$Z
f=f.in$Z
dist=vector(length=length(c))
threshold <- 0.01

id <- 1:length(c)
for (i in length(id)) {
  # First, find all rows within the threshold
  t <- id[abs(f-c[i])<threshold]
  # Second, find the distance to the closest row
  dist[i] <- round(sqrt(min((x[t]-x[i])^2+(y[t]-y[i])^2)))
}

library(raster)
dist.rast <- rasterFromXYZ(x,y,dist)

You could set the values crossing the threshold to NA and then compute the 'as-the-crow-flies' shortest distance using the distance() function alongside the direction() function within the raster package. You'd then have two raster layers or matrices specifying the exact location of each cell (distance and direction), after a small Pythagorean calculation incorporating the raster resolution. For the sake of simplicity, you may want to remove the spatial projection from the rasters beforehand to remove the ellipsoidal component of the distance and direction calculations. They are easily added back later. If all this is too slow, I recommend trying sparseMatrix() or coding it in IDL or MATLAB. If you use RRO, which ships with MKL optimization, that should improve the performance of matrix operations.

You do excellent research by the way :) Say hi to Andreas for me.