isolines do not pass through the grid points that bear the requested value

[jldupouey]

Isoband is a nice anf fast package. But I wonder why segments generated by the isolines function do not pass through the grid points bearing the resquested value (given in levels parameter). In the following example:

v<-c(1,2,0,
     1,1,2,
     2,2,2)
m<-matrix(v,nrow=3,byrow=TRUE)

(iso=isolines(1:3,3:1,m,levels=1))
# $`1`
# $`1`$x
# [1] 2.5 3.0
#
# $`1`$y
# [1] 3.0 2.5
# 
# $`1`$id
# [1] 1 1
# 
# attr(,"class")
# [1] "isolines" "iso"

the "1" values in position [1,3], [1,2] and [2,2] of the matrix do not belong to any segment, although they have the requested value. In mathematical applications, or in spatial analysis, a point bearing the value normally belongs to the isoline of this value. The present algorithm needs additional post-processing to add the missing segments. Of course, these segments could in some cases be of zero length (e.g., when there is only one point in the grid bearing the value of the isoline), but that is better than no information at all. Or perhaps I missed something in the package?

[clauswilke]

I don't remember the details, but it's quite possible that I implemented it the other way round, that the value itself is excluded. If you shift the level by a tiny amount I believe you get the result you were looking for.

A simple workaround would be to shift the levels by a tiny amount, smaller than the accuracy of your data.

library(isoband)
library(grid)

make_plot <- function(m, level) {
  df_lines <- isolines((1:ncol(m))/(ncol(m)+1), (nrow(m):1)/(nrow(m)+1), m, level)[[1]]
  g <- expand.grid(x = (1:ncol(m))/(ncol(m)+1), y = (nrow(m):1)/(nrow(m)+1))
  grid.newpage()
  grid.points(g$x, g$y, default.units = "npc", pch = 19, size = unit(0.5, "char"))
  grid.polyline(df_lines$x, df_lines$y, df_lines$id)
}

v <- c(1, 2, 0,
       1, 1, 2,
       2, 2, 2)
m <- matrix(v, nrow=3, byrow=TRUE)
make_plot(m, 1)

make_plot(m, 1.0001)

^{Created on 2021-12-28 by the reprex package (v2.0.0)}

[clauswilke]

Actually, no, I definitely use a greater than or equals sign in the code: https://github.com/wilkelab/isoband/blob/2478d9768007ce40f51ea8a6788837041ac5a748/src/isoband.cpp#L1458

It only seems to be a problem at the edges of the grid:

library(isoband)
library(grid)

make_plot <- function(m, level) {
  df_lines <- isolines((1:ncol(m))/(ncol(m)+1), (nrow(m):1)/(nrow(m)+1), m, level)[[1]]
  g <- expand.grid(x = (1:ncol(m))/(ncol(m)+1), y = (nrow(m):1)/(nrow(m)+1))
  grid.newpage()
  grid.points(g$x, g$y, default.units = "npc", pch = 19, size = unit(0.5, "char"))
  grid.polyline(df_lines$x, df_lines$y, df_lines$id)
}

v <- c(0, 1, 2, 0,
       0, 1, 1, 2,
       0, 2, 2, 2)
m <- matrix(v, nrow=3, byrow=TRUE)
make_plot(m, 1)

isolines(1:ncol(m), nrow(m):1, m, 1)[[1]]
#> $x
#> [1] 2.0 2.0 1.5 3.5 4.0
#> 
#> $y
#> [1] 3.0 2.0 1.0 3.0 2.5
#> 
#> $id
#> [1] 1 1 1 2 2

^{Created on 2021-12-28 by the reprex package (v2.0.0)}

[jldupouey]

Thanks for the answer. Shifting the value of the isoline level by a tiny amount does not work. It depends on the context of the isolines.

library(isoband)  

v<-c(1,2,0,
     1,1,2,
     2,2,2)
m<-matrix(v,nrow=3,byrow=TRUE)

# works with small amount added:

# (iso=isolines(1:3,3:1,m,1+1e-5))
# $`1.00001`
# $`1.00001`$x
# [1] 1.000010 2.000000 2.000010 2.000000 1.000000 2.499995 3.000000

# $`1.00001`$y
# [1] 3.000000 2.000010 2.000000 1.999990 1.999990 3.000000 2.499995

# $`1.00001`$id
# [1] 1 1 1 1 1 2 2

# attr(,"class")
# [1] "isolines" "iso"     

# reverse 0 and 2 in the m matrix
# => no longer works when adding a small amount:

v<-c(1,0,2,
     1,1,0,
     0,0,0)
m<-matrix(v,nrow=3,byrow=TRUE)

(iso=isolines(1:3,3:1,m,1+1e-5))
# $`1.00001`
# $`1.00001`$x
# [1] 2.500005 3.000000

# $`1.00001`$y
# [1] 3.000000 2.500005

# $`1.00001`$id
# [1] 1 1

# attr(,"class")
# [1] "isolines" "iso" 

One should have to visually inspect or programmatically analyze the grid before deciding wether an epsilon should be added or not. And this can vary within a grid. If you try with:

v<-c(1,0,2,
     0,2,2,
     0,2,1)

at levels=1, adding or not adding a small amount to level, you will get either of the two extreme "1"s, but never both, which is an inconsistant behaviour for a simple search of isolines.

It is not only a problem at the edge of the grid. The following grid does not identify the central "1" at levels=1 either:

v<-c(2,2,2,
     2,1,2,
     2,2,2)

I work with large temperature grids, which take a limited number of values, and where it is common to look for isolines passing in between but also exactly on the grid points.

[clauswilke]

I think it's important to be precise about what the problem is. In https://github.com/wilkelab/isoband/issues/29#issuecomment-1002406726, I mentioned that I think the problem occurs only at the boundaries of the grid. So let's first establish whether you agree with this or not. Do you have a case where the problem occurs in the interior?

[jldupouey]

Thanks. Here is a simple case (same as in my previous comment) where an interior point is not identified as belonging to the isoline:

v<-c(2,2,2,
     2,1,2,
     2,2,2)
m<-matrix(v,nrow=3,byrow=TRUE)

(iso=isolines(1:3,3:1,m,1))

# $`1`
# $`1`$x
# numeric(0)
# 
# $`1`$y
# numeric(0)
# 
# $`1`$id
# integer(0)
# 
# attr(,"class")
# [1] "isolines" "iso"

Probably I do not clearly understand what you mean, because it was already the case in the previous example you provided. In:

library(isoband)
library(grid)

make_plot <- function(m, level) {
  df_lines <- isolines((1:ncol(m))/(ncol(m)+1), (nrow(m):1)/(nrow(m)+1), m, level)[[1]]
  g <- expand.grid(x = (1:ncol(m))/(ncol(m)+1), y = (nrow(m):1)/(nrow(m)+1))
  grid.newpage()
  grid.points(g$x, g$y, default.units = "npc", pch = 19, size = unit(0.5, "char"))
  grid.polyline(df_lines$x, df_lines$y, df_lines$id)
}

v <- c(0, 1, 2, 0,
       0, 1, 1, 2,
       0, 2, 2, 2)
m <- matrix(v, nrow=3, byrow=TRUE)
make_plot(m, 1)

the "1" in row 2 and column 3 of the matrix m is missed by isolines, not recognized as belonging to the isoline levels=1, while the two other "1" are correctly included.

[clauswilke]

This case is definitely being handled correctly:

v <- c(0, 1, 2, 0,
       0, 1, 1, 2,
       0, 2, 2, 2)

The isoline at level 1 is the boundary between values < 1 and values >= 1. The middle 1 is not at the boundary. It is in the interior.

Whether degenerate cases like this one should create any output is debatable.

v<-c(2,2,2,
     2,1,2,
     2,2,2)

It would be an isoline without any extent. I may even have written code to remove those because there may have been issues with downstream tools. I don't fully remember. You'd have to audit the code line-by-line to figure this out.

I also think I've figured out now what happens in a case like this one:

v <- c(1, 2, 0,
       1, 1, 2,
       2, 2, 2)

You may think there should be a vertical line at position 1, but that assumes that values to the left of those 1s are lower, and we don't know that because there is no data. Isoband only places lines when it is certain that there actually is a boundary. So I'm now convinced that this case is handled correctly also.

Note that when you calculate isobands rather than isolines, this case correctly shows an open boundary on the left. The points are included in the isosurface but we can't know whether there's a boundary or not.

library(isoband)

v <- c(1, 2, 0,
       1, 1, 2,
       2, 2, 2)
m <- matrix(v, nrow=3, byrow=TRUE)

plot_iso(m, 1, 1.4)

^{Created on 2021-12-29 by the reprex package (v2.0.0)}

[clauswilke]

Actually, the degenerate case with a single value at your level is also handled correctly. The case with the 1 surrounded by 2s is again a case where there isn't actually any boundary. It's entirely interior. When you surround the 1 by 0s you get an isoline without extent, as expected. The code to remove those that I vaguely remember writing is probably further downstream, in the conversion to sf objects.

library(isoband)

v <- c(0, 0, 0,
       0, 1, 0,
       0, 0, 0)
m <- matrix(v, nrow=3, byrow=TRUE)
isolines(1:ncol(m), nrow(m):1, m, 1)[[1]]
#> $x
#> [1] 2 2 2 2 2
#> 
#> $y
#> [1] 2 2 2 2 2
#> 
#> $id
#> [1] 1 1 1 1 1

v <- c(2, 2, 2,
       2, 1, 2,
       2, 2, 2)
m <- matrix(v, nrow=3, byrow=TRUE)
isolines(1:ncol(m), nrow(m):1, m, 1)[[1]]
#> $x
#> numeric(0)
#> 
#> $y
#> numeric(0)
#> 
#> $id
#> integer(0)

^{Created on 2021-12-29 by the reprex package (v2.0.0)}

[jldupouey]

Thank you, I now understand the choices you made for your package, which are the classic choices of the marching squares algorithm. But you are right, we have to be precise, starting with the definition of isolines. According to Wenger (Isosurfaces. Geometry, topology & algorithms, CRC Press, 2013), given a scalar (and continuous) field φ : ℝ² → ℝ and a constant σ ∈ ℝ, the isoline at level σ is the set of points {x : φ(x) = σ}, i.e. φ^-1(σ). Referring to this definition, with which I agree, it is difficult to justify that:

v<-c(1,0,
     0,1)
m<-matrix(v,nrow=2,byrow=TRUE)

(iso=isolines(1:2,2:1,m,1))

correctly identifies the diagonal 1-1 as an isoline, but that :

v<-c(1,2,
     2,1)
m<-matrix(v,nrow=2,byrow=TRUE)

(iso=isolines(1:2,2:1,m,1))

rejects it.

I understand clearly that it is the basic marching squares algorithm that leads to this inconsistent and rather troublesome result for some isoline search applications. This algorithm is widely used, but it seems to me a bit unfortunate that it takes so little account of this problem. As it is currently most often implemented, it looks more adapted to the search of regions, bands, than of set of points. Maybe there are implementations that do? If there are not and if I find some time, I will try to modify the algorithm to take into account these cases of "degeneracy". As I mentioned above, it is not so much "degeneracy" as in some integer fields of values, it is very common to have isolines passing through a significant proportion of the grid points.

In fact, we don't have 16 types of squares to take into account anymore, as in the basic marching squares algorithm, but 3⁴ = 81 cases (a value at a point of the grid can be less than, equal to or greater than σ). And there are several possible ways to build the isolines for some of these 81 cases. 65 of these cases include points equal to σ, which are new cases to examine compared to the 16 classical cases of the marching squares algorithm. But because of the symmetries, there are finally much less new configurations to analyse, 9 at maximum.

If anyone has thought about this, or knows of an algorithm that implements a solution, I'm interested.