MRaster examples 21.0.0.0
Image Processing Library
Loading...
Searching...
No Matches
mandelbrot_distance.cpp
Go to the documentation of this file.
1// -*- Mode:C++; Coding:us-ascii-unix; fill-column:158 -*-
2/*******************************************************************************************************************************************************.H.S.**/
3/**
4 @file mandelbrot_distance.cpp
5 @author Mitch Richling <https://www.mitchr.me>
6 @brief This program draws a mandelbrot set using the "distance".@EOL
7 @std C++20
8 @copyright
9 @parblock
10 Copyright (c) 1988-2022, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
11
12 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
13
14 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
15
16 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
17 and/or other materials provided with the distribution.
18
19 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
20 without specific prior written permission.
21
22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
24 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
27 DAMAGE.
28 @endparblock
29 @filedetails
30
31 The normal Mandelbrot raster scan iteration fails to illustrate the thin filaments projecting off the Mandelbrot set -- because the complex number at the
32 center of the pixel misses the set, and thus the entire rectangular region of the complex plane covered by the pixel is marked as "escaped". In order to
33 provide some graphical representation of these filaments, we can use a the Milnor and Thurston distance estimator as presented in 'The Science of Fractal
34 Images' to draw pixels that have a center "close" to the Mandelbrot set. This program uses a simple color scheme to label points as described in the
35 comment for the pixelStateEnum enumeration. All of the drawing code is confined to a switch statement at the end of the x loop to make it easy to play
36 around with various color schemes.
37
38 Reference: Peitgen, Heinz-Otto and Saupe, Dietmar; 1988; The Science of Fractal Images; ISBN: 0-387-96608-0; pp 192-196
39*/
40/*******************************************************************************************************************************************************.H.E.**/
41/** @cond exj */
42
43//--------------------------------------------------------------------------------------------------------------------------------------------------------------
44#include "ramCanvas.hpp"
45
46enum class pixelStateEnum {UNKNOWN, // Intermediate state -> ESCAPED || MAXITR1
47 P1BULB, // In the period 1 bulb ....................................................................... COL: RED
48 P2BULB, // In the period 2 bulb ....................................................................... COL: MAGENTA
49 MAXITR1, // Final state: Assumed in the set because it didn't escape aftter MAXITR1 iterations ......... COL: CYAN
50 ESCAPED, // Intermediate state -> ESCAPED_BUT_CLOSE || ESCAPED_BUT_FAR || ESCAPED_UNKNOWN || MAXITR2
51 ESCAPED_BUT_CLOSE, // Final state: escaped (>2), but the point is close to the set ............................... COL: BLUE
52 ESCAPED_BUT_FAR, // Final state: escaped (>2), and the point is NOT close to the set ........................... COL: YELLOW
53 ESCAPED_UNKNOWN, // Final state: escaped (>2), but we could not compute distance (can't happen) ................ COL: WHITE
54 MAXITR2}; // Final state: escaped (>2), but we never got to MAXZSQ2 (shouldn't happen) .................. COL: GREEN
55
56//--------------------------------------------------------------------------------------------------------------------------------------------------------------
57int main(void) {
58 std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now();
59 constexpr int IMXSIZ = 7680;
60 constexpr int IMYSIZ = 7680;
61 constexpr double DISTTH = 0.0002;
62 constexpr int MAXITR2 = 2048;
63 constexpr double MAXZSQ2 = 10000.0;
64 constexpr int MAXITR1 = 2048;
65 constexpr double MAXZSQ1 = 4.0;
66 constexpr double ZROEPS = 1.0e-5;
67 mjr::ramCanvas3c8b theRamCanvas(IMXSIZ, IMYSIZ, -1.9, 0.5, -1.2, 1.2);
68
69# pragma omp parallel for schedule(static,1)
70 for(int y=0; y<theRamCanvas.getNumPixY(); y++) {
71 for(int x=0; x<theRamCanvas.getNumPixX(); x++) {
72 pixelStateEnum pixelState = pixelStateEnum::UNKNOWN;
73 std::complex<double> c = theRamCanvas.int2real(x, y);
74 double p = std::abs(c-0.25);
75 if(std::abs(c+1.0) < 0.25) {
76 pixelState = pixelStateEnum::P2BULB;
77 } else if(std::real(c) < p-2.0*p*p+0.25) {
78 pixelState = pixelStateEnum::P1BULB;
79 } else {
80 std::complex<double> der(0, 0), z(0, 0);
81 int count = 0;
82 // First we do at most MAXITR1 standard iterations to see if the point seems to be in the set.
83 while(pixelState == pixelStateEnum::UNKNOWN) {
84 der = 2.0 * z * der + 1.0;
85 z = z * z + c;
86 if (std::norm(z) > MAXZSQ1) {
87 pixelState = pixelStateEnum::ESCAPED;
88 } else if (count > MAXITR1) {
89 pixelState = pixelStateEnum::MAXITR1;
90 }
91 count++;
92 }
93 // If it escaped, we see if it was close.
94 while(pixelState == pixelStateEnum::ESCAPED) {
95 der = 2.0 * z * der + 1.0;
96 z = z * z + c;
97 double absz = std::abs(z);
98 if (absz > MAXZSQ2) {
99 double der_mag = std::abs(der);
100 if(der_mag > ZROEPS) {
101 double dist = 2.0*std::log(absz)*absz/der_mag;
102 if(dist < DISTTH) {
103 pixelState = pixelStateEnum::ESCAPED_BUT_CLOSE;
104 } else {
105 pixelState = pixelStateEnum::ESCAPED_BUT_FAR;
106 }
107 } else {
108 pixelState = pixelStateEnum::ESCAPED_UNKNOWN;
109 }
110 } else if (count > MAXITR2) {
111 pixelState = pixelStateEnum::MAXITR2;
112 }
113 count++;
114 }
115 }
116 // Now we should have a good pixelState, so we color the pixel.
117 switch(pixelState) {
118 case pixelStateEnum::P2BULB: theRamCanvas.drawPoint(x, y, "magenta"); break;
119 case pixelStateEnum::P1BULB: theRamCanvas.drawPoint(x, y, "red"); break;
120 case pixelStateEnum::MAXITR1: theRamCanvas.drawPoint(x, y, "cyan"); break;
121 case pixelStateEnum::ESCAPED_BUT_CLOSE: theRamCanvas.drawPoint(x, y, "blue"); break;
122 case pixelStateEnum::ESCAPED_BUT_FAR: theRamCanvas.drawPoint(x, y, "yellow"); break;
123 case pixelStateEnum::ESCAPED_UNKNOWN: theRamCanvas.drawPoint(x, y, "white"); break;
124 case pixelStateEnum::MAXITR2: theRamCanvas.drawPoint(x, y, "green"); break;
125 default: theRamCanvas.drawPoint(x, y, "black"); break;
126 }
127 }
128 }
129 theRamCanvas.writeTIFFfile("mandelbrot_distance.tiff");
130 std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime;
131 std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl;
132}
133/** @endcond */
main
int main(int argc, char *argv[])
Definition convertRawToTIFF.cpp:37
Generated by doxygen 1.11.0