# Otsu threshold (no extra packages)
otsu <- function(x, nbins=256){
  h <- hist(as.numeric(x), breaks=seq(0,1,length.out=nbins+1), plot=FALSE)
  p <- h$counts / sum(h$counts)
  m <- h$mids
  omega <- cumsum(p)
  mu <- cumsum(p*m)
  mu_t <- mu[length(mu)]
  # guard against 0/0 at ends
  sb2 <- (mu_t*omega - mu)^2 / pmax(omega*(1-omega), .Machine$double.eps)
  m[which.max(sb2)]
}

library(depstats)
library(imager)

##3
#Correlated Laplace
t3 <- function(n, rho) {
  X <- LaplacesDemon::rmvl(n, c(0, 0),
                           matrix(c(1,   rho,
                                    rho, 1  ), nrow = 2, byrow = TRUE))
  return(X)
}
##5
#Ishigami Style
t5 <- function(n, T = 1) {
  U <- runif(n)
  V <- runif(n)
  W <- runif(n)
  x <- U
  y <- sin(T * U) + 4 * sin(V) ^ 2 + 0.5 * W ^ 4 * sin(T * U)
  X <- cbind(x, y)
  return(X)
}
##6
#Tree Ring
t6 <- function(n, sigma = 0.1, rings = 10) {
  L <- sample(1:rings, n, replace = TRUE)
  theta <- runif(n, 0, 2 * pi)
  epsilonx <- rnorm(n, 0, sigma)
  epsilony <- rnorm(n, 0, sigma)
  x <- L * cos(theta) + epsilonx / 4
  y <- L * sin(theta) + epsilony / 4
  X <- cbind(x, y)
  return(X)
}
##7
#Change in variance
t7 <- function(n, p) {
  x <- runif(n)
  e <- rnorm(n)
  y <- abs(x) ^ p * e
  X <- cbind(x, y)
  return(X)
}



n <- 400

op <- par(
  mfrow = c(2, 3),
  mar   = c(2, 2, 1.5, 1),
  oma   = c(0, 0, 0, 0),
  pty   = "s"    # make each panel a square
)


##3
#Correlated Laplace
set.seed(1)
X1 <- depgen(1, n, 't3(n, runif(1, 0.5, 0.6))',randrotate = FALSE, workers = 1)
plot(X1, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title("Correlated Laplace")
#abline(a = 0, b = 0.55, lwd = 3, lty = 2)
##5
#Ishigami Style
set.seed(1)
X1 <- depgen(1, n, 't5(n, runif(1, pi, 3 * pi))',randrotate = FALSE, workers = 1)
plot(X1, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title("Ishigami style")
myf <- function(x) {y <- sin(2 * pi * x) + 1.092 + sin(2 * pi * x) / 10; return(y)}
#curve(myf, lwd = 3, lty = 2, add = TRUE)
##6
#Tree Ring
set.seed(1)
X1 <- depgen(1, n, 't6(n, runif(1, 0, 0.5), sample(2:10, 1))',randrotate = FALSE, workers = 1)
plot(X1, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title("Tree Ring")
myf1 <- function(x) {y <- 3.5 * sin(acos(x / 3.5)); return(y)}
myf2 <- function(x) {y <- 3.5 * sin(acos(x / 3.5)); return(-y)}
#curve(myf1, lwd = 3, lty = 2, add = TRUE, n = 100000)
#curve(myf2, lwd = 3, lty = 2, add = TRUE, n = 100000)
##7
#Change in variance
set.seed(1)
X1 <- depgen(1, n, 't7(n, runif(1, 0.5, 0.6))',randrotate = FALSE, workers = 1)
plot(X1, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title("Change in variance")
myf1 <- function(x) {y <- rep(0, length(x)); return(y)}
myf2 <- function(x) {y <- sqrt(2 / pi) * x ^ 0.55; return(y)}
myf3 <- function(x) {y <- -sqrt(2 / pi) * x ^ 0.55; return(y)}
#curve(myf1, lwd = 3, lty = 2, add = TRUE)
#curve(myf2, lwd = 3, lty = 2, add = TRUE)
#curve(myf3, lwd = 3, lty = 2, add = TRUE)


I1  <- load.image("I1.jpg")
Z   <- as.array(grayscale(I1))[,,1,1]     # x × y, values in [0,1]
thr <- otsu(Z)
BW <- ifelse(Z >= thr, 1, 0)              # light → white, dark → black
nx <- nrow(BW); ny <- ncol(BW)
x  <- seq(-0.5, 0.5, length.out = nx)
y  <- seq(-0.5, 0.5, length.out = ny)
image(x, y, BW[, ny:1], col = NA, #c("black","white"),
      useRaster=TRUE, asp=1, xaxs="i", yaxs="i",
      xlab="", ylab="", xaxt="n", yaxt="n")
set.seed(1)
XI1 <- depgen(1, n, 'normnoise(randimage(n, I1), runif(1, 0, 0.1))', randrotate = FALSE, print = TRUE, worker = 1)
points(XI1, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title(expression(bold(infinity)), cex.main = 3)

I2  <- load.image("I2.jpg")
Z   <- as.array(grayscale(I2))[,,1,1]     # x × y, values in [0,1]
thr <- otsu(Z)
BW <- ifelse(Z >= thr, 1, 0)              # light → white, dark → black
nx <- nrow(BW); ny <- ncol(BW)
x  <- seq(-0.5, 0.5, length.out = nx)
y  <- seq(-0.5, 0.5, length.out = ny)
image(x, y, BW[, ny:1], col = NA, #c("black","white"),
      useRaster=TRUE, asp=1, xaxs="i", yaxs="i",
      xlab="", ylab="", xaxt="n", yaxt="n")
set.seed(1)
XI2 <- depgen(1, n, 'normnoise(randimage(n, I2), runif(1, 0, 1))', randrotate = FALSE, print = TRUE, worker = 1)
points(XI2, pch = 19, col = "darkblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex = 0.6, xlim = range(X1[, 1]), ylim = range(X1[, 2]))
title(expression(bold(pi)), cex.main = 3)

par(op)


dev.copy(
  jpeg,
  filename = "ADD.jpg",
  width    = 8,      # inches
  height   = 4,      # inches (adjust to your panel aspect)
  units    = "in",
  res      = 300,    # dots per inch
  quality  = 100     # max JPEG quality
)
dev.off()           # close the file device