每日算法 — 使用java实现曼德勃罗集：复数迭代与逃逸时间算法

概述

曼德勃罗集（Mandelbrot Set）是数学中最著名的分形之一，由复数平面上的点通过简单迭代公式生成。本文用Java实现曼德勃罗集的渲染算法，重点讲解复数迭代公式、逃逸时间算法，以及基于逃逸时间着色的可视化方案。该算法仅需数十行代码即可生成无限复杂的边界图案，是计算美学与数学的完美结合。

核心算法设计

1. 复数运算基础

曼德勃罗集的核心迭代公式为 z = z² + c，其中 z 和 c 均为复数。我们先实现一个轻量级的复数类。

public class Complex {
    public double real;
    public double imag;

    public Complex(double real, double imag) {
        this.real = real;
        this.imag = imag;
    }

    /** 复数加法：(a+bi) + (c+di) = (a+c) + (b+d)i */
    public Complex add(Complex other) {
        return new Complex(this.real + other.real, this.imag + other.imag);
    }

    /** 复数乘法：(a+bi)(c+di) = (ac-bd) + (ad+bc)i */
    public Complex multiply(Complex other) {
        double r = this.real * other.real - this.imag * other.imag;
        double i = this.real * other.imag + this.imag * other.real;
        return new Complex(r, i);
    }

    /** 复数模的平方：|z|² = a² + b² */
    public double magnitudeSquared() {
        return real * real + imag * imag;
    }
}

2. 逃逸时间算法

对于复平面上的每个点 c，从 z=0 开始迭代 z = z² + c。如果在限定迭代次数内 |z| 超过逃逸半径（通常取2），则称 c 逃逸，其迭代次数即为逃逸时间。若始终不逃逸，则 c 属于曼德勃罗集。

public class MandelbrotCalculator {
    private static final int MAX_ITERATIONS = 1000; // 最大迭代次数
    private static final double ESCAPE_RADIUS = 2.0;
    private static final double ESCAPE_RADIUS_SQUARED = 4.0;

    /**
     * 计算点c的逃逸时间
     * @return 逃逸迭代次数；若未逃逸则返回MAX_ITERATIONS
     */
    public int computeEscapeTime(double cReal, double cImag) {
        double zReal = 0.0;
        double zImag = 0.0;

        for (int iter = 0; iter < MAX_ITERATIONS; iter++) {
            // z = z² + c，展开为实部和虚部分开计算
            double zRealSq = zReal * zReal;
            double zImagSq = zImag * zImag;

            // 提前检测逃逸：|z|² > 4
            if (zRealSq + zImagSq > ESCAPE_RADIUS_SQUARED) {
                return iter;
            }

            // zImag = 2*zReal*zImag + cImag
            zImag = 2 * zReal * zImag + cImag;
            // zReal = zReal² - zImag² + cReal
            zReal = zRealSq - zImagSq + cReal;
        }
        return MAX_ITERATIONS; // 未逃逸，属于集合
    }

    /**
     * 优化版本：利用周期检测加速计算
     * 若z值在迭代中重复出现，则必定不会逃逸
     */
    public int computeEscapeTimeOptimized(double cReal, double cImag) {
        double zReal = 0.0, zImag = 0.0;
        double periodReal = 0.0, periodImag = 0.0;
        int period = 0, maxPeriod = 10;

        for (int iter = 0; iter < MAX_ITERATIONS; iter++) {
            double zRealSq = zReal * zReal;
            double zImagSq = zImag * zImag;

            if (zRealSq + zImagSq > ESCAPE_RADIUS_SQUARED) {
                return iter;
            }

            zImag = 2 * zReal * zImag + cImag;
            zReal = zRealSq - zImagSq + cReal;

            // 周期检测：若z回到period值附近，增加检测窗口
            if (zReal == periodReal && zImag == periodImag) {
                return MAX_ITERATIONS;
            }
            period++;
            if (period >= maxPeriod) {
                period = 0;
                maxPeriod *= 2;
                periodReal = zReal;
                periodImag = zImag;
            }
        }
        return MAX_ITERATIONS;
    }
}

3. 屏幕坐标到复平面的映射

将像素坐标 (px, py) 映射到复平面上的点 (re, im)，需要定义可视区域的中心、缩放比例和宽高比。

public class Viewport {
    public double centerX = -0.5;   // 复平面中心实部
    public double centerY = 0.0;    // 复平面中心虚部
    public double zoom = 1.0;       // 缩放因子（越大显示范围越小）
    public int width = 800;
    public int height = 600;

    /**
     * 将像素坐标转换为复平面坐标
     */
    public double[] pixelToComplex(int px, int py) {
        // 曼德勃罗集经典可视范围：实部[-2.5, 1.0]，虚部[-1.25, 1.25]
        double baseWidth = 3.5;
        double baseHeight = 2.5;

        double aspectRatio = (double) width / height;
        double viewWidth = baseWidth / zoom;
        double viewHeight = baseHeight / zoom;

        if (aspectRatio > 1.4) {
            viewWidth = viewHeight * aspectRatio;
        } else {
            viewHeight = viewWidth / aspectRatio;
        }

        double re = centerX - viewWidth / 2 + (px / (double) width) * viewWidth;
        double im = centerY + viewHeight / 2 - (py / (double) height) * viewHeight;
        return new double[]{re, im};
    }
}

4. 基于逃逸时间的着色算法

逃逸时间不仅用于判断集合成员，更是着色的核心依据。平滑着色公式 mu = iter + 1 - log(log(|z|)) / log(2) 可消除色带边界。

public class ColorMapper {
    /**
     * 平滑逃逸时间计算
     * 使用对数迭代计数消除色带
     */
    public double smoothEscapeTime(int iter, double zReal, double zImag) {
        if (iter >= 1000) return 1000;
        double magnitude = Math.sqrt(zReal * zReal + zImag * zImag);
        return iter + 1 - Math.log(Math.log(magnitude)) / Math.log(2);
    }

    /**
     * 将逃逸时间映射为RGB颜色（霓虹渐变）
     */
    public int mapToColor(double escapeTime, int maxIter) {
        if (escapeTime >= maxIter) {
            return 0x000000; // 集合内部为黑色
        }

        // 使用正弦波生成周期性颜色
        double t = escapeTime / 50.0;
        int r = (int)(128 + 127 * Math.sin(t * 0.8 + 0.0));
        int g = (int)(128 + 127 * Math.sin(t * 0.8 + 2.09));
        int b = (int)(128 + 127 * Math.sin(t * 0.8 + 4.18));
        return (r << 16) | (g << 8) | b;
    }

    /**
     * 火焰配色方案
     */
    public int mapToFireColor(double escapeTime, int maxIter) {
        if (escapeTime >= maxIter) return 0x000000;
        double ratio = escapeTime / maxIter;
        int r = (int)(255 * Math.min(1.0, ratio * 3));
        int g = (int)(255 * Math.min(1.0, Math.max(0, ratio * 3 - 1)));
        int b = (int)(255 * Math.max(0, ratio * 3 - 2));
        return (r << 16) | (g << 8) | b;
    }
}

5. 并行渲染加速

曼德勃罗集渲染是 embarrassingly parallel （完全可并行）的计算密集型任务，非常适合多线程加速。

import java.awt.image.BufferedImage;
import java.util.concurrent.*;

public class ParallelRenderer {
    private MandelbrotCalculator calculator = new MandelbrotCalculator();
    private ColorMapper colors = new ColorMapper();
    private Viewport viewport = new Viewport();
    private ExecutorService executor;

    public ParallelRenderer(int threads) {
        executor = Executors.newFixedThreadPool(threads);
    }

    /**
     * 并行渲染曼德勃罗集图像
     * 每行作为一个独立任务提交到线程池
     */
    public BufferedImage render() {
        BufferedImage image = new BufferedImage(viewport.width, viewport.height, BufferedImage.TYPE_INT_RGB);
        CountDownLatch latch = new CountDownLatch(viewport.height);

        for (int py = 0; py < viewport.height; py++) {
            final int row = py;
            executor.submit(() -> {
                for (int px = 0; px < viewport.width; px++) {
                    double[] c = viewport.pixelToComplex(px, row);
                    int iter = calculator.computeEscapeTime(c[0], c[1]);
                    int rgb = colors.mapToColor(iter, 1000);
                    image.setRGB(px, row, rgb);
                }
                latch.countDown();
            });
        }

        try {
            latch.await();
        } catch (InterruptedException e) {
            Thread.currentThread().interrupt();
        }
        return image;
    }

    /**
     * 分块渲染：将图像划分为块，减少线程同步开销
     */
    public BufferedImage renderTiled(int tileSize) {
        BufferedImage image = new BufferedImage(viewport.width, viewport.height, BufferedImage.TYPE_INT_RGB);
        int tilesX = (viewport.width + tileSize - 1) / tileSize;
        int tilesY = (viewport.height + tileSize - 1) / tileSize;
        CountDownLatch latch = new CountDownLatch(tilesX * tilesY);

        for (int ty = 0; ty < tilesY; ty++) {
            for (int tx = 0; tx < tilesX; tx++) {
                final int tileX = tx, tileY = ty;
                executor.submit(() -> {
                    int startX = tileX * tileSize;
                    int startY = tileY * tileSize;
                    int endX = Math.min(startX + tileSize, viewport.width);
                    int endY = Math.min(startY + tileSize, viewport.height);

                    for (int py = startY; py < endY; py++) {
                        for (int px = startX; px < endX; px++) {
                            double[] c = viewport.pixelToComplex(px, py);
                            int iter = calculator.computeEscapeTime(c[0], c[1]);
                            image.setRGB(px, py, colors.mapToColor(iter, 1000));
                        }
                    }
                    latch.countDown();
                });
            }
        }

        try {
            latch.await();
        } catch (InterruptedException e) {
            Thread.currentThread().interrupt();
        }
        return image;
    }
}

6. 完整主程序

import javax.imageio.ImageIO;
import java.awt.image.BufferedImage;
import java.io.File;

public class MandelbrotMain {
    public static void main(String[] args) throws Exception {
        ParallelRenderer renderer = new ParallelRenderer(Runtime.getRuntime().availableProcessors());
        Viewport vp = new Viewport();
        vp.width = 1920;
        vp.height = 1080;
        vp.centerX = -0.743643887037151;
        vp.centerY = 0.13182590420533;
        vp.zoom = 1000; // 放大1000倍，展现局部细节

        long start = System.currentTimeMillis();
        BufferedImage image = renderer.renderTiled(64);
        long elapsed = System.currentTimeMillis() - start;
        System.out.println("渲染完成，耗时: " + elapsed + "ms");

        ImageIO.write(image, "png", new File("mandelbrot.png"));
    }
}

复杂度分析

总结

曼德勃罗集的神奇之处在于：一个极简的复数迭代公式 z = z² + c，经过足够多次迭代后，能在复平面上勾勒出无限精细的边界。逃逸时间算法是该分形的计算基石，而平滑着色与并行渲染则将其从数学概念转化为视觉艺术。对于Java开发者而言，该项目也是练习多线程并行计算和图像处理的绝佳案例。