How can I calculate the number of local maxima in a 3D matrix in R?

Question

I am currently working on a number of animal teeth stored as surface meshes in R. Some of the variables I am looking at may be influenced by the number of cusps on the tooth and I would like to try and find a way to dynamically and (relatively) non-subjectively calculate the approximate number of cusps on each mesh.

It occurred to me that it should be possible to roughly estimate the number of cusps by identifying the number of local maxima on each tooth, particularly on the z-axis. Each local maxima in three-dimensional space should correspond to the apex of a cusp.

For example, let's take the Hills sample file from the molaR package.

library(molaR)
DNE3d(DNE(Hills))

Based on this mesh, approximately four distinct peaks are present.

From this, it should be possible to identify the n (4) local minima present on the mesh using the matrix of vertices from Hills$vb, which shows the location of points on the surface of the mesh in the x, y, and z axes. The desired output would be a vector with the positions of the four three-dimensional points a, b, c, and d (or something like this). From there I can take the length of this vector to return the number of local maxima.

However, when I've tried identifying the number of peaks using ggpmisc::find_peaks or pracma::findpeaks, it does not work. I've also tried using the find_peaks function presented here, but the results are also highly variable.

What it looks like is these functions are tailored for finding local optima in a time series or 2D array, but not a 3D array/matrix. Perhaps there is some way to calculate maxima by treating x and y as positional information and z as height, but I am not sure.

Trying to look up this question in Google produced the following code.

find_local_maxima <- function(matrix) {
  local_maxima <- array(FALSE, dim = dim(matrix))
  for (i in 1:dim(matrix)[1]) {
    for (j in 1:dim(matrix)[2]) {
      for (k in 1:dim(matrix)[3]) {
        # Get the current value
        current_value <- matrix[i, j, k]
        
        # Get neighbors (adjust neighborhood as needed)
        neighbors <- c(
          matrix[i-1, j, k], matrix[i+1, j, k],
          matrix[i, j-1, k], matrix[i, j+1, k],
          matrix[i, j, k-1], matrix[i, j, k+1]
        )
        
        # Check if current value is a local maximum
        if (all(current_value >= neighbors)) {
          local_maxima[i, j, k] <- TRUE
        }
      }
    }
  }
  return(local_maxima)
}

Perhaps there is something of value in this code, but I am very skeptical of any result provided by generative AI. When I try applying it to Hills$vb, I get the following error message, even if I restrict the matrix to the first three vectors (Hills$vb[1:3,]).

Error in 1:dim(matrix)[3] : NA/NaN argument

Is there any way to identify the number of maxima in this type of data using R?

@TimG This sort of works, though the number of points recovered is heavily dependent on the k value input. For example when I tried applying it to a real mesh with two cusps I had to change k to 100, whereas 11 returned 66 local maxima. Is there any way to figure out an appropriate value beyond trial and error, especially if this is going to be used in a loop? — user2352714
– user2352714, Commented May 12 at 22:48
This question is too broad for SO, maybe you can ask it at another StackExchange site — Cornelius Roemer
– Cornelius Roemer, Commented May 29 at 10:09

Allan Cameron · Accepted Answer · 2025-05-13 12:55:51Z

The question implies data is in the format of a 3D mesh, like molaR::Hills.

It seems to me that the simplest way forward is to convert the mesh into a 2D raster. We can do this by using the vertex co-ordinates of the mesh, though since they are not in a regular grid we have to interpolate the x, y, z values of the vertices using akima.

Once we have done this, it is straightforward to find local maxima using the focal function from terra. We can then get a data frame of their x, y, z locations:

local_maxima <- function(mesh, window_size = 19) {
  
  df <- setNames(as.data.frame(t(mesh$vb[1:3, ])), c("x","y","z"))
  xi <- seq(min(df$x), max(df$x), length = 200)
  yi <- seq(min(df$y), max(df$y), length = 200)
  
  r <- with(df, akima::interp(x, y, z, xo = xi, yo = yi)) |>
    suppressWarnings() |>
    akima::interp2xyz() |>
    as.data.frame() |>
    terra::rast()
  
  w <- matrix(1, window_size, window_size)
  fmax <- terra::focal(r, w = w, fun = "max", na.rm = TRUE)  
  local_maxima <- r == fmax
  d <- as.data.frame(terra::xyFromCell(r, which(local_maxima[])))
  d$z <- unlist(r[terra::cellFromXY(r, d)])
  d
}

So, using Hills from your example, we can do:

result <- local_maxima(molaR::Hills)

result
#>             x        y        z
#> 1 -13.4188631 14.89494 5.984194
#> 2   0.9787878 14.89494 5.985423
#> 3 -13.4188631  3.84411 5.987416
#> 4   0.9787878  3.84411 5.988645

Of course, the number of cusps is simply

nrow(result)
#> [1] 4

To show these are the correct co-ordinates, we can draw the original mesh and add in the maxima points as red balls:

rgl::shade3d(molaR::Hills, color = "lightblue")
rgl::spheres3d(result, radius = 0.5, color = "red")

Although this works well on the given example, if it doesn't give you the results you want on your actual data, you can try playing around with window_size parameter of local_maxima. A larger window size will identify only the biggest bumps but might miss small cusps, a smaller window will identify more local maxima but runs the risk of false positives if there are irregularities on the tooth surface.

ThomasIsCoding · Accepted Answer · 2025-05-12 22:45:10Z

2

I think you are almost there, but you need to be more cautious when dealing with boundaries (see the indices for listing the neighbors).

The code below is nothing different from yours but just takes additional care of the edges of the 3D matrix, e.g.,

find_local_maxima <- function(matrix) {
  local_maxima <- array(FALSE, dim = dim(matrix))
  for (i in 1:dim(matrix)[1]) {
    for (j in 1:dim(matrix)[2]) {
      for (k in 1:dim(matrix)[3]) {
        # Get the current value
        current_value <- matrix[i, j, k]

        # Get neighbors (adjust neighborhood as needed)
        neighbors <- c(
          matrix[max(1, i - 1), j, k], matrix[min(i + 1, dim(matrix)[1]), j, k],
          matrix[i, max(j - 1, 1), k], matrix[i, min(j + 1, dim(matrix)[2]), k],
          matrix[i, j, max(k - 1, 1)], matrix[i, j, min(k + 1, dim(matrix)[3])]
        )

        # Check if current value is a local maximum
        if (all(current_value >= neighbors)) {
          local_maxima[i, j, k] <- TRUE
        }
      }
    }
  }
  return(local_maxima)
}

Example

Given a dummy 3-D matrix like below

> set.seed(0)

> mat <- array(runif(2 * 3 * 4), dim = c(2, 3, 4))

> print(mat)
, , 1

          [,1]      [,2]      [,3]
[1,] 0.8966972 0.3721239 0.9082078
[2,] 0.2655087 0.5728534 0.2016819

, , 2

          [,1]      [,2]       [,3]
[1,] 0.8983897 0.6607978 0.06178627
[2,] 0.9446753 0.6291140 0.20597457

, , 3

          [,1]      [,2]      [,3]
[1,] 0.1765568 0.3841037 0.4976992
[2,] 0.6870228 0.7698414 0.7176185

, , 4

          [,1]      [,2]      [,3]
[1,] 0.9919061 0.7774452 0.2121425
[2,] 0.3800352 0.9347052 0.6516738

we will see

> find_local_maxima(mat)
, , 1

      [,1]  [,2]  [,3]
[1,] FALSE FALSE  TRUE
[2,] FALSE FALSE FALSE

, , 2

      [,1]  [,2]  [,3]
[1,] FALSE FALSE FALSE
[2,]  TRUE FALSE FALSE

, , 3

      [,1]  [,2]  [,3]
[1,] FALSE FALSE FALSE
[2,] FALSE FALSE FALSE

, , 4

      [,1]  [,2]  [,3]
[1,]  TRUE FALSE FALSE
[2,] FALSE  TRUE FALSE

and

> mat[find_local_maxima(mat)]
[1] 0.9082078 0.9446753 0.9919061 0.9347052

edited May 12 at 22:45

answered May 12 at 22:11

ThomasIsCoding

106k9 gold badges38 silver badges110 bronze badges

3 Comments

user2352714 May 12 at 22:33

I am not sure how to interpret the results here. It seems like the code is reporting if a point is a local minima along the x, y, or z axis, but not in total three-dimensional space. In my case, finding the vertices that represent three-dimensional local maxima are of interest because they are likely apices of cusps on teeth. Am I misunderstanding something?

ThomasIsCoding May 12 at 22:41

@user2352714 it is local maxima, say, comparing one vertex with its surrounding points from 3 dimensions. You can see the output of mat[find_local_maxima(mat)]

user2352714 May 12 at 23:11

I am not sure if this will work for the present purposes. The data is in a matrix of three-dimensional coordinates (x, y, z) and I am trying to find the identity of the coordinates that represent local maxima rather than a singular value. I tried applying this to the Hills mesh above and got the error Error in 1:dim(matrix)[3] : NA/NaN argument, even when I tried transposing the code as @TimG did.

Collectives™ on Stack Overflow

How can I calculate the number of local maxima in a 3D matrix in R?

2 Answers 2

Comments

Example

3 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Example

3 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related