[R] Enumerating unit cell traversals Re: Interpolate? a line

Rainer M Krug r.m.krug at gmail.com
Wed Sep 15 14:58:46 CEST 2010


-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 15/09/10 14:46, Michael Bedward wrote:
> Nice posts David. I like the way that 'simple' problems such as this
> one give rise to an interesting assortment of approaches and gotchas
> :)
> 
> Michael
> 
> On 15 September 2010 22:41, David Winsemius <dwinsemius at comcast.net> wrote:
>>
>> On Sep 15, 2010, at 7:55 AM, David Winsemius wrote:
>>
>>>
>>> On Sep 15, 2010, at 7:24 AM, David Winsemius wrote:
>>>
>>>> Replacing context:
>>>>
>>>>>> Hello everyone.
>>>>>> I have created a 100*100 matrix in R.
>>>>>> Let's now say that I have a line that starts from (2,3) point and ends
>>>>>> to the
>>>>>> (62,34) point. In other words this line starts at cell (2,3) and ends
>>>>>> at cell
>>>>>> (62,34).
>>>>>>
>>>>>> Is it possible to get by some R function all the matrix's cells that
>>>>>> this line
>>>>>> transverses?

This is a known problem in computer graphics and the algorythm used to
identify these cells is the "Bresenham's Line Drawing Algorithm"

You can check out one version at

http://www.falloutsoftware.com/tutorials/dd/dd4.htm

or

http://en.wikipedia.org/wiki/Bresenham's_line_algorithm

Cheers,

Rainer

>>>>>>
>>>>>> I would like to thank you for your feedback.
>>>>>>
>>>>>> Best Regards
>>>>>> Alex
>>>>
>>>> On Sep 15, 2010, at 6:52 AM, Michael Bedward wrote:
>>>>
>>>>> Hello Alex,
>>>>>
>>>>> Here is one way to do it. It works but it's not pretty :)
>>>>
>>>> If you want an alternative, consider that produces the Y cell indices
>>>> (since the x cell indices are already 2:62):
>>>>
>>>>> linefn <- function(x) 3+((34-3)/(62-2)) *(x-2)
>>>>> findInterval(linefn(2:62), 3:34)
>>>> [1]  1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  9  9 10 10 11 11 12
>>>> 12 13 13 14
>>>> [28] 14 15 15 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26
>>>> 26 27 27 28
>>>> [55] 28 29 29 30 30 31 32
>>>> # that seems "off" by two
>>>>> linefn(62)
>>>> [1] 34
>>>>> linefn(2)
>>>> [1] 3 # but that checks out and I realized those were just indices for
>>>> the 3:34 findInterval vector
>>>>
>>>>> (3:34)[findInterval(linefn(2:62), 3:34)]
>>>> [1]  3  3  4  4  5  5  6  6  7  7  8  8  9  9 10 10 11 11 12 12 13 13 14
>>>> 14 15 15 16
>>>> [28] 16 17 17 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28
>>>> 28 29 29 30
>>>> [55] 30 31 31 32 32 33 34
>>>>
>>>> ( no rounding and I think the logic is clearer.)
>>>
>>> But I also realized it didn't enumerate all the the cells were crossed
>>> either, only indicating which cell was associated with an integer value of
>>> x. Also would have even more serious problems if the slope were greater than
>>> unity. To enumerate the cell indices that were crossed, try:
>>>
>>> unique( floor( cbind( seq(2,62, by=0.1), linefn(seq(2,62, by=0.1)) ) )  )
>>>     [,1] [,2]
>>> [1,]    2    3
>>> [2,]    3    3
>>> [3,]    4    4
>>> [4,]    5    4
>>> [5,]    5    5
>>> [6,]    6    5
>>> [7,]    7    5
>>> [8,]    7    6
>>> snipping interior results
>>> [83,]   58   32
>>> [84,]   59   32
>>> [85,]   60   32
>>> [86,]   60   33
>>> [87,]   61   33
>>> [88,]   62   34
>>>
>>> That could probably be passed to rect() to illustrate (and check logic):
>>>
>>> rect(cellidxs[,1], cellidxs[,2], cellidxs[,1]+1, cellidxs[,2]+1,
>>> col="red")
>>>
>>> #redraw line :
>>> lines(2:62, 3+(34-3)/(62-2)*(0:60))
>>
>> And then I got to wondering if every 0.1 was sufficient to catch all the
>> corners and discovered I could identify 3 more cell traversals with by=0.01
>>
>>> cellid2 <-unique( floor(cbind(seq(2,62, by=0.01), linefn(seq(2,62,
>>> by=0.01) )) ) )
>>> NROW(cellid2) # 91 cells
>>> rect(cellid2[,1], cellid2[,2], cellid2[,1]+1, cellid2[,2]+1, col="blue")
>>> rect(cellidxs[,1], cellidxs[,2], cellidxs[,1]+1, cellidxs[,2]+1,
>>> col="red")
>>> lines(2:62, 3+(34-3)/(62-2)*(0:60))
>>
>> (Trying by=0.001 did not change the count, but did take longer)
>>
>> --
>> David.
>>
>>>> -
>>>> David.
>>>>
>>>>>
>>>>> interp <- approx(c(2, 62), c(3, 34), method="linear", xout=2:62)
>>>>> m <- matrix(c(interp$x, round(interp$y)), ncol=2)
>>>>> tie <- m[,2] == c(-Inf, m[-nrow(m),2])
>>>>> m <- m[ !tie, ]
>>>>>
>>>>> You might want to examine the result like this...
>>>>>
>>>>> plot(m)  # plots points
>>>>> lines(c(2,26), c(3, 34))  # overlay line for comparison
>>>>
>>>> you can add a grid with
>>>> abline(v=2:62, h=3:34)
>>>>>
>>>>> Michael
>>>
>>
>>
>> David Winsemius, MD
>> West Hartford, CT
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


- -- 
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Tel:        +33 - (0)9 53 10 27 44
Cell:       +27 - (0)8 39 47 90 42
Fax (SA):   +27 - (0)8 65 16 27 82
Fax (D) :   +49 - (0)3 21 21 25 22 44
Fax (FR):   +33 - (0)9 58 10 27 44
email:      Rainer at krugs.de

Skype:      RMkrug
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAkyQwwUACgkQoYgNqgF2egrDSACeJPgbW5tQtqh7RyUMIezeNNL9
1+QAn2Vgow6lgeCuzZ2N2ysgxVTxmNm4
=NFGN
-----END PGP SIGNATURE-----



More information about the R-help mailing list