I want to use skimage.draw.line to draw an aliased line through two pixels. The issue is when I overlay the original line, there is often a discrepancy due to rounding.
Example:
Attempting line through points (22,30) and (25,39) in a 50x50 image. The full line is the green channel and the overlay line is the red channel. Yellow is overlap. The result I want would yield no visible red, only yellow and green. 
I tried generating a bigger line and then cropping, but still getting issues:
y0x0, ymxm=(np.array([0,0]), np.array([50,50]))
y1x1, y2x2=(np.array([22,30]), np.array([25,39]))
#rawzero3 is essentially trunc() but rounding up
#https://stackoverflow.com/q/46572699
yx_fin1=(y2x2-y1x1)*np.amin(\
rawzero3((ymxm-(y2x2-y1x1))/(y2x2-y1x1)).astype(int))+y1x1
yx_fin2=(y0x0-1)-((y2x2-y1x1)*np.amin(\
rawzero3((y0x0-(y2x2-y1x1))/(y2x2-y1x1)).astype(int))+y1x1)
zimg[tuple(map(tuple,np.array(draw.line(*yx_fin1,*yx_fin2)).T\
[np.all((y0x0<=np.array(draw.line(*yx_fin1,*yx_fin2)).T)\
&(np.array(draw.line(*yx_fin1,*yx_fin2)).T<ymxm), axis=1)].T))]=255
How can my desired result be achieved?
Edit
I updated my code to deal with issues outside of octant 0, similar to Bresenham's Line Algorithm:
def rawzero3(a):
_a = np.abs(a)
np.ceil(_a, _a) # inplace
np.copysign(_a, a, out = _a)
return _a
def fullline_pos(r0,c0,r1,c1,s):
y0x0, ymxm=(np.array([0,0]), np.array(s))
y1x1, y2x2=(np.array([r0,c0]), np.array([r1,c1]))
mdyx=(y2x2-y1x1)//np.gcd.reduce(y2x2-y1x1)#this will always be a whole number, the // makes it an int
yx_fin1=mdyx*np.amin(rawzero3((ymxm-y1x1)/mdyx).astype(int))+y1x1 #should be slightly more efficient
yx_fin2=(mdyx*np.amin(rawzero3((y0x0-y1x1)/mdyx).astype(int))+y1x1) #should be slightly more efficient
if v: print(f'y0x0 {y0x0}, ymxm {ymxm}, y1x1 {y1x1}, y2x2 {y2x2}, (y2x2-y1x1) {(y2x2-y1x1)}, yx_fin1 {yx_fin1}, yx_fin2 {yx_fin2}')
return tuple(map(tuple,np.array(draw.line(*yx_fin1,*yx_fin2)).T[np.all((y0x0<=np.array(draw.line(*yx_fin1,*yx_fin2)).T)&(np.array(draw.line(*yx_fin1,*yx_fin2)).T<ymxm), axis=1)].T))
def fullline_neg(r0,c0,r1,c1,s):
sr,sc=np.array(s)-1
r0_,c0_,r1_,c1_=sr-r0,c0,sr-r1,c1
coords=np.array(fullline_pos(r0_,c0_,r1_,c1_,s))
coords[0]=sr-coords[0]
return tuple(map(tuple,coords))
def fullline(r0,c0,r1,c1,s):#returns a 2xN array of valid yx coordinates
#r0,c0,r1,c1 are int, s is a 2d shape
assert (r0,c0)!=(r1,c1)
if r0==r1: return draw.line(r0,0,r0,s[1]-1)
if c0==c1: return draw.line(0,c0,s[0]-1,c1)
if (r0<r1)==(c0<c1):
return fullline_pos(r0,c0,r1,c1,s,v=v)
else:
return fullline_neg(r0,c0,r1,c1,s,v=v)
For visualization purposes, I ran fullline(y,x,25,25,(50,50)) for all x and y values 0 to 49. The following image shows the relative amount of errors (red pixels) for each y,x location: 
