用显卡画Mandelbrot集

[ColdMooon]

学编程的大都写过画Mandelbrot集的程序。我很早以前写过一个。但是速度很慢。我速龙双核4400的CPU，画一幅1000x1000x256色的图，最多要16秒。前几天我想到，图中每个点都是独立的，也就是高度并行计算。GPU不就是干这个的？我知道DX11有compute shader,但是我不想装WIN7。于是我用pixel shader写了一个。
程序代码如下：
[code=C/C++]#include "stdafx.h"
#include "resource.h"
#include <d3dx9.h>

CString App = _T("Fractal");

struct Vertex
{
float x, y, z;
float zx, zy;
float cx, cy;
};

class Fractal : public CDialogImpl<Fractal>, public CAutoMessageFilter
{
CComPtr<IDirect3DDevice9> Device;
CComPtr<IDirect3DPixelShader9> Shader;
POINTFLOAT z, c;
float zoom;
int wx, wy;
BOOL bDrag;
CPoint HoldPos;
CRect SelRect;

public:
enum {IDD = ID_APP};

Fractal()
{
bDrag = false;

Create(0);
}

~Fractal()
{
DestroyWindow();
}

BOOL PreTranslateMessage(MSG* pMsg)
{
return IsDialogMessage(pMsg);
}

private:
int OnInitDialog(CWindow, LPARAM)
{
SetIcon(AtlLoadIcon(ID_APP));

IDirect3D9* d3d9 = Direct3DCreate9(D3D_SDK_VERSION);
D3DADAPTER_IDENTIFIER9 Adapter;
d3d9->GetAdapterIdentifier(D3DADAPTER_DEFAULT, 0, &Adapter);
::SetDlgItemTextA(*this, ID_GPU, Adapter.Description);

D3DCAPS9 Caps;
d3d9->GetDeviceCaps(D3DADAPTER_DEFAULT, D3DDEVTYPE_HAL, &Caps);
if (Caps.PixelShaderVersion < D3DPS_VERSION(3, 0))
{
::MessageBoxA(*this, "Pixel Shader v3.0 is not Supported.",
Adapter.Description, MB_ICONSTOP);
d3d9->Release();
exit(1);
}

D3DPRESENT_PARAMETERS d3dpp = {0};
d3dpp.Windowed = true;
d3dpp.SwapEffect = D3DSWAPEFFECT_COPY;
d3d9->CreateDevice(D3DADAPTER_DEFAULT, D3DDEVTYPE_HAL,
*this, D3DCREATE_HARDWARE_VERTEXPROCESSING, &d3dpp, &Device);
d3d9->Release();

CComboBox CtlFuns = GetDlgItem(ID_FUNCTION);
CFindFile ff;
if (ff.FindFile(_T("*.pso")))
{
do CtlFuns.AddString(ff.GetFileTitle());
while (ff.FindNextFile());
CtlFuns.SelectString(0, _T("Mand2"));
OnFunction(CBN_SELCHANGE, ID_FUNCTION, CtlFuns);
}
return 0;
}

void OnFunction(UINT, int, CComboBox Ctl)
{
Shader.Release();
CString FileName;
Ctl.GetLBText(Ctl.GetCurSel(), FileName);
HANDLE hFile = CreateFile(FileName + _T(".pso"),
GENERIC_READ, FILE_SHARE_READ, 0, OPEN_EXISTING, 0, 0);
if (hFile == INVALID_HANDLE_VALUE)
{
MessageBox(_T("Can not load shader from ") + FileName + _T(".pso\n\n") +
GetErrorMessage(), App, MB_ICONEXCLAMATION);
return;
}
DWORD size = GetFileSize(hFile, 0);
void* code = _alloca(size);
ReadFile(hFile, code, size, &size, 0);
CloseHandle(hFile);
Device->CreatePixelShader((DWORD*)code, &Shader);
Device->SetPixelShader(Shader);
zoom = 100;
Draw();
}

void OnSize(UINT Type, CPoint p)
{
wx = p.x - 142;
wy = p.y;
if (wx > 0)
{
D3DPRESENT_PARAMETERS d3dpp = {wx, wy};
d3dpp.Windowed = true;
d3dpp.hDeviceWindow = *this;
d3dpp.SwapEffect = D3DSWAPEFFECT_COPY;
Device->Reset(&d3dpp);
Device->SetPixelShader(Shader);
Device->SetFVF(D3DFVF_XYZ|D3DFVF_TEX2);
Draw();
}
::MoveWindow(GetDlgItem(ID_START ), p.x - 130, 10, 120, 30, true);
::MoveWindow(GetDlgItem(ID_FUNCTION), p.x - 130, 50, 120, 200, true);
::MoveWindow(GetDlgItem(ID_LOCATION), p.x - 136, 90, 132, 112, true);
::MoveWindow(GetDlgItem(ID_PX ), p.x - 130, 110, 120, 24, true);
::MoveWindow(GetDlgItem(ID_PY ), p.x - 130, 140, 120, 24, true);
::MoveWindow(GetDlgItem(ID_PZ ), p.x - 130, 170, 120, 24, true);
::MoveWindow(GetDlgItem(ID_JULIA ), p.x - 130, 210, 120, 24, true);
::MoveWindow(GetDlgItem(ID_PARAMS ), p.x - 136, 240, 132, 82, true);
::MoveWindow(GetDlgItem(ID_PJ ), p.x - 130, 260, 120, 24, true);
::MoveWindow(GetDlgItem(ID_PK ), p.x - 130, 290, 120, 24, true);
::MoveWindow(GetDlgItem(ID_GPU ), p.x - 130, 340, 120, 30, true);
::MoveWindow(GetDlgItem(ID_SAVE ), p.x - 130, 400, 120, 30, true);
::MoveWindow(GetDlgItem(ID_ABOUT ), p.x - 130, 430, 120, 30, true);
::MoveWindow(GetDlgItem(ID_CLOSE ), p.x - 130, 460, 120, 30, true);
}

void Draw()
{
Device->BeginScene();
if (IsDlgButtonChecked(ID_JULIA))
{
Vertex v[3] = {-2, -2, 0, z.x - wx / zoom, z.y - wy / zoom, c.x, c.y,
0, 4, 0, z.x, z.y + 2 * wy / zoom, c.x, c.y,
4, 0, 0, z.x + 2 * wx / zoom, z.y, c.x, c.y};
Device->DrawPrimitiveUP(D3DPT_TRIANGLELIST, 1, v, sizeof Vertex);
}
else
{
Vertex v[3] = {-2, -2, 0, c.x, c.y, z.x - wx / zoom, z.y - wy / zoom,
0, 4, 0, c.x, c.y, z.x, z.y + 2 * wy / zoom,
4, 0, 0, c.x, c.y, z.x + 2 * wx / zoom, z.y};
Device->DrawPrimitiveUP(D3DPT_TRIANGLELIST, 1, v, sizeof Vertex);
}
Device->EndScene();

TCHAR ss[20];
int Precision = _sntprintf(ss, 20, _T("%.f"), zoom) - 1;
SetDlgItemText(ID_PZ, ss);
for (int i = 0; i < 4; i++)
{
_sntprintf(ss, 20, _T("%.*f"), Precision, i[&z.x]);
SetDlgItemText(ID_PX + i, ss);
}
Invalidate(false);
}

void OnPaint()
{
CPaintDC(*this);
Device->Present(0, CRect(0, 0, wx, wy), 0, 0);
}

[ColdMooon]

void OnStart()
{
for (int i = 0; i < 5; i++)
{
TCHAR ss[20];
GetDlgItemText(ID_PX + i, ss, 20);
i[&z.x] = (float)_tcstod(ss, 0);
}
Draw();
}

void OnJulia()
{
POINTFLOAT t;
t = z;
z = c;
c = t;
Draw();
}

void OnLButtonDown(UINT Flags, CPoint p)
{
if (p.x < wx)
{
SetCapture();
bDrag = true;
HoldPos = p;
SelRect = CRect(p, p);
}
}

void OnMouseMove(UINT Flags, CPoint p)
{
if (bDrag)
{
CClientDC DC(*this);
DC.DrawFocusRect(SelRect);
SelRect = CRect(HoldPos, p);
SelRect.NormalizeRect();
SelRect &= CRect(0, 0, wx, wy);
DC.DrawFocusRect(SelRect);
}
}

void OnLButtonUp(UINT Flags, CPoint p)
{
if (bDrag)
{
ReleaseCapture();
bDrag = false;
z.x += (SelRect.left + SelRect.right - wx) / zoom / 2;
z.y -= (SelRect.top + SelRect.bottom - wy) / zoom / 2;
zoom *= min(min((float)wx / SelRect.Width(), (float)wy / SelRect.Height()), 10);
Draw();
}
}

void OnRButtonUp(UINT Flags, CPoint p)
{
if (bDrag)
{
DrawFocusRect(CClientDC(*this), SelRect);
ReleaseCapture();
bDrag = false;
}
else
{
zoom = max(zoom / 10, 100);
if (zoom == 100)
z.x = z.y = 0;
Draw();
}
}

void OnSave()
{
TCHAR FileName[MAX_PATH] = _T("");
OPENFILENAME ofn = {sizeof ofn};
ofn.hwndOwner = *this;
ofn.lpstrFilter = _T("Png(*.png)\0*.png\0Bitmap(*.bmp)\0*.bmp\0Jpeg(*.jpg)\0*.jpg\0");
ofn.lpstrFile = FileName;
ofn.nMaxFile = MAX_PATH;
ofn.Flags = OFN_OVERWRITEPROMPT|OFN_ENABLESIZING|OFN_HIDEREADONLY|OFN_NOCHANGEDIR;
ofn.lpstrDefExt = _T("png");

if (GetSaveFileName(&ofn))
{
CComPtr<IDirect3DSurface9> BackBuffer;
Device->GetBackBuffer(0, 0, D3DBACKBUFFER_TYPE_MONO, &BackBuffer);
if (FAILED(D3DXSaveSurfaceToFile(FileName,
D3DXIMAGE_FILEFORMAT(ofn.nFilterIndex + 2 & 3), BackBuffer, 0, 0)))
{
MessageBox(_T("Can not save image to ") + CString(FileName) + _T("\n\n") +
GetErrorMessage(), App, MB_ICONEXCLAMATION);
}
}
}

void OnAbout()
{
MSGBOXPARAMS mb = {sizeof mb};
mb.hwndOwner = *this;
mb.hInstance = GetModuleHandle(NULL);
mb.lpszText = (LPCTSTR)ID_ABOUT;
mb.lpszCaption = App;
mb.dwStyle = MB_USERICON;
mb.lpszIcon = (LPCTSTR)ID_APP;

MessageBoxIndirect(&mb);
}

void OnClose()
{
PostQuitMessage(0);
}

BEGIN_MSG_MAP(?)
MSG_WM_INITDIALOG( OnInitDialog)
MSG_WM_PAINT( OnPaint)
MSG_WM_SIZE( OnSize)
MSG_WM_LBUTTONDOWN( OnLButtonDown)
MSG_WM_MOUSEMOVE( OnMouseMove)
MSG_WM_LBUTTONUP( OnLButtonUp)
MSG_WM_RBUTTONUP( OnRButtonUp)
MSG_WM_CLOSE( OnClose)
COMMAND_CODE_HANDLER_EX(CBN_SELCHANGE, OnFunction)
COMMAND_ID(ID_START, OnStart)
COMMAND_ID(ID_JULIA, OnJulia)
COMMAND_ID(ID_SAVE, OnSave)
COMMAND_ID(ID_ABOUT, OnAbout)
COMMAND_ID(ID_CLOSE, OnClose)
END_MSG_MAP()
}
theFractal;
[/code]
程序用的是WTL框架。

[ColdMooon]

pixel shader代码如下：

#define MAX_ITER 255
#define EPSILON 1e-6

float2 cmul(float2 a, float2 b)
{
    float2 r;
    r.x = a.x * b.x - a.y * b.y;
    r.y = a.y * b.x + a.x * b.y;
    return r;
}

float2 cdiv(float2 a, float2 b)
{
    float2 r;
    r.x = a.x * b.x + a.y * b.y;
    r.y = a.y * b.x - a.x * b.y;
    return r / dot(b, b);
}

vector color(int n)
{
    vector c = 0;
    if (n > 0 && n < MAX_ITER)
    {
        c = (vector(14, 9, 4, 0) + n) * .4;
        c = abs(c % 6 - 3) - 1;
    }
    return c;
}


vector Mand2(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        z = cmul(z, z) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mand3(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        z = cmul(cmul(z, z), z) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mand4(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        z = cmul(z, z);
        z = cmul(z, z) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mand5(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 t = cmul(z, z);
        t = cmul(t, t);
        z = cmul(t, z) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mand6(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 t = cmul(z, z);
        t = cmul(t, z);
        z = cmul(t, t) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mandn1(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        z = cdiv(1, z) + c;
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Mandn2(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        z = cdiv(1, cmul(z, z) + c);
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector Newton3(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 s = cmul(z, z);
        float2 t = cmul(s, z) - c;
        if (dot(t, t) < EPSILON) break;
        z -= cdiv(t, s) / 3;
    }
    return color(n);
}

vector Newton4(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 s = cmul(z, z);
               s = cmul(s, z);
        float2 t = cmul(s, z) - c;
        if (dot(t, t) < EPSILON) break;
        z -= cdiv(t, s) / 4;
    }
    return color(n);
}

vector Newton5(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 s = cmul(z, z);
               s = cmul(s, s);
        float2 t = cmul(s, z) - c;
        if (dot(t, t) < EPSILON) break;
        z -= cdiv(t, s) / 5;
    }
    return color(n);
}

vector Newton6(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        float2 s = cmul(z, z);
               s = cmul(s, s);
               s = cmul(s, z);
        float2 t = cmul(s, z) - c;
        if (dot(t, t) < EPSILON) break;
        z -= cdiv(t, s) / 6;
    }
    return color(n);
}

vector barnsleyj1(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        if (z.x >= 0)
            z.x--;
        else
            z.x++;
        z = cmul(z, c);
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector barnsleyj2(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        if (z.x * c.y + c.x * z.y >= 0)
            z.x--;
        else
            z.x++;
        z = cmul(z, c);
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector barnsleyj3(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    for (int n = 0; n < MAX_ITER; n++)
    {
        if (z.x >= 0)
        {
            z = cmul(z, z);
            z.x--;
        }
        else
        {
            z = cmul(z, z) + z.x * c;
            z.x--;
        }
        if (dot(z, z) > 4) break;
    }
    return color(n);
}

vector henon(float2 z : TEXCOORD0, float2 c : TEXCOORD1) : COLOR
{
    float tx = z.x * z.x, ty;
    for (int n = 0; n < MAX_ITER; n++)
    {
        ty = z.x * c.y;
        z.x = 1 + z.y + tx * c.x;
        z.y = ty;
        if (tx + z.y * z.y > 4) break;
    }
    return color(n);
}

用DX SDK中的FXC编译，如fxc /Tps_3_0 /EMand2 /FoMand2.pso ps.txt
把编译出来的一堆pso文件放在程序所在目录运行程序即可。

有一个问题：pixel shader中for循环次数最多255次。就算MAX_ITER大于255也不行。这是怎么回事？期待DX高手指点。

[wmesci]

"速龙双核4400的CPU，画一幅1000x1000x256色的图，最多要16秒"
算法有问题
我用C#写过一个，AMD Athlon(tm) X2 245 2.91Ghz，分辨率800*600，32位色，速度不到50~60ms，C++应该更快

[fronz]

DX11的就这么慢么？
没概念

[zhuangzhuang1988]

这边有一个mandelbrot的例子 http://www.thecodeway.com/blog/?p=409 直接使用的DC但是他使用了TBB进行加速.挺好玩的.