volfire.br

用于GPU通用计算的编程语言BrookGPU 0.4

/* This code is heavily based on Yury Uralsky's Volumetric Fire shader
on www.cgshaders.org.  The noise function has been changed slightly
and a few other changes have been made during the port to Brook.  The
largest of these changes being that we are not actually rendering, but
only doing the simulation, so we don't have a vertex program.
Therefore we need to force the user to provide some of the information
and we need to generate some of it ourselves in the code below. */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "main.h"
#include "volfire.h"

#define B 256

float g[B + B + 2][4];	// 4d gradient vectors
float p[B + B + 2];		// permutations

extern unsigned char fire_image[64 * 128 * 3 +1];

kernel float noise4(float4 vec, float p[], float4 g[])
{
    float4 b0, b1, r0, r1;
    float4 u, v, u0, v0, u1, v1;
    float i0, i1, i00, i10, i01, i11;
    float i000, i100, i010, i110, i001, i101, i011, i111;
    float4 temp1, temp2, temp3, temp4;
    float4 s;
    float2 a;

    //r0 = modf(vec, b0);
    b0 = floor(vec);
    r0 = vec - b0;
    b1 = b0 + 1;
    r1 = r0 - 1;
    b0 = fmod(b0, 256.);  //hardcoded texture size probably not the best idea...
    b1 = fmod(b1, 256.);

    // our streams (texture) are 1d, so we only need certain parts of
    // the coordinates space

    i0 = p[b0.x];
    i1 = p[b1.x];

    i00 = p[i0 + b0.y];
    i10 = p[i1 + b0.y];
    i01 = p[i0 + b1.y];
    i11 = p[i1 + b1.y];

    i000 = p[i00 + b0.z];
    i100 = p[i10 + b0.z];
    i010 = p[i01 + b0.z];
    i110 = p[i11 + b0.z];
    i001 = p[i00 + b1.z];
    i101 = p[i10 + b1.z];
    i011 = p[i01 + b1.z];
    i111 = p[i11 + b1.z];

    // sample 16 neighbors and calculate 16 dot-products

    // sample at x = 0
    temp1 = float4(0., 1., 0., 0.);
    temp1 = temp1 > 0. ? r1:r0;
    temp2 = float4(0., 0., 0., 1.);
    temp2 = temp2 > 0. ? r1:r0;
    temp3 = float4(0., 1., 0., 1.);
    temp3 = temp3 > 0. ? r1:r0;
    u0 = float4(dot(g[i000 + b0.w], r0),    dot(g[i010 + b0.w], temp1),  
		dot(g[i000 + b1.w], temp2), dot(g[i010 + b1.w], temp3));  

    temp1 = float4(0., 0., 1., 0.);
    temp1 = temp1 > 0. ? r1:r0;
    temp2 = float4(0., 1., 1., 0.);
    temp2 = temp2 > 0. ? r1:r0;
    temp3 = float4(0., 0., 1., 1.);
    temp3 = temp3 > 0. ? r1:r0;
    temp4 = float4(0., 1., 1., 1.);
    temp4 = temp4 > 0. ? r1:r0;
    v0 = float4( dot(g[i001 + b0.w], temp1), dot(g[i011 + b0.w], temp2), 
		 dot(g[i001 + b1.w], temp3), dot(g[i011 + b1.w], temp4));

    // sample at x = 1
    temp1 = float4(1., 0., 0., 0.);
    temp1 = temp1 > 0. ? r1:r0;
    temp2 = float4(1., 1., 0., 0.);
    temp2 = temp2 > 0. ? r1:r0;
    temp3 = float4(1., 0., 0., 1.);
    temp3 = temp3 > 0. ? r1:r0;
    temp4 = float4(1., 1., 0., 1.);
    temp4 = temp4 > 0. ? r1:r0;
    u1 = float4(dot(g[i100 + b0.w], temp1), dot(g[i110 + b0.w], temp2), 
		dot(g[i100 + b1.w], temp3), dot(g[i110 + b1.w], temp4)); 

    temp1 = float4(1., 0., 1., 0.);
    temp1 = temp1 > 0. ? r1:r0;
    temp2 = float4(1., 1., 1., 0.);
    temp2 = temp2 > 0. ? r1:r0;
    temp3 = float4(1., 0., 1., 1.);
    temp3 = temp3 > 0. ? r1:r0;     
    v1 = float4( dot(g[i101 + b0.w], temp1), dot(g[i111 + b0.w], temp2), 
		 dot(g[i101 + b1.w], temp3), dot(g[i111 + b1.w], r1));

    // now, perform quadro-cubic interpolation between these 16 noise samples

    s = r0*r0*(3-2*r0);

    // interpolation in 4d-space requires 15 1d interpolations,
    // but thanks to parallel nature of lerp instruction and some swizzling
    // we can do the same with only 5 lerp instructions

    // lerp between two pairs of 8 values along x axis

    u = lerp(u0, u1, s.x);
    v = lerp(v0, v1, s.x);

    // lerp between two pairs of 4 values (resulted from the previous lerp) along z axis

    u = lerp(u, v, s.z);

    // now, lerp between two pairs of 2 values along y axis

    a = lerp(u.xz, u.yw, s.y);

    // finally, lerp along w axis

    return abs(lerp(a.x, a.y, s.w));
}

// turbulence in 4 dimensions - add 5 octaves of gradient noise at different
// frequencies

kernel float turbulence4(float4 vec, float p[], float4 g[])
{
    float n, i;

    n = noise4(vec, p, g) +
	 0.5 * noise4(vec * 2.0, p, g) +
	 0.25 * noise4(vec * 4.0, p, g) +
	 0.125 * noise4(vec * 8.0, p, g) +
	 0.0625 * noise4(vec * 16.0, p, g);
    return n;
}


kernel void fireFP(float  p[],         // permutations texture
                   float4 g[],         // random gradient 4d vectors texture
                   float3 firetex[][], // color texture
                   float turb_scale,   // turbulence scale
                   float z_pos,        // what slice do you want?
                   float time,         // the time factor 
                   iter float2 it<>,   // hold position information
                   out float3 fire<>)
{

    //normally, we would be dealing with actual vertex data...
    //but we are going to hack it.

    float2 uv;
    float4 v;
    float4 frag_pos;
    float2 up, diag, right;
    up   = float2(0.0f,1.0f);
    diag = float2(1.0f,1.0f);
    right = float2(1.0f,0.0f);
	
    frag_pos = float4( it.x, it.y, z_pos, time);
    uv = float2( frag_pos.x - 0.5, z_pos - 0.5 );
    uv = uv+uv;
    uv = uv*uv;

    v = float4( frag_pos.x, frag_pos.y+(10000.0f - time), z_pos, time);

    uv = float2( sqrt(uv.x + uv.y), frag_pos.y + sqrt(frag_pos.y) * turb_scale * turbulence4(v * 1.5, p, g));   

    fire = firetex[uv];
}

void create_noise_lookups()
{ 
     int i = 0, j = 0;
     float x,y,z,w;
     float length;

     srand(rand());     
     while (i < B)
     {
	  x = 2 * rand()/(float)RAND_MAX - 1.;
	  y = 2 * rand()/(float)RAND_MAX - 1.;
	  z = 2 * rand()/(float)RAND_MAX - 1.;
	  w = 2 * rand()/(float)RAND_MAX - 1.;

	  length = sqrtf( x*x+y*y+z*z+w*w );
	  g[i][0] = x/length;
	  g[i][1] = y/length;
	  g[i][2] = z/length;
	  g[i][3] = w/length;
	  ++i;
     }
     
     // Generate permutations     
     for (i = 0; i<B; i++) p[i] = (float)i;
     
     while (--i) {
	  float k = p[i];
	  int j;
	  p[i] = p[j = rand() % B];
	  p[j] = k;
     }
     
     for (i = 0; i<B + 2; i++) {
	  p[B + i] = p[i];
	  for (j = 0; j<4; j++) g[B + i][j] = g[i][j];
     }
}

int volfire_main(int argc, char* argv[]) {
    unsigned char* input = NULL;
    unsigned char* output = NULL;
    float* outputf = NULL;
    float* inputf = NULL;
    int i, j, m;
    int size;
    float turbscale = 1.0;
    float zpos = 0.5;
    float time = 1000.;
    char filename[255];
    char base_file[255];
    char zeros[8];
    float time_in_ms = 0.0f;
    int iters;
    int timer = 0;
    FILE* outfile;

    if(argc < 6){
	fprintf(stderr, "Usage: %s <size> <turbscale> <zpos> <outfile> <iterations>\n", argv[0]);
	exit(1);
    }

    if(argc == 7 && !strcmp("timer", argv[6]))
	 timer = 1;
 
    size = atoi(argv[1]);
    turbscale = (float)atof(argv[2]);
    zpos = (float)atof(argv[3]);
    strncpy(base_file, argv[4], 254);
    iters = atoi(argv[5]);

    create_noise_lookups();
    {	
	 float* fire_image_float;
	 float3 fire<size,size>;
	 float p_stream<256>;
	 float4 g_stream<256>;
	 float3 firetex<128,64>;
	 
	 fire_image_float = (float*)malloc(3*128*64*sizeof(float));
	 
	 for(i=0; i<3*128*64; i++){
	      fire_image_float[i] = fire_image[i];
	 }
	 
	 output  = (unsigned char*)malloc(size*size*3*sizeof(unsigned char));
	 outputf = (float*)malloc(size*size*3*sizeof(float));
	  //Start timing
	 start = GetTime();

	 streamRead(p_stream, p);
	 streamRead(g_stream, g);
	 streamRead(firetex, fire_image_float);
	 
	 for(m=0; m<iters; m++){
	      
	      iter float2 it<size,size> = iter( float2(-32.0f, 0.0f), float2(32.0f, 128.0f));
	      fireFP(p_stream, g_stream, firetex, turbscale, zpos, m*2., it, fire );
	      
	      streamWrite(fire, outputf);
	      
	      if(!timer){
		   // convert output back to bytes for (ppm) output
		   for(i=0; i<size; i++){
			for(j=0; j<size; j++){
			     output[i*size*3+j*3+0] = (unsigned char)outputf[i*size*3+j*3+0];
			     output[i*size*3+j*3+1] = (unsigned char)outputf[i*size*3+j*3+1];
			     output[i*size*3+j*3+2] = (unsigned char)outputf[i*size*3+j*3+2];
			}
		   }
		   
		   if(m<10)
			sprintf(zeros, "000");
		   else if(m<100)
			sprintf(zeros, "00");
		   else if(m<1000)
			sprintf(zeros, "0");
		   else
			sprintf(zeros, "\0");
		   
		   sprintf(filename, "%s%s%d.ppm", base_file, zeros, m);
		   outfile = fopen(filename, "wb");
		   fprintf(outfile, "P6\n%d\n%d\n255\n", size, size);
		   for(i=0; i<size; i++){
			for(j=0; j<size; j++){
			     fwrite(&output[i*size*3+j*3], sizeof(unsigned char), 3, outfile); 
			}
		   }
		   fclose(outfile);
	      }
	 }
	 if(time)
		streamWrite(fire, outputf);
	 stop = GetTime();       
	 time_in_ms = (float)(stop-start)/1000.;
	 fprintf(stderr, "That took %.2fms total\n", time_in_ms);
    }
    return 0;
}