Mathematica, 166 137 bytes
l:={i,j};s=Sign;f[p_,q_,h_,w_]:=Grid@Table[(1-Max[s[p-l]s[q-l],0])Boole[Abs@Mean[s@Det@{p-l+#,p-q}&/@Tuples[.5{1,-1},2]]<.6],{i,h},{j,w}]
More readable version:
l := {i, j}; s = Sign;
f[p_, q_, h_, w_] :=
Grid@Table[(1 - Max[s[p - l] s[q - l], 0]) Boole[
Abs@Mean[
s@Det@{p - l + #, p - q} & /@
Tuples[.5 {1, -1}, 2]] < .6], {i, h}, {j, w}]
This defines a function called f. I interpreted the input and output specifications fairly liberally. The function f takes input in the format f[{x0, y0}, {x1, y1}, height, width], and the grid is 1-indexed, starting in the top left. Outputs look like
with the line displayed as 1s and the background as 0s (shown here for f[{2, 6}, {4, 2}, 5, 7]). The task of turning a Mathematica matrix of 1s and 0s into a string of #s and .s has been golfed in many other challenges before, so I could just use a standard method, but I don't think that adds anything interesting.
Explanation:
The general idea is that if the line passes through some pixel, then at least one of the four corners of the pixel is above the line, and at least one is below. We check if a corner is above or below the line by examining the angle between the vectors ({x0,y0} to corner) and ({x0,y0} to {x1,y1}): if this angle is positive, the corner is above, and if the angle is negative, the corner is below.
If we have two vectors {a1,b1} and {a2,b2}, we can check if the angle between them is positive or negative by finding the sign of the determinant of the matrix {{a1,b1},{a2,b2}}. (My old method of doing this used arithmetic of complex numbers, which was way too…well, complex.)
The way this works in the code is as follows:
{p-l+#,p-q}&/@Tuples[.5{1,-1},2] gets the four vectors from {x0,y0} and the four corners of the pixel (with l:={i,j}, the coordinates of the pixel, defined earlier), and also the vector between {x0,y0} and {x1,y1}.s@Det@... finds the signs of the angles between the line and the four corners (using s=Sign). These will equal -1, 0 or 1.Abs@Mean[...]<.6 checks that some of the angles are positive and some negative. The 4-tuples of signs that have this property all have means between -0.5 and 0.5 (inclusive), so we compare to 0.6 to save a byte by using < instead of <=.
There's still a problem: this code assumes that the line extends forever in both directions. We therefore need to crop the line by multiplying by 1-Max[s[p-l]s[q-l],0] (found by trial and error), which is 1 inside the rectangle defined by the endpoints of the line, and 0 outside it.
The rest of the code makes a grid of these pixels.
(As a bonus, here's an earlier attempt with a completely different method, for 181 bytes:)
Quiet@Grid@Table[(1-Max[Sign[{i,j}-#3]Sign[{i,j}-#4],0])Boole[#3==#4=={i,j}||2Abs@Tr[Cross@@Thread@{{i,j},#3,#4}]/Norm[d=#3-#4]<2^.5Cos@Abs[Pi/4-Mod[ArcTan@@d,Pi/2]]],{i,#},{j,#2}]&
