按素数绘画（2018）

两个星期前，我偶然发现了主要肖像的概念。简而言之，原始人像是为每种颜色分配了编号的图片。当我们对齐每个像素时，结果应该是质数。因为找不到本文中使用的代码，所以我决定用Python重新创建这些主要肖像-彩色！以下是Jupyter笔记本，其中包含我用来生成《蒙娜丽莎》的代码（也可在Github上找到）：

素数肖像是素数，格式为矩阵，每行X位数字。当我们为每个数字选择颜色时，我们可以生成图像。

在找到与已知图像相似的东西之前，我没有针对许多质数和配色方案进行此操作，而是改过了这个过程。我拍摄了标志性的图像，例如《蒙娜丽莎》和《星夜》，并将它们转换为仅10种颜色的图像。我为每种颜色分配了一个数字。然后，我生成了许多类似的图像，并添加了一些“噪音”。噪音稍微改变了图像中的颜色，从而改变了数字。如果图像中的数字形成质数，我发现了质数！

导入matplotlib.pyplot作为plt ％matplotlib内联 将numpy导入为np 导入skimage.transform 导入skimage.io 从sklearn.cluster导入KMeans 导入scipy.misc 从PIL导入Image，ImageFont，ImageDraw 将numpy导入为np 从多处理导入池 导入时间 随机导入 进口警告

当使用极长的素数工作时，测试每个除数的幼稚方法会花费太多时间。质数很多，检查每个可能除数的数字需要很长时间才能测试许多质数肖像。

相反，我们对素数采取概率方法。 Miller Rabin检验使用启发式方法确定数字是否（可能）是质数。我从Rosetta代码下载了以下Miller-Rabin素数测试代码，以检查素数：

def try_composite（a，d，n，s）： ＃测试基数a以查看它是否是n的复合性的见证 如果pow（a，d，n）== 1： 返回False 对于范围（s）中的i： 如果pow（a，2 ** i * d，n）== n-1： 返回False 返回True＃n绝对是复合的 def is_probable_prime（n，mrpt_num_trials = 5）： ＆＃34;＆＃34;＆＃34; Miller-Rabin素数检验。 返回值为False意味着n肯定不是素数。返回值为 True表示n很可能是素数。 ＆＃34;＆＃34;＆＃34; 断言n＆gt; = 2 ＃特例2 如果n == 2： 返回True ＃确保n为奇数 如果n％2 == 0： 返回False ＃将n-1写为2 ** s * d ＃反复尝试将n-1除以2 s = 0 d = n-1 而True： 商，余数= divmod（d，2） 如果余数== 1： 打破 s + = 1 d =商 断言（2 ** s * d == n-1） 对于我在范围内（mrpt_num_trials）： a =随机的。 randrange（2，n） 如果try_composite（a，d，n，s）： 返回False 返回True＃没有测试的基础显示n为复合 断言is_probable_prime（2） 断言is_probable_prime（3） 断言不是is_probable_prime（4） 断言is_probable_prime（5） 断言不是is_probable_prime（123456789） primes_under_1000 = [如果is_probable_prime（i）表示i在范围（2，1000）中的i] 断言len（primes_under_1000）== 168 断言primes_under_1000 [-10：] == [937，941，947，953，967，971，977，983，991，997] 断言is_probable_prime（6438080068035535544392301298549614926991513861075340134 \ 3291807343952413826484237063006136971539473913409092293733259059038472039 \ 7133335969549256322620979036686633213903952966175107096769180017646161 \ 851573147596390153） 断言不是is_probable_prime（7438080068035535592392301298549614926991513861075340134 \ 3291807343952413826484237063006136971539473913409092293733259059038472039 \ 7133335969549256322620979036686633213903952966175107096769180017646161 \ 851573147596390153）

该算法的基础是生成仅具有10种可能颜色值的图像。我不知道有什么通用方法，所以我认为应用k-means算法会行得通。为此，我对输入图像进行了整形（因此每个像素都是3D空间中的一个点），对其进行了聚类，然后将创建的10个聚类用作颜色值。旁注：当您的锤子是机器学习工具时，一切看起来都像钉子一样……

对于素数生成而言，重要的是要向原始图像添加一点噪点，以创建看起来相同但略有不同的图像。这意味着您创建了一个不同的数字，这是我们找到质数的新机会。

def get_k_means（图片）： pointcloud = np。重塑（image，（-1，3）） kmeans = KMeans（n_clusters = 10）。适合（pointcloud） 返回kmeans def create_numbered_image（image，kmeans）： ＆＃34;＆＃34;＆＃34; 将RGB图像转换为带有提供的kmeans分类器的簇标签的图像。 ＆＃34;＆＃34;＆＃34; ＃应用噪音 stdev = np。标准（图片） random_noise = np。随机的。 random_sample（大小=图片。形状）*（stdev / 3） 图片=图片+ random_noise orig_shape =图片。形状 图片= np。重塑（image，（-1，3）） numbered_image = kmeans。预测（图片） numbered_image = np。重塑（numbered_image，orig_shape [：2]） ＃确保结尾不均匀 如果numbered_image [-1，-1]％2 == 0： numbered_image [-1，-1] + = 1 返回numbered_image def numbered_image_to_normal_image（numbered_image，kmeans）： ＆＃34;＆＃34;＆＃34; 通过使用聚类中心将只有0到9之间的值的图像转换为彩色图像 提供的kmeans分类器。 ＆＃34;＆＃34;＆＃34; shape =（numbered_image。shape [0]，numbered_image。shape [1]，3） 图片= np。零（形状） 对于标签，以zip格式显示的颜色（范围（10），kmeans。cluster_centers_）： 图片[numbered_image ==标签] =颜色 返回图片

在这里，我创建了一些蒙娜丽莎小数字图像的示例。如果仔细观察，您会发现每张图像都略有不同，但是您仍然可以识别原始图片。这对我们至关重要，因为这使我们能够生成可以测试素数的多个数字化图像。

def load_and_resize_image（filename，resize_factor = 18）： image = skimage。 io。读（文件名） 图片=图片。 astype（np。float64） 图片=图片/ 255。 oldshape =图片。形状 带有警告。 catch_warnings（）： 警告。 simplefilter（＆＃34; ignore＆＃34;） resized_image = skimage。转变 。调整大小（image，（oldshape [0] // resize_factor，oldshape [1] // resize_factor）） 返回resized_image resize_image = load_and_resize_image（＆＃39; input_pictures / monalisa.jpg＆＃39;） kmeans = get_k_means（resized_image） N_IMAGES = 3 f，axarr = plt。子图（1，N_IMAGES） 对于范围（N_IMAGES）中的索引： n_image = create_numbered_image（resize_image，kmeans） normal_image = numbered_image_to_normal_image（n_image，kmeans） 轴[索引]。 imshow（normal_image）

我必须创建一些辅助函数来返回字符串，整数以及加载和显示图像。

def image_to_number（numbered_image）： to_be_number = numbered_image。重塑（-1） as_string =＆＃39;＆＃39; 。加入（[to_be_number中x的[str（int（x））]） as_int = int（as_string） 返回as_int，as_string def show_and_save_image（image，n_image，filename，fontsize = 16）： oldshape =图片。形状 resized_image = scipy。杂项imresize（image，（oldshape [0] * fontsize，oldshape [1] * fontsize），interp =＆＃39; nearest＆＃39;） img =图片。 fromarray（resized_image）。转换（＆＃34; RGBA＆＃34;） txt =图片。新（＆＃39; RGBA＆＃39;，img。size，（255，255，255，0）） 画= ImageDraw。画（txt） 字体= ImageFont。 truetype（＆＃34; pirulen rg.ttf＆＃34;，fontsize） 对于y_i，输入枚举（n_image）： 对于x_i，枚举（totype）中的字母： xpos = x_i *字体大小+ 1 ypos = y_i *字体大小 如果字母== 1： xpos + = 4 画 。文本（（xpos，ypos），str（字母），（255，255，255，128），字体=字体） img =图片。 alpha_composite（img，txt） img。保存（文件名） 请图（figsize =（20，20）） 请imshow（img） 请节目 （） def result_filename（filename）： 返回文件名。分割（＆＃39;。＆＃39;）[0] +＆＃34; -prime.png＆＃34;

为了验证某物是否是好的基本画像，我进行了以下辅助功能：

def is_good_prime_portrait（n_image）： 整数，字符串= image_to_number（n_image） 如果is_probable_prime（integer）： 返回integer，string，n_image 其他： 不返回 def print_result（string，n_image）： 打印（字符串） 打印（＆＃34;-＆＃34; * 100） 对于n_image中的行： 打印（＆＃39;＆＃39;。join（[x的行的str（x）]））

最初，我创建了一个单线程解决方案，该解决方案起作用了。但是，我的笔记本电脑有多个核心，可以同时处理多个数字，因此我以多线程的方式重写了最初的方法：

def print_result（string，n_image）： 打印（＆＃34;找到结果：＆＃34; +＆＃34;-＆＃34; * 100） 打印（字符串） 打印（＆＃34;代表肖像：＆＃34; +＆＃34;-＆＃34; * 100） 对于n_image中的行： print（＆＃39;＆＃39;。join（[x的行的str（x）]）） def multi_threaded_prime_generator（resized_image，kmeans，filename，线程= 4，log_process = True）： image_generator =（范围（1000000）中的_的create_numbered_image（resize_image，kmeans）） 开始=时间。时间 （） 与池（线程）作为池： 结果=池。 imap_unordered（is_good_prime_portrait，image_generator） total_results = 0 为了结果： total_results + = 1 ＃可能会花时间搜索此素数 如果log_process和total_results％30 == 0： 经过的时间=时间。时间 （） 过去=过去-开始 打印（＆＃34;在（函数名）中花费的秒数是{}每个结果的时间：{}＆＃34;。format（str（elapsed），str（elapsed / total_results））） 如果结果！=无： ＃找到素数，打印并保存！ 整数，字符串，n_image =结果 print_result（string，n_image） normal_image = numbered_image_to_normal_image（n_image，kmeans） 请imshow（normal_image） 请节目 （） show_and_save_image（normal_image，n_image，result_filename（filename）） 打破 def search_prime_portrait（filename，resize_factor = 16，log_process = True，线程= 4）： resize_image = load_and_resize_image（filename，resize_factor = resize_factor） 打印（＆＃34;使用大小＆＃34; + str（resized_image。shape）） kmeans = get_k_means（resized_image） multi_threaded_prime_generator（resize_image，kmeans，文件名，log_process = log_process，线程=线程）

现在，我们定义了所有功能，终于可以生成基本肖像了！试一试，尝试输入自己的图像。请注意，您可以将resize_factor变量传递给函数“ search prime portrait”。 resize_factor越大，找到基本画像的速度就越快。有关生成更大的素像所需的指数时间的更多信息，请参见下文。

文件夹=＆＃39; input_pictures /＆＃39; 文件名= [ 文件夹+＆＃39monalisa.jpg＆＃39; ， #folder +＆＃39; sunflowers.jpg＆＃39 ;， #folder +＆＃39; starrynight.jpg＆＃39; ] 文件名=文件名[0] 用于文件名中的文件名： search_prime_portrait（文件名，resize_factor = 20，log_process = False，线程= 4）

使用尺寸（19、13、3） 发现结果：---------------------------------------------- -------------------------------------------------- ---- 666666666666266626626222222233122222222239514622206668573106228000175818000800085381000000668877180203338899538673333159957818287188370138838188188188188188818118118888311814414101055555514411181111757388411111333881144814118814114114444 代表肖像：---------------------------------------------- -------------------------------------------------- ---- 6666666666662 6666262662622 2222233122222 2222395146222 0666857310622 8000175818000 8000853810000 0066887718020 3338899538673 3331599578182 8718837011483 8788188113188 0888181181188 8831181441410 1185551441118 1111757388411 1113338811448 1411881444411 4111444444441

显然，最大的已知质数长为22338618数字。这将使我们能够创建4726 * 4726原始图像，这将是4K中的原始图像。我已经很喜欢“星空之夜”的大图了，我绝对希望拥有这张画的4K黄金肖像。

但是：我注意到找到质数并不能真正扩展。质数的长度与图像的分辨率成平方增长，这会引起很大的问题。尽管找到质数的概率仅随质数的长度线性降低，但是当您获得较大的图像分辨率时，其概率呈二次方降低。我在这张图中绘制了它，以显示高达300×300的效果：

但更糟糕的是，计算时间随着质数的长度呈二次方增长。我还为笔记本电脑绘制了高达30×30的图像分辨率的图像：

def time_expected_img（w，h）： 阶数= w * h 开始=时间。时间 （） 对于_范围（100）： posprim =随机的。 randrange（10 **（订单-1），10 **订单） a = is_probable_prime（posprim） 结束=时间。时间 （） 返回结束-开始 test_range =范围（2，30，4） time_needed = [test_range中x的time_expected_img（x，x）] 请情节（test_range，time_needed） 请标题（＆＃34;预计第二步，以计算XxY图像是否是基本肖像＆＃34;） 请xlabel（＆＃34; Width image＆＃34;） 请ylabel（＆＃34; Expected seconds＆＃34;）

结合两条线性曲线，您可以轻松估算找到一定尺寸的肖像所需的预期时间：

total_needed = [对于test_range中的x，需要trials_expected_img（x，x）* time_expected_img（x，x）] 请情节（test_range，total_needed） 请标题（＆＃34;预计第二次会找到XxY素像&& 34;） 请xlabel（＆＃34; Width image＆＃34;） 请ylabel（＆＃34; Expected seconds＆＃34;）

如果您可以找到一个更高效的概率素数检查器Python实现，或者以更适合于高要求计算的语言来实现它，则可以减少生成这些肖像所需的时间。每个额外的计算机内核也可以减少此时间，但是目前我们不能期望很快有4K素像。

在过去的几周里，我真的很喜欢让我的计算机查找基本肖像。这是我的计算机上找到的一些艺术品：梵高的向日葵Hokusai Girl的《大波浪》和Vermeer的珍珠耳环

尽管现有艺术品的图片通常会使我们感觉很好，但是要创建美观的主要人物肖像却相当困难。最终的图像通常对比度不足，仅使用十种颜色。我决定为我喜欢在网上关注的人尝试一些肖像。使我们感到满意的图像是：

每个星期六在荷兰报纸上写有趣的文章，不久前在荷兰漫画《唐老鸭》中登载：

在Numberphile Youtube频道上制作有关数学概念的有趣视频。另外：有一只可爱的狗😉

当我将代码发布到网上时，您可以尝试为自己生成一些基本的肖像！如果您真的认为我应该尝试某些艺术品，请随时与我联系。另外：如果您有个好主意如何处理这些结果，请随时与我们联系。我知道我个人将在墙上挂满黄金星夜！