my Cpp&R image convolution

cpp代码实现黑白照片的处理

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
int mirrorIndex_c(int fetchI, int length){
  // code for mirror index at boundary
  if(fetchI<0)
    fetchI=-1-fetchI;
  if(fetchI>=length)
    fetchI=length-1-(fetchI-length);
  return fetchI;
}
// [[Rcpp::export]]
mat imageConv_c(mat img,mat kernel){
  // code for image convolution
  int m,n,x,y;
  m=img.n_rows;
  n=img.n_cols;
  mat new_img(m, n);
  double fin;
  for(int i=1;i<m-1;i++){
    for(int j=1;j<n-1;j++){
      fin=0;
      for(int p=-1;p<=1;p++){
        for(int q=-1;q<=1;q++){
          x=mirrorIndex_c(i+p,m);
          y=mirrorIndex_c(j+q,n);
          fin += img(x,y)*kernel(p+1,q+1);
        }
      }
    if(fin<0)fin=0;
    if(fin>1)fin=1;
    new_img(i,j)=fin;
    }
  }
  return new_img;
}

R代码实现黑白照片的处理

mirrorIndex_R<-function(fetchI, length){
  if(fetchI<1)
    fetchI=1-fetchI
  if(fetchI>=length)
    fetchI=length-1-(fetchI-length)
  return(fetchI)
}

imageConv_R<-function(img,kernel){
  m<-dim(img)[1]
  n<-dim(img)[2]
  new_img<-matrix(data=0,nrow=m,ncol=n)
  x<-0
  y<-0
  for(i in 1:m){
    for(j in 1:n){
      fin<-0
      for(p in -1:1)
        for(q in -1:1){
          x<-mirrorIndex_R(i+p,m)
          y<-mirrorIndex_R(j+q,n)
          fin = fin + img[x,y]*kernel[p+2,q+2]
        }
      if(fin<0){
        fin=0}else if(fin>1){fin=1}
      new_img[i,j]=fin
    }
  }
  return(new_img)
}

R中调用上述cpp code

Rcpp::sourceCpp('C:/Users/tianj/mycon.cpp')
smoothM = matrix(c(1,2,1,2,4,2,1,2,1)/16, 3,3)
sharpenM = matrix(c(0,-1,0,-1,5,-1,0,-1,0),3,3)
edgeM = matrix(c(-1,-1,-1,-1,8,-1,-1,-1,-1),3,3)
library(jpeg)
img = readJPEG("C:/Users/tianj/Desktop/cat.jpg")
img = img[,,1]
img2C = imageConv_c(img, smoothM) 
img3C = imageConv_c(img, sharpenM) 
img4C = imageConv_c(img, edgeM) 
par(mar = c(0,0,0,0))
plot(0,0,type = "n")
rasterImage(img, -1, 0, 0, 1)
rasterImage(img2C, 0, 0,1,1)
rasterImage(img3C,-1,-1,0,0)
rasterImage(img4C,0,-1,1,0)

R中调用上述R code

source('~/.active-rstudio-document')
img2R = imageConv_R(img, smoothM) 
img3R = imageConv_R(img, sharpenM) 
img4R = imageConv_R(img, edgeM) 
par(mar = c(0,0,0,0))
plot(0,0,type = "n")
rasterImage(img, -1, 0, 0, 1)
rasterImage(img2R, 0, 0,1,1)
rasterImage(img3R,-1,-1,0,0)
rasterImage(img4R,0,-1,1,0)

测试cpp卷积函数和R卷积函数的速度

library(microbenchmark)
testConv = microbenchmark(imageConv_c(img, smoothM),imageConv_R(img, smoothM),times = 100)
testConv