options(warn=-1) # turns off warnings, to turn on: "options(warn=0)" library(imager) library(png) library(misc3d) library(SynchWave) 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) M <- read_bin("nt_toolbox/data/vessels.bin", ndims=3) M <- rescale(M) n <- dim(M)[2] options(repr.plot.width=6, repr.plot.height=6) slices <- seq(10,n-10, length = 4) for (i in 1:length(slices)){ s <- slices[i] s <- as.integer(s) imageplot(M[,,s], paste("Z = ",s), c(2,2,i)) } options(repr.plot.width=8, repr.plot.height=6) isosurface(M,0.5,3) MW <- M a <- (MW[seq(1,n,2),,] + MW[seq(2,n,2),,])/sqrt(2) b <- (MW[seq(1,n,2),,] - MW[seq(2,n,2),,])/sqrt(2) c <- array(0, c(n, n, n)) c[1:(n/2),,] <- a c[(n/2+1):n,,] <- b MW <- c a <- (MW[,seq(1,n,2),] + MW[,seq(2,n,2),])/sqrt(2) b <- (MW[,seq(1,n,2),] - MW[,seq(2,n,2),])/sqrt(2) c <- array(0, c(n, n, n)) c[,1:(n/2),] <- a c[,(n/2+1):n,] <- b MW <- c a <- (MW[,,seq(1,n,2)] + MW[,,seq(2,n,2)])/sqrt(2) b <- (MW[,,seq(1,n,2)] - MW[,,seq(2,n,2)])/sqrt(2) c <- array(0, c(n, n, n)) c[,,1:(n/2)] <- a c[,,(n/2+1):n] <- b MW <- c options(repr.plot.width=8, repr.plot.height=6) imageplot(MW[,,30], "Horizontal slice", c(1,2,1)) imageplot((MW[,30,]), "Vertical slice", c(1,2,2)) source("nt_solutions/multidim_2_volumetric/exo1.R") ## Insert your code here. m <- round(0.01*n**3) MWT <- perform_thresholding(MW, m, type="largest") source("nt_solutions/multidim_2_volumetric/exo2.R") s <- 30 imageplot(M[, ,s], 'Original', c(1,2,1)) imageplot(clamp(M1[, ,s]), 'Approximation', c(1,2,2)) isosurface(M1,.5,2) sigma <- 0.06 Mnoisy <- M + array(rnorm(n**3, sd = sigma), c(n,n,n)) imageplot(Mnoisy[,,n/2],"X slice",c(1,2,1)) imageplot(Mnoisy[,n/2,],"Y slice",c(1,2,2)) x <- (-round(n/2):(round(n/2)-1)) grid <- meshgrid_3d(x, x, x) X <- grid$X ; Y <- grid$Y ; Z <-grid$Z s <- 2 #width h <- exp(-(X**2 + Y**2 + Z**2)/(2*(s**2))) h <- h/sum(h) Mh <- Re(fft( fft(Mnoisy) * fft(fftshift_3d(h)), inverse=T)/length(h)) options(repr.plot.width=8, repr.plot.height=6) i <- 40 imageplot(Mnoisy[,,i], "Noisy", c(1,2,1)) imageplot(Mh[,,i], "Denoised", c(1,2,2)) isosurface(rescale(Mh),.5,3) options(repr.plot.width=5, repr.plot.height=5) source("nt_solutions/multidim_2_volumetric/exo3.R") ## Insert your code here. options(repr.plot.width=8, repr.plot.height=6) isosurface(Mblur,.5,2) paste("Filtering, SNR = ", snr(M, Mblur), "dB") options(repr.plot.width=5, repr.plot.height=5) source("nt_solutions/multidim_2_volumetric/exo4.R") ## Insert your code here. options(repr.plot.width=8, repr.plot.height=6) isosurface(Mwav,.5,2) paste("Soft thresholding, SNR = ", snr(M, Mwav), "dB") w <- 4 grid <- meshgrid_3d(0:(w-1), 0:(w-1), 0:(w-1)) dZ <- grid$X ; dY <- grid$Y ; dX <- grid$Z Mspin <- array(0, c(n,n,n)) circshift <- function(x,v){ x <- roll(x,v[1], axis = 1) x <- roll(x,v[2], axis = 2) x <- roll(x,v[3], axis = 3) return(x) } for (i in 1:w**3){ # shift the image MnoisyC <- circshift(Mnoisy, c(dX[i],dY[i],dZ[i])) # denoise the image to get a result M1 M1 <- MnoisyC; # replace this line by some denoising # shift inverse M1 <- circshift(M1, c(-dX[i],-dY[i],-dZ[i])) # average the result Mspin <- Mspin*(i-1)/(i) + M1/(i) } source("nt_solutions/multidim_2_volumetric/exo5.R") ## Insert your code here. options(repr.plot.width=8, repr.plot.height=6) isosurface(Mspin,.5,2) paste("Cycle spinning, SNR = ", snr(M, Mspin), "dB")