Straight line drawing by matrix-Inventor's original research of Python image processing-
A story about drawing an antialiasing straight line ** in your own way ** by just matrix operations without relying on an image processing library. Also possible with Pythonista
** Click here for basics ** ** Click here for drawing **
Independent research !?
[Bresenham's algorithm](https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%AC%E3%82%BC%E3%83%B3%E3%83%8F%E3 Read% 83% A0% E3% 81% AE% E3% 82% A2% E3% 83% AB% E3% 82% B4% E3% 83% AA% E3% 82% BA% E3% 83% A0) If so, there is something called Xiaolin Wu's algorithm, and there is also Gupta-Sproull's algorithm. There was .ne.jp/ousttrue/20090224/1235499257). However, I thought I would do my best if I used matrix calculation. I called it my own research because I couldn't find it, but I don't think it is unique. (In short, reinventing the wheel)
concept
From here Borrowed
The figure above shows that a straight line with a thickness of 1 spans various pixels. I want to give each pixel a density according to this straddling area.
Situation settings
Excuse me for handwriting
The figure above shows how a straight line is drawn in the pixel. That is, the squares indicate pixels, and the three diagonally drawn straight lines indicate one thick straight line T (the middle straight line M is the center line). Consider a pixel whose center is the point A $ (n_x, n_y) $ and a straight line T $ ax + by + c = 0 $ whose thickness is $ w / 2 $. However, $ a \ geq-b \ geq0, a> 0 $. Also, the slope of the straight line is $ \ theta =-\ arctan \ frac {b} {a} $.
This pixel has four corner points on the upper left $ A_ {ul} (n_x-0.5, n_y-0.5) $, Upper right $ A_ {ur} (n_x-0.5, n_y + 0.5) $, Bottom left $ A_ {dl} (n_x + 0.5, n_y-0.5) $, Lower right $ A_ {dr} (n_x + 0.5, n_y + 0.5) $.
On the other hand, for the thick straight line T, three lines, U indicating the upper limit, D indicating the lower limit, and M running in the center, are shown in the figure. Both the upper and lower lines are separated by $ w / 2 $ from the center line, and the total thickness is $ w $.
The straight lines that coincide with the left and right edges of the pixel at the center of A are L and R, respectively, and the subscripts of the letters L and R change depending on which straight line they intersect.
Problem setting
In the above situation, find the area S of the intersection (light blue) of pixel A and the parallelogram B $ L_UR_UR_DL_D $ (point B is a point on the straight line M in the positive direction of A's x-axis).
Definition of S_X
When $ t = (x coordinate of B)-(x coordinate of A_ {ul} A_ {ur}) $, below the line segment $ A_ {ul} A_ {ur} $ (x-axis positive direction) Let the area of a certain figure $ X $ be $ S_X (t) $ (image of translating the figure from top to bottom)
S=S_B(t)-S_B(t-1)
Will be. The second term is the area of the part below the line segment $ A_ {dl} A_ {dr} $. Therefore, $ S_B (t) $ should be calculated.
Decomposition of parallelogram B
As shown in the figure, the parallelogram B $ L_UR_UR_DL_D $ can be divided into three parts: a red right-angled triangle $ \ triangle Red $, a blue right-angled triangle $ \ triangle Blue $, and a white rectangle $ \ Box White $. .. (In the figure below, the white rectangle has a negative area)
Then
S_B(t) = S_{\triangle Red}(t) +S_{\triangle Blue}(t) +S_{\Box White}(t)
Will be.
Therefore, consider $ S_X (t) $ of each figure. In addition, it should be noted
\tau_\triangle=\frac{1}{2}\left(w/\cos\theta +tan\theta\right)\\
\tau_\Box=\frac{1}{2}\left(w/\cos\theta -tan\theta\right)
And. ($ \ Tau_ \ triangle $ corresponds to the $ x $ coordinate difference between B and $ R_D $ and $ L_U $, and $ \ tau_ \ Box $ corresponds to the $ x $ coordinate difference between B and $ R_U $ and $ L_D $. Yes)
First, the total area of these three parts is
\triangle Red=\triangle Blue = 0.5\tan\theta = \frac{\tau_\triangle-\tau_\Box}{2}\\
\Box White = 2\tau_\Box
\\
Next, consider each $ S_X (t) $. $ \ Triangle Red $ exceeds $ A_ {ul} A_ {ur} $ when $ t =-\ tau_ \ triangle $, and the entire figure fits when $ t =-\ tau_ \ Box $. Also, the increase in between is a downwardly convex quadratic function.
S_{\triangle Red}(t) =
\begin{cases}
0 &
(t<-\tau_\triangle)\\
\frac{1}{2}\frac{1}{\tau_\triangle-\tau_\Box}\left(t+\tau_\triangle\right)^2 &
(-\tau_\triangle\leq t\leq-\tau_\Box)\\
\frac{\tau_\triangle-\tau_\Box}{2} &
(-\tau_\Box<t)
\end{cases}
$ \ Box White $ exceeds $ A_ {ul} A_ {ur} $ when $ t = \ min (-\ tau_ \ Box, \ tau_ \ Box) $, and $ t = \ max (-\ tau_ \ Box) , \ tau_ \ Box) When $, the whole figure fits. Also, since the increase in the meantime is linear,
S_{\Box White}(t) =
\begin{cases}
0 &
(t<-\tau_\Box)\\
t+ \tau_\Box &
(-\tau_\Box\leq t)\\
\end{cases}
+
\begin{cases}
0 &
(t<\tau_\Box)\\
-t+ \tau_\Box &
(\tau_\Box\leq t))\\
\end{cases}
$ \ Triangle Blue $ exceeds $ A_ {ul} A_ {ur} $ when $ t = \ tau_ \ Box $, and when $ t = \ tau_ \ triangle) $, the entire figure fits. Also, the increase in between is a quadratic function that is convex upwards. 
Therefore, this is equivalent to the function of moving $ S_ {\ triangle Red} (t) $ point-symmetrically to the origin target and then adding $ \ frac {\ tau_ \ triangle- \ tau_ \ Box} {2} $. ..
S_{\triangle Blue}(t) = \frac{\tau_\triangle-\tau_\Box}{2}-S_{\triangle Red}(-t)
code
The initial conditions are as follows.
import numpy as np
import matplotlib.pyplot as plt
x, y = np.mgrid[:10,:10]
a, b, c = 2, -1, -1
w = 1
calculation of t
This is easy. Just add 0.5 to the distance between the pixel and the straight line.
import numpy as np
import matplotlib.pyplot as plt
x, y = np.mgrid[:10,:10]
a, b, c = 2, -1, -1
w = 1
t = - b/a * y - c/a - x + 0.5
Calculation of τ
Calculate without going through $ \ theta $.
tau_triangle = (w*(a**2 + b**2)**0.5/a - b/a)/2
tau_square = (w*(a**2 + b**2)**0.5/a + b/a)/2
Calculation of S_X
#S_Because it is also used in blue, S_Declare red as a function
S_red_func = lambda t :\
np.where(t<-tau_triangle,0,
np.where(t<= -tau_square,(t+tau_triangle)**2/(2*(tau_triangle - tau_square)),
(tau_triangle - tau_square)/2))
S_red = S_red_func(t)
S_blue = (tau_triangle - tau_square)/2 - S_red_func(-t)
S_white = np.where(t<-tau_square,0, t+tau_square) + \
np.where(t<tau_square,0, -t+tau_square)
#S_Create B
S_B = S_red + S_blue + S_white
Calculate the difference in S_X
Use np.diff. To calculate $ S_B (t)-S_B (t-1) $, -np.diff (S_B) to arrange before and after subtraction
use.
It also performs conversion to make a black straight line on a white background.
img_how_long_did_it_take_to_get_this = -np.diff(S_B,axis = 0)
#Maximum value is 1 pixel area=1
img_black = 1-img_how_long_did_it_take_to_get_this
plt.imshow(img_black, cmap = 'gray')
Note that the vertical pixels are decremented by 1 due to the diff.
Summarize as a function
Create a function after removing the condition of $ a \ geq-b \ geq0, a> 0 $. (In addition, exception handling is not performed when a and b are 0)
def my_antialias(height, width, a, b, c, w=1):
'''Straight line ax with thickness w on a plane with height height and width width+ by + c =Draw 0 with antialiasing.'''
#Add the amount that decreases when you take the diff first
x, y = np.mgrid[:height + 1,:width]
#Calculate t(a times)
t = 0.5*a - (a*x + b*y + c)
#Calculate τ(a times)
tau_triangle = (w*(a**2 + b**2)**0.5 +abs(b))/2
tau_square = (w*(a**2 + b**2)**0.5 - abs(b))/2
#Calculation of S
S_red_func = lambda t :\
np.where(t<-tau_triangle,0,
np.where(t<= -tau_square,(t+tau_triangle)**2/(2*(tau_triangle - tau_square)),
(tau_triangle - tau_square)/2))
S_red = S_red_func(t)
S_blue = (tau_triangle - tau_square)/2 - S_red_func(-t)
S_white = np.where(t<-tau_square,0, t+tau_square) + \
np.where(t<tau_square,0, -t+tau_square)
#t,Since τ is multiplied by a, divide by a
S_B = (S_red + S_blue + S_white)*a
return -np.diff(S_B,axis = 0)
plt.imshow(my_antialias(20,20,-3,-7,60,1),cmap = 'gray',vmin = 0,vmax=1)
plt.show()
By the way, the following is [Bresenham's algorithm imitation](http://qiita.com/secang0/items/b45272a6f1212510ef99#%E3%83%96%E3%83%AC%E3%82%BC%E3%83%B3% E3% 83% 8F% E3% 83% A0% E3% 81% AE% E3% 82% A2% E3% 83% AB% E3% 82% B4% E3% 83% AA% E3% 82% BA% E3% 83% A0% E6% 93% AC% E3% 81% 8D% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E7% 9B% B4% E7% B7% 9A% E6% 8F% It is a straight line generated by 8F% E7% 94% BB). 
In my execution environment, using Bresenham's algorithm was 3-4 times faster.
Recommended Posts