MRaster examples 23.0.0.0
Image Processing Library
Loading...
Searching...
No Matches
poly_root_cloud.cpp
Go to the documentation of this file.
1// -*- Mode:C++; Coding:us-ascii-unix; fill-column:158 -*-
2/*******************************************************************************************************************************************************.H.S.**/
3/**
4 @file poly_root_cloud.cpp
5 @author Mitch Richling http://www.mitchr.me/
6 @brief Draw the roots of paramaterized polynomials. @EOL
7 @keywords polynomial root
8 @std C++23
9 @see https://github.com/richmit/mraster/
10 @copyright
11 @parblock
12 Copyright (c) 2025, Mitchell Jay Richling <http://www.mitchr.me/> All rights reserved.
13
14 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
15
16 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
17
18 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
19 and/or other materials provided with the distribution.
20
21 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
22 without specific prior written permission.
23
24 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
26 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28 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
29 DAMAGE.
30 @endparblock
31 @filedetails
32
33 Plot roots of the polynomial p. Two of the coefficients of p (c1 & c2) are parameterized by polynomials (p2 & p2) over the unit circle.
34
35 In pseudo-code:
36 Loop
37 Select sample points on the unit circle for t1 & t2
38 Plug t1 & t2 into p1 & p2
39 Set coefficients c1 & c2 in p to the value of p1 & p2
40 Find the roots of p and plot them
41 Repeat
42
43 I made a prototype for this problem in Matlab. In the Matlab code I placed semitransparent filled circles at each root, and allowed color to build up with
44 repeated layers of roots. This results in a "fuzzy" look spreading color more smoothly over the envelope of the root cloud. That code may be found here:
45 https://github.com/richmit/PolyRootCloud/
46
47 In this code I wanted to emphasize more directly the way the circles in the parametrization map to the root cloud space. So I'm actually creating a
48 histogram of root hits in the plane with a very fine grid, and using many more sample points on one parameter. This produces a Spirograph'esq effect.
49*/
50/*******************************************************************************************************************************************************.H.E.**/
51/** @cond exj */
52
53//--------------------------------------------------------------------------------------------------------------------------------------------------------------
54#include "ramCanvas.hpp"
55#include "MRMathUPLY.hpp"
56
57//--------------------------------------------------------------------------------------------------------------------------------------------------------------
58// Overide the complex types in lapacke.h.
59#define LAPACK_COMPLEX_CUSTOM 1
60#define lapack_complex_float std::complex<float>
61#define lapack_complex_double std::complex<double>
62
63#include <lapacke.h> /* LAPACK C Interface LAPACKE */
64
65typedef lapack_complex_double cplx;
66
67//--------------------------------------------------------------------------------------------------------------------------------------------------------------
68
69typedef mjr::ramCanvas1c16b drc_t;
70typedef mjr::color3c8b oc_t;
71typedef mjr::ramCanvasPixelFilter::FuncHomoTransform<drc_t, oc_t> hpf_t;
72
73//--------------------------------------------------------------------------------------------------------------------------------------------------------------
74int main(void) {
75 std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now();
76
77 const int n1 = 1000;
78 const int n2 = 100;
79 const int c1 = 8;
80 const int c2 = 12;
81 const int XSIZ = 2048*8;
82 const int YSIZ = 2048*8;
83 drc_t theRamCanvas(XSIZ, YSIZ, -2.5, 2.5, -2.5, 2.5);
84
85# pragma omp parallel for schedule(static,1)
86 for(int m1=0; m1<n1; m1++) {
87 std::vector<cplx> p = {{1.0, 0.0}, {-1.0, 0.0}, {0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, {0.0, 0.0},
88 {0.0, 0.0}, { 0.0, 0.0}, {0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, {0.0, 0.0},
89 {0.0, 0.0}, { 0.0, 0.0}, {0.0, 0.0}, {-0.1, 0.0}, { 0.1, 0.0}};
90 std::vector<cplx> p1 = {{0.0, 0.0}, { 0.0, 100.0}, {0.0, -100.0}, { 0.0, 100.0}, { 0.0, -100.0}, {-100.0, 0.0}, {0.0, 100.0}};
91 std::vector<cplx> p2 = {{0.0, -100.0}, { 0.0, -100.0}, {0.0, 100.0}, { 0.0, 100.0}, {100.0, 0.0}};
92 int pdg = (int)p.size() - 1;
93 lapack_int n = (lapack_int)p.size() - 1; // Companion matrix dimensions are the polynomial degree
94 lapack_int lwork = -1; // Set to -1 for work size query
95 cplx *work; // Allocated later after query
96 cplx *a = new cplx[n*n];
97 cplx *w = new cplx[n];
98 double *rwork = new double[n*2];
99 cplx tmp;
100 // Query the optimal value for lwork & allocate it
101 if( 0 != LAPACKE_zgeev_work(LAPACK_COL_MAJOR, 'N', 'N', n, a, n, w, NULL, n, NULL, n, &tmp, lwork, rwork)) {
102 std::cout << "ZGEEV Error (lwork query)" << std::endl;
103 exit(1);
104 }
105 lwork = LAPACK_Z2INT(tmp) + 1;
106 work = new cplx[lwork];
107 // Update p
108 cplx t1 = std::exp(cplx(0.0, 1.0) * static_cast<double>(m1) * 2.0 * std::numbers::pi / static_cast<double>(n1));
109 p[c1] = mjr::math::uply::eval(p1, t1);
110 for(int m2=0; m2<n2; m2++) {
111 // Update p
112 cplx t2 = std::exp(cplx(0.0, 1.0) * static_cast<double>(m2) * 2.0 * std::numbers::pi / static_cast<double>(n2));
113 p[c2] = mjr::math::uply::eval(p2, t2);
114 // Compute the companion matrix and store it in a
115 for(int i=0;i<n*n;i++)
116 a[i] = cplx(0.0, 0.0);
117 for(int i=1;i<n;i++)
118 a[i*n+(i-1)] = cplx(1.0, 0.0);
119 for(int i=1;i<=n;i++)
120 a[(n-i)*n+(n-1)] = - p[i] / p[0];
121 // Now find the eigenvalues.
122 if( 0 != LAPACKE_zgeev_work(LAPACK_COL_MAJOR, 'N', 'N', n, a, n, w, NULL, n, NULL, n, work, lwork, rwork)) {
123 std::cout << "ZGEEV Error (eigenvalue)" << std::endl;
124 exit(1);
125 }
126 // Plot all the roots
127# pragma omp critical
128 for(int j=0; j<pdg; j++)
129 if (theRamCanvas.getPxColor(w[j]).getC0() < 1275)
130 theRamCanvas.incPxChan(w[j]);
131 }
132 // Update on progress
133# pragma omp critical
134 std::cout << "Completed Cycle " << m1 << " of " << n1 << std::endl;
135 // Free the memory for new'ed arrays.
136 delete[] a;
137 delete[] w;
138 delete[] work;
139 delete[] rwork;
140 }
141
142 // Transform for better colors
143 theRamCanvas.applyHomoPixTfrm(&drc_t::colorType::tfrmStdPow, 0.5);
144 theRamCanvas.scaleDownMax(8);
145
146 // Write with color map
147 theRamCanvas.writeTIFFfile("poly_root_cloud.tiff", hpf_t(theRamCanvas,[](auto inColor) { return oc_t::csCCfractal0RYBCW::c(inColor.getC0()); }));
148
149 // Wrap things up
150 std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime;
151 std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl;
152 return 0;
153}
154/** @endcond */
main
int main(int argc, char *argv[])
Definition convertRawToTIFF.cpp:37
Generated by doxygen 1.14.0