How do I line up the color scales of two heat maps (image2D function from the plot3D package)?

Question

EDIT: Turns out it was "just" a misplaced parenthesis in densite_cauchy_norm that skewed the results, everything works fine now.
I have three 2D histograms (heat maps) that I need to compare.
For some reason, one of them (the second on the enclosed image) maps values to colors differently than the two others (e.g a density of 0.020 on the second histogram has the same color that a density of 0.015 on the other two).

These come from Monte Carlo simulation, where I need to compare histograms of 2 samples (X_1, X_2) simulated by rejection sampling with two different density envelopes (plots 1 and 2) to the actual "density" (f_tilde, unnormalized) from which they're sampled (plot 3). So, until they all have the same color mapping, I can't really show how close they both are to the theoretical.

For the first 2 histograms I used the same commands for both samples X_1 and X_2.

Look at the code below; you can find the data at the bottom:

library(plot3D)  
par(mfcol = c(1, 3)); 

n <- 100000;  
X_1 <- simu_f_1(n);  
X_2<- simu_f_2(n);  
x_s <- 50;  y_s <- 50;    

mon_histo(X_1, x_s, y_s, opt_3D = FALSE);  
mon_histo(X_2, x_s, y_s, opt_3D = FALSE);  

z <- outer(seq.int(0, 4, length= x_s), seq.int(0, 2, length = y_s), f_tilde);

image2D(z=z/sum(z), x = seq.int(0, 4, length= x_s),y = seq.int(0, 2, 
                    length = y_s), border="black", contour = FALSE);

The second histogram should approximate the best, but without the same color scales as the others, it seems flat out wrong.

Data and required functions for making a reproducible example:

densite_unif_norm <- function(x,y) { return(dnorm(x, mean = 2, sd = 1) * dunif(y, min = 0, max = 2)); }  
densite_cauchy_norm <- function(x,y) { return( dnorm(x, mean = 2, sd = 1)) * dcauchy(y, location = 1, scale = 0.5); }  
f_tilde <- function(x,y) {  
    return( (x>=0)*(x<=4)*(y>=0)*(y<=2) * exp(-0.5*(x-2)**2)*  
                (cos(x)**2+(2*sin(y)**2)*(cos(x)**4)) / (1+4*(y-1)**2)  );  
}
simu_f_1 <- function(n) { #simulates size n sample X_1 ~ "density" f_tilde by rejection sampling  
    M <- 6*sqrt(2*pi); # M = 15.04  
    f <- function(x) { 
        repeat{  
            Y <- c(rnorm(n = 1, mean = 2, sd = 1), runif(n = 1, min = 0, max = 2));  
            U <- runif(n = 1, min = 0, max = M*densite_unif_norm(Y[1], Y[2]));  
            if(U < f_tilde(Y[1], Y[2])) { break; } #la boucle s'arrête à la dernière valeur de y t.q u < f_tilde(y)  
        }  
        return(rbind(Y[1], Y[2]));  
    }
    return(Vectorize(f)(1:n));  
}

simu_f_2 <- function(n) { 
    M <- 3*(pi^1.5)/sqrt(2); # M = 11.81  
    f <- function(x) {  
        repeat{  
            Y <- c(rnorm(n = 1, mean = 2, sd = 1), rcauchy(n = 1, location = 1, scale = 0.5));  
            U <- runif(n = 1, min = 0, max = M*densite_cauchy_norm(Y[1], Y[2]));  
            if(U < f_tilde(Y[1], Y[2])) { break; }
        }  
        return(rbind(Y[1], Y[2]));  
    }  
    return(Vectorize(f)(1:n));  
}  
mon_histo <- function(X, x_sub, y_sub, opt_3D = FALSE) {  
    x_c <- cut(X[1,], x_sub);  
    y_c <- cut(X[2,], y_sub);  
    z <- table(x_c, y_c);  
    if(opt_3D) {hist3D(z=z, border="black"); }  
    else {   
        u <- seq(from = range(X[1,])[1], to = range(X[1,])[2], length.out = x_sub);  
        v = seq(from = range(X[2,])[1], to = range(X[2,])[2], length.out = y_sub);  
        image2D( z=z/length(X[1,]), border="black", x = u, y = v);  
    }  
}