Precise subpixel line drawing algorithm (rasterization algorithm)

Question

I need an algorithm which can be (a bit) slower than the Bresenham line drawing algorithm but has to be a lot more exact. With 'exact' I mean: every touched pixel should be printed. No more, but also no less! Which means using a more thick line or similar is not an option as too many pixels will be involved. Also I don't need a graphic framework or similar like it was asked before, I need the algorithm! The application is not really in 'graphics' it is in the geography area where pixels are 'tiles'.

The main problem for me is that I need subpixel precision which means that a line could start at 0.75/0.33 and not just at 0/0 like it is the case for integer values. I tried to create a working solution for the last several hours but cannot make it working - there are too many edge cases.

First I thought an anti-aliased version like the algorithm from Wu should make it but it prints too many pixels (especially for start and end points) and in certain cases it still misses some pixels e.g. for very short lines.

Then I tried to make Bresenham working where I replaced the second 'if' with 'else if' as pointed out here, and it is closer but still not there. Then I tried to move the Bresenham from integer- to float-precision which resulted in an endless loop (as the x,y values jumped over the finish condition if (y1 == y2 && x1 == x2)).

I could use the naive line drawing solution but which delta should I use? E.g. if I use 0.1 I will still miss some pixels and using smaller values it will probably take too long (and still miss pixels).

A working solution in C/Java/... would be appreciated. At least it should work for octant 1 but a full blown solution would be even nicer.

Update: I came up with the following idea: using the naive line rasterization and you can calculate 4 pixel-candidates for every point. Then check for those 4 pixels if the line really crosses them. But I'm not sure if line/box intersection can be fast enough.

Answer 1

If you need just constant color (not interpolated by used area of pixel) then use DDA:

void line_DDA_subpixel(int x0,int y0,int x1,int y1,int col) // DDA subpixel -> thick
    {
    int kx,ky,c,i,xx,yy,dx,dy;
    x1-=x0; kx=0; if (x1>0) kx=+1; if (x1<0) { kx=-1; x1=-x1; } x1++;
    y1-=y0; ky=0; if (y1>0) ky=+1; if (y1<0) { ky=-1; y1=-y1; } y1++;
    if (x1>=y1)
     for (c=x1,i=0;i<x1;i++,x0+=kx)
        {
        pnt(x0,y0,col); // this is normal pixel the two below are subpixels 
        c-=y1; if (c<=0) { if (i!=x1-1) pnt(x0+kx,y0,col); c+=x1; y0+=ky; if (i!=x1-1) pnt(x0,y0,col); }
        }
    else
     for (c=y1,i=0;i<y1;i++,y0+=ky)
        {
        pnt(x0,y0,col); // this is normal pixel the two below are subpixels 
        c-=x1; if (c<=0) { if (i!=y1-1) pnt(x0,y0+ky,col); c+=y1; x0+=kx; if (i!=y1-1) pnt(x0,y0,col); }
        }
    }

where:

void pnt(int x,int y,int col);

is routine that rasterize pixel (x,y) with color col The source is in C++

I think it is strait forward but anyway

DDA use parametric line equation y=k*x+q or x=ky+q dependent on the difference (if is bigger x or y difference so there are no holes). The k is dy/dx or dx/dy and the whole division is reduced to substraction+addition inside loop (last line of each loop). This can be easily modified to any number of dimensions (I usually use 7D or more with this). On modern machines is the speed sometimes better then Bresenham (depends on the Platform and usage).

This is how it looks like compared to simple DDA

[edit2] double coordinates // originally [edit1]

OK here is new code:

void line_DDA_subpixel1(double x0,double y0,double x1,double y1,int col)    // DDA subpixel -> thick
    {
    int i,n,x,y,xx,yy;
    // prepare data n-pixels,x1,y1 is line dx,dy step per pixel
    x1-=x0; i=ceil(fabs(x1));
    y1-=y0; n=ceil(fabs(y1));
    if (n<i) n=i; if (!n) n=1;
    x1/=double(n);
    y1/=double(n); n++;
    // rasterize DDA line
    for (xx=x0,yy=y0,i=0;i<=n;i++,x0+=x1,y0+=y1)
        {
        // direct pixel
        pnt(x,y,col);
        // subpixels on change in both axises
        x=x0; y=y0;
        if ((i<n)&&(x!=xx)&&(y!=yy)) { pnt(xx,y,col); pnt(x,yy,col); }
        xx=x; yy=y;
        }
    }

And this is how it looks like:

Angle should be in double precision now but pnt(x,y,col) is still on integers !!!

[edit3] pixel grid crossing

void DDAf_line_subpixel(float x0,float y0,float x1,float y1,int col)    // DDA subpixel -> thick
    {
    int i,n; float a,a0,a1,aa,b,d;
    // end-points
    pnt(x0,y0,col);
    pnt(x1,y1,col);
    // x-axis pixel cross
    a0=1; a1=0; n=0;
    if (x0<x1) { a0=ceil(x0); a1=floor(x1); d=(y1-y0)/(x1-x0); a=a0; b=y0+(a0-x0)*d; n=fabs(a1-a0); } else
    if (x0>x1) { a0=ceil(x1); a1=floor(x0); d=(y1-y0)/(x1-x0); a=a0; b=y1+(a0-x1)*d; n=fabs(a1-a0); }
    if (a0<=a1) for (aa=a,i=0;i<=n;i++,aa=a,a++,b+=d) { pnt(aa,b,col); pnt( a,b,col); }
    // y-axis pixel cross
    a0=1; a1=0; n=0;
    if (y0<y1) { a0=ceil(y0); a1=floor(y1); d=(x1-x0)/(y1-y0); a=a0; b=x0+(a0-y0)*d; n=fabs(a1-a0); } else
    if (y0>y1) { a0=ceil(y1); a1=floor(y0); d=(x1-x0)/(y1-y0); a=a0; b=x1+(a0-y1)*d; n=fabs(a1-a0); }
    if (a0<=a1) for (aa=a,i=0;i<=n;i++,aa=a,a++,b+=d) { pnt(b,aa,col); pnt(b, a,col); }
    }

Finally had some time for this so I tweaked DDA a little but id lead to many ifs so I change rasterization quite a bit. Now all pixel grid crossing (intersections) are computed and then for each the right sub-pixel is added. This is how it looks like (no wrong sub-pixels):

For each x or y grid lines is the first cross point computed (a,b) and step is in one axis 1 pixel and in second the rest according to dy/dx or dx/dy. After this the for loop fill the sub-pixels ...

Answer 2

If your line is thin and pixels are rectangular (square):

then consider using of voxel grid traversal algorithms, for example, see article "Fast Voxel Traversal Algorithm..." by Woo and Amanatides.

Practical implementation (in grid traversal section)

Answer to comment:
Proper initialization for X-coordinate variables (the same for Y)

  DX = X2 - X1
  tDeltaX = GridCellWidth / DX
  tMaxX = tDeltaX * (1.0 - Frac(X1 / GridCellWidth)) 
  //Frac if fractional part of float, for example, Frac(1.3) = 0.3

example in my answer here