画像っぽい素数をつくる

# 経緯

「２が現れる素数」という記事がありました。
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`
です。

# 処理の方針

1. 画像のリサイズ
2. 画像の減色
3. 数字のあてはめ、素数判定
4. 素数が見つかれば終了、結果表示

# 画像のリサイズ

特に難しいことは無いです。
ただ普通は縮小処理だと思いますが、ごくまれに拡大することがあるかもしれません。このとき、線形補間ではなく最近傍を用いた方が補間時に色が増えず、素数画像の見た目が良くなるかと思います。

ここで縦横ともにピクセル数を素数 $p$、$q$ にすると、$pq$ 桁の素数を渡されたときに$p$ × $q$ に並べてみるという発想に自然に行き着くかなと思います。$p = q$ だとさらに都合がいいです。

トリミングは考えてなかった。前もって適当なレタッチソフトで余白を削っている前提とします。

# 画像の減色処理

ピクセルごとに0～9の数字しかあてはめられないので、最大でも10色しか使っていない画像にする必要があります。

ここでは2通りの方法を紹介します。

* グレースケールののちLook up tableを用いる

```python
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


```python
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色に落としてみました。

元画像：![Penguins.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258019/812cb6a3-e557-ab5a-cabd-2a4f2ae2926e.jpeg)

減色：![Penguin_post.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258019/b22d990a-e702-3497-92c8-33b8fd3eb03b.png)

黄色の部分を残すのに６色必要でした。

ただ、今回の話は数字を当てはめるという関係上ピクセルごとの色のラベルがあればいいので、`center` から画像を作ることをせず、途中で関数から抜けます。

```python
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))
```

これでこの関数はピクセルごとのラベルを返します。同じラベルなら同じ色（数字）。

## 減色画像が得られたら

さて、ラベルが得られたら各ラベルが元の画像ではどんな色かを調べないといけません。

```python
# 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`にちょうどいい関数がありますので、それを使います。

```python
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`という文字列があるかどうかを確かめる関数とします。

```python
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`オプションでより厳しくする、素数と判定されたものに対しては別の判定プログラムと組み合わせるなどが考えられます。

## ちょうどいい数字割りあてのサーチ

単純なループです。

```python
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

# 実行例

次の画像を使ってみます。

![qiita-favicon.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258019/18dcdbaa-507c-327a-6322-bad8da3fde12.png)

これを縦横サイズを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)
```

以下の画像のようなのが得られます。

![prime_qiita-favicon.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258019/41086802-935e-0ad2-d94e-32e124935149.png)

無事Qiita素数が得られました。

# その他

```
7777777777777730000669999994606000097777777777777777777777777777777900666699999946660600477777777777777777777777777777774006666999999966660000377777777777777777777777777777700066669999999666660009777777777777777777777777777777006666699999996666000097777777777777777777777777777770006666999999906666000977777777777777777777777777777700666669999999066666669777777777777777777777777777777906666619999919666666097777777777777777777777777777779006399999999999166660177777777777777777777777777777790599992999299999996006777777777777777777777777777777619999229929299999981907777777777777777777777777777773999992922392999999998937777777777777777777777777777099999999229932999999999997777777777777777777777777779999918999999999999999911997777777777777777777777777999999899999999999999999999997777777777777777777777709999918989999999919999999989997777777777777777777777999999189899999999439199999899937777777777777777777779999191898747699894499849999999977777777777777777777799991811889979998144198999991999977777777777777777773999811318899999984444988999999999577777777777777777799998334918899999344449889999999999777777777777777777999885144881999995444498888999999995777777777777777759988141444488999944444158889999999997777777777777777999881434222418981222445588899999999977777777777777779998884427622488807722444888899999999177777777777777681988854772264188722223441888999999999777777777777771899884147222345847222274418881999988997777777777777788918853471385445472891744188189999819977777777777777189988154729874444628917448811899998199777777777777773899918447722644444722664485488999981997777777777777788999985444444444444444444348889998819977777777777777889999884444444444444444444888899988999777777777777778189999844444444444444444888888999889997777777777777738889998844444411154444418888899988898677777777777777788889988814444444444458888888999887987777777777777777987589188815444444438888873889988779777777777777777777577777788891144434311187778188677777777777777777777777777777788199453499199777788777777777777777777777777777777777999899119999997777777777777777777777777777777777777779999989999919977777777777777777777777777777777777777199441899986664643777777777777777777777777777777777754666153995636665444777777777777777777777777777777777446663666366566614444777777777777777777777777777777744466666660666666044447777777777777777777777777777777444666666606666666444457777777777777777777777777777764446666666066666664444477777777777777777777777777777544566666660666666344444577777777777777777777777777774440666666666666666344444767777777777777777777777773344406666666066666661444366977777777777777777777777763444666666663666666634366669777777777777777777777777666046666666636666666666666647777777777777777777777776666666666666366666666666666677777777777777777777777766666666666666666666663666666777777777777
```

![prime_sarval.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258019/d462c904-f903-e251-8765-67b727ddf59c.jpeg)