MRaster examples 21.0.0.0
Image Processing Library
Loading...
Searching...
No Matches
mandelbrot_precomp.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_precomp.cpp
5 @author Mitch Richling <https://www.mitchr.me>
6 @brief Produce several images related to the period/cycle structure of the Mandelbrot set.@EOL
7 @std C++20
8 @copyright
9 @parblock
10 Copyright (c) 1988-2015,2021 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 Produce several images related to the period/cycle structure of the Mandelbrot set:
32
33 - mandelbrot_precompPER.tiff -- Period of the point. 0 means no period known.
34 - mandelbrot_precompSTB.tiff -- Stability of the period. 0 means stability unknown.
35 - mandelbrot_precomps.tiff -- Number of iterations required to escape. 0 means it didn't escape.
36 - mandelbrot_precompNOE.tiff -- Points that didn't escape -- the mandelbrot set. 255
37 - mandelbrot_precomp.tiff -- A composite of the above with a few notable period regions labeled.
38
39 On my 2022 vintage Intel i7, this takes about 30min to run. The runtime is directly proportional to the NUMITR, so lower that number if you want it to go
40 faster. Lowering NUMITR will have cause more non-escaping points to not have a known period -- the green points in
41
42*/
43/*******************************************************************************************************************************************************.H.E.**/
44/** @cond exj */
45
46//--------------------------------------------------------------------------------------------------------------------------------------------------------------
47#include "ramCanvas.hpp"
48
49//--------------------------------------------------------------------------------------------------------------------------------------------------------------
50typedef mjr::ramCanvas1c16b rcCNT; // for counts!
51typedef mjr::ramCanvas1c8b rcM8;
52typedef mjr::ramCanvas3c8b rcC8;
53
54//--------------------------------------------------------------------------------------------------------------------------------------------------------------
55int main(void) {
56 std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now();
57 rcC8::colorType aColor;
58
59 const rcCNT::colorChanType MAXITR = rcCNT::colorType::maxChanVal/1 - 2;
60 const int MCSIZE = 7680;
61 const int CSIZEF = 1;
62 const int CSIZE = MCSIZE/CSIZEF;
63 const double MAXZSQ = 4.0;
64 const double MAXZSQU = MAXZSQ + 0.1;
65 rcC8 theRamCanvas(CSIZE, CSIZE, -2.1, 0.75, -1.4, 1.4);
66 rcCNT perRamCanvas(CSIZE, CSIZE, -2.1, 0.75, -1.4, 1.4); // Period -- 0 => not a periodic point
67 rcCNT stbRamCanvas(CSIZE, CSIZE, -2.1, 0.75, -1.4, 1.4); // Number of iterations period structure was stable
68 rcM8 noeRamCanvas(CSIZE, CSIZE, -2.1, 0.75, -1.4, 1.4); // No escape: Exceeded iteration count
69 rcCNT escRamCanvas(CSIZE, CSIZE, -2.1, 0.75, 1.4, 1.4); // Excaped: Iteration count for cases where |z|>2
70
71#pragma omp parallel for schedule(static,1)
72 for(int y=0;y<perRamCanvas.getNumPixY();y++) {
73 rcCNT::colorChanType maxCheckpointITR = 0;
74 double meanCheckpointITR = 0;
75 std::vector<std::complex<double>> lastZs(MAXITR);
76 std::chrono::time_point<std::chrono::system_clock> rowStartTime = std::chrono::system_clock::now();
77
78 for(int x=0;x<perRamCanvas.getNumPixX();x++) {
79 std::complex<double> c(perRamCanvas.int2realX(x), perRamCanvas.int2realY(y));
80 std::complex<double> z(0.0, 0.0);
81 rcCNT::colorChanType count = 1;
82 bool seekingUnderstanding = true;
83 rcCNT::colorChanType checkpointITR = 1024*8;
84
85 while(seekingUnderstanding) {
86 if ((MAXITR - checkpointITR) < checkpointITR)
87 checkpointITR = MAXITR;
88 else
89 checkpointITR *= 2;
90
91 if (maxCheckpointITR < checkpointITR)
92 maxCheckpointITR = checkpointITR;
93
94 while((std::norm(z)<MAXZSQU) && (count<checkpointITR)) {
95 z=std::pow(z, 2) + c;
96 lastZs[count] = z;
97 count++;
98 }
99
100 if (std::norm(z)>MAXZSQ) { // Escaped
101 escRamCanvas.drawPoint(x, y, count);
102 seekingUnderstanding = false;
103 } else { // no escape. Perhaps a cycle?
104 for(rcCNT::colorChanType period=1; period<(checkpointITR-2); period++) {
105 if(std::abs(z-lastZs[checkpointITR-1-period])<1e-7) { // Found an identical point 'period' away.
106 rcCNT::colorChanType stab;
107 for(stab=0; stab<(checkpointITR-period); stab++) {
108 if(std::abs(lastZs[checkpointITR-1-stab]-lastZs[checkpointITR-1-period-stab])>1e-7) {
109 break;
110 }
111 }
112 if (stab > period) { // Definatly found a cycle
113 stbRamCanvas.drawPoint(x, y, checkpointITR-stab);
114 perRamCanvas.drawPoint(x, y, period);
115 noeRamCanvas.drawPoint(x, y, "white");
116 seekingUnderstanding = false;
117 break;
118 }
119 }
120 }
121 if (seekingUnderstanding && (checkpointITR == MAXITR)) { // Didn't escape. No cycle. No more time.
122 noeRamCanvas.drawPoint(x, y, "white");
123 seekingUnderstanding = false;
124 break;
125 }
126 }
127 }
128 meanCheckpointITR += checkpointITR / static_cast<double>(CSIZE);
129
130 }
131 std::chrono::duration<double> rowRunTime = std::chrono::system_clock::now() - rowStartTime;
132 std::cout << "my: " << CSIZE << " y: " << y << " max: " << maxCheckpointITR << " mean: " << meanCheckpointITR << " secs: " << rowRunTime.count() << std::endl;
133 }
134
135 perRamCanvas.writeTIFFfile("mandelbrot_precompPER.tiff");
136 stbRamCanvas.writeTIFFfile("mandelbrot_precompSTB.tiff");
137 escRamCanvas.writeTIFFfile("mandelbrot_precompESC.tiff");
138 noeRamCanvas.writeTIFFfile("mandelbrot_precompNOE.tiff");
139
140 int numConNoCyc = 0;
141 int numCyc = 0;
142 int numEsc = 0;
143 int numPts = 0;
144 int maxPer = 0;
145 for(int y=0;y<theRamCanvas.getNumPixY();y++) {
146 for(int x=0;x<theRamCanvas.getNumPixX();x++) {
147 auto period = perRamCanvas.getPxColorNC(x, y).getC0();
148 if (period > maxPer)
149 maxPer = period;
150 if (period > 0) {
151 if (period > (rcC8::colorType::csCBDark2::maxNumC-1)) {
152 theRamCanvas.drawPoint(x, y, "red");
153 } else {
154 theRamCanvas.drawPoint(x, y, rcC8::colorType::csCBDark2::c(period));
155 }
156 numCyc++;
157 } else {
158 if (noeRamCanvas.getPxColorNC(x, y).getC0() > 0) {
159 theRamCanvas.drawPoint(x, y, "green");
160 numConNoCyc++;
161 } else {
162 rcC8::csFltType c = static_cast<rcC8::csFltType>(escRamCanvas.getPxColorNC(x, y).getC0()) / MAXITR;
163 theRamCanvas.drawPoint(x, y, rcC8::colorType::csCCdiag01::c(c*30));
164 numEsc++;
165 }
166 }
167 numPts++;
168 }
169 }
170 theRamCanvas.drawString("1", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.1500, 0.0000, "black", 6.0/CSIZEF, 20);
171 theRamCanvas.drawString("2", mjr::hershey::font::ROMAN_SL_SANSERIF, -1.0000, 0.0000, "black", 6.0/CSIZEF, 20);
172
173 theRamCanvas.drawString("3", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.1200, 0.7400, "black", 6.0/CSIZEF, 20);
174 theRamCanvas.drawString("3", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.1200, -0.7400, "black", 6.0/CSIZEF, 20);
175 theRamCanvas.drawString("3", mjr::hershey::font::ROMAN_SL_SANSERIF, -1.758, 0.0000, "black", 1.0/CSIZEF, 20);
176
177 theRamCanvas.drawString("4", mjr::hershey::font::ROMAN_SL_SANSERIF, -1.3100, 0.0000, "black", 6.0/CSIZEF, 20);
178 theRamCanvas.drawString("4", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.2800, 0.5320, "black", 6.0/CSIZEF, 20);
179 theRamCanvas.drawString("4", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.2800, -0.5320, "black", 6.0/CSIZEF, 20);
180
181 theRamCanvas.drawString("5", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.5060, -0.5620, "black", 4.0/CSIZEF, 20);
182 theRamCanvas.drawString("5", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3780, -0.3370, "black", 4.0/CSIZEF, 20);
183 theRamCanvas.drawString("5", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.5060, 0.5620, "black", 4.0/CSIZEF, 20);
184 theRamCanvas.drawString("5", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3780, 0.3370, "black", 4.0/CSIZEF, 20);
185
186 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3890, -0.2150, "black", 2.0/CSIZEF, 20);
187 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.1140, -0.8620, "black", 4.0/CSIZEF, 20);
188 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, -1.1370, -0.2390, "black", 4.0/CSIZEF, 20);
189 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3890, 0.2150, "black", 2.0/CSIZEF, 20);
190 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.1140, 0.8620, "black", 4.0/CSIZEF, 20);
191 theRamCanvas.drawString("6", mjr::hershey::font::ROMAN_SL_SANSERIF, -1.1370, 0.2390, "black", 4.0/CSIZEF, 20);
192
193 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.6220, -0.4240, "black", 2.0/CSIZEF, 20);
194 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.1210, -0.6100, "black", 2.0/CSIZEF, 20);
195 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3760, -0.1440, "black", 1.0/CSIZEF, 20);
196 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, -0.6220, 0.4240, "black", 2.0/CSIZEF, 20);
197 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.1210, 0.6100, "black", 2.0/CSIZEF, 20);
198 theRamCanvas.drawString("7", mjr::hershey::font::ROMAN_SL_SANSERIF, 0.3760, 0.1440, "black", 1.0/CSIZEF, 20);
199
200 theRamCanvas.writeTIFFfile("mandelbrot_precomp.tiff");
201
202 std::cout << "numConNoCyc ... " << numConNoCyc << std::endl;
203 std::cout << "numCyc ........ " << numCyc << std::endl;
204 std::cout << "numInSet ...... " << numConNoCyc+numCyc << std::endl;
205 std::cout << "numEsc ........ " << numEsc << std::endl;
206 std::cout << "numPts ........ " << numPts << std::endl;
207 std::cout << "maxPer ........ " << maxPer << std::endl;
208
209 std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime;
210 std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl;
211}
212/** @endcond */
int main(int argc, char *argv[])