while loop may never exit (solution would be to check tMaxX, tMaxY,tMaxZ to surpass T1)

[SkybuckFlying]

while (current_X_index != end_X_index || current_Y_index != end_Y_index || current_Z_index != end_Z_index)

if for whatever reason these indexes never happen the loop will hang forever. I am not yet sure why these indexes never happen, it might be because there are off by 1 or perhaps floating errors make them go off by 1.

It might be saver to check for tMaxX, tMaxY, tMaxZ to be smaller than the (real space) line length.

And then the code should go from t0 to t1.

while (tMaxX < RayLength) or (tMaxY < RayLength) or (tMaxZ < RayLength)

This could and would most probably make the code a lot saver.

[cgyurgyik]

Hi, thanks for pointing this out. Feel free to submit a pull request! I haven't touched this code in quite some time, but I'm sure others would appreciate it.

[SkybuckFlying]

Hi, my original comment/issue was slightly wrong, I updated it to replace T1 with RayLength in this comment.

Also I cannot perform a pull request for you, because I re-wrote the code heavily to be more readable and to compile with older c++ compilers.

Basically everything with just double variables, instead of special classes and routines, everything embedded into a single function and all () removed which make the code hard to read, and some further re-formatting of the code.

I see many new c++ features/keywords and don't really know what they do. I do have visual studio 2019 and I might be able to test it in a seperate test project. However this repository does not come with test code. I would be nice if there was a test project, that would save me some time.

Especially with extereme test cases. Perhaps you could add some test code ?

[cgyurgyik]

From the README.md:

Note, this has not been tested, and is not guaranteed to be bug-free.

I don't intend on up-keeping this myself, but will happily accept pull requests. This compiled at one point, but I have no idea if it compiles today :-). I was using this for prototyping rather than a golden reference model.

[SkybuckFlying]

My strategy for now is to try out different code and find one that works. This code of yours uses some weird/new compiler things, don't see why that is necessary.

For now the difficult with this algorithm is with setting things up and ending the loop. So me will try some different implementations and see which one comes out on top =D

I even tried to create a version which could later be parallized, it could theoretically use three for loop to loop over the x,y,z integers between the two points, inspired by somebody elses posting, but couldn't get the math quite right. It was close but not perfect.

I also took a look at my very old version in same game, but couldn't get the code ported/working well in a simple example, might have to spent more time debugging it, not worth it, for now, so I try other peoples implementations first to see if I can cut some time and also for better math or inspiration and also consensus on best way to calculate things, though I see plenty of problems with divisions, and possible floating points errors.

The parallel version would strangle enough be more exact, it does not use iterative processing/adding. It calculates each point individually/directly which is a bit interesting. However since it's integers and not floating points the math gets a bit strange...

Also I am starting to doubt this entire algorithm. It might be much more complex than the original paper makes it seem to be.

What if lines are completely vertical or completely horizontal. In that case the other dimension is "inactive" and perhaps this algorithm will fail and never complete, thus the actual implementation that is flawless will become much more complex and needs a lot more logic to actually work perfect.

So for now I consider this paper half-baked and perhaps even misleading, especially since setup and end phase seem to be shady, plus not all corner cases examined, but it remains interesting none the less.

Oh that reminds me, I did do another implementation a while ago in cuda, now I have two to look at one from myself and one from somebody else... My c cuda one was actually used for light casting/shading algorithm and it worked well, so that is why I am convinced there is value in this algorithm.

I knew I had an algorithm somewhere that was good, just couldn't remember where it was lol, but now I remember, wrote it in Delphi and ported to Cuda C, so it was usefull discussing this with you ! =D

I tried to do this light casting idea in 1995, with turbo pascal, 64 kb memory, ms-dos, 4 mb ram, vga, etc, but I was total noobish programmer and cpu power wasn't there... I got as far as shading 1 little ray/line lol... but in 2013... I did it ! =D

[SkybuckFlying]

This code was converted to Pascal/Delphi for testing and I can confirm it works flawlessly, it does need some slight adjustments:

The one that worked for me in 3D for both positive and negative directions (in CUDA C):

#define SIGN(x) (x > 0 ? 1 : (x < 0 ? -1 : 0))
#define FRAC0(x) (x - floorf(x))
#define FRAC1(x) (1 - x + floorf(x))

float tMaxX, tMaxY, tMaxZ, tDeltaX, tDeltaY, tDeltaZ;
int3 voxel;

float x1, y1, z1; // start point
float x2, y2, z2; // end point

int dx = SIGN(x2 - x1);
if (dx != 0) tDeltaX = fmin(dx / (x2 - x1), 10000000.0f); else tDeltaX = 10000000.0f;
if (dx > 0) tMaxX = tDeltaX * FRAC1(x1); else tMaxX = tDeltaX * FRAC0(x1);
voxel.x = (int) x1;

int dy = SIGN(y2 - y1);
if (dy != 0) tDeltaY = fmin(dy / (y2 - y1), 10000000.0f); else tDeltaY = 10000000.0f;
if (dy > 0) tMaxY = tDeltaY * FRAC1(y1); else tMaxY = tDeltaY * FRAC0(y1);
voxel.y = (int) y1;

int dz = SIGN(z2 - z1);
if (dz != 0) tDeltaZ = fmin(dz / (z2 - z1), 10000000.0f); else tDeltaZ = 10000000.0f;
if (dz > 0) tMaxZ = tDeltaZ * FRAC1(z1); else tMaxZ = tDeltaZ * FRAC0(z1);
voxel.z = (int) z1;

while (true) {
// process voxel
if (tMaxX < tMaxY) {
if (tMaxX < tMaxZ) {
voxel.x += dx;
tMaxX += tDeltaX;
} else {
voxel.z += dz;
tMaxZ += tDeltaZ;
}
} else {
if (tMaxY < tMaxZ) {
voxel.y += dy;
tMaxY += tDeltaY;
} else {
voxel.z += dz;
tMaxZ += tDeltaZ;
}
}
if (tMaxX > 1 && tMaxY > 1 && tMaxZ > 1)
{
// process voxel here
break;
}
}

Code above slightly adjusted to make it flawless, immediately after loop statement process first voxel, and before break process end voxel.

The way to call this code is to divide coordinates by the cell width, height, dpeth of the grid without any truncate or rounding... just leave the floating points as they are and this algorithm will then automatically traverse the grid in integer style.

The division should be floating point style, example of setup of x,y,z:

x1 := point1.x / cell_width;
y1 := point1.y / cell_height;
z1 := point1.z / cell_depth;

x2 := point2.x / cell_width;
y2 := point2.y /cell_height;
z2 := point2.z / cell_depth;

^ make sure the above division is floating point and then algorithm works flawlessly ! =D

[SkybuckFlying]

So what is noticeable for the above posted code, it is very clean. It does not incorporate any grid. It does not need to, grids always have the same indexing internally. Strange but true.

All that is necessary is to convert the x,y,z input coordinates in such a way that they lie on a perfect 1-sized-bucket-grid which is symmetrical... the divisions by cell width and such take care of that.

It would be nice to have some seperate clipping code though ;)

[SkybuckFlying]

Apperently this algorithm works on the basis of ratios of 0..1 for the cell width, height, depth... not sure but this might have been hard to understand in the original document. I just checked the document and it's indeed a bit misleading because the document does not use symmetrical boxes for the visualizations of the lines and it's probably not mentioned anywhere.

These ratios can ofcourse exceed 1 but that means it will enter next cell...

Example suppose cell width is 160 and x coordinate is 160

160/160 = 1 means, it just enter the cell 1 instead of 0.

Other example:

159/160 = 0.999 or something, still in cell 0

;)

[SkybuckFlying]

So having to scale any line or coordinates into a grid which has cells of 1 by 1 by 1... might actually be a down side of this algorithm... at least the traversing is fast. So basically this algorithm in it's current form would require 6 divisions per ray... unless the grid is already cells of 1x1x1.... for gaming it's more interesting to have bigger cells...

[SkybuckFlying]

The posted C algorithm does have one major pain in the ass issue, it seems to mostly overshoot the grid by 1 voxel in certain cases, possibly caused by floating point instability. When replacing algorithm by direct computation of each integer between the coordinates it also overshoots, from strange, somehow the starting point must be adjusted depending on which direction one is travelling.