R: Mandelbrot (Rcpp 0.12.6)

ふとおもいついて、マンデルブロー集合を描く昔のコードを最近のR 3.3.1とRcp 0.12.6で うごくようになおしてみました。

[2016-07-30] バグ修正

library(Rcpp)
library(inline)

## C++ source
mandelbrot.src <- "
Rcpp::NumericVector x(xx);
Rcpp::NumericVector y(yy);
unsigned short t = Rcpp::as<unsigned short>(threshold);
int nx = x.size(), ny = y.size();
Rcpp::NumericVector u(nx);

if (nx != ny) {
  return Rcpp::wrap(-1);
} else {
  for (int i = 0; i < nx; i++) {
    std::complex<double> z(0.0, 0.0);
    std::complex<double> c(x[i], y[i]);
    unsigned short k = 0;

    while (k < t) {
      z = z * z + c;
      if (abs(z) > 2.0) break;
      k++;
    }
    u[i] = k;
  }
  return Rcpp::wrap(u);
}
"

## Call C++
mandelbrot.c <- rcpp(signature(xx = "numeric",
                               yy = "numeric",
                               threshold = "numeric"),
                     mandelbrot.src, 
                     includes = "#include <complex>")


## R function
mandelbrot <- function(min.x, max.x, min.y, max.y,
                       px = 256,
                       py = round((max.y - min.y) / (max.x - min.x) * px),
                       threshold = 255) {
  if (min.x > max.x) {
    a <- min.x
    min.x <- max.x
    max.x <- a
  }
  if (min.y > max.y) {
    a <- min.y
    min.y <- max.y
    max.y <- a
  }
  if (px < 0) px <- -px
  if (py < 0) py <- -py

  rx <- seq(min.x, max.x, length = px)
  ry <- seq(min.y, max.y, length = py)
  xx <- rep(rx, py)
  yy <- rep(ry, each = px)

  z <- mandelbrot.c(xx, yy, threshold)
  matz <- matrix(z, nrow = px)
}

## sample parameters
px <- 900
py <- 900
threshold <- 1023
min.x <- -0.7625
max.x <- -0.7605
min.y <- -0.0860
max.y <- -0.0840
matz <- mandelbrot(min.x, max.x, min.y, max.y, px, py,
                   threshold = threshold)

png(file = "Mandelbrot%03d.png", width = 900, height = 900)
image(seq(min.x, max.x, length = nrow(matz)),
      seq(min.y, max.y, length = ncol(matz)),
      matz,
      col = rep(rainbow(256), (threshold + 1) / 256),
      xlab = "x", ylab = "y")
dev.off()

サンプル画像

タグ：R