C:/Users/thebigMuh/Documents/Visual Studio 2005/Projects/Volumizer/Volumizer/RayEngine.cpp

Go to the documentation of this file.

#include "Volumizer.h"

unsigned short *rawdat, rawmax;         
float *voldat;                                          
unsigned int volx, voly, volz;          

extern HWND wnd_main;
extern camera cam;
extern HDC rayHdc;
extern HGLRC rayHglrc;
extern vect ldir;

INT_PTR CALLBACK ProgressProc(HWND hwndDlg, UINT uMsg, WPARAM wParam, LPARAM lParam)
{
        return FALSE;
}

void ShovelData()
{
        if (rawdat)
        {

                HWND wnd_progress = CreateDialog(hInst, 
                        MAKEINTRESOURCE(IDD_PROGRESS), 
                        wnd_main, 
                        ProgressProc);

                ShowWindow(wnd_progress, SW_SHOW);

                HWND progressbar = GetDlgItem(wnd_progress, IDC_PROGRESS);

                SendMessage(progressbar, PBM_SETRANGE, 0, MAKELPARAM(0, volx));
        
                unsigned int i = 0;

                for (unsigned int x = 0; x < volx; x++)
                {
                        for (unsigned int y = 0; y < voly; y++)
                        {
                                for (unsigned int z = 0; z < volz; z++)
                                {
                                        voldat[i] = OogleVal((float)rawdat[i] / (float)rawmax);
                                        i++;
                                }
                        }
                        SendMessage(progressbar, PBM_SETPOS, x, 0);
                }

                BuildOrthoSlice(wnd_front);
                BuildOrthoSlice(wnd_left);
                BuildOrthoSlice(wnd_top);

                DestroyWindow(wnd_progress);

                DrawOrtho(wnd_front);
                DrawOrtho(wnd_left);
                DrawOrtho(wnd_top);

                CalcHistogram();

                DrawHistogram(wnd_histogram);

                RayRender();
        }
}

void LoadDataFile(wchar_t *fileName)
{
        int infile;
        _wsopen_s(&infile, fileName, _O_BINARY | _O_SEQUENTIAL, _SH_DENYNO, 0);

        _read(infile, &volx, 2);
        _read(infile, &voly, 2);
        _read(infile, &volz, 2);

        delete [] voldat;

        voldat = new float[volx*voly*volz];

        delete [] rawdat;

        rawdat = new unsigned short[volx*voly*volz];

        _read(infile, rawdat, volx*voly*volz*2);

        rawmax = 0;

        for (unsigned int i = 0; i < volx * voly * volz; i++)
        {
                if (rawdat[i] > rawmax)
                        rawmax = rawdat[i];
        }

        _close(infile);

        xslice = volx / 2;
        yslice = voly / 2;
        zslice = volz / 2;

        InitCam();
        InitLight();

        ShovelData();
}

float bias(float _val, float _bias)
{
        return _val / ((1.0f / _bias - 2.0f) * (1.0f - _val) + 1.0f);
}

float gain(float _val, float _gain)
{
        if (_val < 0.5f)
        {
                return _val / ((1.0f / _gain - 2.0f) * (1.0f - 2.0f * _val) + 1.0f);
        }
        else
        {
                return ((1.0f / _gain - 2.0f) * (1.0f - 2.0f * _val) - _val) / ((1.0f / _gain - 2.0f) * (1.0f - 2.0f * _val) - 1.0f);
        }
}

static __inline float OogleVal(float inp)
{
        float tmp = bias(inp, brightness);

        tmp = gain(tmp, contrast);

        return (gamma == 0.0f)?0.0f:pow(tmp, 1.0f/gamma);
}

static __inline float GetSample(int x, int y, int z)
{
        return voldat[x + volx * y + voly * volx * z];
}

static __inline float GetInterpolatedSample(vect pos)
{
        float s1, s2, s3, s4;

        float intx, inty, intz;

        intx = fmod(pos.x, 1.0f);
        inty = fmod(pos.y, 1.0f);
        intz = fmod(pos.z, 1.0f);

        unsigned int basex, basey, basez;

        basex = (unsigned int)pos.x;
        basey = (unsigned int)pos.y;
        basez = (unsigned int)pos.z;

        if (basex > volx - 2)
                basex = volx - 2;

        if (basey > voly - 2)
                basey = voly - 2;

        if (basez > volz - 2)
                basez = volz - 2;

        s1 = GetSample(basex, basey, basez) * (1.0f - intz) + GetSample(basex, basey, basez + 1) * intz;
        s2 = GetSample(basex, basey+1, basez) * (1.0f - intz) + GetSample(basex, basey+1, basez + 1) * intz;
        s3 = GetSample(basex+1, basey, basez) * (1.0f - intz) + GetSample(basex+1, basey, basez + 1) * intz;
        s4 = GetSample(basex+1, basey+1, basez) * (1.0f - intz) + GetSample(basex+1, basey+1, basez + 1) * intz;

        s1 = s1 * (1.0f - inty) + s2 * inty;
        s2 = s3 * (1.0f - inty) + s4 * inty;

        s1 = s1 * (1.0f - intx) + s2 * intx;

        return s1;
}

static __inline
int GetBoxIntersection(vect &vnear, 
                                           vect &vfar,  
                                           vect &org,   
                                           vect &dir    
                                           )
{
        vect sect;
        float dist_near, dist_far, dist;
        int nsect = 0;
        int vx = volx - 1;
        int vy = voly - 1;
        int vz = volz - 1;

        dist_near = dist_far = 100000000.0f;

        if (fabs(dir.x) > 0.001f)
        {
                // 0-plane
                dist = -org.x / dir.x;

                sect.x = 0.0f;
                sect.y = org.y + dir.y * dist;
                sect.z = org.z + dir.z * dist;

                if (sect.y >= 0.0f && sect.y <= vy &&
                        sect.z >= 0.0f && sect.z <= vz &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }

                // volx-plane
                dist = (vx - org.x) / dir.x;

                sect.x = (float)vx;
                sect.y = org.y + dir.y * dist;
                sect.z = org.z + dir.z * dist;

                if (sect.y >= 0.0f && sect.y <= vy &&
                        sect.z >= 0.0f && sect.z <= vz &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }
        }

        // intersect with y-planes
        if (fabs(dir.y) > 0.001f)
        {
                // 0-plane
                dist = -org.y / dir.y;

                sect.x = org.x + dir.x * dist;
                sect.y = 0.0f;
                sect.z = org.z + dir.z * dist;

                if (sect.x >= 0.0f && sect.x <= vx &&
                        sect.z >= 0.0f && sect.z <= vz &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }

                // voly-plane
                dist = (vy - org.y) / dir.y;

                sect.x = org.x + dir.x * dist;
                sect.y = (float)vy;
                sect.z = org.z + dir.z * dist;

                if (sect.x >= 0.0f && sect.x <= vx &&
                        sect.z >= 0.0f && sect.z <= vz &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }
        }

        // intersect with z-planes
        if (fabs(dir.z) > 0.001f)
        {
                // 0-plane
                dist = -org.z / dir.z;

                sect.x = org.x + dir.x * dist;
                sect.y = org.y + dir.y * dist;
                sect.z = 0.0f;

                if (sect.x >= 0.0f && sect.x <= vx &&
                        sect.y >= 0.0f && sect.y <= vy &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }

                // volz-plane
                dist = (vz - org.z) / dir.z;

                sect.x = org.x + dir.x * dist;
                sect.y = org.y + dir.y * dist;
                sect.z = (float)vz;

                if (sect.x >= 0.0f && sect.x <= vx &&
                        sect.y >= 0.0f && sect.y <= vy &&
                        dist > 0.0f)
                {
                        // valid intersection - find if closer or farther than existing
                        nsect++;

                        if (dist < dist_near)
                        {
                                vfar = vnear;
                                dist_far = dist_near;

                                vnear = sect;
                                dist_near = dist;
                        }
                        else
                        {
                                vfar = sect;
                                dist_far = dist;
                        }
                }
        }

        return nsect;
}

static __inline float GetRand01()
{
        return (float)rand() / (float)RAND_MAX;
}

static __inline vect GetGradient(vect pos,              
                                                                 float base             
                                                                 )
{
        vect gradient;

        pos.x += vol_opts.gradstep;
        gradient.x = GetInterpolatedSample(pos) - base;

        pos.x -= vol_opts.gradstep;
        pos.y += vol_opts.gradstep;
        gradient.y = GetInterpolatedSample(pos) - base;

        pos.y -= vol_opts.gradstep;
        pos.z += vol_opts.gradstep;
        gradient.z = GetInterpolatedSample(pos) - base;

        return gradient;
}

bool ShadeStep(vect &dir,                       
                           vect &pos,                   
                           color &col,                  
                           float steplen                
                           )
{
        if (col.a > 0.99f)
                return false;

        float val = GetInterpolatedSample(pos);

        vect grad = GetGradient(pos, val);

        float len = norm(grad);

        if (len > 0.0f)
        {
                grad.x /= len;
                grad.y /= len;
                grad.z /= len;
        }

        histo_handle    hh;

        CreateInterpolatedHandle(hh, val);

        len = 1.0f * (1.0f - hh.gradient) + len * hh.gradient;
        hh.out_alpha *= len;

        hh.out_alpha *= steplen;

        hh.out_alpha = clamp(hh.out_alpha, 0.0f, 1.0f);

        color end;

        end.r = hh.col_ambient[0];
        end.g = hh.col_ambient[1];
        end.b = hh.col_ambient[2];

        float kd = dot(ldir, grad);

        if (kd < 0.0f)
                kd = 0.0f;

        end.r += kd * hh.col_diffuse[0];
        end.g += kd * hh.col_diffuse[1];
        end.b += kd * hh.col_diffuse[2];

        vect half;

        half.x = dir.x + ldir.x;
        half.y = dir.y + ldir.y;
        half.z = dir.z + ldir.z;

        float half_len = norm(half);

        if (half_len > 0.0f)
        {
                half.x /= half_len;
                half.y /= half_len;
                half.z /= half_len;
        }

        float ks = dot(half, grad);

        if (ks < 0.0f)
                ks = 0.0f;

        ks = pow(ks, 1.0f / hh.spec_exponent);

        end.r += ks * hh.col_specular[0];
        end.g += ks * hh.col_specular[1];
        end.b += ks * hh.col_specular[2];


        float newweight = (1.0f - col.a) * hh.out_alpha;

        col.r += end.r * newweight;
        col.g += end.g * newweight;
        col.b += end.b * newweight;
        col.a += newweight;

        return true;
}

color ShootRay(vect org,                
                           vect dir                     
                           )
{
        color col;
        col.r = col.g = col.b = col.a = 0.0f;

        vect sect1, sect2, pos;

        int nsects = GetBoxIntersection(sect1, sect2, org, dir);

        if (nsects > 0)
        {
                if (nsects == 1)
                {
                        sect2 = sect1;
                        sect1 = org;
                }
                
                vect dist;

                dist.x = sect1.x - sect2.x;
                dist.y = sect1.y - sect2.y;
                dist.z = sect1.z - sect2.z;

                float range = norm(dist);
                float travelled = 0.0f;

                pos = sect1;

                if (vol_opts.rendermode == MODE_FULL)
                {
                        while (travelled < range)
                        {
                                float rnd = (GetRand01() * 2.0f - 1.0f) * vol_opts.jitter * vol_opts.stepsize + vol_opts.stepsize;

                                if (!ShadeStep(dir, pos, col, rnd))
                                {
                                        travelled = range;
                                }

                                pos.x += dir.x * rnd;
                                pos.y += dir.y * rnd;
                                pos.z += dir.z * rnd;

                                travelled += rnd;
                        }
                }
                else
                {
                        float maxval = 0.0f;

                        while (travelled < range)
                        {
                                float rnd = (GetRand01() * 2.0f - 1.0f) * vol_opts.jitter * vol_opts.stepsize + vol_opts.stepsize;

                                float val = GetInterpolatedSample(pos);

                                if (val > maxval)
                                {
                                        maxval = val;
                                }

                                if (maxval >= 0.999f)
                                        break;

                                pos.x += dir.x * rnd;
                                pos.y += dir.y * rnd;
                                pos.z += dir.z * rnd;

                                travelled += rnd;
                        }

                        histo_handle hh;
                        CreateInterpolatedHandle(hh, maxval);

                        col.r = hh.col_diffuse[0];
                        col.g = hh.col_diffuse[1];
                        col.b = hh.col_diffuse[2];
                        col.a = 1.0f;
                }
        }

        return col;
}

void DrawRect(int x, int y,                             
                          int height, int width,        
                          color c, int size,            
                          color *fbuffer)                       

{
        for (int j = y; j < (y+size >= height?height-1:y+size); j++)
        {
                for (int i = x; i < (x+size >= width?width-1:x+size); i++)
                {
                        fbuffer[j * width + i] = c;
                }
        }
}

void RayRender()
{
        color *fbuffer;
        RECT crect;
        MSG msg;

        GetClientRect(wnd_3d, &crect);

        int width = crect.right;
        int height = crect.bottom;

        fbuffer = new color[width * height];

        vect dir;
        dir.x = cam.org.x - cam.pos.x;
        dir.y = cam.org.y - cam.pos.y;
        dir.z = cam.org.z - cam.pos.z;

        normalize(dir);

        vect up = cam.up;

        vect right = cross(dir, up);
        normalize(right);

        vect raydir;
        
        float aspect = (float)width / (float)height;
        float vert_open = atan((float)M_PI / 6.0f); // vertical FOV of 60 deg
        float horz_open = vert_open * aspect;

        wglMakeCurrent(rayHdc, rayHglrc);

        RECT clientRect;
        GetClientRect(wnd_3d, &clientRect);

        glViewport(0, 0, clientRect.right, clientRect.bottom);

        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();

        gluOrtho2D(0, clientRect.right, clientRect.bottom, 0.0f);

        glMatrixMode(GL_MODELVIEW);

        glDisable(GL_BLEND);
        glDisable(GL_DEPTH_TEST);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        for (int n = MAX_ITERATION; n >= 0; n--)
        {
                for (int y = 0; y < height; y += (1 << n))
                {
                        for (int x = 0; x < width; x += (1 << n))
                        {
                                if (n == MAX_ITERATION || 
                                        !(x % (1 << (n+1)) == 0 &&
                                         y % (1 << (n+1)) == 0))
                                {
                                        float xpart = (float)(x * 2 - width) / (float)width * horz_open;
                                        float ypart = (float)(y * 2 - height) / (float)height * vert_open;

                                        raydir.x = dir.x + right.x * xpart;
                                        raydir.y = dir.y + right.y * xpart;
                                        raydir.z = dir.z + right.z * xpart;

                                        raydir.x = raydir.x + up.x * ypart;
                                        raydir.y = raydir.y + up.y * ypart;
                                        raydir.z = raydir.z + up.z * ypart;

                                        normalize(raydir);

                                        DrawRect(x, y, height, width, 
                                                ShootRay(cam.pos, raydir), 1 << n, fbuffer);
                                }
                        }
                        
                        if (y % (1 << (n + 5)) == 0)
                        {
                                glDrawPixels(width, height, GL_RGBA, GL_FLOAT, fbuffer);

                                glBegin(GL_LINES);
                                        glColor3f(0.0f, 0.0f, 1.0f);
                                        glVertex2i(0, height - y);
                                        glVertex2i(width, height - y);
                                glEnd();

                                glFinish();

                                SwapBuffers(rayHdc);
                        }

                        while (PeekMessage(&msg, wnd_main, WM_MOUSEFIRST, WM_MOUSELAST, PM_REMOVE))
                        {
                                switch(msg.message)
                                {
                                case WM_LBUTTONDOWN:
                                case WM_MBUTTONDOWN:
                                case WM_RBUTTONDOWN:
                                        glEnable(GL_DEPTH_TEST);
                                        glEnable(GL_BLEND);

                                        delete [] fbuffer;

                                        PostMessage(msg.hwnd, msg.message, msg.wParam, msg.lParam);

                                        return;
                                }
                        }

                        if (PeekMessage(&msg, wnd_main, WM_KEYFIRST, WM_KEYLAST, PM_NOREMOVE))
                        {
                                switch(msg.message)
                                {
                                case WM_KEYDOWN:
                                        glEnable(GL_DEPTH_TEST);
                                        glEnable(GL_BLEND);

                                        delete [] fbuffer;

                                        return;
                                }
                        }
                }
        }

        glDrawPixels(width, height, GL_RGBA, GL_FLOAT, fbuffer);

        glFinish();

        SwapBuffers(rayHdc);

        glEnable(GL_DEPTH_TEST);
        glEnable(GL_BLEND);

        delete [] fbuffer;
}

---