R/main.R

In tcie: Topologically Correct Isosurface Extraction

#============================================================================================
#============================================================================================
#' Returns the dataset size and the isovalue range of a given dataset.
#'
#' @param file_path A string: the path to the nhdr file
#' @return A vector with the dataset size and the minimum and maximum values of a given dataset.
#' @description
#' This function returns the size and isovalue range of a given dataset.
#' @examples
#' DataInfo(system.file("extdata", "f3.nhdr", package = "tcie"))
#' DataInfo(system.file("extdata", "f3.nhdr", package = "tcie"))
#' @export
DataInfo<-function(file_path)
{

  if(nat::is.nrrd(file_path) == FALSE){stop("The input file must to be in nrrd format")}
  f<-nat::read.nrrd(file_path, ReadByteAsRaw = "none")

  #Fills data structure
  size_x <- attr(f,"header")$sizes[1]
  size_y <- attr(f,"header")$sizes[2]
  size_z <- attr(f,"header")$sizes[3]

  min_coord<- which(f == min(f,na.rm=T), arr.ind = T)
  max_coord<- which(f == max(f,na.rm=T), arr.ind = T)

  min<-f[min_coord[1,1],min_coord[1,2],min_coord[1,3]]
  max<-f[max_coord[1,1],max_coord[1,2],max_coord[1,3]]

  cat("Grid size:", size_x, size_y, size_z)
  cat("\n")
  cat("Data range:", min, max)
  cat("\n")

  return(c(size_x, size_y, size_z,min, max))
}
#============================================================================================
#============================================================================================
#' Builds and plots  the mesh representing the desired isosurface
#'
#' @param file_path A string: the path to the nhdr file
#' @param isovalue A number: the value corresponding to the desired isosurface
#' @param verification A boolean: determines whether the grid verification will be performed
#' @param color_mesh A string: the color to rendirind the resulting mesh
#' @param opacity A number: the opacity-level
#' @param new_window A boolean: determines whether a new view window will be opened
#' @return The visualization of the generated mesh.
#' @description
#' This function returns the visualization (rendered by the rgl package) of the mesh generated by the Marching Cubes 33 algotithm.
#' Optionally, the original dataset grid is preprocessed (subdivided in specifics points) to eliminate some configuration
#' which results in non-manifold edges.
#' @examples
#' ManifoldContour(system.file("extdata", "f3.nhdr", package = "tcie"), 0.0, FALSE,"red",1.0, TRUE)
#' ManifoldContour(system.file("extdata", "f3.nhdr", package = "tcie"), 0.0, TRUE,"red",1.0, TRUE)
#' ManifoldContour(system.file("extdata", "f9.nhdr", package = "tcie"), 0.0, TRUE, "blue",1.0, TRUE)
#' @export
ManifoldContour <- function(file_path,isovalue,verification,color_mesh,opacity,new_window)
{
  if(nat::is.nrrd(file_path) == FALSE){stop("The input file must to be in nrrd format")}
  grid_data<-nat::read.nrrd(file_path, ReadByteAsRaw = "none")

  #Fills data structure
  size_x <-attr(grid_data,"header")$sizes[1]
  size_y <-attr(grid_data,"header")$sizes[2]
  size_z <-attr(grid_data,"header")$sizes[3]

  sx <-attr(grid_data,"header")$spacings[1]
  sy <-attr(grid_data,"header")$spacings[2]
  sz <-attr(grid_data,"header")$spacings[3]


  if(verification == TRUE)
  {

  cat("Analizing grid...")

  NewGrid<-AnalizeGrid(grid_data, size_x, size_y, size_z,isovalue)

  divide_info <- NewGrid$Sub_info

  conty <-NewGrid$conty
  contz <-NewGrid$contz

  cat("done!\n")

  if((conty + contz)>0)
  {
    cat(conty + contz, "non-manifold edge would be generated. \n")
    cat("Fixing grid...")

    if((conty != 0)&&(contz !=0))
    {
      divide_info_y <- divide_info[1:conty,]
      divide_info_z <- divide_info[(conty+1):(conty+contz),]
    }
    if(conty == 0)
    {
      divide_info_y <- NULL
      divide_info_z <- divide_info[1:contz,]
    }
    if(contz == 0)
    {
      divide_info_y <- divide_info[1:conty,]
      divide_info_z <- NULL
    }

    if(contz == 1)
      divide_info_z <-matrix(divide_info_z, nrow = 1, ncol = 2)

    if(conty == 1)
      divide_info_y <-matrix(divide_info_y, nrow = 1, ncol = 2)

    Sub_grid<-SubdivideGrid(conty, contz,size_x, size_y, size_z, divide_info_y,divide_info_z,grid_data)

    grid_data <- Sub_grid$grid

    size_y <- size_y + Sub_grid$Conty
    size_z <- size_z + Sub_grid$Contz


    cat("done!\n")
  }
  else
    cat("The grid does not need to be fixed. \n")
}

  MC <-ExtendedRunMC(grid_data, isovalue, size_x, size_y, size_z, sx,sy,sz)

  n_trigs <-MC$ntrig
  if(n_trigs == 0)
  {
    warning("There is no surface associated with the chosen isovalue.")
  }
  else
  {
    verts <- MC$verts
    trigs <- MC$trigs

  mesh<- rgl::tmesh3d(as.vector(t(verts)),as.vector(t(trigs)),homogeneous = FALSE, material = NULL,  normals = NULL, texcoords = NULL)

  if(new_window == TRUE) rgl::open3d()

  rgl::rgl.material(color = c(color_mesh),lit = TRUE, alpha = opacity, ambient = "black",specular = "white",emission = "black",
                   shininess = 40.0,smooth = TRUE,texture = NULL,textype = "rgb",texmipmap = FALSE,texminfilter = "linear",
                    texmagfilter = "linear",texenvmap = FALSE,front = "fill",back = "fill",size = 3.0,lwd = 1.0,fog = TRUE,
                    point_antialias = FALSE,line_antialias = FALSE,depth_mask = TRUE,depth_test = "less")

  rgl::shade3d(mesh)
  }
}

#============================================================================================
#============================================================================================
#' Builds and exports the mesh representing the desired isosurface
#'
#' @param file_path A string: the path to the nhdr file
#' @param isovalue A number: the value corresponding to the desired isosurface
#' @param verification A boolean: determines whether the grid verification will be performed
#' @param export_path A string: the path to a directory to put the mesh generated, following by the file name.
#' @return The file of the mesh in the ply format.
#' @description
#' This function returns the file (ply format) of the mesh generated by the Marching Cubes 33 algorithm.
#' Optionally, the original dataset grid is preprocessed (subdivided in specifics points) to eliminate the configuration that
#' results in non-manifold edges.
#' @examples
#' ExpManifoldContour(system.file("extdata", "f3.nhdr", package = "tcie"), 0.0, FALSE, "mesh.ply")
#' ExpManifoldContour(system.file("extdata", "f3.nhdr", package = "tcie"), 0.0, TRUE, "mesh.ply")
#' ExpManifoldContour(system.file("extdata", "f9.nhdr", package = "tcie"), 0.0, TRUE, "mesh.ply")
#' @export
ExpManifoldContour <- function(file_path,isovalue,verification, export_path)
{

  if(nat::is.nrrd(file_path) == FALSE){stop("The input file must to be in nrrd format")}
  grid_data<-nat::read.nrrd(file_path, ReadByteAsRaw = "none")

  #Fills data structure
  size_x <-attr(grid_data,"header")$sizes[1]
  size_y <-attr(grid_data,"header")$sizes[2]
  size_z <-attr(grid_data,"header")$sizes[3]

  sx <-attr(grid_data,"header")$spacings[1]
  sy <-attr(grid_data,"header")$spacings[2]
  sz <-attr(grid_data,"header")$spacings[3]

  if(verification == TRUE)
  {

    cat("Analizing grid...")

    NewGrid<-AnalizeGrid(grid_data, size_x, size_y, size_z,isovalue)

    divide_info <- NewGrid$Sub_info

    conty <-NewGrid$conty
    contz <-NewGrid$contz

    cat("done!\n")

    if((conty + contz)>0)
    {
      cat(conty + contz, "non-manifold edge would be generated. \n")
      cat("Fixing grid...")

      if((conty != 0)&&(contz !=0))
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- divide_info[(conty+1):(conty+contz),]
      }
      if(conty == 0)
      {
        divide_info_y <- NULL
        divide_info_z <- divide_info[1:contz,]
      }
      if(contz == 0)
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- NULL
      }

      if(contz == 1)
        divide_info_z <-matrix(divide_info_z, nrow = 1, ncol = 2)

      if(conty == 1)
        divide_info_y <-matrix(divide_info_y, nrow = 1, ncol = 2)

      Sub_grid<-SubdivideGrid(conty, contz,size_x, size_y, size_z, divide_info_y,divide_info_z,grid_data)

      grid_data <- Sub_grid$grid

      size_y <- size_y + Sub_grid$Conty
      size_z <- size_z + Sub_grid$Contz


      cat("done!\n")
    }
    else
      cat("The grid does not need to be fixed. \n")
  }

  MC <-ExtendedRunMC(grid_data, isovalue, size_x, size_y, size_z, sx,sy,sz)

  n_trigs <- MC$ntrig

  if(n_trigs == 0)
  {
    warning("There is no surface associated with the chosen isovalue or data subset.")
  }
  else
  {
    verts <- MC$verts
    trigs <- MC$trigs
    mesh<- rgl::tmesh3d(as.vector(t(verts)),as.vector(t(trigs)),homogeneous = FALSE, material = NULL,  normals = NULL, texcoords = NULL)
    Rvcg::vcgClean(mesh, sel = 0, tol = 0, silent = TRUE, iterate = FALSE)
    Rvcg::vcgPlyWrite(mesh,export_path, binary = FALSE,addNormals = TRUE)
  }
}

#============================================================================================
#============================================================================================
#' Builds and plots the mesh representing the desired isosurface of an specific region of the dataset
#' @param file_path A string: the path to the nhdr file
#' @param isovalue A number: the value corresponding to the desired isosurface
#' @param verification A boolean: determines whether the grid verification will be performed
#' @param color_mesh A string: the color to rendirind the resulting mesh
#' @param opacity A number: the opacity-level
#' @param new_window A boolean: determines whether a new view window will open
#' @param range_x A vector contaning the clipping limits of the dataset in the x axis.
#' @param range_y A vector contaning the clipping limits of the dataset in the y axis.
#' @param range_z A vector contaning the clipping limits of the dataset in the z axis.
#' @return The visualization of the generated mesh.
#' @description
#' This function returns the visualization of an specific region of the dataset.
#' @examples
#' Clipping(system.file("extdata","f3.nhdr",package ="tcie"),0,FALSE,"red",1,TRUE,c(2,3),c(2,4),c(2,4))
#' Clipping(system.file("extdata","f3.nhdr",package ="tcie"),0,TRUE,"red",1,TRUE,c(2,3),c(2,4),c(2,4))
#' Clipping(system.file("extdata","f9.nhdr",package ="tcie"),0,TRUE,"blue",1,TRUE,c(1,5),c(2,4),c(3,5))
#' @export
Clipping <- function(file_path, isovalue,verification,color_mesh,opacity, new_window, range_x, range_y, range_z)
{
  if(nat::is.nrrd(file_path) == FALSE){stop("The input file must to be in nrrd format")}
  grid_data<-nat::read.nrrd(file_path, ReadByteAsRaw = "none")

  initial_size_x <-attr(grid_data,"header")$sizes[1]
  initial_size_y <-attr(grid_data,"header")$sizes[2]
  initial_size_z <-attr(grid_data,"header")$sizes[3]

  if(!is.null(attr(grid_data,"header")$spacings[1])) sx <-attr(grid_data,"header")$spacings[1]
  else sx<-1

  if(!is.null(attr(grid_data,"header")$spacings[2])) sy <-attr(grid_data,"header")$spacings[2]
  else sy<-1

  if(!is.null(attr(grid_data,"header")$spacings[3])) sz <-attr(grid_data,"header")$spacings[3]
  else sz<-1

  if((min(range_x) < 1)||(max(range_x))>initial_size_x) {stop(sprintf("The input range_x must a interval contained in [1, %d ]",initial_size_x))}
  if((min(range_y) < 1)||(max(range_y))>initial_size_y) {stop(sprintf("The input range_y must a interval contained in [1, %d]", initial_size_y))}
  if((min(range_z) < 1)||(max(range_z))>initial_size_z) {stop(sprintf("The input range_z must a interval contained in [1, %d]", initial_size_z))}

  grid_data <- grid_data[range_x[1]:range_x[2],range_y[1]:range_y[2],range_z[1]:range_z[2]]

  size_x <-range_x[2]-range_x[1] + 1
  size_y <-range_y[2]-range_y[1] + 1
  size_z <-range_z[2]-range_z[1] + 1

  if(verification == TRUE)
  {

    cat("Analizing grid...")

    NewGrid<-AnalizeGrid(grid_data, size_x, size_y, size_z,isovalue)

    divide_info <- NewGrid$Sub_info

    conty <-NewGrid$conty
    contz <-NewGrid$contz

    cat("done!\n")

    if((conty + contz)>0)
    {
      cat(conty + contz, "non-manifold edge would be generated. \n")
      cat("Fixing grid...")

      if((conty != 0)&&(contz !=0))
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- divide_info[(conty+1):(conty+contz),]
      }
      if(conty == 0)
      {
        divide_info_y <- NULL
        divide_info_z <- divide_info[1:contz,]
      }
      if(contz == 0)
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- NULL
      }

      if(contz == 1)
        divide_info_z <-matrix(divide_info_z, nrow = 1, ncol = 2)

      if(conty == 1)
        divide_info_y <-matrix(divide_info_y, nrow = 1, ncol = 2)

      Sub_grid<-SubdivideGrid(conty, contz,size_x, size_y, size_z, divide_info_y,divide_info_z,grid_data)

      grid_data <- Sub_grid$grid

      size_y <- size_y + Sub_grid$Conty
      size_z <- size_z + Sub_grid$Contz


      cat("done!\n")
    }
    else
      cat("The grid does not need to be fixed. \n")
  }

  MC <-ExtendedRunMC(grid_data, isovalue, size_x, size_y, size_z, sx,sy,sz)

  n_trigs <- MC$ntrig

  if(n_trigs == 0)
  {
    warning("There is no surface associated with the chosen isovalue or data subset.")
  }
  else
  {

    verts <- MC$verts
    trigs <- MC$trigs
  mesh<- rgl::tmesh3d(as.vector(t(verts)),as.vector(t(trigs)),homogeneous = FALSE, material = NULL,  normals = NULL, texcoords = NULL)

  if(new_window == TRUE) rgl::open3d()
  rgl::rgl.material(color = c(color_mesh),lit = TRUE, alpha = opacity, ambient = "black",specular = "white",emission = "black",
                    shininess = 100.0,smooth = TRUE,texture = NULL,textype = "rgb",texmipmap = FALSE,texminfilter = "linear",
                    texmagfilter = "linear",texenvmap = FALSE,front = "fill",back = "fill",size = 3.0,lwd = 1.0,fog = TRUE,
                    point_antialias = FALSE,line_antialias = FALSE,depth_mask = TRUE,depth_test = "less")
  rgl::shade3d(mesh)
  }
}

#============================================================================================
#============================================================================================
#' Builds and exports  the mesh representing the desired isosurface of an specific region of the dataset
#'
#' @param file_path A string: the path to the nhdr file
#' @param isovalue A number: the value corresponding to the desired isosurface
#' @param verification A boolean: determines whether the grid verification will be performed
#' @param range_x A vector contaning the clipping limits of the dataset in the x axis.
#' @param range_y A vector contaning the clipping limits of the dataset in the y axis.
#' @param range_z A vector contaning the clipping limits of the dataset in the z axis.
#' @param export_path A string: the path to a directory to put the mesh generated, following by the file name.
#' @return The file of the mesh in the ply format.
#' @description
#' This function returns the file (ply format) of the mesh of an specific region of the dataset.
#' @examples
#' ExpClipping(system.file("extdata","f3.nhdr",package ="tcie"),0,FALSE,c(2,3),c(2,4),c(2,4),"m.ply")
#' ExpClipping(system.file("extdata","f3.nhdr",package ="tcie"),0,TRUE,c(2,3),c(2,4),c(2,4),"m.ply")
#' ExpClipping(system.file("extdata","f9.nhdr",package ="tcie"),0,TRUE,c(1,5),c(2,4),c(3,5),"m.ply")
#' @export
ExpClipping <- function(file_path,isovalue, verification, range_x, range_y, range_z,export_path)
{
  if(nat::is.nrrd(file_path) == FALSE){stop("The input file must to be in nrrd format")}
  grid_data<-nat::read.nrrd(file_path, ReadByteAsRaw = "none")

  initial_size_x <-attr(grid_data,"header")$sizes[1]
  initial_size_y <-attr(grid_data,"header")$sizes[2]
  initial_size_z <-attr(grid_data,"header")$sizes[3]

  if(!is.null(attr(grid_data,"header")$spacings[1])) sx <-attr(grid_data,"header")$spacings[1]
  else sx<-1

  if(!is.null(attr(grid_data,"header")$spacings[2])) sy <-attr(grid_data,"header")$spacings[2]
  else sy<-1

  if(!is.null(attr(grid_data,"header")$spacings[3])) sz <-attr(grid_data,"header")$spacings[3]
  else sz<-1

  if((min(range_x) < 1)||(max(range_x))>initial_size_x) {stop(sprintf("The input range_x must a interval contained in [1, %d ]",initial_size_x))}
  if((min(range_y) < 1)||(max(range_y))>initial_size_y) {stop(sprintf("The input range_y must a interval contained in [1, %d]", initial_size_y))}
  if((min(range_z) < 1)||(max(range_z))>initial_size_z) {stop(sprintf("The input range_z must a interval contained in [1, %d]", initial_size_z))}

  grid_data <- grid_data[range_x[1]:range_x[2],range_y[1]:range_y[2],range_z[1]:range_z[2]]

  size_x <-range_x[2]-range_x[1] + 1
  size_y <-range_y[2]-range_y[1] + 1
  size_z <-range_z[2]-range_z[1] + 1

  if(verification == TRUE)
  {

    cat("Analizing grid...")

    NewGrid<-AnalizeGrid(grid_data, size_x, size_y, size_z,isovalue)

    divide_info <- NewGrid$Sub_info

    conty <-NewGrid$conty
    contz <-NewGrid$contz

    cat("done!\n")

    if((conty + contz)>0)
    {
      cat(conty + contz, "non-manifold edge would be generated. \n")
      cat("Fixing grid...")

      if((conty != 0)&&(contz !=0))
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- divide_info[(conty+1):(conty+contz),]
      }
      if(conty == 0)
      {
        divide_info_y <- NULL
        divide_info_z <- divide_info[1:contz,]
      }
      if(contz == 0)
      {
        divide_info_y <- divide_info[1:conty,]
        divide_info_z <- NULL
      }

      if(contz == 1)
        divide_info_z <-matrix(divide_info_z, nrow = 1, ncol = 2)

      if(conty == 1)
        divide_info_y <-matrix(divide_info_y, nrow = 1, ncol = 2)

      Sub_grid<-SubdivideGrid(conty, contz,size_x, size_y, size_z, divide_info_y,divide_info_z,grid_data)

      grid_data <- Sub_grid$grid

      size_y <- size_y + Sub_grid$Conty
      size_z <- size_z + Sub_grid$Contz


      cat("done!\n")
    }
    else
      cat("The grid does not need to be fixed. \n")
  }

  MC <-ExtendedRunMC(grid_data, isovalue, size_x, size_y, size_z, sx,sy,sz)

  n_trigs <- MC$ntrig

  if(n_trigs == 0)
  {
    warning("There is no surface associated with the chosen isovalue or data subset.")
  }
  else
  {
    verts <- MC$verts
    trigs <- MC$trigs
    mesh<- rgl::tmesh3d(as.vector(t(verts)),as.vector(t(trigs)),homogeneous = FALSE, material = NULL,  normals = NULL, texcoords = NULL)
    Rvcg::vcgPlyWrite(mesh,export_path, binary = FALSE,addNormals = TRUE)
  }

}

#============================================================================================
#============================================================================================
# A utility function to find the subset of an element i
find_<-function(parent,i)
{
  if (parent[i] == -1) return (i)
  else(find_(parent, parent[i]))
}

# A utility function to do union of two subsets
union_<-function(parent, x, y)
{
  xset <- find_(parent, x)
  yset <- find_(parent, y)
  parent[xset] <- yset

  return(parent)
}
#' Returns the  Betti numbers of a given mesh.
#'
#' @param file_path A string: the path to a triangular mesh file in the ply format
#' @description
#' This function returns the Betti numbers b0, b1 and b2, which represents the number of connected components,
#' the number of independent tunnelsand the number of closed regions in space, respectively.
#' The function implementation followsthe algorithms described by Konkle, Moran, Hamann, and Joy in the work
#' Fast Methods for Computing Isosur-face Topology with Betti Numbers.
#' @examples
#' BettiNumbers(system.file("extdata", "m3.ply", package = "tcie"))
#' @export
BettiNumbers<-function(file_path)
{
  mesh <- geomorph::read.ply(file_path,ShowSpecimen=F, addNormals = FALSE)

  n_v<- Rvcg::nverts(mesh)
  n_f<- Rvcg::nfaces(mesh)
  edges <-Rvcg::vcgGetEdge(mesh, unique = TRUE)

  parent <- rep(-1, n_v)

  b0 <- n_v

  n_e_b <- 0
  board_verts <-c()

  for(t in 1: n_f)
  {
    v <- mesh$it[,t]

    cont_face_neighbor <- 0
    for(i in 1:3)
    {
      j <- i + 1
      if(i == 3) j <- 1

      x <- find_(parent, v[i])
      y <- find_(parent, v[j])

      if(x != y)
      {
        b0 <-b0 -1
        parent <- union_(parent, x, y)
      }

      v1 <- v[i]
      v2 <- v[j]

      cont_edge_neighbor <-0

      for(l in 1:3)
      {
        faces_v1 <- which(mesh$it[l,] == v1)
        l_v1_v2  <- which(mesh$it[,faces_v1] == v2)

        if((length(faces_v1)!=0)&&(length(l_v1_v2)!=0))
        {
          for(i in 1:length(l_v1_v2))
          {
            id_neighbor<- ceiling(l_v1_v2[i]/3)
            if(faces_v1[id_neighbor] != t)
            {
              cont_face_neighbor <- cont_face_neighbor + 1
              cont_edge_neighbor <- 1
            }
          }
        }
      }

      if(cont_edge_neighbor == 0) board_verts <- c(board_verts, v1)

      if((i == 3)&&(cont_face_neighbor == 0)) n_e_b<- n_e_b + 3
      if((i == 3)&&(cont_face_neighbor == 1)) n_e_b<- n_e_b + 2
      if((i == 3)&&(cont_face_neighbor == 2)) n_e_b<- n_e_b + 1
    }
  }

  n_b <- 0
  n_b_verts <- length(board_verts)
  parent_board_vert <- c()

  if(n_b_verts == 1)  n_b <- 1
  if(n_b_verts > 1)
  {
    for(i in 1: n_b_verts)
    {
      options(expressions = 100000)
      parent_board_vert[i]<-find_(parent, board_verts[i])
    }

    n_b <-length(unique(parent_board_vert))
  }
  n_e <- length(edges[,1])
  e_c <- n_v - n_e + n_f

  b2 <- b0 - n_b
  b1 <- b0 + b2 - e_c

  return(c(b0, b1, b2))
}
#============================================================================================
#============================================================================================
SubdivideGrid<-function(conty, contz, size_x, size_y, size_z, subdivide_info_y, subdivide_info_z,f)
{

  if(conty != 0)
  {
    aux <-subdivide_info_y[,1]
    index_y <- unique(sort(aux, decreasing = TRUE))

    conty <- length(index_y)
  }

  if(contz != 0)
  {
    aux <-subdivide_info_z[,1]

    index_z <- unique(sort(aux, decreasing = TRUE))

    contz <- length(index_z)
  }

  f_ <- array(dim=c(size_x,size_y+conty,size_z+contz))


  aux_conty <- conty
  aux_contz <- contz

  if(contz !=0)
  {

    for(c in 1 : contz)
    {

      aux_contz <- aux_contz -1
      k <- index_z[c] + 1
      id<-which(subdivide_info_z[,1] == index_z[c])
      h <- subdivide_info_z[id[1],2]

      for (i in 1:(size_x))
        for (j in 1:(size_y))
        {
          a <-f[i,j,k]
          b <-f[i,j,(k+1)]
          x <- a + h*(b-a)

          data <-f[i,j,]
          f_[i,j,]<-c(data[1:k],x,data[(k+1):size_z],rep(0, aux_contz))

        }
    }
  }

  if(conty !=0)
  {
    for(c in 1 : conty)
    {
      aux_conty <- aux_conty -1
      j <- index_y[c] + 1
      id<-which(subdivide_info_y[,1] == index_y[c])
      h <- subdivide_info_y[id[1],2]

      for (k in 1:(size_z))
        for (i in 1:(size_x))
        {
          a <-f[i,j,k]
          b <-f[i,j + 1,k]
          x <- a + h*(b-a)

          if(c == 1)
            {
              data <- f[i,,k]
              f_[i,,k]<-c(data[1:j],x,data[(j+1):size_y],rep(0, aux_conty))
            }
          else
            {
              data <-f_[i,,k]
              f_[i,,k]<-c(data[1:j],x,data[(j+1):(size_y + conty - 1)])
            }
        }
    }
  }
  return(list(grid = f_, Conty = conty, Contz = contz))
}