How to create correlated Latin Hypercube samples?

Question

I have been trying to generate 2-d correlated Latin Hypercube samples using a 1-factor Gaussian copula model. However, the resulting correlated quantiles do not necessarily have 1 sample point in each row-column in the 2-d Latin cube.

Below is an example R code using a 2-d case. I've coded a simple LHS random number generator and have not used an R package for that, but I believe this should not matter. The plot shows the 2-d latin cube for uncorrelated (black points) and correlated (red points) samples. You can see that although this seems correct, the correlated samples do not occur at each row/column as it should be for in LHS.

Simple LHS random number generator:

rnlhs <- function(n) {
  # define strata
  brks <- seq(0, 1, length.out=(n + 1))
  strata <- data.frame(start=brks[1:n], end=brks[2:(n + 1)])
  s <- runif(n, strata$start, strata$end) 
  # reshuffle
  sample(s)
}

rho=0.5; n=10

### Random case
set.seed(9)
q1 <- rnlhs(n)
q2 <- rnlhs(n)

plot(q1, q2, xlim=c(0, 1), ylim=c(0, 1))
abline(h=seq(0, 1, 0.1), v=seq(0, 1, 0.1), col="gray", lty=3)

### Correlated samples (factor model)
q1 <- rnlhs(n)
q2 <- rnlhs(n)
q0 <- rnlhs(n)
X0 <- qnorm(q0, 0, 1)
X1 <- qnorm(q1, 0, 1)
X2 <- qnorm(q2, 0, 1)
Y1 <- sqrt(rho)*X0 + sqrt(1 - rho)*X1
Y2 <- sqrt(rho)*X0 + sqrt(1 - rho)*X2

# correlated quantiles using the normal CDF
q11 <- pnorm(Y1, 0, 1)
q22 <- pnorm(Y2, 0, 1)

points(q11, q22, col="red", pch=16)

The resulting correlated quantiles, q11 and q22, cannot be really defined as Latin Hypercube samples since – for example – both q11, q22 do not have a sample point within [0.1, 0.2]. How is it possible to insure that one correlated samples point occurs in each row and column?