基于区间算法的像素函数绘图方法（附matlab代码）（仿GrafEq）

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

0 前言

这个文章写作目的是很久之前看到了Matrix67的一篇博客，介绍一款软件叫GrafEq。
http://www.matrix67.com/blog/archives/4447#more-4447
和常规的数学软件不同，它不是给出一些离散点，然后进行连线。它是基于像素级别的绘制定义域内的图像。因此，它可以做到只要能写出函数表达式，就一定能绘制出图像，无论是隐函数还是什么其它复杂函数。

现在它的官网依然存在，可以体验一下。
http://www.peda.com/grafeq/
下面是官网中一个示例函数，还有很多神奇的图像也是根据函数绘制的，这里就不再介绍了。

当然，最近突然想起这个茬了，就大概看了看基本原理，自己仿照着算法做了一款基于matlab的图像绘制程序。

参考文献：
Reliable Two-Dimensional Graphing Methods for Mathematical Formulae with Two Free Variables (作者：Jeff Tupper)

1 算法原理

1.1 判断该像素是否是解

首先把每一个像素都考虑其长和宽，以一个区间的形式来描述。
这种算法是基于区间分析，当区间内存在解的时候，或者不确定区间内是否存在解的时候，将像素所代表的区间标为解。当区间内不存在解的时候，将像素所代表的区间标为背景。

其中像素所代表的区间，可以举例说明：
比如像素在x轴上的坐标为0,1,2,3,4…
则每一个像素对应的区间为[-0.5,0.5]，[0.5,1.5]，[1.5,2.5]…

对于输入的函数，我们最终可以把函数表示为下面的形式：
f(x,y)>g(x,y)或f(x,y)<g(x,y)或f(x,y)=g(x,y)
比如y=x^ 2-1，其中f(x,y)=y，g(x,y)=x^2-1

根据区间算法，可以得出2个逻辑值T和F。
以小于<为例，[a,b]<[A,B]（这里默认a<b，A<B）
如果b<A，也就是区间ab全部小于区间AB，则返回[T,T]；
如果B<a，也就是区间ab全部大于区间AB，则返回[F,F]；
如果两个区间存在交集，则返回[T,F]。

之后，所有返回[F,F]的标记为确定区间不存在解。所有返回[T,T]的标记为确定区间存在解。所有返回[T,F]的标记为不确定是否存在解。

1.2 算法示例

以y<x-2的函数为例，总的区间为x[-4,4]，y[-4,4]

如果像素是1×1像素，则这个像素点x区间为[-4,4]，y区间为[-4,4]
计算式左侧为[-4,4]，计算式右侧为[-6,2]
进行判断[-4,4]<[-6,2]，区间存在交集，返回[T,F]。
说明有可能存在解。

如果像素是2×2像素，则遍历循环。
1号像素x区间为[-4,0]，y区间为[-4,0]
进行判断[-4,0]<[-6,-2]，区间存在交集，返回[T,F]。
2号像素x区间[0,4]，y区间为[-4,0]
进行判断[-4,0]<[-2,2]，区间存在交集，返回[T,F]。
3号像素x区间[-4,0]，y区间为[0,4]
进行判断[0,4]<[-6,-2]，区间完全大于另一个区间，返回[F,F]。
4号像素x区间[0,4]，y区间为[0,4]
进行判断[0,4]<[-2,2]，区间存在交集，返回[T,F]。

因此，我们可以绘制出2*2像素的图像：

同理，如果取更多的像素，则可以看到更精细的图像：

可以看到，当把像素取到1024×1024时，则看不到明显的像素值了。

1.3 区间算法

前面的是以y<x-2为例，算法较为简单，如果碰上其它计算，则会复杂一些。比如y=xx，xy=x+sin(x)之类的。

比如f=x*y，对于区间运算，就是找到f可能的所有可能范围，以x[-1,2]，y[1,3]为例。x*y所能取到的最小值为-1*3=-3，x*y所能取到的最大值为2*3=6，因此[-1,2] *[1,3]=[-3,6]。

比如f=sin(x)，对于区间运算，同样是找到f可能的所有可能范围。以x[0,2]为例，最小值是0，最大值为1。因此，sin([0,2])=[0,1]。

其它区间算法就不再一一列举了。好像matlab内没有相关的计算，所以需要自己进行编程计算。

2 计算示例

以下面的函数
f(x,y) =cos(cos(min(sinx+y,x+siny))) − cos(sin(max(siny+x,y+sinx))) > 0
为例：

绘图区间x为[-10,10]，y为[-10,10]，总共划分为1024份。

绘图matlab代码如下：

clear
clc
close all
%定义求解区间
X=[-10,10];
Y=[-10,10];


%求解域的长与宽
W=X(2)-X(1);
H=Y(2)-Y(1);


N_Split=10;%分割10次，采用2分法
D_Split=2;%采用2分法

dX=W/D_Split^N_Split;
dY=H/D_Split^N_Split;
x=X(1):dX:X(2);
y=Y(1):dY:Y(2);
%定义本次切割的所有网格格点
Interval_k_Ind=Mesh2IntervalInd(x,y);
%0是未定义，1是确定存在解，2是确定不存在解
Solve=zeros(length(y)-1,length(x)-1);
%遍历循环
for m=1:size(Interval_k_Ind,1)
    x_I_New=Interval_k_Ind(m,1);
    y_I_New=Interval_k_Ind(m,2);
    X_I_k=x(x_I_New:x_I_New+1);
    Y_I_k=y(y_I_New:y_I_New+1);
    I_TF=I_Formula(X_I_k,Y_I_k);
    %然后是对结果的判断
    if xor(I_TF(1),I_TF(2))
        Solve(y_I_New,x_I_New)=0;
    elseif I_TF(1)==true && I_TF(2)==true
        Solve(y_I_New,x_I_New)=1;
    elseif I_TF(1)==false && I_TF(2)==false
        Solve(y_I_New,x_I_New)=2;
    end
    
end
Solve_Old=Solve;

    
%将所有等于0的替换为1
Solve(Solve==0)=1;
Solve=2-Solve;
%绘图
[X_Draw,Y_Draw]=meshgrid(0.5*(x(1:end-1)+x(2:end)),0.5*(y(1:end-1)+y(2:end)));
figure()
imagesc(0.5*(x(1:end-1)+x(2:end)),0.5*(y(1:end-1)+y(2:end)),Solve)
set(gca,'YDir','normal')
colormap([[1,1,1];[0,0,0]])%白色的是背景，解为黑色
shading flat
axis equal
xlim(X);ylim(Y);




function I_TF=I_Formula(I_x,I_y)
% %示例等式1：y < x-2
% I_TF=Less(I_y,I_x-2);

% %示例等式2：y > x*x-1
% I_TF=Greater(I_y,I_Multi(I_x,I_x)-1);

%示例等式3：f(x,y) =cos(cos(min(sinx+y,x+siny))) − cos(sin(max(siny+x,y+sinx))) > 0
F1=I_Cos(I_Cos(I_Min( I_PLUS(I_Sin(I_x),I_y) , I_PLUS(I_Sin(I_y),I_x))));
F2=I_Cos(I_Sin(I_Max( I_PLUS(I_Sin(I_y),I_x) , I_PLUS(I_Sin(I_x),I_y))));
I_TF=Greater(F1,F2);

end

function I_Save=Mesh2Interval(x,y)
[x2,y2]=meshgrid(x,y);
I_Save=zeros((length(x)-1)*(length(y)-1),4);
k=1;
for kx=1:length(x)-1
    for ky=1:length(y)-1
        %储存为x1,x2,y1,y2的形式
        I_Save(k,:)