経緯
「２が現れる素数」という記事がありました。
http://integers.hatenablog.com/entry/2017/11/29/082604
こちらでは
700000000000000007
000000222222000000
000022222222220000
000222000002222000
000000000022220000
000000000222200000
000000002222000000
000000222200000000
000022220000000000
000222222222222000
000222222222222000
700000000000000003
という素数を紹介していました。
同様に、マリオが現れる素数というのも紹介されています。
https://akiyah.hatenablog.com/entry/2017/12/05/144440
一方で、2が現れる素数が奇跡だという人に物申す
https://qiita.com/elzup/items/1d882f3af040506aec8b
という記事にもあるように、桁数の多い数字でも素数の確率は意外と大きく、コメント欄では $\pi(x)/x \sim 1/\ln(x)$ だと見積もられています。
さて、マリオ素数のように画像を素数化する場合だと多いときは１０色を１０個の数字に割り当てるのですが、その場合の数は $10! = 3628800$ 通りです。左上は０を使えないとか、右下は5以外の奇数だとして割り当て方を制限しても相当の場合の数があります。この場合の数だと相当大きい画像を素数化しない限り、素数はほぼ必ず見つかるのではと考え、やってみました。
僕は数学はあまり詳しくないので、愚直に画像を１０色以下に減色、数字を全パターン当てはめ、素数になるまで繰り返すこととしました。
ということで、この記事の個人的なポイントは、
1. 減色のアルゴリズム
2. 素数判定のコマンドssl prime
です。
処理の方針
- 画像のリサイズ
- 画像の減色
- 数字のあてはめ、素数判定
- 素数が見つかれば終了、結果表示
画像のリサイズ
特に難しいことは無いです。
ただ普通は縮小処理だと思いますが、ごくまれに拡大することがあるかもしれません。このとき、線形補間ではなく最近傍を用いた方が補間時に色が増えず、素数画像の見た目が良くなるかと思います。
ここで縦横ともにピクセル数を素数 $p$、$q$ にすると、$pq$ 桁の素数を渡されたときに$p$ × $q$ に並べてみるという発想に自然に行き着くかなと思います。$p = q$ だとさらに都合がいいです。
トリミングは考えてなかった。前もって適当なレタッチソフトで余白を削っている前提とします。
画像の減色処理
ピクセルごとに0～9の数字しかあてはめられないので、最大でも10色しか使っていない画像にする必要があります。
ここでは2通りの方法を紹介します。
- グレースケールののちLook up tableを用いる
def quantize_img_LUT(img, ncolor):
img_p = None
img_g = img
if len(img.shape) == 3:
img_g = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
max_bin = img_g.max()
min_bin = img_g.min()
n = min(ncolor, max_bin - min_bin + 1)
bins = np.linspace(min_bin, max_bin + 1, n + 1)
y = np.array([bins[i - 1] for i in np.digitize(np.arange(256), bins)]).astype(int)
img_p = np.array(cv2.LUT(img_g, y), dtype=np.uint8)
return img_p
画像をグレースケールにし、一番暗い画素値と一番明るい画素値を取ってきて ncolor 分割、あとは digitize で代表値を設定して cv2.LUT で画像を減色します。
ただし、カラーだと違う色でもグレースケール観点では同じ色になってしまうことがあり、そうなると同じ数字を当てはめてしまうので、今回は違う方法を考えます。
- k-means法による色のクラスタ
画像中の各ピクセルは $(b, g, r)$ の値を持っていますが、これを $(b, g, r)$ 空間中にピクセルが散らばっているとして、クラスタ分けするという発想です。$x$、 $y$ 座標といった空間的な情報は今回は捨てています。
参考：https://docs.opencv.org/4.0.0/d1/d5c/tutorial_py_kmeans_opencv.html
def quantize_img_kmeans(img, ncolor):
h, w, c = img.shape
Z = np.float32(img.reshape((-1, 3)))
crit = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.1)
ret, label, center = cv2.kmeans(Z, ncolor, None, crit, 10, cv2.KMEANS_RANDOM_CENTERS)
center = np.uint8(center)
result = center[label.flatten()]
img_p = result.reshape(img.shape)
return img_p
これでこの関数は ncolor 色に減色した画像を返します。
コードの解説を少し。
k-means法の実装はcv2.kmeansを使います。OpenCV、ほんとにいろいろあるな。
center = np.uint8(center) は各クラスタの重心座標を表していて、ここでは各点が $(r, g, b)$ の3次元なので、 center.shapeは [ncolor, 3] です。
laeblはサンプルの各点がどのクラスタに属するかのラベルが並んでいて、全ピクセル数をnpixelとして、成分が npixelあります。 なんか 1 の成分があるかもしれないので、 shapeはちょっと断言できません。
以上を合わせると result = center[label.flatten()] はnumpyのインデクシングにより、各点の $(r, g, b)$ 値が並んだ配列になります。 result.shapeは[npixel, 3]となります。
あとは img_p = result.reshape(img.shape) として画像としてのサイズの情報を持たせ、完成です。
たとえばWindowsのサンプルピクチャにあるペンギンの写真を6色に落としてみました。
元画像：
減色：
黄色の部分を残すのに６色必要でした。
ただ、今回の話は数字を当てはめるという関係上ピクセルごとの色のラベルがあればいいので、center から画像を作ることをせず、途中で関数から抜けます。
def quantize_img_kmeans(img, ncolor):
h, w, c = img.shape
Z = np.float32(img.reshape((-1, 3)))
crit = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.1)
ret, label, center = cv2.kmeans(Z, ncolor, None, crit, 10, cv2.KMEANS_RANDOM_CENTERS)
return label.reshape((h, w))
これでこの関数はピクセルごとのラベルを返します。同じラベルなら同じ色（数字）。
減色画像が得られたら
さて、ラベルが得られたら各ラベルが元の画像ではどんな色かを調べないといけません。
# img_ico: 入力画像をリサイズした画像
img = quantize_img_kmeans(img_ico, ncolor)
pixels = list(set(img.flatten()))
avgcolor = {}
for ipix, pix in enumerate(pixels):
avgcolor[ipix] = np.array(np.mean(img_ico[img == pix].reshape(-1, 3), axis=0), dtype=np.uint8)
なにか配列などiterableなもの hoge があったとき、set(hoge) とすれば重複したものが取り除けます。もともとは set は集合をあらわすものですが、その性質をうまく使ったかんじですね。
これで画像中にある画素値のリスト pixels が得られます。LUTでつくっていればy、k-means法で作っていれば label に当たる変数です。
avgcolorというdictに、ipix 番目の色が実際どんな色を持っているかを保存します。
img_ico[img == pix].reshape(-1, 3)は減色画像中で画素値 pix をもつ画素を抽出する処理になります。元はカラー画像だから reshape に 3 が入っています。また、そのため平均するのも画素間であって同じ画素中の $(r, g, b)$ までは平均しないように、np.mean の引数に axis=0 を指定しています。これでavgcolorはkmeans法ではcenterと同等の量が入ります。
k-means法だとcenterなど途中に出てきた変数を使えば必要なかった処理ですが、一応k-means法以外の方法で減色処理をするかもしれないので、二度手間でもこの処理を書いておきます。
素数判定
素数判定アルゴリズム
Pythonでこういうアルゴリズムを実装するとすごい重そうなので、おとなしくLinuxコマンドにあるものを使うことにしました。
RSA暗号では巨大な桁数の素数を使います。その一環でopensslには素数かどうかを判定するコマンド openssl prime がありますので、これを用います。
まず画像のピクセルに0-9の数字が割り当てられていて、それを並べた２次元行列があったとして、これを数字にします。といっても2次元のndarrayを単純に１行の文字列にするだけです。numpyにちょうどいい関数がありますので、それを使います。
def array_to_str(arr):
# arr: array of 0-9 integers
ret = np.array2string(arr.flatten(), max_line_width=arr.size+2, threshold=arr.size, separator="")[1:-1]
return ret
そしてその文字列が素数かどうかのチェック関数です。
$ openssl prime xxxx
とした時、このコマンドはxxxxが素数だとxxxx is prime、 合成数だとxxxx is not primeと標準出力に出力します。なのでsubprocessあたりで実行して標準出力をとり、 is primeという文字列があるかどうかを確かめる関数とします。
def check_prime(string):
cmdret = subprocess.run(["openssl", "prime", string], stdout=subprocess.PIPE)
return "is prime" in cmdret.stdout.decode()
subprocess.runのcapture_output引数はpython 3.7からです。
確認してないこと
factor コマンドは因数分解ですが、 openssl prime は判定だけするコマンドかなと思い、openssl を使ってみました。 factor コマンドを使うとどれくらい実行時間が変わるかは確認してません。
ただちょっと確かめた感じでは、 openssl prime だと１秒程度だった素数に対し、それから1/3程度切り出した数字にたいして factor を使ってみると数時間たっても終わりませんでした。
また、 openssl prime はミラーラビン法を使っているかもしれません。この場合、素数でないものを素数だと判定する可能性があります。対策としては、-checksオプションでより厳しくする、素数と判定されたものに対しては別の判定プログラムと組み合わせるなどが考えられます。
ちょうどいい数字割りあてのサーチ
単純なループです。
npix = len(pixels)
idx_mat = np.zeros_like(img)
for i, px in enumerate(pixels):
idx_mat[img == px] = i
lefttop = idx_mat[0, 0]
rightbottom = idx_mat[-1, -1]
last_img = np.array(idx_mat)
for newpx in itertools.permutations(range(10), npix):
if newpx[lefttop] == 0:
continue
if newpx[rightbottom] in (0, 2, 4, 5, 6, 8):
continue
for i in range(npix):
last_img[idx_mat == i] = newpx[i]
numstring = array_to_str(last_img)
if check_prime(numstring):
# 適当な出力
break
new_pxは色ラベル番号と数字の対応をもちます。例えばnew_px[i]が3だったら、色ラベルiには3を当てはめるという事になります。
idx_matが色ラベル番号を保持する行列です。あるピクセル [y,x] が元の画像でどんな色に対応しているかを知りたい時、idx_mat[y, x]の値を取り出し、これがkだったらこのピクセルは avgcolor[k] だという事になります。
last_imgは実際に置かれる数字を並べた画像です。last_img[y, x]がkだった場合、このピクセルには数字kが置かれます。
その他の処理
あとは途中経過だったり結果画像を作って表示・保存するコードを足してやってひとまずは動きます。
コマンドライン引数を一通り整備したスクリプトをgithubに上げました。
https://github.com/sage-git/prime_image_maker
実行例
次の画像を使ってみます。
これを縦横サイズを37ピクセルに縮小し、数字を当てはめて画像になる素数を作ってみたく思います。
checking (1, 3)
checking (1, 7)
checking (1, 9)
checking (2, 1)
checking (2, 3)
checking (2, 7)
checking (2, 9)
checking (3, 1)
found good permutation
----
3333333333111111113333333333331111111333333311111111111111333333333111111133333311111111111111113333333331111113333311111111111111111113333333311111333311111111111111111111333333331111133311111111111111111111113333333311113311111111111111111111111133333331111331111111111111111111111111333333111131111111111111111111111111133333331113111111111111111111111111113333333111111111111111111111111111111133333311111111111111111111111111111113333331111111111111111111111111111111333333111111111111111111111111111111133333311111111111111111111111111111113333331111111111111111111111111111111333333111111111111111111111111111111133333311111111111111111111111111111113333331113111111111111111111111111113333333111311111111111111111111111111333333311133111111111111111111111111333333311113311111111111111111111111133333331111333111111111111111111111133333333111133331111111111111111111133333333111113333311111111111111111133333333311111333333111111111111111133333333311111133333333111111111111333333333331111113333333333311111133333333333333311111333333333333333333333333333333333111113333333333333333333333333333333331111133333333333333333333333333333333311111333333333333333333333311333333333111111333333333333333333111113333333331111111133333333333311111111133333331111111111111111111111111111111333331111111111111111111111111111111113331111111111111111111111111111111111131111
----
3333333333111111113333333333331111111
3333333111111111111113333333331111111
3333331111111111111111333333333111111
3333311111111111111111113333333311111
3333111111111111111111113333333311111
3331111111111111111111111333333331111
3311111111111111111111111133333331111
3311111111111111111111111113333331111
3111111111111111111111111113333333111
3111111111111111111111111113333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
1111111111111111111111111111333333111
3111111111111111111111111113333333111
3111111111111111111111111113333333111
3311111111111111111111111133333331111
3311111111111111111111111133333331111
3331111111111111111111111333333331111
3333111111111111111111113333333311111
3333311111111111111111133333333311111
3333331111111111111111333333333111111
3333333311111111111133333333333111111
3333333333311111133333333333333311111
3333333333333333333333333333333331111
1333333333333333333333333333333333111
1133333333333333333333333333333333311
1113333333333333333333333113333333331
1111133333333333333333311111333333333
1111111133333333333311111111133333331
1111111111111111111111111111113333311
1111111111111111111111111111111333111
1111111111111111111111111111111131111
Avg color:
3 (253, 254, 253)
1 (13, 187, 89)
以下の画像のようなのが得られます。
無事Qiita素数が得られました。
その他
7777777777777730000669999994606000097777777777777777777777777777777900666699999946660600477777777777777777777777777777774006666999999966660000377777777777777777777777777777700066669999999666660009777777777777777777777777777777006666699999996666000097777777777777777777777777777770006666999999906666000977777777777777777777777777777700666669999999066666669777777777777777777777777777777906666619999919666666097777777777777777777777777777779006399999999999166660177777777777777777777777777777790599992999299999996006777777777777777777777777777777619999229929299999981907777777777777777777777777777773999992922392999999998937777777777777777777777777777099999999229932999999999997777777777777777777777777779999918999999999999999911997777777777777777777777777999999899999999999999999999997777777777777777777777709999918989999999919999999989997777777777777777777777999999189899999999439199999899937777777777777777777779999191898747699894499849999999977777777777777777777799991811889979998144198999991999977777777777777777773999811318899999984444988999999999577777777777777777799998334918899999344449889999999999777777777777777777999885144881999995444498888999999995777777777777777759988141444488999944444158889999999997777777777777777999881434222418981222445588899999999977777777777777779998884427622488807722444888899999999177777777777777681988854772264188722223441888999999999777777777777771899884147222345847222274418881999988997777777777777788918853471385445472891744188189999819977777777777777189988154729874444628917448811899998199777777777777773899918447722644444722664485488999981997777777777777788999985444444444444444444348889998819977777777777777889999884444444444444444444888899988999777777777777778189999844444444444444444888888999889997777777777777738889998844444411154444418888899988898677777777777777788889988814444444444458888888999887987777777777777777987589188815444444438888873889988779777777777777777777577777788891144434311187778188677777777777777777777777777777788199453499199777788777777777777777777777777777777777999899119999997777777777777777777777777777777777777779999989999919977777777777777777777777777777777777777199441899986664643777777777777777777777777777777777754666153995636665444777777777777777777777777777777777446663666366566614444777777777777777777777777777777744466666660666666044447777777777777777777777777777777444666666606666666444457777777777777777777777777777764446666666066666664444477777777777777777777777777777544566666660666666344444577777777777777777777777777774440666666666666666344444767777777777777777777777773344406666666066666661444366977777777777777777777777763444666666663666666634366669777777777777777777777777666046666666636666666666666647777777777777777777777776666666666666366666666666666677777777777777777777777766666666666666666666663666666777777777777
