R 语言生成分形图像的几种方法

• 下载新版本代码和示例 - 2.6 KB
• 下载完整项目源码 - 23 KB

引言

在 R 语言中开发了一组非常简单但功能强大的函数，用于生成分形图像。两种主要方法是： Kronecker 积方法 (KPB) 和 IFS (迭代函数系统) 方法 (IFSB)。此外，还开发了几种算法来生成不同的布朗树。

关于分形以及上述 3 类分形的描述，请参见 [1-6]。

我的动机是：填补互联网上 R 语言分形图像（尤其是这 3 类）的空白。令我惊讶的是，我在 The R Journal [7] 上发现了几十篇关于布朗运动的文章/图片，但没有一张布朗树的图片。我认为 CodeProject (CP) 的读者对此会感兴趣，因为我通过 Google 和 R 网站的广泛搜索确认：在互联网或其他来源中，关于 R 语言中的这些类型分形，没有重要的实现/发布。同样，在 C、C++、C#、VB、Python 等语言中，有许多示例代码，包括一些 JavaScript 在线生成器，但在 R 语言中却什么也没有。

值得一提的是，这 3 种方法都可以生成几乎无限数量的分形。唯一的限制将是所用计算机的性能和用户的创造力。

在在线文章 [8] 中，介绍了一种使用 Kronecker 积应用于仅包含零和一的简单矩阵来构建分形的简单方法。作者在 SAS/IML 程序中实现了它，并构建并展示了两种分形：Vicsek 分形和“La Guardia 机场降落”分形。本文将展示在 R 语言中实现的第一个分形以及许多其他分形，包括许多原创的分形。

我已经用 PARI/GP 和 JavaScript [10-12] 实现了一个 Barnsley 蕨分形。但我偶然发现了 [13] 中这篇关于在线 IFS 分形生成器的文章，该生成器能够生成 30 多种不同的分形，包括 Barnsley 蕨。这个在线工具基于处理 IFS 表的通用方法 [10]，而这种方法已在 R 语言中实现。再次强调，我只使用了想法。源代码并未提供下载，但 IFS 表是免费提供的。有关表格示例和描述，请参见 [4,10,13]。

为了在 R 语言中实现 KPB 分形和布朗树 [9]，首先需要创建一组绘图辅助函数。这组函数最初在 PARI/GP 中创建，后来被翻译成 R。用于构建和绘制布朗树和其他分形而填充 0/1 矩阵的想法在 PARI/GP [11, 12] 中被独立“重新发明”，并很容易地将其应用于 R 语言中的函数。

绘图辅助函数

需要强调的是，尽管任何辅助函数都非常有用，但它们在许多情况下可能无法应用，或者不如直接绘图有效。例如，请看下面的 `pIfsFractal()` 函数，或者 [14] 中的 `pVoronoiD()` 函数。

只有 2 个绘图辅助函数：`plotmat()` 和 `plotv2()`，它们是真正通用的绘图函数，已经被使用多次，并证明了它们的有用性和可靠性 [12]。

让我们仔细看看这些函数。

## R Helper Functions
#
## HFR#1 plotmat(): Simple plotting using matrix mat (filled with 0/nonzero int).
#  Where: mat - matrix; fn - file name (no extension); clr - color;
#         ttl - plot title; dflg - writing dump file flag (0-no/1-yes):
#         psz - picture size; cx - cex or scale.
plotmat <- function(mat, fn, clr, ttl, dflg=0, psz=600, cx=1.0) {
  m = nrow(mat); d = 0; X=NULL; Y=NULL;
  pf = paste0(fn, ".png"); df = paste0(fn, ".dmp");
  # Building X and Y arrays for plotting from not equal to zero values in mat.
  for (i in 1:m) {
    for (j in 1:m) {if(mat[i,j]==0){next} else {d=d+1; X[d]=i; Y[d]=j} }
  };
  cat(" *** Matrix(", m,"x",m,")", d, "DOTS\n");
  # Dumping if requested (dflg=1).
  if (dflg==1) {dump(c("X","Y"), df); cat(" *** Dump file:", df, "\n")};
  # Plotting
  if (ttl!="") {
    plot(X,Y, main=ttl, axes=FALSE, xlab="", ylab="", col=clr, pch=20, cex=cx)}
  else {par(mar=c(0,0,0,0));
    plot(X,Y, axes=FALSE, xlab=NULL, ylab=NULL, col=clr, pch=20, cex=cx)};
  # Writing png-file
  dev.copy(png, filename=pf, width=psz, height=psz);
  # Cleaning
  dev.off(); graphics.off();
}

## HFR#2 plotv2(): Simple plotting using 2 vectors (dumped into ".dmp" file).
# Where: fn - file name; clr - color; ttl - plot title; psz - picture size;
#        cx - cex or scale.
plotv2 <- function(fn, clr, ttl, psz=600, cx=1.0) {
  cat(" *** START:", date(), "clr=", clr, "psz=", psz, "\n");
  pf = paste0(fn, ".png"); df = paste0(fn, ".dmp");
  source(df); d = length(X);
  cat(" *** Plot file -", pf, "Source dump-file:", df, d, "DOTS\n");
  # Plotting
  if (ttl!="") {
    plot(X,Y, main=ttl, axes=FALSE, xlab="", ylab="", col=clr, pch=20, cex=cx)}
  else {par(mar=c(0,0,0,0));
    plot(X,Y, axes=FALSE, xlab=NULL, ylab=NULL, col=clr, pch=20, cex=cx)};
  # Writing png-file
  dev.copy(png, filename=pf, width=psz, height=psz);
  # Cleaning
  dev.off(); graphics.off();
  cat(" *** END:", date(), "\n");
}

首先，它们非常小巧、简单且清晰。每个函数（不含注释）的代码都少于 15 行。因此，很容易将它们翻译成其他语言。

主要函数是 `plotmat()`。它旨在处理用 `0` 和 `1` 填充的矩阵的绘图。实际上，在当前版本中，矩阵可以填充零和任何其他整数。

此函数包含几个简单的步骤（参见上面的代码）：

• 使用双重“for”循环从 `mat` 矩阵中选择值不等于零的点/点，并为选定的点构建 X 和 Y 数组。这些数组稍后用于转储和绘图。
• 如果请求，则编写一个转储文件。
• 绘制 X,Y 点。
• 将图保存为 png 文件。

另一个函数 `plotv2()` 甚至更简单。它绘制由 `plotmat()` 创建的转储文件中的 X,Y 点。此外，它还可以是任何函数创建的任何类似的转储文件。

如果生成时间较长，请求 `plotmat()` 的转储文件会很有用。拥有一个转储文件可以轻松快速地使用 `plotv2()` 以不同的颜色、标题和大小重新绘制。

通常，`pIfsFractal()` 函数和与布朗树相关的函数生成时间较长。如果生成的点数量巨大，则绘图时间会很长。

重要说明

• 所有演示的生成/绘图函数（`pIfsFractal()` 除外）都使用 `plotmat()` 来处理 `mat` 矩阵。
• 如果使用 `plotv2()` 进行重新绘图：将使用由 `plotmat()` 创建的转储文件中的 2 个向量 X,Y。
• `plotmat()` 和 `plotv2()` 中使用的文件名不带扩展名（在需要时将添加“.png”和“.dmp”）。

基于 Kronecker 积的分形生成方法

基于 Kronecker 积的分形的起源和性质在 [3] 中有更详细的解释。

让我们仔细看看以下 2 个函数。

## Kronecker power of a matrix.
## Where: m - initial matrix, n - power.
matkronpow <- function(m, n) {
  if (n<2) {return (m)};
  r = m; n = n-1;
  for(i in 1:n) {r = r%x%m};
  return (r);
}

## Generate and plot Kronecker product based fractals.
## gpKronFractal(m, n, pf, clr, ttl, dflg, psz, cx):
## Where: m - initial matrix (filled with 0/int); n - order of the fractal;
## fn - plot file name (without extension); clr - color; ttl - plot title;
## dflg - writing dump file flag (0/1); psz - picture size; cx - cex.
gpKronFractal <- function(m, n, fn, clr, ttl, dflg=0, psz=640, cx=1.0) {
  fign="Kpbf";
  cat(" *** START:", date(), "n=", n, "clr=", clr, "psz=", psz, "\n");
  if(fn=="") {fn=paste0(fign,"o", n)} else {fn=paste0(fn)};
  if(ttl!="") {ttl=paste0(ttl,", order ", n)};
  cat(" *** Plot file -", fn, "title:", ttl, "\n");
  r = matkronpow(m, n);
  plotmat(r, fn, clr, ttl, dflg, psz, cx);
  cat(" *** END:", date(), "\n");
}

正如你所见，它们非常简单明了。`matkronpow(m, n)` 返回 m x m x m ... (n 次积)。第二个函数实际上只有 2 行代码用于生成和绘图，所有其他语句仅用于日志记录。

提示：创建和使用任何类型的辅助函数。这将有助于保持其他函数的代码简洁、清晰和稳定。

一个只有 5-10 行 R 代码的函数可能需要解释什么？

基于 Kronecker 积分形的设计/绘图示例

使用简单的文本表示初始矩阵，可以轻松设计基于 Kronecker 积的分形。在这种情况下，设计好的矩阵代表什么非常清楚。当然，不总是如此。特别是，如果初始矩阵是通过应用许多不同矩阵进行操作的几个步骤创建的。

注意：[3] 中的 `KPBFdesign.html` 页面可能非常有用。

为了展示如何设计/绘制 KPB 分形，我们选择了新的“O”分形。这个“O”分形有以下 3 个简单的步骤。

1. 基本矩阵设计

使用简单的文本表示，如下所示：

   Prime matrix   or   Inverted matrix (0 and 1 are swaped)
   1111111		             0000000
   1000001		             0111110
   1011101		             0100010
   1011101		             0100010
   1011101		             0100010
   1000001		             0111110
   1111111		             0000000

2. 将其“翻译”为如下的 1-2 行 R 代码：

  M <- matrix(c(1,1,1,1,1,1,1, 1,0,0,0,0,0,1, 1,0,1,1,1,0,1, 1,0,1,1,1,0,1,
   1,0,1,1,1,0,1, 1,0,0,0,0,0,1, 1,1,1,1,1,1,1), ncol=7, nrow=7, byrow=TRUE);

3. 在 R GUI 窗口中执行它

  gpKronFractal(M, 3, "OF3t", "maroon", "'O' fractal");

另一种设计方法是使用两个不同矩阵的 Kronecker 积创建初始矩阵，并将其应用于 `gpKronFractal()`。

下面是一个此类示例：

M <- matrix(c(0,1,0, 1,1,1, 0,1,0), ncol=3, nrow=3, byrow=TRUE);
M1 <- matrix(c(1,0,1, 0,1,0, 1,0,1), ncol=3, nrow=3, byrow=TRUE);
R=M%x%M1;
gpKronFractal(R, 2, "Crosses2F", "maroon", "2 crosses based fractal", 1);

重要说明

• 生成函数 `gpKronFractal()` 使用初始矩阵 M 来构建结果矩阵并绘制它。
• 请先尝试设置一个较低的 n，因为较大的 n 需要大量内存和时间。

以下是测试 Vicsek、Sierpinski 地毯、Rug 和“T”分形的代码。

M <- matrix(c(0,1,0, 1,1,1, 0,1,0), ncol=3, nrow=3, byrow=TRUE);
gpKronFractal(M, 4, "VicsekFractal1","red", "Vicsek Fractal");

M <- matrix(c(1,1,1, 1,0,1, 1,1,1), ncol=3, nrow=3, byrow=TRUE);
gpKronFractal(M, 4, "SierpinskiCF1", "maroon", "Sierpinski carpet fractal");

M <- matrix(c(1,1,1,1,1, 1,0,0,0,1, 1,0,0,0,1, 1,0,0,0,1, 1,1,1,1,1),
 ncol=5, nrow=5, byrow=TRUE);
gpKronFractal(M, 4, "RugF", "brown", "Rug fractal", 1);

M <- matrix(c(1,1,1,1,1, 1,1,1,0,1, 1,0,0,0,1, 1,1,1,0,1, 1,1,1,1,1),
 ncol=5, nrow=5, byrow=TRUE);
gpKronFractal(M,4,"GPFRT","darkviolet","",1,600);

注意：如果你的计算机不是超快的，那么这可能需要很长时间。使用转储标志 = 1，以保存生成的图形以供重新绘制。

说明：另一方面，在旧笔记本电脑 (WinXP) 上测试“T”分形时，我得到了非常奇怪的结果：使用 `gpKronFractal()` 比使用 `plotv2()` 快 6 倍。我认为这是由于内存不足导致了内存交换过程。

在下面的代码片段及其输出日志中自行查看。

M = matrix(c(1,1,1,1,1, 1,1,1,0,1, 1,0,0,0,1, 1,1,1,0,1, 1,1,1,1,1),ncol=5,nrow=5,byrow=T);
gpKronFractal(M,4,"GPFRT","darkviolet","",1,600,0.5);
 *** START: Thu Jun 08 16:57:12 2017 n= 4 clr= darkviolet psz= 600
 *** Plot file - GPFRT title:
 *** Matrix( 625 x 625 ) 160000 DOTS
 *** Dump file: GPFRT.dmp
 *** END: Thu Jun 08 16:59:08 2017

plotv2("GPFRT", "maroon", "'T' fractal", 640, 0.5);
 *** START: Thu Jun 08 17:00:31 2017 clr= maroon psz= 640
 *** Plot file - GPFRT.png Source dump-file: GPFRT.dmp 160000 DOTS
 *** END: Thu Jun 08 17:12:58 2017

我需要再次强调：分形的应用，以及执行时间的评估或测量，不属于我的兴趣范围。

提示：在你的计算机上测试生成/绘图的速度，并为你选择最佳函数。

你可以使用上面显示的 2 个示例“按原样”进行此类测试。

下面图 1 - 图 4 展示了所有四个生成的分形。

提示：右键单击并保存上方的任何图片，即可查看更大的图像。

技巧：“T”分形上方看起来是 3 阶的，因为只有 3 种大小的“T”。要看到它实际上是 4 阶分形，应该以不同的方式绘制。

M <- matrix(c(1,1,1,1,1, 1,1,1,0,1, 1,0,0,0,1, 1,1,1,0,1, 1,1,1,1,1),
 ncol=5,nrow=5,byrow=T);
gpKronFractal(M,4,"GPFRT2","darkviolet","",0, 1280, 0.5);
# Note: bigger size + scaling are requested.

此外，它应该被放大，例如使用 Microsoft Office Picture Manager。现在才能看到所有 4 种大小的“T”。

R 脚本的新版本

R 脚本的新版本（主要与 KPBF 绘图相关）在许多方面都得到了升级。现在 `plotmat()` 函数执行速度更快，消耗的内存更少，支持旋转，支持背景颜色，并有其他一些改进。

测试表明，在许多情况下，现在使用 `plotmat()` 函数重新绘制分形比使用转储文件更合适。

如果生成时间很长，`plotv2()` 函数仍然有用。例如，用于绘制布朗树。

基于 IFS 的分形生成方法

只有一个函数实现了基于 IFS 的方法。让我们仔细看看这个函数。

## Plotting fractals using IFS style
## Plotting is based on already calculated M x 7 table of coefficients
## in the input file.
## Note: 1. Input ifs-file should be dat-file; output is set as png-file.
##       2. Swap 2nd and 3rd column if you've got data used in Java,
##          JavaScript, etc.

## pIfsFractal(fn, n, clr, ttl, psz, cx): Plot fractals using IFS style.
## Where: fn - file name; n - number of dots; clr - color; ttl - plot title,
##        psz - plot size, cx - cex.

pIfsFractal <- function(fn, n, clr, ttl, psz=600, cx=0.5) {
  # pf - plot file name; df - data/ifs file name;
  pf=paste0(fn,".png"); df=paste0(fn,".dat");
  cat(" *** IFSSTART:", date(), "n=", n, "clr=", clr, "psz=", psz, "\n");

  # Reading a complete data table from the file: space delimited, no header.
  # Table has any number of rows, but always 7 columns is a must.
  (Tb = as.matrix(read.table(df, header=FALSE)))
  tr = nrow(Tb)
  # Creating matrix M1 from 1st 4 columns of each row.
  M1 = vector("list",tr);
  for (i in 1:tr) {M1[[i]] = matrix(c(Tb[i,1:4]),nrow=2)}
  # Creating matrix M2 from columns 5,6 of each row.
  M2 = vector("list",tr);
  for (i in 1:tr) {M2[[i]] = matrix(c(Tb[i,5:6]),nrow=2)}
  ## Creating matrix M3 (actually a vector) from column 7 of each row.
  M3 = c(Tb[1:tr,7])

  x = numeric(n); y = numeric(n);
  x[1] = y[1] = 0;

  # Main loop
  for (i in 1:(n-1)) {
    k = sample(1:tr, prob=M3, size=1);
    M = as.matrix(M1[[k]]);
    z = M%*%c(x[i],y[i]) + M2[[k]];
    x[i+1] = z[1];
    y[i+1] = z[2];
  }
  # Plotting
  if (ttl!="") {
    plot(x,y, main=ttl, axes=FALSE, xlab="", ylab="", col=clr, pch=20, cex=cx)}
  else {par(mar=c(0,0,0,0));
    plot(x,y, axes=FALSE, xlab=NULL, ylab=NULL, col=clr, pch=20, cex=cx)};
  # Writing png-file
  dev.copy(png, filename=pf,width=psz,height=psz);
  # Cleaning
  dev.off(); graphics.off();
  cat(" *** IFS  END:", date(), "\n");
}

上面的代码注释得很清楚，不需要额外的解释，但请参阅 [10] 以获取算法、表格和计算机生成的描述。

重要说明

• 生成基于输入文件中的一个已计算的 m x 7 系数表 [10]。
• 输入 ifs 文件应为 dat 文件（在下载的 zip 文件中查找示例）；输出设置为 png 文件。
• 如果你有 Java、JavaScript 等中使用的数据，请交换 Ifs 表的第 2 列和第 3 列。

测试几种分形的生成/绘图

注意：相关数据文件在 `GPFRdatafiles.txt` 文件和 `GPRFDATA` 文件夹中提供。

pIfsFractal("BarnsleyFern", 100000, "dark green", "Barnsley fern fractal", ,0.25);

pIfsFractal("Pentaflake", 100000, "blue", "Pentaflake fractal");

pIfsFractal("Sierpinski3", 100000, "red", "Sierpinski triangle fractal");

pIfsFractal("TriangleF", 100000, "maroon", "Triangle fractal");

下面图 5 - 图 8 展示了所有四个生成的分形。

生成布朗树分形

所有四个生成布朗树分形的函数都从 PARI/GP [11, 12] 翻译而来。但反过来，PARI/GP 函数是从其他语言翻译的。这个事实证明了每个函数的基本算法是正确的。

让我们仔细看看其中一个函数。

# ALL 4 versions are in GFIRALLfuncs.R
# translation of PARI/GP: http://rosettacode.org/wiki/Brownian_tree#PARI.2FGP

# Generate and plot Brownian tree. Version #1.
# gpBrownianTree1(m, n, clr, fn, ttl, dflg)
# Where: m - defines matrix m x m; n - limit of the number of moves;
#   fn - file name (.ext will be added); ttl - plot title; dflg - 0-no dump,
#   1-dump.
gpBrownianTree1 <- function(m, n, clr, fn, ttl, dflg=0) {
  cat(" *** START:", date(),"m=",m,"n=",n,"clr=",clr,"\n");
  M = matrix(c(0),ncol=m,nrow=m,byrow=T);
  # Seed in center
  x = m%/%2; y = m%/%2;
  M[x,y]=1;
  pf=paste0(fn,".png");
  cat(" *** Plot file -",pf,"\n");
  # Main loops: Generating matrix M
  for (i in 1:n) {
    if(i>1) {
      x = sample(1:m, 1, replace=F)
      y = sample(1:m, 1, replace=F)}
    while(1) {
      ox=x; oy=y;
      x = x + sample(-1:1, 1, replace=F);
      y = y + sample(-1:1, 1, replace=F);
      if(x<=m && y<=m && x>0 && y>0 && M[x,y])
        {if(ox<=m && oy<=m && ox>0 && oy>0){M[ox,oy]=1; break}}
      if(!(x<=m && y<=m && x>0 && y>0)) {break}
    }
  }
  plotmat(M, fn, clr, ttl, dflg); ## Plotting matrix M
  cat(" *** END:",date(),"\n");
}
#gpBrownianTree1(400,15000,"red", "BT13", "Brownian Tree v.1-1",1);   ## Dump (Seed in center alwys now)

其他 3 个函数可以在 `GFIRALLfuncs.R` 文件中找到。

所有四个函数都有类似的 3 个步骤：

1. 设置初始“种子”点/点（预定义或随机）。
2. 使用不同但始终只有 2 个循环：“for”和“while”来模拟“随机游走”。此步骤中的重要部分是：
   • 朝着最近的可用空间进行随机“行走”步骤
   • 确保每个“行走”步骤都在矩阵范围内
   • 将成功的“行走”步骤位置保存在矩阵中
3. 绘制代表布朗树的结果矩阵。

重要说明

• 所有绘图函数都使用 `plotmat()`。
• 因为使用“随机游走”来填充结果矩阵，- 每次所有 4 个函数都会产生一个外观不同的树（即使请求的点数和矩阵大小相同）。

测试布朗树 v.1 - v.4 的生成/绘图

gpBrownianTree1(400,15000,"red", "BTv1", "Brownian Tree v.1", 1);

gpBrownianTree2(400,5000,"brown", "BTv2", "Brownian Tree v.2", 1);

gpBrownianTree3(400,5000,"dark green", "BTv31", "Brownian Tree v.3-1", 1);

gpBrownianTree4(400,15000,"navy", "BTv41", "Brownian Tree v.4-1", 1);

下面图 9 - 图 12 展示了所有四个版本的布朗树分形。

如果你需要其他语言，请在 [9] 中查找 47 种语言的布朗树示例。

结论

再次强调：本项目演示了 R 语言中 3 种类型分形的绘图技术。在许多其他语言中也可以找到相同类型的分形 [3,8,9,11,12]，但在 R 语言中却不行。

演示的 3 种分形生成方法（连同支持函数集）可以生成几乎无限数量的分形。实际上，有无限数量的矩阵（可用于 KPB 方法），也有无限数量的初始数据表（可用于 IFSB 方法）。更不用说无限的布朗树了。

真正的限制仅仅是所用计算机的性能，当然，还有新分形设计者的创造力。

在文件 `GPFRsamples.txt` 和 `v2Samples.txt` 中，为所有 3 种类型分形提供了许多额外的测试样本，供您尝试并享受新的分形。

注意：`GPFRsamples.txt` 文件中的所有 KPBF 样本都可以在新版本中“按原样”使用。

参考文献

1. 分形，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Fractal
2. Kronecker 积，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Kronecker_product。
3. Voevudko, A.E. (2017) 生成基于 Kronecker 积的分形. Code Project。
4. 迭代函数系统，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Iterated_function_system
5. 布朗树，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Brownian_tree
6. 扩散限制聚集，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Diffusion-limited_aggregation
7. The R Journal，网址：https://journal.r-project.org/
8. Wicklin R. (2014)，从 Kronecker 积中自相似的结构，网址：http://blogs.sas.com/content/iml/2014/12/17/self-similar-structures-from-kronecker-products.html
9. 布朗树，Rosetta Code Wiki，网址：http://rosettacode.org/wiki/Brownian_tree
10. Barnsley 蕨，维基百科，自由的百科全书，网址：https://en.wikipedia.org/wiki/Barnsley_fern
11. Voevudko, A.E.，用户页面，OEIS Wiki，网址：http://oeis.org/wiki/User:Anatoly_E._Voevudko
12. Voevudko, A.E.，用户页面，Rosetta Code Wiki，网址：http://rosettacode.org/wiki/User:AnatolV
13. Toal R.，迭代函数系统，网址：http://cs.lmu.edu/~ray/notes/ifs/
14. Voevudko, A.E. (2017) 生成随机 Voronoi 图. Code Project

历史

• 2018/04/17 -- 添加包含 `plotmat()`、`plotv2()` 等升级版本的代码文件，以及包含新示例的文件。许多小编辑：主要改变“外观”。修正了一些拼写错误。添加了“历史”章节和“R 脚本新版本”子章节。
• 2017/07/05 -- 首次发布