First commit!

.Rbuildignore

@@ -0,0 +1,3 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^\.travis\.yml$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

.travis.yml

@@ -0,0 +1,5 @@
+# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r
+
+language: R
+sudo: false
+cache: packages

DESCRIPTION

@@ -0,0 +1,13 @@
+Package: RED
+Type: Package
+Title: REgularization by Denoising (RED)
+Version: 0.1.0
+Author: person("Adriano", "Passos", email="adriano.utfpr@gmail.com", role=c("aut","cre"))
+Maintainer: The package maintainer <yourself@somewhere.net>
+Description: More about what it does (maybe more than one line)
+    Use four spaces when indenting paragraphs within the Description.
+Depends: R (>= 3.3.0), imager
+License: GPL-3
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 6.0.1

NAMESPACE

@@ -0,0 +1,9 @@
+# Generated by roxygen2: do not edit by hand
+
+export(MAE)
+export(MSE)
+export(RED)
+export(degrade)
+export(resample)
+export(shift)
+import(imager)

R/RED_pkg.R

@@ -0,0 +1,10 @@
+#' SRI: Super-resolution imaging package
+#'
+#' Super-resolution imaging (SR) is a class of techniques that enhance the resolution of an imaging system.
+#'
+#' @import imager
+#'
+#' @docType package
+#' @name SRI
+NULL
+

R/data.R

@@ -0,0 +1,19 @@
+#' Photograp of a cameraman
+#'
+#' This image is uslay used as benchmar in SR problems
+#'
+#' @format an image of class \code{cimg}
+'cameraman'
+
+# cameraman <- grayscale(cameraman)
+# devtools::use_data(cameraman)
+
+#' Photograp of Lenna
+#'
+#' The Lenna (or Lena) picture is one of the most widely used standard test
+#' images used for compression algorithms
+#'
+#' @format an image of class \code{cimg}
+'lenna'
+
+# devtools::use_data(lenna)

R/degrade.R

@@ -0,0 +1,49 @@
+#' Degradation of an image
+#'
+#' This function degrades a high resolution image into a low resolution
+#' image. The degradation folows this: bla bla bla
+#'
+#' @export
+#' @param z a cimg object containing the high resolution image
+#' @inheritParams resample
+#' @inheritParams shift
+#' @param noise numeric indicating the standard deviation of the noise or an cimg object that will be added to the resampled z
+#' @param blur numeric indicating the blur range (for uniform blur) or an cimg object with the blur kernel to be convolved with z
+#'  if nothing is provided an default kernel will be used.
+#'
+#' @return A degraded cimg object
+#'
+#' @examples
+#' degraded.lenna <- degrade(lenna, L = 4, noise = 0.05, blur = 9)
+#' par(mfrow = c(1,2), mar = c(0,0,1,0)+0.1)
+#' plot(lenna, axes = FALSE, interp = FALSE, main = 'Original Lenna')
+#' plot(degraded.lenna, axes = FALSE, interp = FALSE, main = 'Degraded Lenna')
+degrade <- function(z, L = 1, s = c(0,0), noise = 0, blur = 1, L1 = L, L2 = L){
+
+  N <- dim(z)
+  M <- round(c(N[1]/L1, N[2]/L2, 1, 1))
+
+  if (is.cimg(noise))
+    noise <- noise
+  else if (noise == 0)
+    noise <- noise
+  else if(noise > 0)
+    noise <- imager::as.cimg(array(stats::rnorm(prod(M), mean = 0, sd = noise), M))
+  else
+    stop("noise must be non-negative scalar or cimg object")
+
+  if (is.cimg(blur))
+    blur <- blur
+  else if (length(blur) == 1)
+    blur <- imager::imfill(blur, blur, 1, 1/blur^2)
+  else if (is.null(blur)) # change condition to a name
+    blur <- imager::imfill(5, 5, 1, 1/5^2) #change to PSF #incoherent psf finite detector
+
+  y <- z
+  y <- shift(y, s)
+  y <- convolve(y, blur)
+  y <- resample(y, L1 = 1/L1, L2 = 1/L2)
+  y <- y + noise
+
+  return(y)
+}

R/error.R

@@ -0,0 +1,35 @@
+#' Error measurements of images
+#'
+#' This function calculates error between two images
+#'
+#' @param x,y \code{cimg} objects
+#'
+#' @examples
+#' degraded.lenna <- degrade(lenna, noise = 0.05)
+#' MSE(lenna, degraded.lenna)
+#' MAE(lenna, degraded.lenna)
+#' #alternatively it can be done like:
+#' MSE(lenna - degraded.lenna)
+#' MAE(lenna - degraded.lenna)
+#' @name error
+NULL
+
+#' @export
+#' @describeIn error Mean Squared Error
+MSE <- function(x, y = NULL){
+
+  if(is.null(y))
+    return(mean(x^2))
+  else
+    return(mean((x - y)^2))
+}
+
+#' @export
+#' @describeIn error Mean Absolute Error
+MAE <- function(x, y = NULL){
+
+  if(is.null(y))
+    return(mean(abs(x)))
+  else
+    return(mean(abs(x - y)))
+}

R/red.R

@@ -0,0 +1,238 @@
+#' REgularization by Denoising
+#'
+#' @param y cimg object with the low resolution frame(s)
+#' @param mu numeric indicating the step size
+#' @param lambda,sigma numeric indicating the regularization parameters
+#' @param functional character with the optimization task or function with the functional to be used
+#' @param engine character indicating the denoised engine or function with the denoiser engine to be used
+#' @param niter numeric indicating the maximum number of iterations
+#' @param args arguments to be passed implicity to H HT and f
+#'
+#' @export
+#' @examples
+#'
+#' y <- degrade(lenna, noise = 0.05)
+#' x <- RED(y, sigma = 1, lambda = 5, mu = 0.1, 'DN', niter = 10)
+#' par(mfrow = c(1,2), mar = c(0,0,2,0)+0.1)
+#' plot(y, interp = FALSE, axes = FALSE, main = 'Degraded Lenna')
+#' mtext(MSE(y, lenna), side = 1, line = -2)
+#' plot(x, interp = FALSE, axes = FALSE, main = 'Restored Lenna')
+#' mtext(MSE(x, lenna), side = 1, line = -2)
+#'
+RED <- function(y, lambda, sigma, mu = NULL, functional = 'SR', engine = 'MF', niter = 50, args = NULL){
+
+#  sigma <- sigma * 255
+#  y <- y * 255
+
+  f <- NULL
+
+  if (engine == 'MF')
+    f$dn <- function(x) medianblur(x, n = 3, threshold = 0)
+  else if(is.function(engine))
+    f$dn <- engine
+  else
+    stop('Unsupported denoise engine')
+
+  if (functional == 'SR'){
+    f$H <- function(im){ #transform HR_FRAME to LR_FRAME
+      im <- lapply(1:nrow(args$par), function(i){
+        res <- im
+        #res <- imager::imshift(res, par[i,1], par[i,2])
+        res <- shift(res, args$par[i,])
+        res <- imager::isoblur(res, args$sigma_blur)
+        #res <- imager::resize(res, ncol(im)/args$scale, nrow(im)/args$scale, interpolation_type = 2, boundary_conditions = 1)
+        res <- imager::resize(res, ncol(im)/args$scale, nrow(im)/args$scale, interpolation_type = 5)
+        return(res)
+      })
+      im <- imappend(as.imlist(im), 'z')
+      return(im)
+    }
+    f$HT <- function(im){ #transform LR_FRAME to HR_FRAME
+      im <- imsplit(im, 'z')
+      im <- lapply(1:length(im), function(i){
+        res <- im[[i]]
+        #res <- imager::resize(res, ncol(im[[i]])*args$scale, nrow(im[[i]])*args$scale, interpolation_type = 2, boundary_conditions = 1)
+        res <- imager::resize(res, ncol(im[[i]])*args$scale, nrow(im[[i]])*args$scale, interpolation_type = 1)
+        res <- imager::imsharpen(res, args$amplitude, type = 'shock', alpha = 5, sigma = args$sigma_blur)
+        #res <- imager::imshift(res, -par[i,1], -par[i,2])
+        res <- shift(res, -args$parpar[i,])
+        return(res)
+      })
+      im <- average(as.imlist(im))
+      return(im)
+    }
+    x <- f$HT(imsplit(y, 'z')[[1]])
+  }
+
+  if (functional == 'DN'){
+    f$H <- function(im){
+      return(im)
+    }
+    f$HT <- function(im){
+      return(im)
+    }
+    x <- y
+  }
+
+  N <- prod(dim(x))
+
+  if(is.null(mu))
+    mu <- 2/(1/(sigma^2) + lambda)
+  p <- dim(y)[3]
+
+  ## stepest descent
+  for(i in 1:niter){
+    dif <- f$H(x) - y
+    difn <- x - f$dn(x)
+
+    grad <- (1/(sigma^2))*f$HT(dif) + lambda*difn
+    loss <- (1/(2*(sigma^2)))*sum(dif^2)/p + (lambda/2)*sum(x*difn)
+    #x <- x - mu*zero.outliers(grad)
+    x <- x - mu*grad
+    cat(loss/N,'\n')
+#    plotOR(x, q = c(0,1), interp = F)
+
+  }
+
+  return(x)
+}
+
+
+### checking H and HT stability
+if(F){
+
+  ###
+  args <- list(sigma_blur = 1.6, amplitude = 1.6)
+  f <- NULL
+  f$H <- function(im){ #transform HR_FRAME to LR_FRAME
+    im <- isoblur(im, args$sigma_blur, gaussian = TRUE)
+    return(im)
+  }
+  f$HT <- function(im){ #transform LR_FRAME to HR_FRAME
+    im <- imsharpen(im, args$amplitude, type = 'shock', alpha = 5, sigma = args$sigma_blur)
+    return(im)
+  }
+
+  im <- im0 <- lenna
+  loss <- NULL
+
+  for (i in 1:20){
+    im <- f$HT(f$H(im))
+    loss <- c(loss, MSE(im, im0))
+  }
+  par(mfrow = c(1,2))
+  plot(loss)
+  plot(im, axes = F)
+
+  ####
+  args <- list(scale = 4)
+  f <- NULL
+  f$H <- function(im){ #transform HR_FRAME to LR_FRAME
+    im <- resize(im, ncol(im)/args$scale, nrow(im)/args$scale, interpolation_type = 5)
+    return(im)
+  }
+  f$HT <- function(im){ #transform LR_FRAME to HR_FRAME
+    im <- resize(im, ncol(im)*args$scale, nrow(im)*args$scale, interpolation_type = 1)
+    return(im)
+  }
+
+  im <- im0 <- lenna
+  loss <- NULL
+
+  for (i in 1:20){
+    im <- f$HT(f$H(im))
+    loss <- c(loss, MSE(im, im0))
+  }
+  par(mfrow = c(1,2))
+  plot(loss)
+  plot(im, axes = F)
+
+  ####
+  args <- list(par = c(1,1))
+  f <- NULL
+  f$H <- function(im){ #transform HR_FRAME to LR_FRAME
+    im <- shift(im, par = args$par)
+    return(im)
+  }
+  f$HT <- function(im){ #transform LR_FRAME to HR_FRAME
+    im <- shift(im, par = - args$par)
+    return(im)
+  }
+
+  im <- im0 <- lenna
+  loss <- NULL
+
+  for (i in 1:20){
+    im <- f$HT(f$H(im))
+    loss <- c(loss, MSE(im, im0))
+  }
+  par(mfrow = c(1,2))
+  plot(loss)
+  plot(im, axes = F)
+
+  ####
+  args <- list(par = c(0.5,0.5))
+  f <- NULL
+  f$H <- function(im){ #transform HR_FRAME to LR_FRAME
+    im <- shift(im, par = args$par)
+    return(im)
+  }
+  f$HT <- function(im){ #transform LR_FRAME to HR_FRAME
+    im <- shift(im, par = - args$par)
+    return(im)
+  }
+
+  im <- im0 <- lenna
+  loss <- NULL
+
+  for (i in 1:20){
+    im <- f$HT(f$H(im))
+    loss <- c(loss, MSE(im, im0))
+  }
+  par(mfrow = c(1,2))
+  plot(loss)
+  plot(im, axes = F)
+
+}
+
+if(F){
+  Lambda <- 10^(seq(-2,1,0.25))
+  y <- degrade(lenna, sd_noise = 0.05, sd_blur = 0)
+  x <- NULL
+  mse <- MSE(y, lenna)
+  for(lambda in Lambda){
+    im <- RED(y, sigma = 1, lambda = lambda, mu = 0.1, 'DN', niter = 10)
+    mse <- c(mse, MSE(im, lenna))
+    x <- c(x, imsplit(im, 'z'))
+  }
+  plot(c(0,Lambda), mse)
+  display(as.imlist(x))
+}
+
+### deconvolução!!!! PORRAAA FILHA DA PUTAAA FOIII!!!!!
+if(F){
+  filter <- imfill(9,9,val = 1)
+  filter <- filter/sum(filter)
+  im <- im0 <- lenna
+
+  ffilter <- pad(filter, 1, 'xy', 1)
+  ffilter <-pad(ffilter, nrow(im)-10, 'xy', 0)
+  fft.filter <- FFT(ffilter)
+  fft.filter <- fft.filter$real + fft.filter$imag*1i
+  fft.im <- FFT(im)
+  fft.im <- fft.im$real + fft.im$imag*1i
+  fft.conv <- fft.filter * fft.im
+  conv <- FFT(Re(fft.conv), Im(fft.conv), inverse = TRUE)$real
+  conv <- imappend(imsplit(imappend(imsplit(conv, 'x', 2)[2:1], 'x'), 'y', 2)[2:1], 'y')
+  plot(conv, interp = FALSE)
+  plot(convolve(im, filter), interp = FALSE)
+
+  im <- conv
+  fft.im <- FFT(im)
+  fft.im <- fft.im$real + fft.im$imag*1i
+  fft.deconv <- fft.im / fft.filter
+  deconv <- FFT(Re(fft.deconv), Im(fft.deconv), inverse = TRUE)$real
+  deconv <- imappend(imsplit(imappend(imsplit(deconv, 'x', 2)[2:1], 'x'), 'y', 2)[2:1], 'y')
+  plot(deconv, interp = FALSE)
+  plot(im0, interp = FALSE)
+}