A shortcut to Bresenham's algo?

General Usage
1

Hi,

I have a set mycurves of broken lines in 3d. I have a grid specified by 3 ranges rx,ry,rz.

For each curve in mycurve, I want to tag the cells in the grid by which the curve passes. Hence, the cells hold the value 0 or 1. In short, this can be done using the Bresenham algorithm.

Now, I want to do this for all curves so I increment the value of the cell when a curve passes through it. The max value hold by a cell is thus length(mycurves).

Before I embark into coding the Bresenham algo in 3d, is there a shortcut for doing this (like using Fhist.jl for example)?

Thank you for your advices!

Here is the set of curves for example:

using GeometryBasics

function mycurve()
    lines = Point3{Float64}[]
    X0 = Point3(1.,0,0)
    lines = [X0]
    dt = 0.001
    for i=1:100
        X0 = @. X0 + 10dt * (-X0) + 0.1 * sqrt(dt) * randn()
        push!(lines, X0)
    end

    lines
end

# lines
mycurves = [myline() for i=1:100]
2

Which one do you want?

As far as I understand your helpful wiki link, Bresenham does something entirely different than tagging all cells that have nonempty intersection with your curves. Especially see the examples in your link.

3

Bresenham’s is for drawing a line by coloring pixels, isn’t it? Perfect for tagging intersecting cells, I’d say.

4

No. Tagging intersecting cells is a well-defined task – it begins and ends in the platonic realm of mathematics.

Drawing a line by coloring pixels is an open-ended issue that targets human perception of visual fidelity. It is not judged as “correct” or not, it is judged by somebody doing a study on how some mix of experts and laypeople perceive the colored pixels as a line, on what kind of display technology.

Consider this example on the helpfully linked wikipedia page. You see that not all pixels that intersect the line are colored, because that would look like shit.

1 Like
5

Good point.

6

The algorithm is already implemented in Meshes.jl for various 2D geometries. You can write

indices(grid, geom)

to retrieve the indices efficiently where the grid is CartesianGrid or a RegularGrid with “orthogonal” CRS.

Check the source here for inspiration:

Contributions are welcome in the 3D case.

4 Likes
7

I want the “true” interesections!

1 Like
8

Can you give a little bit of context for why you want what you want? Such that most questions answer themselves.

Next clarifying question: Suppose you have a single curve, which is a large figure-eight in 2d. This curve has a self-intersection. Are you sure that you only want to count the grid square with the self-intersection as multiplicity one? As opposed to counting it with multiplicity two?

Next question: What about curves that don’t intersect the interior of a grid cell, but that do have nonempty intersection with the boundary. I can imagine 3 reasonable answers: (1) tag the cell!, (2) don’t tag the cell!, and (3) the curves are assumed in general position, so it doesn’t matter (so throwing an exception is fine).

As an example, consider the 2d grid (1:1:100) x (1:1:100) and the line segment (2.5,2.5)-> (50,50). Is the grid cell 50:51 x 50:51 tagged? Is the grid cell 30:31 x 29:30 tagged?

Is your grid uniform, i.e. each axis is a linearly spaced grid? Are the spacings of each axis the same? Is the number of axis ticks on each axis the same? Is the total size of each axis the same (i.e. a cube) or do you rather have a rectangular grid?

In the following I am assuming the following partially clarified problem statement: The curves are considered as closed sets, i.e. including endpoints; grid cells are considered as closed sets, i.e. including boundaries, and are therefore grid cells are not mutually disjoint; a curve consisting of a single point can tag 8 cells if it sits on a grid point. Then you want to count, for each grid cell, the number of curves that have nonempty intersection with that grid cell.

Now we have a clarified real analysis problem statement. But we’re in the puny finitary realm of floating point numbers. So you need to tell us how you want to handle the approximations involved.

Do you really need an exact count? Your desired function is only upper semi-continuous. This suggests that you really only need to compute an upper bound (which you want to be reasonably tight).

To give a handwavy example: Suppose you have a real configuration where the curve does not intersect a grid-cell in Float64, but does intersect the grid-cell when viewed in Float32. Upper semi-continuity suggests that you should be OK with tagging the cell, even if you decide to work in the Float64 world. On the other hand, suppose your configuration is such that the cell doesn’t intersect the curve in Float64, but it does intersect in Float128. Upper semi-continuity suggests that you want the cell tagged, even if you work in Float64 exclusively.

To go further on that, what does it mean for a curve to intersect a grid cell “in Float64” / “in Float32” / etc? There are two reasonablish definitions:

  1. The whole Float64 discretization scheme applies to all grid and curve parameters. Then you view it as a real analysis problem.
  2. For the two points p1, p2 you ask whether there exist a floating point number 0 <= t1 <= 1 and 0 <= t2 <= 1 such that t1+t2 == 1 and p1*t1 + t2*p2 lies in the grid cell, when everything is viewed as floating point. (note that it would be something different and very weird asymmetric to ask about 0 <= t <= 1 such that p1*t + (1-t)*p2 lies in the grid cell! Because values 0 < t2 < eps cannot be expressed as t2 = (1-t1) in float).

All this subtleties suggest that you really care about a reasonably tight upper bound, not about exact counts.

4 Likes
9

The curves are numerical solutions of an SDE sampled on a GPU with Float32 precision.

I’d say two indeed.

Next question: What about curves that don’t intersect the interior of a grid cell, but that do have nonempty intersection with the boundary. I can imagine 3 reasonable answers: (1) tag the cell!, (2) don’t tag the cell!, and (3) the curves are assumed in general position, so it doesn’t matter (so throwing an exception is fine).

I’d say(1) but that situation should occur with zero probability.

Is your grid uniform

yes. It has a rectangular shape.

1 Like