1#' Mandelbrot convergence counts
2#'
3#' @param Z               A complex matrix for which convergence
4#'                        counts should be calculated.
5#' @param xmid,ymid,side,resolution Alternative specification of
6#'                        the complex plane `Z`, where
7#'                        `mean(Re(Z)) == xmid`,
8#'                        `mean(Im(Z)) == ymid`,
9#'                        `diff(range(Re(Z))) == side`,
10#'                        `diff(range(Im(Z))) == side`, and
11#'                        `dim(Z) == c(resolution, resolution)`.
12#' @param maxIter         Maximum number of iterations per bin.
13#' @param tau             A threshold; the radius when calling
14#'                        divergence (Mod(z) > tau).
15#'
16#' @return Returns an integer matrix (of class Mandelbrot) with
17#' non-negative counts.
18#'
19#' @examples
20#' counts <- mandelbrot(xmid = -0.75, ymid = 0, side = 3)
21#' str(counts)
22#' \dontrun{
23#' plot(counts)
24#' }
25#'
26#' \dontrun{
27#' demo("mandelbrot", package = "future", ask = FALSE)
28#' }
29#'
30#' @author The internal Mandelbrot algorithm was inspired by and
31#' adopted from similar GPL code of Martin Maechler available
32#' from ftp://stat.ethz.ch/U/maechler/R/ on 2005-02-18 (sic!).
33#'
34#' @aliases as.raster.Mandelbrot plot.Mandelbrot mandelbrot_tiles
35#' @export
36#'
37#' @keywords internal
38mandelbrot <- function(...) UseMethod("mandelbrot")
39
40#' @export
41mandelbrot.matrix <- function(Z, maxIter = 200L, tau = 2.0, ...) {
42  stop_if_not(is.matrix(Z), mode(Z) == "complex")
43
44  ## By default, assume none of the elements will converge
45  counts <- matrix(maxIter, nrow = nrow(Z), ncol = ncol(Z))
46
47  ## But as a start, mark all be non-diverged
48  idx_of_non_diverged <- seq_along(Z)
49
50  ## SPEEDUP: The Mandelbrot sequence will only be calculated on the
51  ## "remaining set" of complex numbers that yet hasn't diverged.
52  sZ <- Z ## The Mandelbrot sequence of the "remaining" set
53  Zr <- Z ## The original complex number of the "remaining" set
54
55  for (ii in seq_len(maxIter - 1L)) {
56    sZ <- sZ * sZ + Zr
57
58    ## Did any of the "remaining" points diverge?
59    diverged <- (Mod(sZ) > tau)
60    if (any(diverged)) {
61      ## Record at what iteration divergence occurred
62      counts[idx_of_non_diverged[diverged]] <- ii
63
64      ## Early stopping?
65      keep <- which(!diverged)
66      if (length(keep) == 0) break
67
68      ## Drop from remain calculations
69      idx_of_non_diverged <- idx_of_non_diverged[keep]
70
71      ## Update the "remaining" set of complex numbers
72      sZ <- sZ[keep]
73      Zr <- Zr[keep]
74    }
75  }
76
77  attr(counts, "params") <- list(Z = Z, maxIter = maxIter, tau = tau)
78
79  class(counts) <- c("Mandelbrot", class(counts))
80
81  counts
82}
83
84
85#' @export
86mandelbrot.numeric <- function(xmid = -0.75, ymid = 0.0, side = 3.0,
87                               resolution = 400L, maxIter = 200L,
88                               tau = 2.0, ...) {
89  ## Validate arguments
90  stop_if_not(side > 0)
91  resolution <- as.integer(resolution)
92  stop_if_not(resolution > 0)
93
94  maxIter <- as.integer(maxIter)
95  stop_if_not(maxIter > 0)
96
97  ## The nx-by-ny bins
98  nx <- ny <- resolution
99
100  ## Setup (x, y) bins
101  xrange <- xmid + c(-1, 1) * side / 2
102  yrange <- ymid + c(-1, 1) * side / 2
103  x <- seq(from = xrange[1], to = xrange[2], length.out = nx)
104  y <- seq(from = yrange[1], to = yrange[2], length.out = ny)
105
106  ## Set of complex numbers to be investigated
107  Z <- outer(y, x, FUN = function(y, x) complex(real = x, imaginary = y))
108
109  mandelbrot(Z, maxIter = maxIter, tau = tau)
110}
111
112
113#' @export
114#' @importFrom grDevices as.raster hsv
115#' @keywords internal
116as.raster.Mandelbrot <- function(x, ...) {
117  maxIter <- attr(x, "params", exact = TRUE)$maxIter
118  img <- hsv(h = x / maxIter, s = 1, v = 1)
119  img[x == maxIter] <- "#000000"
120  dim(img) <- dim(x)
121  img <- t(img)
122  img <- structure(img, class = "raster")
123  img
124}
125
126
127#' @export
128#' @importFrom grDevices as.raster
129#' @importFrom graphics par plot
130#' @keywords internal
131plot.Mandelbrot <- function(x, y, ..., mar = c(0, 0, 0, 0)) {
132  if (!is.null(mar)) {
133    opar <- par(mar = c(0, 0, 0, 0))
134    on.exit(par(opar))
135  }
136  plot(as.raster(x), ...)
137}
138
139
140#' @export
141mandelbrot_tiles <- function(xmid = -0.75, ymid = 0.0, side = 3.0,
142                             nrow = 2L, ncol = nrow,
143                             resolution = 400L, truncate = TRUE) {
144  ## Validate arguments
145  stop_if_not(side > 0)
146  resolution <- as.integer(resolution)
147  stop_if_not(resolution > 0)
148
149  ## The nx-by-ny bins
150  nx <- ny <- resolution
151
152  ## Bins per tile
153  dx <- ceiling(nx / ncol)
154  dy <- ceiling(ny / nrow)
155  stop_if_not(dx > 0, dy > 0)
156
157  ## Truncate so all tiles have identical dimensions?
158  if (truncate) {
159    nx <- ncol * dx
160    ny <- nrow * dy
161  }
162
163  ## Setup (x, y) bins
164  xrange <- xmid + c(-1, 1) * side / 2
165  yrange <- ymid + c(-1, 1) * side / 2
166  x <- seq(from = xrange[1], to = xrange[2], length.out = nx)
167  y <- seq(from = yrange[1], to = yrange[2], length.out = ny)
168
169
170  ## Generate tiles row by row
171  res <- list()
172  for (rr in seq_len(nrow)) {
173    yrr <- if (rr < nrow) y[1:dy] else y
174    y <- y[-(1:dy)]
175
176    xrr <- x
177    for (cc in seq_len(ncol)) {
178      xcc <- if (cc < ncol) xrr[1:dx] else xrr
179      xrr <- xrr[-(1:dx)]
180
181      Ccc <- outer(yrr, xcc, FUN = function(y, x) {
182        complex(real = x, imaginary = y)
183      })
184      attr(Ccc, "region") <- list(xrange = range(xcc), yrange = range(yrr))
185      attr(Ccc, "tile") <- c(rr, cc)
186      res <- c(res, list(Ccc))
187    }
188  }
189  dim(res) <- c(nrow, ncol)
190
191  res
192}
193