diff --git a/libid/engine/PertEngine.cpp b/libid/engine/PertEngine.cpp index 9336b3e55..3de1e1f6c 100644 --- a/libid/engine/PertEngine.cpp +++ b/libid/engine/PertEngine.cpp @@ -9,9 +9,11 @@ // #include "engine/PertEngine.h" +#include "engine/perturbation.h" #include "engine/calcfrac.h" #include "engine/id_data.h" +#include "engine/bailout_formula.h" #include "fractals/fractalp.h" #include "fractals/pickover_mandelbrot.h" #include "math/biginit.h" @@ -28,7 +30,7 @@ // Raising this number makes more calculations, but less variation between each calculation (less chance // of mis-identifying a glitched point). -constexpr double GLITCH_TOLERANCE{1e-6}; +double GLITCH_TOLERANCE = g_perturbation_tolerance; void PertEngine::initialize_frame( const BFComplex ¢er_bf, const std::complex ¢er, double zoom_radius) @@ -42,6 +44,8 @@ void PertEngine::initialize_frame( m_center_bf.y = alloc_stack(g_bf_length + 2); copy_bf(m_center_bf.x, center_bf.x); copy_bf(m_center_bf.y, center_bf.y); + m_c_bf.x = alloc_stack(g_r_bf_length + 2); + m_c_bf.y = alloc_stack(g_r_bf_length + 2); } else { @@ -58,27 +62,12 @@ int PertEngine::calculate_one_frame() std::complex reference_coordinate; m_reference_points = 0; - m_glitch_point_count = 0L; - m_remaining_point_count = 0L; - - m_points_remaining.resize(g_screen_x_dots * g_screen_y_dots); - m_glitch_points.resize(g_screen_x_dots * g_screen_y_dots); m_perturbation_tolerance_check.resize(g_max_iterations * 2); m_xn.resize(g_max_iterations + 1); // calculate the pascal's triangle coefficients for powers > 3 pascal_triangle(); - // Fill the list of points with all points in the image. - for (long y = 0; y < g_screen_y_dots; y++) - { - for (long x = 0; x < g_screen_x_dots; x++) - { - Point pt(x, g_screen_y_dots - 1 - y); - m_points_remaining[y * g_screen_x_dots + x] = pt; - m_remaining_point_count++; - } - } - + double magnified_radius = m_zoom_radius; int window_radius = std::min(g_screen_x_dots, g_screen_y_dots); BigStackSaver saved; @@ -93,117 +82,35 @@ int PertEngine::calculate_one_frame() tmp_bf = alloc_stack(g_r_bf_length + 2); } - while (m_remaining_point_count > (g_screen_x_dots * g_screen_y_dots) * (m_percent_glitch_tolerance / 100)) - { - m_reference_points++; + m_reference_points++; // Determine the reference point to calculate // Check whether this is the first time running the loop. - if (m_reference_points == 1) - { - if (g_bf_math != BFMathType::NONE) - { - copy_bf(c_bf.x, m_center_bf.x); - copy_bf(c_bf.y, m_center_bf.y); - copy_bf(reference_coordinate_bf.x, c_bf.x); - copy_bf(reference_coordinate_bf.y, c_bf.y); - } - else - { - c = m_center; - reference_coordinate = m_center; - } - - m_delta_real = 0; - m_delta_imag = 0; - - if (g_bf_math != BFMathType::NONE) - { - reference_zoom_point(reference_coordinate_bf, g_max_iterations); - } - else - { - reference_zoom_point(reference_coordinate, g_max_iterations); - } - } - else - { - if (!m_calculate_glitches) - { - break; - } - - std::srand(g_random_seed); - if (!g_random_seed_flag) - { - ++g_random_seed; - } - const int index{(int) ((double) std::rand() / RAND_MAX * m_remaining_point_count)}; - Point pt{m_points_remaining[index]}; - // Get the complex point at the chosen reference point - double delta_real = magnified_radius * (2 * pt.get_x() - g_screen_x_dots) / window_radius; - double delta_imag = -magnified_radius * (2 * pt.get_y() - g_screen_y_dots) / window_radius; - - // We need to store this offset because the formula we use to convert pixels into a complex point - // does so relative to the center of the image. We need to offset that calculation when our - // reference point isn't in the center. The actual offsetting is done in calculate point. - - m_delta_real = delta_real; - m_delta_imag = delta_imag; - - if (g_bf_math != BFMathType::NONE) - { - float_to_bf(tmp_bf, delta_real); - add_bf(reference_coordinate_bf.x, c_bf.x, tmp_bf); - float_to_bf(tmp_bf, delta_imag); - add_bf(reference_coordinate_bf.y, c_bf.y, tmp_bf); - } - else - { - reference_coordinate.real(c.real() + delta_real); - reference_coordinate.imag(c.imag() + delta_imag); - } - - if (g_bf_math != BFMathType::NONE) - { - reference_zoom_point(reference_coordinate_bf, g_max_iterations); - } - else - { - reference_zoom_point(reference_coordinate, g_max_iterations); - } - } + if (g_bf_math != BFMathType::NONE) + { + copy_bf(c_bf.x, m_center_bf.x); + copy_bf(c_bf.y, m_center_bf.y); + copy_bf(reference_coordinate_bf.x, c_bf.x); + copy_bf(reference_coordinate_bf.y, c_bf.y); + } + else + { + c = m_center; + reference_coordinate = m_center; + } - int last_checked{-1}; - m_glitch_point_count = 0; - for (int i = 0; i < m_remaining_point_count; i++) - { - if (i % 1000 == 0 && driver_key_pressed()) - { - return -1; - } - Point pt{m_points_remaining[i]}; - if (calculate_point(pt, magnified_radius, window_radius) < 0) - { - return -1; - } - // Everything else in this loop is just for updating the progress counter. - double progress = (double) i / m_remaining_point_count; - if (int(progress * 100) != last_checked) - { - last_checked = int(progress * 100); - m_status = "Pass: " + std::to_string(m_reference_points) + ", Ref (" + - std::to_string(int(progress * 100)) + "%)"; - } - } + m_delta_real = 0; + m_delta_imag = 0; - // These points are glitched, so we need to mark them for recalculation. We need to recalculate them - // using Pauldelbrot's glitch fixing method (see calculate point). - std::memcpy(m_points_remaining.data(), m_glitch_points.data(), sizeof(Point) * m_glitch_point_count); - m_remaining_point_count = m_glitch_point_count; + if (g_bf_math != BFMathType::NONE) + { + reference_zoom_point(reference_coordinate_bf, g_max_iterations); } - - cleanup(); + else + { + reference_zoom_point(reference_coordinate, g_max_iterations); + } + restore_stack(m_saved_stack); return 0; } @@ -213,235 +120,10 @@ void PertEngine::cleanup() { restore_stack(m_saved_stack); } - m_points_remaining.clear(); - m_glitch_points.clear(); m_perturbation_tolerance_check.clear(); m_xn.clear(); } -int PertEngine::calculate_point(const Point &pt, double magnified_radius, int window_radius) -{ - // Get the complex number at this pixel. - // This calculates the number relative to the reference point, so we need to translate that to the center - // when the reference point isn't in the center. That's why for the first reference, - // m_calculated_real_delta and m_calculated_imaginary_delta are 0: it's calculating relative to the - // center. - const double delta_real = - ((magnified_radius * (2 * pt.get_x() - g_screen_x_dots)) / window_radius) - m_delta_real; - const double delta_imaginary = - ((-magnified_radius * (2 * pt.get_y() - g_screen_y_dots)) / window_radius) - m_delta_imag; - std::complex delta_sub_0{delta_real, delta_imaginary}; - std::complex delta_sub_n{delta_real, delta_imaginary}; - int iteration{}; - bool glitched{}; - - double min_orbit{1e5}; // orbit value closest to origin - long min_index{}; // iteration of min_orbit - double magnitude; - do - { - if (g_cur_fractal_specific->pert_pt == nullptr) - { - throw std::runtime_error("No perturbation point function defined for fractal type (" + - std::string{g_cur_fractal_specific->name} + ")"); - } - g_cur_fractal_specific->pert_pt(m_xn[iteration], delta_sub_n, delta_sub_0); - iteration++; - magnitude = mag_squared(m_xn[iteration] + delta_sub_n); - - if (g_inside_color == BOF60 || g_inside_color == BOF61) - { - if (magnitude < min_orbit) - { - min_orbit = magnitude; - min_index = iteration + 1L; - } - } - - // This is Pauldelbrot's glitch detection method. You can see it here: - // http://www.fractalforums.com/announcements-and-news/pertubation-theory-glitches-improvement/. As - // for why it looks so weird, it's because I've squared both sides of his equation and moved the - // |ZsubN| to the other side to be precalculated. For more information, look at where the reference - // point is calculated. I also only want to store this point once. - if (m_calculate_glitches && !glitched && - magnitude < m_perturbation_tolerance_check[iteration]) - { - m_glitch_points[m_glitch_point_count] = Point(pt.get_x(), pt.get_y(), iteration); - m_glitch_point_count++; - glitched = true; - break; - } - } while (magnitude < g_magnitude_limit && iteration < g_max_iterations); - - if (!glitched) - { - int index; - const double rq_lim2{std::sqrt(g_magnitude_limit)}; - const std::complex w{m_xn[iteration] + delta_sub_n}; - - if (g_biomorph >= 0) - { - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - if (std::abs(w.real()) < rq_lim2 || std::abs(w.imag()) < rq_lim2) - { - index = g_biomorph; - } - else - { - index = iteration % 256; - } - } - } - else - { - switch (g_outside_color) - { - case 0: // no filter - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = iteration % 256; - } - break; - - case ZMAG: - if (iteration == g_max_iterations) - { - index = (int) ((w.real() * w.real() + w.imag() + w.imag()) * (g_max_iterations >> 1) + 1); - } - else - { - index = iteration % 256; - } - break; - - case REAL: - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = iteration + (long) w.real() + 7; - } - break; - - case IMAG: - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = iteration + (long) w.imag() + 7; - } - break; - - case MULT: - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else if (w.imag()) - { - index = (long) ((double) iteration * (w.real() / w.imag())); - } - else - { - index = iteration; - } - break; - - case SUM: - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = iteration + (long) (w.real() + w.imag()); - } - break; - - case ATAN: - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = (long) std::abs(atan2(w.imag(), w.real()) * 180.0 / PI); - } - break; - - default: - if (g_potential_flag) - { - index = potential(mag_squared(w), iteration); - } - else // no filter - { - if (iteration == g_max_iterations) - { - index = g_max_iterations; - } - else - { - index = iteration % 256; - } - } - break; - } - - if (g_inside_color >= 0) // no filter - { - if (iteration == g_max_iterations) - { - index = g_inside_color; - } - else - { - index = iteration % 256; - } - } - else - { - switch (g_inside_color) - { - case ZMAG: - if (iteration == g_max_iterations) - { - index = (int) (mag_squared(w) * (g_max_iterations >> 1) + 1); - } - break; - case BOF60: - if (iteration == g_max_iterations) - { - index = (int) (std::sqrt(min_orbit) * 75.0); - } - break; - case BOF61: - if (iteration == g_max_iterations) - { - index = min_index; - } - break; - } - } - } - g_plot(pt.get_x(), g_screen_y_dots - 1 - pt.get_y(), index); - } - return 0; -} - void PertEngine::reference_zoom_point(const BFComplex ¢er, int max_iteration) { BigStackSaver saved; @@ -529,3 +211,171 @@ void PertEngine::reference_zoom_point(const std::complex ¢er, int ma g_cur_fractal_specific->pert_ref(center, z); } } + +int PertEngine::perturbation_per_pixel(int x, int y, double bailout) +{ + double magnified_radius = m_zoom_radius; + int window_radius = std::min(g_screen_x_dots, g_screen_y_dots); + +// if (m_glitches[x + y * g_screen_x_dots] == 0) // processed and not glitched +// return -2; + + const double delta_real = ((magnified_radius * (2 * x - g_screen_x_dots)) / window_radius) - m_delta_real; + const double delta_imaginary = ((-magnified_radius * (2 * y - g_screen_y_dots)) / window_radius) - m_delta_imag; + m_delta_sub_0 = {delta_real, -delta_imaginary}; + m_delta_sub_n = m_delta_sub_0; +// m_points_count++; +// m_glitched = false; + return 0; +} + +int PertEngine::calculate_orbit(int x, int y, long iteration) +{ + // Get the complex number at this pixel. + // This calculates the number relative to the reference point, so we need to translate that to the center + // when the reference point isn't in the center. That's why for the first reference, calculatedRealDelta + // and calculatedImaginaryDelta are 0: it's calculating relative to the center. + + std::complex temp, temp1; + double magnitude; + + temp1 = m_xn[iteration - 1] + m_delta_sub_n; + g_cur_fractal_specific->pert_pt(m_xn[iteration - 1], m_delta_sub_n, m_delta_sub_0); + + if (g_cur_fractal_specific->pert_pt == nullptr) + { + throw std::runtime_error("No perturbation point function defined for fractal type (" + + std::string{g_cur_fractal_specific->name} + ")"); + } + temp = m_xn[iteration] + m_delta_sub_n; + magnitude = mag_squared(temp); + + // This is Pauldelbrot's glitch detection method. You can see it here: + // http://www.fractalforums.com/announcements-and-news/pertubation-theory-glitches-improvement/. As for + // why it looks so weird, it's because I've squared both sides of his equation and moved the |ZsubN| to + // the other side to be precalculated. For more information, look at where the reference point is + // calculated. I also only want to store this point once. + // if (magnitude < m_perturbation_tolerance_check[iteration]) + { + // here is where the magic happens... eventually + +// calculate_reference(x, y); +/* + perturbation_per_pixel(x, y, g_magnitude_limit); + for (long i = 0; i < g_max_iterations; i++) + { + temp1 = m_xn[i] + m_delta_sub_n; + g_cur_fractal_specific->pert_pt(m_xn[i], m_delta_sub_n, m_delta_sub_0); + + if (g_cur_fractal_specific->pert_pt == nullptr) + { + throw std::runtime_error("No perturbation point function defined for fractal type (" + + std::string{g_cur_fractal_specific->name} + ")"); + } + temp = m_xn[iteration] + m_delta_sub_n; + int status = g_bailout_float(); + if (status == 1) + { + g_color_iter = 24; + break; + } + } +*/ +// return g_bailout_float(); +// m_reference_points++; + // return true; + } + g_new_z.x = temp.real(); + g_new_z.y = temp.imag(); + // the following are needed because although perturbation operates in doubles, the math type may not be + if (g_bf_math == BFMathType::BIG_NUM) + { + float_to_bn(g_new_z_bn.x, temp.real()); + float_to_bn(g_new_z_bn.y, temp.imag()); + } + else if (g_bf_math == BFMathType::BIG_FLT) + { + float_to_bf(g_new_z_bf.x, temp.real()); + float_to_bf(g_new_z_bf.y, temp.imag()); + } + return g_bailout_float(); +} + +int PertEngine::calculate_reference(int x, int y) +{ + double magnified_radius = m_zoom_radius; + int window_radius = std::min(g_screen_x_dots, g_screen_y_dots); + BFComplex reference_coordinate_bf{}; + BigFloat tmp_bf{}; + std::complex reference_coordinate; + int saved = save_stack(); + + m_reference_points++; + + // if (m_calculate_glitches == false) +// return 0; + + if (g_bf_math != BFMathType::NONE) + { + reference_coordinate_bf.x = alloc_stack(g_r_bf_length + 2); + reference_coordinate_bf.y = alloc_stack(g_r_bf_length + 2); + tmp_bf = alloc_stack(g_r_bf_length + 2); + } +/* + std::srand(g_random_seed); + if (!g_random_seed_flag) + { + ++g_random_seed; + } + + + const int index{(int) ((double) std::rand() / RAND_MAX * m_remaining_point_count)}; + Point pt{m_points_remaining[index]}; +*/ + // Get the complex point at the chosen reference point +// double delta_real = magnified_radius * (2 * pt.get_x() - g_screen_x_dots) / window_radius; +// double delta_imag = -magnified_radius * (2 * pt.get_y() - g_screen_y_dots) / window_radius; + double delta_real = magnified_radius * (2 * x - g_screen_x_dots) / window_radius; + double delta_imag = -magnified_radius * (2 * y - g_screen_y_dots) / window_radius; + + // We need to store this offset because the formula we use to convert pixels into a complex point + // does so relative to the center of the image. We need to offset that calculation when our + // reference point isn't in the center. The actual offsetting is done in calculate point. + + m_delta_real = delta_real; + m_delta_imag = delta_imag; + + if (g_bf_math != BFMathType::NONE) + { + float_to_bf(tmp_bf, delta_real); + add_bf(reference_coordinate_bf.x, m_c_bf.x, tmp_bf); + float_to_bf(tmp_bf, delta_imag); + sub_bf(reference_coordinate_bf.y, m_c_bf.y, tmp_bf); + } + else + { + reference_coordinate.real(m_c.real() + delta_real); + reference_coordinate.imag(m_c.imag() - delta_imag); + } + + if (g_bf_math != BFMathType::NONE) + { + reference_zoom_point(reference_coordinate_bf, g_max_iterations); + } + else + { + reference_zoom_point(reference_coordinate, g_max_iterations); + } + + if (g_bf_math != BFMathType::NONE) + { + restore_stack(saved); + } + + return 0; +} + +int PertEngine::get_number_references() +{ + return m_reference_points; +} diff --git a/libid/engine/calcfrac.cpp b/libid/engine/calcfrac.cpp index 143a61b8b..970edb093 100644 --- a/libid/engine/calcfrac.cpp +++ b/libid/engine/calcfrac.cpp @@ -35,6 +35,7 @@ #include "engine/tesseral.h" #include "engine/wait_until.h" #include "engine/work_list.h" +#include "engine/perturbation.h" #include "fractals/fractalp.h" #include "fractals/frothy_basin.h" #include "fractals/lyapunov.h" @@ -881,6 +882,10 @@ static void perform_work_list() bool (*sv_per_image)() = nullptr; // once-per-image setup int alt = find_alternate_math(g_fractal_type, g_bf_math); + if (g_use_perturbation) // setup per-image perturbation + { + mandel_perturbation_setup(); + } if (alt > -1) { sv_orbit_calc = g_cur_fractal_specific->orbit_calc; @@ -1149,10 +1154,10 @@ static void perform_work_list() case 'g': // TODO: fix this // horrible cludge preventing crash when coming back from perturbation and math = bignum/bigflt - if (g_calc_status != CalcStatus::COMPLETED) - { +// if (g_calc_status != CalcStatus::COMPLETED) +// { solid_guess(); - } +// } break; case 'd': @@ -1431,6 +1436,24 @@ int standard_fractal() // per pixel 1/2/b/g, called with row & col set } g_overflow = false; // reset integer math overflow flag + if (g_use_perturbation) + { +// int temp_color = get_color(g_col, g_row); + perturbation_per_pixel(); // initialize the pert calculations +// if (perturbation_per_pixel() == -2) // initialize the calculations -2 means not glitched +// { +// g_color_iter = temp_color; // we have done this pixel in an earlier pass and it's not glitched// +// g_color = std::abs((int) g_color_iter); +// g_pixel_is_complete = true; +// return g_color; +// } + } + else + { + g_cur_fractal_specific->per_pixel(); // initialize the calculations + } + + g_cur_fractal_specific->per_pixel(); // initialize the calculations attracted = false; @@ -1531,9 +1554,22 @@ int standard_fractal() // per pixel 1/2/b/g, called with row & col set } // the usual case - else if ((g_cur_fractal_specific->orbit_calc() && g_inside_color != STAR_TRAIL) || g_overflow) + else { - break; + if (g_use_perturbation) + { + if ((perturbation_per_orbit() && g_inside_color != STAR_TRAIL) || g_overflow) + { + break; + } + } + else + { + if ((g_cur_fractal_specific->orbit_calc() && g_inside_color != STAR_TRAIL) || g_overflow) + { + break; + } + } } if (g_show_orbit) { diff --git a/libid/engine/fractalb.cpp b/libid/engine/fractalb.cpp index a4a0275f8..89a6660ff 100644 --- a/libid/engine/fractalb.cpp +++ b/libid/engine/fractalb.cpp @@ -9,6 +9,7 @@ #include "engine/calcfrac.h" #include "engine/fractals.h" #include "engine/id_data.h" +#include "engine/perturbation.h" #include "engine/type_has_param.h" #include "fractals/divide_brot.h" #include "fractals/fractalp.h" @@ -356,15 +357,6 @@ bool mandel_bn_setup() } } - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - mandel_perturbation_setup(); - // TODO: figure out crash if we don't do this - g_std_calc_mode ='g'; - g_calc_status = CalcStatus::COMPLETED; - return true; - } - g_c_exponent = (int) g_params[2]; switch (g_fractal_type) { @@ -475,10 +467,6 @@ bool mandel_bf_setup() { case FractalType::MANDEL: case FractalType::BURNING_SHIP: - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - return mandel_perturbation_setup(); - } break; case FractalType::JULIA: @@ -487,19 +475,6 @@ bool mandel_bf_setup() break; case FractalType::MANDEL_Z_POWER: - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - // only allow integer values of real part - if (const int degree = (int) g_params[2]; degree > 2) - { - return mandel_z_power_perturbation_setup(); - } - else if (degree == 2) - { - return mandel_perturbation_setup(); - } - } - init_big_pi(); if ((double) g_c_exponent == g_params[2] && (g_c_exponent & 1)) // odd exponents { @@ -536,7 +511,7 @@ bool mandel_bf_setup() int mandel_bn_per_pixel() { - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) + if (g_use_perturbation && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) { return true; } diff --git a/libid/engine/perturbation.cpp b/libid/engine/perturbation.cpp index 753dbdcea..78bcbb914 100644 --- a/libid/engine/perturbation.cpp +++ b/libid/engine/perturbation.cpp @@ -7,6 +7,7 @@ #include "engine/PertEngine.h" #include "engine/convert_center_mag.h" #include "engine/id_data.h" +#include "fractals/fractalp.h" #include "math/biginit.h" #include @@ -16,6 +17,11 @@ static PertEngine s_pert_engine; +PerturbationMode g_perturbation{PerturbationMode::AUTO}; +double g_perturbation_tolerance{1e-6}; +bool g_use_perturbation{}; // select perturbation code +int g_number_referrences{}; // number of references used + bool perturbation() { BigStackSaver saved; @@ -59,3 +65,44 @@ bool perturbation() g_calc_status = CalcStatus::COMPLETED; return false; } + +int perturbation_per_pixel() +{ + int result = s_pert_engine.perturbation_per_pixel(g_col, g_row, g_magnitude_limit); + return result; +} + +int perturbation_per_orbit() +{ + int status = s_pert_engine.calculate_orbit(g_col, g_row, g_color_iter); + return status; +} + +int get_number_references() +{ + return s_pert_engine.get_number_references(); +} + +bool is_perturbation() +{ + if (bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) + { + if (g_perturbation == PerturbationMode::AUTO && g_bf_math != BFMathType::NONE) + { + g_use_perturbation = true; + } + else if (g_perturbation == PerturbationMode::YES) + { + g_use_perturbation = true; + } + else + { + g_use_perturbation = false; + } + } + else + { + g_use_perturbation = false; + } + return g_use_perturbation; +} diff --git a/libid/fractals/frasetup.cpp b/libid/fractals/frasetup.cpp index a18c83b9c..3fb1562cb 100644 --- a/libid/fractals/frasetup.cpp +++ b/libid/fractals/frasetup.cpp @@ -30,10 +30,6 @@ // Mandelbrot Routine bool mandel_setup() { - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - return mandel_perturbation_setup(); - } // use the main processing loop g_calc_type = standard_fractal; return true; @@ -48,13 +44,7 @@ standalone_setup() bool mandel_perturbation_setup() { - return perturbation(); -} - -bool mandel_z_power_perturbation_setup() -{ - constexpr int MAX_POWER{28}; - g_c_exponent = std::min(std::max(g_c_exponent, 2), MAX_POWER); + g_calc_type = standard_fractal; return perturbation(); } @@ -91,10 +81,7 @@ mandel_fp_setup() calcmandfp() can currently handle invert, any rqlim, potflag zmag, epsilon cross, and all the current outside options */ - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - return mandel_perturbation_setup(); - } + if (g_debug_flag != DebugFlags::FORCE_STANDARD_FRACTAL && !g_distance_estimator && g_decomp[0] == 0 @@ -104,6 +91,7 @@ mandel_fp_setup() && g_use_init_orbit != InitOrbitMode::VALUE && (g_sound_flag & SOUNDFLAG_ORBIT_MASK) < SOUNDFLAG_X && !g_using_jiim + && !g_use_perturbation && g_bailout_test == Bailout::MOD && (g_orbit_save_flags & OSF_MIDI) == 0) { @@ -118,17 +106,6 @@ mandel_fp_setup() break; case FractalType::MANDEL_Z_POWER: - if (g_std_calc_mode == 'p' && bit_set(g_cur_fractal_specific->flags, FractalFlags::PERTURB)) - { - if (g_c_exponent == 2) - { - return mandel_perturbation_setup(); - } - if (g_c_exponent > 2) - { - return mandel_z_power_perturbation_setup(); - } - } if ((double)g_c_exponent == g_params[2] && (g_c_exponent & 1)) // odd exponents { g_symmetry = SymmetryType::XY_AXIS_NO_PARAM; diff --git a/libid/include/engine/PertEngine.h b/libid/include/engine/PertEngine.h index 76db4664b..9c8b38f4c 100644 --- a/libid/include/engine/PertEngine.h +++ b/libid/include/engine/PertEngine.h @@ -14,9 +14,12 @@ class PertEngine public: void initialize_frame(const BFComplex ¢er_bf, const std::complex ¢er, double zoom_radius); int calculate_one_frame(); + int perturbation_per_pixel(int x, int y, double bailout); + int calculate_orbit(int x, int y, long iteration); + int get_number_references(); private: - int calculate_point(const Point &pt, double magnified_radius, int window_radius); + int calculate_reference(int x, int y); void reference_zoom_point(const BFComplex ¢er, int max_iteration); void reference_zoom_point(const std::complex ¢er, int max_iteration); void cleanup(); @@ -26,15 +29,15 @@ class PertEngine std::vector m_perturbation_tolerance_check; double m_delta_real{}; double m_delta_imag{}; - std::vector m_points_remaining; - std::vector m_glitch_points; - long m_glitch_point_count{}; - long m_remaining_point_count{}; BFComplex m_center_bf{}; std::complex m_center{}; double m_zoom_radius{}; - bool m_calculate_glitches{true}; double m_percent_glitch_tolerance{0.1}; // What percentage of the image is okay to be glitched. int m_reference_points{}; int m_saved_stack{}; + std::complex m_c{}; + BFComplex m_c_bf{}; + + std::complex m_delta_sub_0{}; + std::complex m_delta_sub_n{}; }; diff --git a/libid/include/engine/perturbation.h b/libid/include/engine/perturbation.h index 7faffd008..6e9532aef 100644 --- a/libid/include/engine/perturbation.h +++ b/libid/include/engine/perturbation.h @@ -3,3 +3,28 @@ #pragma once bool perturbation(); +extern bool mandel_perturbation_setup(); +extern bool is_perturbation(); +extern int perturbation_per_pixel(); +extern int perturbation_per_orbit(); +//extern bool is_pixel_finished(int x, int y); +//extern long get_glitch_point_count(); +//extern int calculate_reference(); +//extern void cleanup_perturbation(); +extern int get_number_references(); + +extern bool g_use_perturbation; // select perturbation code + +extern int g_row; +extern int g_col; +extern long g_color_iter; +extern double g_magnitude_limit; + +enum class PerturbationMode +{ + AUTO = 0, // the default + YES = 1, + NO = 2 +}; +extern PerturbationMode g_perturbation; +extern double g_perturbation_tolerance; diff --git a/libid/ui/cmdfiles.cpp b/libid/ui/cmdfiles.cpp index 769fc64d5..0186507c4 100644 --- a/libid/ui/cmdfiles.cpp +++ b/libid/ui/cmdfiles.cpp @@ -16,6 +16,7 @@ #include "engine/get_prec_big_float.h" #include "engine/id_data.h" #include "engine/log_map.h" +#include "engine/perturbation.h" #include "engine/soi.h" #include "engine/sticky_orbits.h" #include "fractals/check_orbit_name.h" @@ -3759,8 +3760,35 @@ static CmdArgFlags cmd_xy_shift(const Command &cmd) return CmdArgFlags::FRACTAL_PARAM | CmdArgFlags::PARAM_3D; } +// tolerance=? +static CmdArgFlags cmd_tolerance(const Command &cmd) +{ + g_perturbation_tolerance = cmd.float_vals[0]; + return CmdArgFlags::NONE; +} + +// perturbation=? +static CmdArgFlags cmd_perturbation(const Command &cmd) +{ + std::string p = cmd.value; + if (p == "auto") + { + g_perturbation = PerturbationMode::AUTO; + } + else if (p == "yes") + { + g_perturbation = PerturbationMode::YES; + } + else + { + g_perturbation = PerturbationMode::NO; + } + g_use_perturbation = is_perturbation(); + return CmdArgFlags::NONE; +} + // Keep this sorted by parameter name for binary search to work correctly. -static std::array s_commands{ +static std::array s_commands{ CommandHandler{"3d", cmd_3d}, // CommandHandler{"3dmode", cmd_3d_mode}, // CommandHandler{"ambient", cmd_ambient}, // @@ -3859,6 +3887,7 @@ static std::array s_commands{ CommandHandler{"passes", cmd_passes}, // CommandHandler{"periodicity", cmd_periodicity}, // CommandHandler{"perspective", cmd_perspective}, // + CommandHandler{"perturbation", cmd_perturbation}, // CommandHandler{"pixelzoom", cmd_pixel_zoom}, // CommandHandler{"plotstyle", cmd_deprecated}, // deprecated print parameters CommandHandler{"polyphony", cmd_polyphony}, // @@ -3900,6 +3929,7 @@ static std::array s_commands{ CommandHandler{"tempdir", cmd_temp_dir}, // CommandHandler{"textcolors", cmd_text_colors}, // CommandHandler{"title", cmd_deprecated}, // deprecated print parameters + CommandHandler{"tolerance", cmd_tolerance}, // perturbation "glitch" error tolerance CommandHandler{"tplus", cmd_tplus}, // CommandHandler{"translate", cmd_deprecated}, // deprecated print parameters CommandHandler{"transparent", cmd_transparent}, // diff --git a/libid/ui/get_toggles.cpp b/libid/ui/get_toggles.cpp index 2a4020842..2cc2dd426 100644 --- a/libid/ui/get_toggles.cpp +++ b/libid/ui/get_toggles.cpp @@ -5,7 +5,9 @@ #include "engine/calcfrac.h" #include "engine/id_data.h" #include "engine/log_map.h" +#include "engine/perturbation.h" #include "fractals/fractype.h" +#include "fractals/fractalp.h" #include "helpdefs.h" #include "io/save_timer.h" #include "ui/cmdfiles.h" @@ -40,24 +42,24 @@ int get_toggles() char prev_save_name[ID_FILE_MAX_DIR + 1]; FullScreenValues values[25]; int old_sound_flag; - const char *calc_modes[] = { - "1", "2", "3", "g", "g1", "g2", "g3", "g4", "g5", "g6", "b", "s", "t", "d", "o", "p"}; - const char *sound_modes[5] = {"off", "beep", "x", "y", "z"}; - const char *inside_modes[] = { - "numb", "maxiter", "zmag", "bof60", "bof61", "epsiloncross", "startrail", "period", "atan", "fmod"}; - const char *outside_modes[] = {"numb", "iter", "real", "imag", "mult", "summ", "atan", "fmod", "tdis"}; + char const *calc_modes[] = { + "1", "2", "3", "g", "g1", "g2", "g3", "g4", "g5", "g6", "b", "s", "t", "d", "o"}; + char const *sound_modes[5] = {"off", "beep", "x", "y", "z"}; + char const *inside_modes[] = { + "numb", "maxiter", "zmag", "bof60", "bof61", "epsiloncross", "startrail", "atan", "fmod"}; + char const *outside_modes[] = {"numb", "iter", "real", "imag", "mult", "summ", "atan", "fmod", "tdis"}; + char const *perturbation_modes[] = {"auto", "yes", "no"}; int k = -1; - choices[++k] = "Passes (1-3, g[es], b[ound], t[ess], d[iff], o[rbit], p[ert])"; + choices[++k] = "Passes (1-3, g[es], b[ound], t[ess], d[iff], o[rbit])"; values[k].type = 'l'; values[k].uval.ch.vlen = 3; values[k].uval.ch.list_len = std::size(calc_modes); values[k].uval.ch.list = calc_modes; - values[k].uval.ch.val = - (g_user_std_calc_mode == '1') ? 0 - : (g_user_std_calc_mode == '2') ? 1 - : (g_user_std_calc_mode == '3') ? 2 + values[k].uval.ch.val = (g_user_std_calc_mode == '1') ? 0 + : (g_user_std_calc_mode == '2') ? 1 + : (g_user_std_calc_mode == '3') ? 2 : (g_user_std_calc_mode == 'g' && g_stop_pass == 0) ? 3 : (g_user_std_calc_mode == 'g' && g_stop_pass == 1) ? 4 : (g_user_std_calc_mode == 'g' && g_stop_pass == 2) ? 5 @@ -65,12 +67,12 @@ int get_toggles() : (g_user_std_calc_mode == 'g' && g_stop_pass == 4) ? 7 : (g_user_std_calc_mode == 'g' && g_stop_pass == 5) ? 8 : (g_user_std_calc_mode == 'g' && g_stop_pass == 6) ? 9 - : (g_user_std_calc_mode == 'b') ? 10 - : (g_user_std_calc_mode == 's') ? 11 - : (g_user_std_calc_mode == 't') ? 12 - : (g_user_std_calc_mode == 'd') ? 13 - : (g_user_std_calc_mode == 'o') ? 14 - : /* "p"erturbation */ 15; + : (g_user_std_calc_mode == 'b') ? 10 + : (g_user_std_calc_mode == 's') ? 11 + : (g_user_std_calc_mode == 't') ? 12 + : (g_user_std_calc_mode == 'd') ? 13 + : /*(g_user_std_calc_mode == 'o') ?*/ 14; + char old_user_std_calc_mode = g_user_std_calc_mode; int old_stop_pass = g_stop_pass; choices[++k] = "Maximum Iterations (2 to 2,147,483,647)"; @@ -228,6 +230,29 @@ int get_toggles() values[k].uval.dval = old_close_proximity; const HelpLabels old_help_mode = g_help_mode; + choices[++k] = "Use Perturbation (yes, no, auto - internal choice)"; + values[k].type = 'l'; + values[k].uval.ch.vlen = 4; + values[k].uval.ch.list_len = std::size(perturbation_modes); + values[k].uval.ch.list = perturbation_modes; + PerturbationMode old_perturbation = g_perturbation; + if (g_perturbation == PerturbationMode::AUTO) + { + values[k].uval.ch.val = 0; + } + else if (g_perturbation == PerturbationMode::YES) + { + values[k].uval.ch.val = 1; + } + else + { + values[k].uval.ch.val = 2; + } + + choices[++k] = "Perturbation tolerance"; + values[k].type = 'd'; + double old_perturbation_tolerance = g_perturbation_tolerance; + values[k].uval.dval = old_perturbation_tolerance; g_help_mode = HelpLabels::HELP_X_OPTIONS; int i = full_screen_prompt( "Basic Options\n(not all combinations make sense)", k + 1, choices, values, 0, nullptr); @@ -414,5 +439,35 @@ int get_toggles() j++; } + int tmp = values[++k].uval.ch.val; + if (tmp >= 0) + { + switch (tmp) + { + case 0: + g_perturbation = PerturbationMode::AUTO; + break; + case 1: + g_perturbation = PerturbationMode::YES; + break; + case 2: + g_perturbation = PerturbationMode::NO; + break; + } + } + + g_use_perturbation = is_perturbation(); + + if (old_perturbation != g_perturbation) + { + j++; + } + + g_perturbation_tolerance = values[++k].uval.dval; + if (old_perturbation_tolerance != g_perturbation_tolerance) + { + j++; + } + return j; } diff --git a/libid/ui/make_batch_file.cpp b/libid/ui/make_batch_file.cpp index 67f46999c..9dd7613fc 100644 --- a/libid/ui/make_batch_file.cpp +++ b/libid/ui/make_batch_file.cpp @@ -11,6 +11,7 @@ #include "engine/get_prec_big_float.h" #include "engine/id_data.h" #include "engine/log_map.h" +#include "engine/perturbation.h" #include "engine/sticky_orbits.h" #include "engine/type_has_param.h" #include "fractals/fractalp.h" @@ -828,6 +829,16 @@ static void write_batch_params(const char *color_inf, bool colors_only, int max_ put_param(" passes=%c%c", g_user_std_calc_mode, (char)g_stop_pass + '0'); } + if (g_perturbation != PerturbationMode::AUTO) + { + put_param(" %s=%s", "perturbation", (g_perturbation == PerturbationMode::YES) ? "yes" : "no"); + } + + if (g_perturbation_tolerance != 1e-6) + { + put_param(" %s=%lg", "tolerance", g_perturbation_tolerance); + } + if (g_use_center_mag) { LDouble magnification; diff --git a/libid/ui/tab_display.cpp b/libid/ui/tab_display.cpp index a69a57242..de7d4202f 100644 --- a/libid/ui/tab_display.cpp +++ b/libid/ui/tab_display.cpp @@ -9,6 +9,7 @@ #include "engine/engine_timer.h" #include "engine/id_data.h" #include "engine/param_not_used.h" +#include "engine/perturbation.h" #include "engine/pixel_grid.h" #include "engine/soi.h" #include "engine/type_has_param.h" @@ -350,6 +351,15 @@ int tab_display() // display the status of the current image } start_row += 1; + if (g_use_perturbation) + { + int ref = get_number_references(); + std::sprintf(msg, (ref == 1) ? " (%d reference)" : " (%d references)", ref); + driver_put_string(start_row, 45, C_GENERAL_HI, "Perturbation"); + driver_put_string(-1, -1, C_GENERAL_HI, msg); + } + start_row += 1; + if (g_calc_status == CalcStatus::IN_PROGRESS || g_calc_status == CalcStatus::RESUMABLE) { if (bit_set(g_cur_fractal_specific->flags, FractalFlags::NO_RESUME))