OpenGLを使ってマンデルブロ

前にCUDAでマンデルブロを描いてビットマップに出力するのを作ったんだけど、毎回ビットマップ見るのめんどいしディスクIOが遅いからショボーンだった。
だから今日1日でGL使って描くように作り直してみた。まぁ手抜きでglut使ってるんだけどね。でもGLの仕様がﾌﾟｷﾞｬｰなのでかなり大変だった。
あとCUDA SDKのソースとか読んだりしてた。

CUDA3.0でグラフィック関係のAPIが変わったりしたので、それもついでに弄ってみた。
3.0だとcudaGraphicsGLRegisterImage()の第3引数にGL_TEXTURE_...とかGL_RENDERBUFFERを指定するといいよってなってるんだけど、GL_RENDERBUFFERが定義されてないとかエラー吐かれるorz。glext.hインクルードしてるんだけどなぁ。誰かわからない？

とりあえずソースを晒してみる。とりあえずやってみる為に書いたソースなので最適化とか全然してないです。ただ読みやすいとは思う。式を素直にコードに落としたからね。
コンパイル時に（まぁいつでもいいけど）USE_CUDA3を定義するとcuda 3.0の関数を呼ぶようになります。

このソースも大学のところに置いておこう。ここ
たまに更新とかしておきます。ちなみにライセンスは3または4条項BSDです。適当に遊んであげてください。

#include <GL/glew.h>
#include <GL/glut.h>

#include <cuda.h>
#include <cuda_gl_interop.h>


__global__ void mandelbrot( uchar4 *buffer, ulong2 size, float2 rrange, float2 irange, unsigned long n )
{
    unsigned long x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned long y = blockIdx.y * blockDim.y + threadIdx.y;
    if ( x >= size.x || y >= size.y )
    {
        return;
    }

    float2 C = make_float2( ( rrange.y - rrange.x ) / size.x * x + rrange.x,
                            ( irange.y - irange.x ) / size.y * y + irange.x );

    uchar4 buf = make_uchar4( 0, 0, 0, 0 );
    if ( C.x * C.x + C.y * C.y <= 4.0f )
    {
        float2 z = make_float2( 0.0f, 0.0f );
        unsigned long cnt = 0;
        while ( cnt < n )
        {
            ++cnt;
            z = make_float2( z.x * z.x - z.y * z.y + C.x, 2.0f * z.x * z.y + C.y );
            if ( z.x * z.x + z.y * z.y > 4.0f )
            {
                unsigned char c = ( 255 - 15 ) * __powf( ( float )cnt / n, 0.8f ) + 15;
                buf = make_uchar4( c, c, c, 0 );
                break;
            }
        }
    }

    buffer[ x + y * size.x ] = buf;
}

namespace
{

uchar4 *h_Src = NULL;

ulong2 size = make_ulong2( 800, 600 );
float2 rrange = make_float2( -2.0f, 2.0f ), irange = make_float2( -2.0f, 2.0f );

#ifdef USE_CUDA3
cudaGraphicsResource *crb = NULL;
#endif
GLuint tex = NULL, rb = NULL;

GLuint shader = NULL;

}

void cleanup( void )
{
    if ( h_Src )
    {
        free( h_Src );
        h_Src = NULL;
    }

    if ( tex )
    {
        glDeleteTextures( 1, &tex );
        tex = NULL;
    }
    if ( rb )
    {
#ifdef USE_CUDA3
        cudaGraphicsUnregisterResource( crb );
        crb = NULL;
#else
        cudaGLUnregisterBufferObject( rb );
#endif
        glDeleteBuffers( 1, &rb );
        rb = NULL;
    }
}

void InitGL( void )
{
    cleanup();

    h_Src = reinterpret_cast< uchar4 * >( malloc( sizeof( uchar4 ) * size.x * size.y ) );

    glEnable( GL_TEXTURE_2D );
    glGenTextures( 1, &tex );
    glBindTexture( GL_TEXTURE_2D, tex );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST );
    glTexImage2D( GL_TEXTURE_2D, 0, GL_RGBA8, size.x, size.y, 0, GL_RGBA, GL_UNSIGNED_BYTE, h_Src );

    glGenBuffers( 1, &rb );
    glBindBuffer( GL_PIXEL_UNPACK_BUFFER_ARB, rb );
    glBufferData( GL_PIXEL_UNPACK_BUFFER_ARB, sizeof( uchar4 ) * size.x * size.y, h_Src, GL_STREAM_COPY );
#ifdef USE_CUDA3
    cudaGraphicsGLRegisterBuffer( &crb, rb, cudaGraphicsMapFlagsNone );
#else
    cudaGLRegisterBufferObject( rb );
#endif

    const char * const code = "!!ARBfp1.0\nTEX result.color, fragment.texcoord, texture[ 0 ], 2D;\nEND";
    glGenProgramsARB( 1, &shader );
    glBindProgramARB( GL_FRAGMENT_PROGRAM_ARB, shader );
    glProgramStringARB( GL_FRAGMENT_PROGRAM_ARB, GL_PROGRAM_FORMAT_ASCII_ARB, static_cast< GLsizei >( strlen( code ) ), ( GLubyte * )code );

    GLint err;
    glGetIntegerv( GL_PROGRAM_ERROR_POSITION_ARB, &err );
    if ( err != -1 )
    {
        shader = NULL;
    }
}

void cmDraw( unsigned long convergence )
{
    uchar4 *d_ptr = NULL;

#ifdef USE_CUDA3
    size_t start;
    cudaGraphicsMapResources( 1, &crb, NULL );
    cudaGraphicsResourceGetMappedPointer( ( void ** )&d_ptr, &start, crb );
#else
    cudaGLMapBufferObject( ( void ** )&d_ptr, rb );
#endif

    dim3 dimBlock( 16, 4 );
    dim3 dimGrid( ( size.x + dimBlock.x - 1 ) / dimBlock.x, ( size.y + dimBlock.y - 1 ) / dimBlock.y );
    mandelbrot<<< dimGrid, dimBlock >>>( d_ptr, size, rrange, irange, convergence );
    cudaThreadSynchronize();

#ifdef USE_CUDA3
    cudaGraphicsUnmapResources( 1, &crb, NULL );
#else
    cudaGLUnmapBufferObject( rb );
#endif
}

void reshape( int w, int h )
{
    glViewport( 0, 0, w, h );
    glMatrixMode( GL_MODELVIEW );
    glLoadIdentity();

    glMatrixMode( GL_PROJECTION );
    glLoadIdentity();
    glOrtho( 0.0f, 1.0f, 0.0f, 1.0f, 0.0f, 1.0f );

    size = make_ulong2( w, h );
    float icenter = ( irange.x + irange.y ) * 0.5f, iwidth = ( rrange.y - rrange.x ) / w * h * 0.5f;
    irange = make_float2( icenter - iwidth, icenter + iwidth );
    InitGL();
}

void display( void )
{
    cmDraw( 1024 );

    glBindTexture( GL_TEXTURE_2D, tex );
    glTexSubImage2D( GL_TEXTURE_2D, 0, 0, 0, size.x, size.y, GL_RGBA, GL_UNSIGNED_BYTE, 0 );

    glBindProgramARB( GL_FRAGMENT_PROGRAM_ARB, shader );
    glEnable( GL_FRAGMENT_PROGRAM_ARB );
    glDisable( GL_DEPTH_TEST );

    glBegin( GL_QUADS );
    glTexCoord2f( 0.0f, 0.0f ); glVertex2f( 0.0f, 0.0f );
    glTexCoord2f( 1.0f, 0.0f ); glVertex2f( 1.0f, 0.0f );
    glTexCoord2f( 1.0f, 1.0f ); glVertex2f( 1.0f, 1.0f );
    glTexCoord2f( 0.0f, 1.0f ); glVertex2f( 0.0f, 1.0f );
    glEnd();

    glBindTexture( GL_TEXTURE_2D, 0 );
    glDisable( GL_FRAGMENT_PROGRAM_ARB );

    glutSwapBuffers();
}

int main( int argc, char **argv )
{
    glutInit( &argc, argv );
    glutInitDisplayMode( GLUT_RGBA | GLUT_DOUBLE );
    glutInitWindowSize( size.x, size.y );
    glutCreateWindow( *argv );

    glewInit();

    glutDisplayFunc( display );
    glutReshapeFunc( reshape );

    cudaGLSetGLDevice( 0 );

    atexit( cleanup );

    glutMainLoop();

    glDeleteProgramsARB( 1, &shader );
    cleanup();
    cudaThreadExit();
    return 0;
}