実質8行の再帰でフラクタル図形（シェルピンスキーのカーペット）を描く

　こういうやつを描きます。

　任意の座標・スケールで出来ます。

コード

import numpy as np
from PIL import Image
import math

def generate_sierpinski_carpet(
   canvas_size, u_start, u_end, v_start, v_end, output_file
):
   eps = 1e-12
   assert abs(u_start + v_end - u_end - v_start) < eps, "must be a square"
   r = u_end - u_start

   canvas = np.ones((canvas_size, canvas_size), dtype=np.uint8) * 255

   def dfs(u, v, L):
       if L * canvas_size < r:
           return 255
       p = math.floor(u)
       q = math.floor(v)
       if p == 1 and q == 1:
           return 0
       else:
           return dfs((u - p) * 3, (v - q) * 3, L / 3)

   for i in range(canvas_size):
       for j in range(canvas_size):
           u = u_start + r * i / canvas_size
           v = v_start + r * j / canvas_size
           canvas[i][j] = dfs(u, v, 3)

   image = Image.fromarray(canvas)
   image.save(output_file)
   
generate_sierpinski_carpet(729,0,1,0,1,'output.png')

ギミック

　重要なのはこの $8$ 行の部分

    def dfs(u, v, L):
        if L * canvas_size < r:
            return 255
        p = math.floor(u)
        q = math.floor(v)
        if p == 1 and q == 1:
            return 0
        else:
            return dfs((u - p) * 3, (v - q) * 3, L / 3)

　長さ $3 \times 3$ の、上のような領域を考えます。
　この中から、実際に描画する領域を赤枠とします。便宜的にこれを UV マップと呼びます。

　まず、この UV マップ上の点が、領域の中のどの部分にあるかを考えます。

　真ん中の領域、つまり、左から $2$ 番目、上から $2$ 番目の領域にある時、ここは空白であることが確定します。ここは黒に塗ってしまって問題ありません。

　問題ははみ出た領域（図で「？」をつけた部分）です。これは現段階では良くわかりません。なので、次の探索に進みます。

　この領域は、更に $3 \times 3$ 分割する必要がありますが、逆に、そこに所属する点が $3$ 倍の座標になったと考えると、全く同じ話ができます。

　左上がわかりやすいと思うので図にすると、こんな感じになります。

　なお、このままだと座標によっては無限にループしてしまうので、どこかで探索を打ち切ります。現在見ている領域のサイズ（L）をキャンパスサイズの数だけ並べても元の UV マップを越えない場合、これ以上細かく見ても意味がないと判断して打ち切ります。