---

data.cpp

Go to the documentation of this file.

#include "stdafx.h"
#include "math.h"
#include "data.h"

gradient_t operator * (gradient_t g, float s)
{
        gradient_t res;
        res.x = g.x * s;
        res.y = g.y * s;
        res.z = g.z * s;
        return res;
}

gradient_t operator + (gradient_t g1, gradient_t g2)
{
        gradient_t res;
        res.x = g1.x + g2.x;
        res.y = g1.y + g2.y;
        res.z = g1.z + g2.z;
        return res;
}

int
Data::GetDensity(int x, int y, int z)
{
        if(x <0) x =0;
        if(y<0) y=0;
        if(z<0) z=0;
        if(x > xDim - 1) x = xDim-1;
        if(y > yDim - 1) y = yDim-1;
        if(z > zDim - 1) z = zDim-1;
        return (int)data[z*xyDim+y*xDim+x];
}

bool
Data::CalcDensityHistogram()
{
        for (int i = 0; i < 4096; i++)
        {
                histogram[i] = 0;
        }

        for (int j = 0; j < size; j++)
        {
                histogram[data[j]] += 1;
        }
        
        return true;
}

int
Data::GetHistogram(int density)
{
        if(density < 0) density = 0;
        if(density >= 4096) density = 4095;
        return histogram[density];
}


gradient_t
Data::CalcGrad(int x, int y, int z) {
        if(precalc && gradients == NULL) CalcGradients();
        if(x < 0) x = 0;
        if(x > xDim-1) x = xDim-1;
        if(z < 0) z = 0;
        if(z > zDim-1) z = zDim-1;
        if(y < 0) y = 0;
        if(y > yDim-1) y = yDim-1;
        gradient_t res;
        if(!precalc) {
                res.x = ((float)(GetDensity(x+1,y,z) - GetDensity(x-1,y,z)))/2.0;
                res.y = ((float)(GetDensity(x,y+1,z) - GetDensity(x,y-1,z)))/2.0;
                res.z = ((float)(GetDensity(x,y,z+1) - GetDensity(x,y,z-1)))/2.0;
        } else {
                res = gradients[z*xyDim+y*xDim+x];
        }
        return res;
}

bool
Data::CalcGradients() {
        if (gradients != NULL) return false;
        gradients = new gradient_t[xDim*yDim*zDim];

        for(int z = 0; z < zDim; z++) {
                for(int y = 0; y < yDim; y++) {
                        for(int x = 0; x < xDim; x++) {
                                gradients[z*xyDim+y*xDim+x].x = ((float)(GetDensity(x+1,y,z) - GetDensity(x-1,y,z)))/2.0;
                                gradients[z*xyDim+y*xDim+x].y = ((float)(GetDensity(x,y+1,z) - GetDensity(x,y-1,z)))/2.0;
                                gradients[z*xyDim+y*xDim+x].z = ((float)(GetDensity(x,y,z+1) - GetDensity(x,y,z-1)))/2.0;
                        }
                }
        }
        
        return true;
}

void
Data::LoadData(char* fname)
{

        char *filename = (char *)fname;
        
        FILE *filePtr;
        
        // open File
        filePtr = fopen(filename, "rb");
        if(!filePtr) return;// false;

        //getting headerinfo
        fread(&xDim, sizeof(short),1,filePtr);
        fread(&yDim, sizeof(short),1,filePtr);
        fread(&zDim, sizeof(short),1,filePtr);


        xyDim = (int)xDim*yDim;


        size = xyDim * zDim;
        data = new short[size];

        
        CWnd* pFrame = AfxGetMainWnd();

        int loadSize = (size) / 100;
        int counter = loadSize;
        int pos = 0;
        int fileSize = 0;

        /* 
        size means the size in data-sets!!
        partially loading if approx filesize is greater than 40MB
        (size * 16) / 8 / 1024 / 1024  =  (size >> 19)
        */
        if ((m_FileThreshold != 0) &&
                (fileSize = DATASET_TO_MB(size)) >= m_FileThreshold) {
                int pieces = size / m_DataPackets;




                pFrame->SendMessage(MYWM_PROGRESS_MIN_MAX, 0, pieces);

                int tempSize = size;

                for (int i = 0; i < pieces; i++) {
                        pFrame->SendMessage(MYWM_PROGRESS, i);
                        fread(data + i * m_DataPackets, sizeof(short), m_DataPackets, filePtr);
                        tempSize -= m_DataPackets;
                }

                pFrame->SendMessage(MYWM_PROGRESS, i);
                if (tempSize > 0) 
                        fread(data + i * m_DataPackets, sizeof(short), tempSize, filePtr);
        }
        else {
                fread(data,sizeof(short),size,filePtr);
        }

        fclose(filePtr);

        pFrame->SendMessage(MYWM_PROGRESS_MIN_MAX, 0, 100);

        //only 12 bit of data are used for each voxel (4096 density values)
        for(int z=0; z < zDim; z++) {
                for(int y=0; y < yDim; y++) {
                        for(int x=0; x < xDim; x++) {
                                data[z*xyDim+y*xDim+x] &= 0x0fff;

                                //progressbar
                                if (counter-- <= 0) {
                                        counter = loadSize;
                                        pFrame->SendMessage(MYWM_PROGRESS, pos++);
                                }
                        }
                }
        }

        pFrame->SendMessage(MYWM_PROGRESS, 0);

        
        this->filename = filename;

        if(!CalcDensityHistogram()) return;// false;
        
//      return true;
}

short
Data::GetXDim()
{
        return xDim;
}

short
Data::GetYDim()
{
        return yDim;
}

short
Data::GetZDim()
{
        return zDim;
}

---