options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"

library(imager)
library(png)

for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}

for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}

source("nt_toolbox/toolbox_wavelet_meshes/meshgrid.R")
options(repr.plot.width=3.5, repr.plot.height=3.5)

gamma0 <- c(.78, .14, .42, .18, .32, .16, .75, .83, .57, .68, .46, .40, .72, .79, .91, .90) + 1i*c(.87, .82, .75, .63, .34, .17, .08, .46, .50, .25, .27, .57, .73, .57, .75, .79)

periodize <- function(gamma){ c(gamma, gamma[1]) }

cplot <- function(gamma, type, pch, lty, col, lw=2, add=FALSE){
    if (add==FALSE){
        plot(Re(periodize(gamma)), Im(periodize(gamma)), type=type, pch=pch, lty=lty, col=col, lwd=lw, xaxt='n', yaxt='n', ann=FALSE, bty="n")
    }
    else{
        lines(Re(periodize(gamma)), Im(periodize(gamma)), type=type, pch=pch, lty=lty, col=col, lwd=lw, xaxt='n', yaxt='n', ann=FALSE, bty="n")
    }
    
}

options(repr.plot.width=3.5, repr.plot.height=4.5)

cplot(gamma0, type="o", pch=1, lty=1, col="blue")

p <- 256

interpc <- function(x, xf, yf){ approx(xf, Re(yf), x)$y + 1i*approx(xf, Im(yf), x)$y }
curvabs <- function(gamma){ c(0, cumsum(1e-5 + abs(gamma[1:(length(gamma)-1)]-gamma[2:length(gamma)]))) }
resample1 <- function(gamma, d){ interpc((1:p)/p, d/d[length(d)], gamma)}
resample <- function(gamma){ resample1( periodize(gamma), curvabs(periodize(gamma)) ) }

gamma1 <- resample(gamma0)

cplot(gamma1, type="l", pch=1, lty=1, col="black")

shiftR <- function(c){ c(c[length(c)],c[1:(length(c)-1)]) }
shiftL <- function(c){ c(c[2:length(c)],c[1]) }
BwdDiff <- function(c){ c - shiftR(c) }
FwdDiff <- function(c){ shiftL(c) - c }

normalize <- function(v){ return( v/pmax( abs(v), 1e-10 ) ) }
tangent <- function(gamma) { normalize( FwdDiff(gamma) ) }
normal <- function(gamma){ -1i*tangent(gamma) }

delta <- .03
gamma2 <- gamma1 + delta * normal(gamma1)
gamma3 <- gamma1 - delta * normal(gamma1)

cplot(gamma1, type="l", pch=1, lty=1, col="black")
cplot(gamma2, type="l", pch=1, lty=2, col="red", add=TRUE)
cplot(gamma3, type="l", pch=1, lty=2, col="blue", add=TRUE)

normalC <- function(gamma){ BwdDiff(tangent(gamma)) / abs( FwdDiff(gamma) ) }

dt <- 0.001 / 100

Tmax <- 3.0 / 100
niter <- round(Tmax/dt)

gamma <- gamma1

gamma <- gamma + dt * normalC(gamma)

gamma <- resample(gamma)

source("nt_solutions/segmentation_2_snakes_param/exo1.R")

n <- 200

nbumps <- 40
theta <- runif(nbumps)*2*pi
r <- .6*n/2
a <- c(.62*n,.6*n)
x <- round(a[1] + r*cos(theta))
y <- round(a[2] + r*sin(theta))
W <- matrix(rep(0, n*n), c(n,n))
for (i in 1:nbumps){
    W[x[i], y[i]] <- 1
}
W <- gaussian_blur(W,6.0)
W <- rescale( -apply(W, c(1,2), min, 0.05), .3,1)

options(repr.plot.width=3.5, repr.plot.height=3.5)

imageplot(W)

G <- grad(W)
G <- G[,,1] + 1i*G[,,2]

imageplot(abs(G))

EvalG <- function(gamma){ bilinear_interpolate(G, Im(gamma), Re(gamma)) }
EvalW <- function(gamma){ bilinear_interpolate(W, Im(gamma), Re(gamma)) }

r <- .98*n/2 # radius
p <- 128 # number of points on the curve
theta <- t( seq(0,2*pi,length=p+1) )
theta <- theta[1:length(theta)-1]
gamma0 <- n/2 * (1 + 1i) +  r*(cos(theta) + 1i*sin(theta))

gamma <- gamma0

dt <- 1

Tmax <- 5000
niter <- round(Tmax/dt)

lw <- 2
imageplot(t(W))
cplot(gamma, type="l", pch=1, lty=1, col="red", add=TRUE)

dotp <- function(c1,c2){ Re(c1)*Re(c2) + Im(c1)*Im(c2) }

N <- normal(gamma)
g <- - EvalW(gamma) * normalC(gamma) + dotp(EvalG(gamma), N) * N
gamma <- gamma - dt*g

gamma <- resample(gamma)

source("nt_solutions/segmentation_2_snakes_param/exo2.R")

n <- 256
name <- 'nt_toolbox/data/cortex.png'
f <- as.matrix(load_image(name, n))

imageplot(f)

G <- grad(f)
d0 <- sqrt(apply(G**2, c(1,2), sum))
imageplot(d0)

a <- 2
d <- gaussian_blur(d0, a)
imageplot(d)

d <- pmin(d, 0.4)
W <- rescale(-d, .8, 1)

imageplot(W)

p <- 128

source("nt_solutions/segmentation_2_snakes_param/exo3.R")

dt <- 2

Tmax <- 9000
niter <- round(Tmax/ dt)

source("nt_solutions/segmentation_2_snakes_param/exo4.R")

n <- 256
f <- as.matrix(load_image(name, n))
f <- f[46:105, 61:120]
n <- dim(f)[1]

imageplot(f)

source("nt_solutions/segmentation_2_snakes_param/exo5.R")

x0 <- 4 + 55i
x1 <- 53 + 4i

p <- 128
t <- t(seq(0, 1, length=p))
gamma0 <- t*x1 + (1-t)*x0

gamma <- gamma0

imageplot(t(W))
cplot(gamma, type="l", pch=1, lw=2, lty=1, col="red", add=TRUE)
points(Re(gamma[1]), Im(gamma[1]), pch=16, cex=1.5, col="blue")
points(Re(gamma[length(gamma)]), Im(gamma[length(gamma)]), pch=16, cex=1.5, col="blue")

curvabs <- function(gamma){ c(0, cumsum(1e-5 + abs(gamma[1:(length(gamma)-1)]-gamma[2:length(gamma)]))) }
resample1 <- function(gamma, d){ interpc((1:p)/p, d/d[length(d)], gamma)}
resample <- function(gamma){ resample1( gamma, curvabs(gamma) ) }

dt <- 1/10

Tmax <- 2000*4/ 7
niter <- round(Tmax/ dt)

source("nt_solutions/segmentation_2_snakes_param/exo6.R")
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>