I was supposed to write a png writing function for our students to output grid data from a diffusion simulation. This is the example code to go along with it... (it's a raymarcher with antialiasing and reflections).

pngraymarcher.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

#include <png.h>
#include <stdio.h>
#include <stdlib.h>

int write_png(const char *file_name, int width, int height, unsigned char *data) {
    int i;
    png_structp png_ptr;
    png_infop info_ptr;
    png_bytep *row_pointers = NULL;
    FILE *fp = fopen(file_name, "wb");
    if(!fp) {
       return 0;
    }

    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if(!png_ptr) {
        fclose(fp);
        return 0;
    }

    info_ptr = png_create_info_struct(png_ptr);
    if(!info_ptr) {
        png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
        fclose(fp);
        return 0;
    }

    if(setjmp(png_jmpbuf(png_ptr))) {
        png_destroy_write_struct(&png_ptr, &info_ptr);
        fclose(fp);
        free(row_pointers);
        return 0;
    }

    png_init_io(png_ptr, fp);
    png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB_ALPHA, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);

    row_pointers = (png_bytep*)malloc(height*sizeof(png_bytep*));
    for(i = 0;i<height;++i) {
        row_pointers[i] = data+4*width*i;
    }
    png_set_rows(png_ptr, info_ptr, row_pointers);
    png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
    free(row_pointers);
    return 1;
}

void normalize(double &x, double &y, double &z) {
    double invsqrt = 1.0/std::sqrt(x*x+y*y+z*z);
    x *= invsqrt; y *= invsqrt; z *= invsqrt;
}

double clamp(double x, double a, double b) {
    return x<a?a:x>b?b:x;
}


double sphere(double x, double y, double z, double cx, double cy, double cz, double r) {
    x-=cx; y-=cy; z-=cz;
    return std::sqrt(x*x+y*y+z*z) - r;
}

double box(double x, double y, double z, double x0, double x1, double y0, double y1, double z0, double z1, double r) {
    x-=clamp(x, x0, x1); y-=clamp(y, y0, y1); z-=clamp(z, z0, z1);
    return std::sqrt(x*x+y*y+z*z) - r;
}

double diffuse(double x, double y, double z) {
    return std::min(
            std::min(
                sphere(x,y,z, 3,0,3, 0.7f),
                sphere(x,y,z, -1.5f,0,0, 1.0)
            ),
            box(x,y,z, -4,-3,-2,4,2,3, 0.2)
        );
}

double mirror(double x, double y, double z) {
    return std::min(
                y+2,
                sphere(x,y,z, 0,3,3, 1.5f)
            );
}

double all(double x, double y, double z) {
    return std::min(
            diffuse(x, y, z),
            mirror(x,y,z)
        );
}

void grad(double x, double y, double z, double &dx, double &dy, double &dz) {
    double eps = 1.e-5f*(std::abs(x)+std::abs(y)+std::abs(z));
    double inv2eps = 0.5f/eps;
    dx = (all(x+eps,y,z)-all(x-eps,y,z))*inv2eps;
    dy = (all(x,y+eps,z)-all(x,y-eps,z))*inv2eps;
    dz = (all(x,y,z+eps)-all(x,y,z-eps))*inv2eps;
}

void reflect(double &dx, double &dy, double &dz, double nx, double ny, double nz) {
    double dot = dx*nx+dy*ny+dz*nz;
    dx -= 2*dot*nx; dy -= 2*dot*ny; dz -= 2*dot*nz;
}

void sky(double dx, double dy, double dz, double &r, double &g, double &b) {
    (void)dx; (void)dy; (void)dz;
    //~ r += 0x87/255.0; g += 0xCE/255.0; b += 0xEB/255.0;
    r += 0; g += 0; b += 0;
}

double lights[][6] = {
    {300,100,-300, 100000,100000,100000},
    {10,10,0, 100,30,5}
};
int lightcount = 2;

double trace(double &x, double &y, double &z, double dx, double dy, double dz, double max_dist = 1.e6) {
    const int max_depth = 16*1024;
    const double precision = 1.e-10f;

    double d_tot = 0;
    double d = all(x,y,z);
    int i = 0;
    for(;i<max_depth && d>precision && d<max_dist;++i) {
        d_tot += d;
        x += d*dx; y += d*dy; z += d*dz;
        d = all(x,y,z);
    }
    return std::min(d_tot, max_dist);
}

void trace(double &x, double &y, double &z, double dx, double dy, double dz, double &r, double &g, double &b, int recursion = 10) {
    if(recursion == 0) return;

    const int max_depth = 16*1024;
    const double max_dist = 1.e6;
    const double precision = 1.e-10f;

    double d = all(x,y,z);
    int i = 0;
    for(;i<max_depth && d>precision && d<max_dist;++i) {
        x += d*dx; y += d*dy; z += d*dz;
        d = all(x,y,z);
    }

    if(i == max_depth || d>=max_dist) {
        sky(dx, dy, dz, r, g, b);
    } else {
        double nx, ny, nz;
        grad(x, y, z, nx, ny, nz);
        normalize(nx, ny, nz);

        for(int i = 0;i<lightcount;++i) {
            double lx = lights[i][0]-x;
            double ly = lights[i][1]-y;
            double lz = lights[i][2]-z;
            double dist2 = lx*lx+ly*ly+lz*lz;
            double dist = std::sqrt(dist2);
            double attenuation = 1.0/dist2;
            normalize(lx, ly, lz);
            double dot = lx*nx+ly*ny+lz*nz;
            double bright = 0;
            bright += 0.1*std::max(0.0, dot); // diffuse

            double dx2 = dx, dy2 = dy, dz2 = dz;
            reflect(dx2, dy2, dz2, nx, ny, nz);
            double dot2 = dx2*lx+dy2*ly+dz2*lz;
            bright += 0.9*std::pow(std::max(0.0, dot2), 64); // specular

            double x1 = x+10.0*precision*lx, y1 = y+10.0*precision*ly, z1 = z+10.0*precision*lz;
            double free_dist = trace(x1, y1, z1, lx, ly, lz, dist);
            bright *= dist-free_dist>1-2*precision?0:1;
            bright *= attenuation;

            r += bright*lights[i][3];
            g += bright*lights[i][4];
            b += bright*lights[i][5];
        }
        if(mirror(x,y,z) < diffuse(x,y,z)) {
            double dx2 = dx, dy2 = dy, dz2 = dz;
            reflect(dx2, dy2, dz2, nx, ny, nz);
            x+=10.0*precision*dx2; y+=10.0*precision*dy2; z+=10.0*precision*dz2;
            trace(x,y,z, dx2,dy2,dz2, r,g,b, recursion-1);
        }
    }
}

double samples[][2] = {
    {0,0},
    {0.33,0.33},{-0.33,0.33},{-0.33,-0.33},{0.33,-0.33},
    {0.33,0},{-0.33,0},{0,0.33},{0,-0.33},
};

double gamma(double x) {
    return std::pow(x, 1.0/2.2);
}

int main() {
    int width = 640, height = 480;
    int samplecount = 9;

    double aspect = double(width)/height;
    std::vector<unsigned char> data(4*width*height);
    #pragma omp parallel for
    for(int x = 0;x<width;++x) {
        for(int y = 0;y<height;++y) {
            double r=0.0, g=0.0, b=0.0;
            for(int s = 0;s<samplecount;++s) {
                double fx = 0, fy = 0, fz = -5;
                double dx = (2.0*(x+samples[s][0])/width-1.0)*aspect;
                double dy = 1.0-2.0*(y+samples[s][1])/height;
                double dz = 1.0;
                normalize(dx, dy, dz);
                trace(fx, fy, fz, dx, dy, dz, r, g, b);
            }
            r/=samplecount;g/=samplecount;b/=samplecount;
            int index = 4*(y*width+x);
            data[index+0] = clamp(gamma(r), 0.0, 1.0)*0xFFu; // red
            data[index+1] = clamp(gamma(g), 0.0, 1.0)*0xFFu; // green
            data[index+2] = clamp(gamma(b), 0.0, 1.0)*0xFFu; // blue
            data[index+3] = 0xFFu; // alpha (opacity)
        }
    }

    write_png("test.png", width, height, &data[0]);

    return 0;
}