[R-sig-Geo] Tessellation error in spatstat

Adrian.Baddeley at csiro.au Adrian.Baddeley at csiro.au
Tue Oct 19 03:25:45 CEST 2010


Richard Mort wrote (on R-help):

...................
> I have confocal images 512x512 from which I have extracted x,y positions of the coordiates 
> of labelled cells exported from ImageJ as a.csv file. I also have images that define an 
> underlying pattern in the tissue defined as areas of different pixel values 0 or 255 
> (also 512x512) I've exported these images as .txt files. I'd like to use spatstat to 
> look at how the individual cells plot onto these patterns.

> #I can import my x,y data and do basic stats with spatstat
>  by doing the following:

> library(spatstat)
> data01 <- read.csv("250810.csv", skip=1) 
> x <- data01[[2]]
> y <- data01[[3]]
> data02 <- ppp(x, y, xrange=c(0,512), yrange=c(0,512)) 
> unitname(data02) <- c("pixel", "pixels") 
> plot(data02)

> #I can also import my text images and convert them to a tessellation
> using the following:

> a <- read.table("FOLLICLES.txt")
> win <- owin(c(0,512), c(0,512))
> unitname(win) <- c("pixel", "pixels")
> b <- as.matrix(a, xrange=c(0,512), yrange=c(0,512)) 

      *** Error here ^^^ 

> unitname(b) <- c("pixel", "pixels")
> c <- im(b, xcol=seq(ncol(b)), yrow=seq(nrow(b)), unitname="pixels") 

      *** Error here ^^^ 

> d <- tess(image = c)
> plot(d)

> #I can then plot my ppp pattern on top of my tessellation using:
> plot(d)
> plot(data02, add=TRUE)

> #So far so good, but when I run for example:
> qb <- quadratcount(data02, tess = d)

> #I get the following error:

> Warning message:
> In quadratcount.ppp(data02, tess = d) :
>   Tessellation does not contain all the points of X

> #When I investigate further the following is returned for each object:

> > d
> Tessellation
> Tessellation is determined by a factor-valued image with 2 levels window: binary image mask
> 512 x 512 pixel array (ny, nx)
> enclosing rectangle: [0.5, 512.5] x [0.5, 512.5] pixels  
> > data02
>  planar point pattern: 254 points
> window: rectangle = [0, 512] x [0, 512] units

> #I don't understand why my tessellation is a different size but even
> taking this into account my ppp should plot inside the window 
> and does when i do plot(data02, add=TRUE). 
> I've spent ages playing around with redefining the windows but I must be missing something (probably obvious).
................................

The output shows that the enclosing rectangle for 'd' is [0.5, 512.5] x [0.5, 512.5] 
while the point pattern 'data02' is enclosed in the rectangle [0, 512] x [0, 512]. 
If there are data points which have x-coordinates between 0 and 0.5 (for example) then these points 
will lie outside the domain of the tessellation 'd', causing an error, although the problem 
would probably not be visible in a plot. This appears to have happened with the pattern 'data02'.

The problem is with the following commands

> b <- as.matrix(a, xrange=c(0,512), yrange=c(0,512)) 
> unitname(b) <- c("pixel", "pixels")
> c <- im(b, xcol=seq(ncol(b)), yrow=seq(nrow(b)), unitname="pixels") 

'as.matrix' is part of the base R system. It doesn't have arguments 'xrange', 'yrange'
and these will be ignored (silently!). So b is just an ordinary matrix.
Then in the call to 'im' you have specified the coordinates of the pixel CENTRES to be
(1,1), ...., (512, 512) and the pixels are squares of side 1, so the pixels cover the
rectangle [0.5, 512.5] x [0.5, 512.5]. 

If I understand correctly you want 'b' to be a 512 x 512 matrix, 
and 'c' to be an image defined on the rectangle [0, 512] x [0,512].
To do that, it's probably easiest to use 'as.im':

     b <- matrix(b, nrow=512)
     c <- as.im(b, owin(c(0,512),c(0,512)))

or more portably
 
     c <- as.im(b, as.owin(data02))

the last line ensures that the windows of 'data02' and 'c' are identical.

It's probably wise to avoid using 'c' as the name of a dataset, 
because 'c' is also a function name

regards
Adrian Baddeley



More information about the R-sig-Geo mailing list