R/LKE_side_cv.R

In LuXiaoEei/image-denosing-by-single-sided-local-linear-kernel-estimation: image-denosing-by-single-sided-local-linear-kernel-estimation

### local kernel estimator 局部线性核估计 针对边缘采用分割支撑集的方式
#requireNamespace('data.table')
LKE_side_cv <- function(m,n,h,
                     data,
                     SearchPoints0,
                     kernel='bivariate kernal',
                     xmin=-1000,
                     xmax=1000,
                     ymin=-1000,
                     ymax=1000,
                     scale){
  #require(data.table)
  data <- as.data.table(data)
  Area <- SearchPoints(m,n,SearchPoints0,xmin,xmax,ymin,ymax)
  Area <-merge(Area,data,by=c("x","y"))
  #Area <- Area[!(Area$x==m&Area$y==n),]
  #if(kernel=='bivariate kernal'){
  Area$x <-  Area$x-m
  Area$y <-  Area$y-n
  Area$K <- BivKernal(x=Area[['x']]/scale,y=Area[['y']]/scale,h=h/scale)
  #}
  Area <- as.data.frame(Area)
  beta <- WLS(Area[,-ncol(Area)],W=diag(Area$K))#拟合值，梯度

  beta <- rbind(beta,WRMS(Z = Area$z,
                          a=beta[1,1],
                          x_grad = beta[2,1],
                          X = Area$x,
                          y_grad = beta[3,1],
                          Y = Area$y,
                          W = Area$K))

  Area$sep <- beta[2,1]*Area[,1]+beta[3,1]*Area[,2]#计算在切线上下
  Area1 <- subset(Area,sep>=0)
  Area2 <- subset(Area,sep<=0)#支撑集分割
  # beta1 <- NA
  # beta2 <- NA
  # if (nrow(Area1)>1){
  #   beta1 <- WLS(Area1[,c('x','y','z')],W=diag(Area1$K))
  # }
  # if (nrow(Area2)>1){
  #   beta2 <- WLS(Area2[,c('x','y','z')],W=diag(Area2$K))
  # }
  # if(nrow(Area1)<=1|nrow(Area2)<=1){
  #   beta1 <- beta
  #   beta2 <- beta
  # }else{
  #   beta1 <- WLS(Area1[,c('x','y','z')],W=diag(Area1$K))
  #   beta2 <- WLS(Area2[,c('x','y','z')],W=diag(Area2$K))
  # }

  #一以下对于报错一方面是某类只有一个点，或者某类的所有点的某个维度的坐标相等导致矩阵不可逆，这里为了计算的效率，采用NA，这不影响最后的结果。
  beta1 <- tryCatch({
    WLS(Area1[,c('x','y','z')],W=diag(Area1$K))
  },error=function(e){
    #WLS(Area2[,c('x','y','z')],W=diag(Area2$K))
    matrix(c(NA,NA,NA))
  })

  beta1 <- rbind(beta1,WRMS(Z = Area1$z,
                            a = beta1[1,1],
                            x_grad = beta1[2,1],
                            X = Area1$x,
                            y_grad = beta1[3,1],
                            Y = Area1$y,
                            W = Area1$K))

  beta2 <-tryCatch({
    WLS(Area2[,c('x','y','z')],W=diag(Area2$K))
  },error=function(e){
    #beta1
    matrix(c(NA,NA,NA))
  })

  beta2 <- rbind(beta2,WRMS(Z = Area2$z,
                            a=beta2[1,1],
                            x_grad = beta2[2,1],
                            X = Area2$x,
                            y_grad = beta2[3,1],
                            Y = Area2$y,
                            W = Area2$K))
  rownames(beta) <- NULL
  rownames(beta1) <- NULL
  rownames(beta2) <- NULL

  result <- data.frame(beta,beta1,beta2)

  return(cbind(result[1,],result[4,]))#返回三个估计值和对应的梯度和WRMS
}