D:/C++/Visualisierung/vislab/src/Volume.h

#pragma once

#include "common.h"
#include <stdlib.h>
#include <vector>

#include <algorithm>
#include <functional>

#include "Matrix.h"


class Volume
{
public:

        class Voxel
        {
        public:

                Voxel()
                {
                        SetValue(0.0f);
                };

                Voxel(const Voxel & datOther)
                {
                        m_fValue = datOther.m_fValue;
                };

                Voxel(const float fValue)
                {
                        SetValue(fValue);
                };


                ~Voxel()
                {
                };


                void SetValue(const float fValue)
                {
                        m_fValue = fValue;
                };

                const float GetValue() const
                {
                        return m_fValue;
                };

                const bool operator==(const Voxel &datOther) const
                {
                        return (GetValue() == datOther.GetValue());
                };

                const bool operator!=(const Voxel &datOther) const
                {
                        return !(*this == datOther);
                };

                const bool operator>(const Voxel &datOther) const
                {
                        return GetValue() > datOther.GetValue();
                };

                const bool operator>=(const Voxel &datOther) const
                {
                        return GetValue() >= datOther.GetValue();
                };

                const bool operator<(const Voxel &datOther) const
                {
                        return GetValue() < datOther.GetValue();
                };

                const bool operator<=(const Voxel &datOther) const
                {
                        return GetValue() <= datOther.GetValue();
                };

                const Voxel & operator+=(const Voxel & datOther)
                {
                        m_fValue += datOther.m_fValue;
                        return *this;
                };

                const Voxel & operator-=(const Voxel & datOther)
                {
                        m_fValue -= datOther.m_fValue;
                        return *this;
                };

                const Voxel & operator*=(const float & fOther)
                {
                        m_fValue *= fOther;
                        return *this;
                };

                const Voxel & operator/=(const float & fOther)
                {
                        m_fValue /= fOther;
                        return *this;
                };

                const Voxel operator+(const Voxel & datOther) const
                {
                        Voxel voxNew = *this;
                        voxNew += datOther;
                        return voxNew;
                };

                const Voxel operator-(const Voxel & datOther) const
                {
                        Voxel voxNew = *this;
                        voxNew -= datOther;
                        return voxNew;
                };

                const Voxel operator*(const float & fOther) const
                {
                        Voxel voxNew = *this;
                        voxNew *= fOther;
                        return voxNew;
                };

                const Voxel operator/(const float & fOther) const
                {
                        Voxel voxNew = *this;
                        voxNew /= fOther;
                        return voxNew;
                };

        private:
                float m_fValue;
        };

public:

        Volume() : m_iWidth(1), m_iHeight(1), m_iDepth(1), m_vecVoxels(1)
        {
        };

        Volume(const std::string &strFilename) : m_iWidth(1), m_iHeight(1), m_iDepth(1), m_vecVoxels(1)
        {
                load(strFilename);
        };


        ~Volume(void)
        {
        };


        const Voxel & Get(const int iX, const int iY, const int iZ) const
        {
                return m_vecVoxels[iX + iY*m_iWidth + iZ*m_iWidth*m_iHeight];
        };

        const Voxel & Get(const int iIndex) const
        {
                return m_vecVoxels[iIndex];
        };


        const Voxel * Get() const
        {
                return &(m_vecVoxels.front());
        };


        const int GetWidth() const
        {
                return m_iWidth;
        };


        const int GetHeight() const
        {
                return m_iHeight;
        };


        const int GetDepth() const
        {
                return m_iDepth;
        };

        const int GetSize() const
        {
                return int(m_vecVoxels.size());
        };
        
        void load(const std::string & strFilename)
        {
                std::cout << "- Loading file '" << strFilename << "' ... " << std::endl;
                FILE *fp = NULL;
                
                fopen_s(&fp,strFilename.c_str(),"rb");

                if (!fp)
                {
                        std::cerr << "+ Error loading file." << std::endl << std::endl;
                }
                else
                {

                        char vcPath[1024];
                        char *pFileName = NULL;
                        GetFullPathName(strFilename.c_str(),1024,vcPath,&pFileName);
                        char vcDrive[1024], vcDirectory[1024], vcFilename[1024], vcExtension[1024];
                        _splitpath_s(vcPath,vcDrive,vcDirectory,vcFilename,vcExtension);
                        const std::string strAdditionalFilename =  std::string(vcDrive)+std::string(vcDirectory)+std::string(vcFilename)+std::string(".ini");

                        char vpSpacingX[1024],vpSpacingY[1024],vpSpacingZ[1024];
                        GetPrivateProfileString("DatFile","oldDat Spacing X","1.0",vpSpacingX,256,strAdditionalFilename.c_str());
                        GetPrivateProfileString("DatFile","oldDat Spacing Y","1.0",vpSpacingY,256,strAdditionalFilename.c_str());
                        GetPrivateProfileString("DatFile","oldDat Spacing Z","1.0",vpSpacingZ,256,strAdditionalFilename.c_str());


                        unsigned short uWidth, uHeight, uDepth;
                        fread(&uWidth,sizeof(unsigned short),1,fp);
                        fread(&uHeight,sizeof(unsigned short),1,fp);
                        fread(&uDepth,sizeof(unsigned short),1,fp);

                        m_iWidth = int(uWidth);
                        m_iHeight = int(uHeight);
                        m_iDepth = int(uDepth);

                        const int iSlice = m_iWidth * m_iHeight;
                        const int iSize = iSlice * m_iDepth;
                        m_vecVoxels.resize(iSize);

                        std::vector<unsigned short> vecData;
                        vecData.resize(iSize);

                        fread((void*)&(vecData.front()),sizeof(unsigned short),iSize,fp);
                        fclose(fp);

                        std::cout << "- File loaded." << std::endl;

                        for (int k=0;k<m_iDepth;k++)
                        {
                                for (int j=0;j<m_iHeight;j++)
                                {
                                        for (int i=0;i<m_iWidth;i++)
                                        {
                                                //we convert the data to float values in an interval of [0..1]
                                                const float fValue = min(1.0f,float(vecData[i + j*m_iWidth + k*iSlice]) / 4095.0f);
                                                m_vecVoxels[i+j*m_iWidth+k*iSlice] = Voxel(fValue);
                                        }
                                }
                                std::cout << "\r- Preparing data (" << (k*100) / (m_iDepth-1) << "%) ...";
                        }
                        std::cout << std::endl << "- Data prepared." << std::endl;
                }
        };

        void computeGradients()
        {
                float *volume = (float*)((void * ) Get());
                int actIndex;
                int xpreIndex;
                int xpostIndex;
                int ypreIndex;
                int ypostIndex;
                int zpreIndex;
                int zpostIndex;
                int gradIndex;

                //allocate enough space for the gradients
                gradients = new float[3 * m_iWidth * m_iHeight * m_iDepth];

                float *gradient = new float[3];
                const int iSlice = m_iWidth * m_iHeight;
                float min = 999.99f;
                float max = -999.99f;

                for(int k = 0, kThird = 0; kThird < m_iDepth; k = k+3, kThird++)
                {
                        for(int j = 0, jThird = 0; jThird < m_iHeight; j = j+3, jThird++)
                        {
                                for(int i = 0, iThird = 0; iThird < m_iWidth; i = i+3, iThird++)
                                {
                                        //Calculate the necessary indices
                                        gradIndex = i + j * m_iWidth + k * iSlice;
                                        actIndex = iThird + jThird * m_iWidth + kThird * iSlice;
                                        xpreIndex = actIndex - 1;
                                        xpostIndex = actIndex + 1;
                                        ypreIndex = actIndex - m_iWidth;
                                        ypostIndex = actIndex + m_iWidth;
                                        zpreIndex = actIndex - iSlice;
                                        zpostIndex = actIndex + iSlice;

                                        //x-direction
                                        if(iThird == 0)
                                                //border left
                                                gradient[0] = volume[xpostIndex] - volume[actIndex];
                                        else if (iThird == m_iWidth - 1)
                                                //right border
                                                gradient[0] = volume[actIndex] - volume[xpreIndex];
                                        else
                                                //somewhere between or in the middle
                                                gradient[0] = 0.5f * (volume[xpostIndex] - volume[xpreIndex]);

                                        gradients[gradIndex] = (gradient[0] + 1.0f) / 2.0f;     //transform gradients into [0...1]


                                        if(gradient[0] < min)
                                                min = gradient[0];
                                        if(gradient[0] > max)
                                                max = gradient[0];


                                        //y-direction
                                        if(jThird == 0)
                                                //bottom border
                                                gradient[1] = volume[ypostIndex] - volume[actIndex];
                                        else if (jThird == m_iHeight - 1)
                                                //top border
                                                gradient[1] = volume[actIndex] - volume[ypreIndex];
                                        else
                                                //somewhere between
                                                gradient[1] = 0.5f * (volume[ypostIndex] - volume[ypreIndex]);

                                        gradients[gradIndex + 1] = (gradient[1] + 1.0f) / 2.0f; //transform gradients into [0...1]


                                        if(gradient[1] < min)
                                                min = gradient[1];
                                        if(gradient[1] > max)
                                                max = gradient[1];


                                        //z-direction
                                        if(kThird == 0)
                                                //front bordern
                                                gradient[2] = volume[zpostIndex] - volume[actIndex];
                                        else if (kThird == m_iDepth - 1)
                                                //back border
                                                gradient[2] = volume[actIndex] - volume[zpreIndex];
                                        else
                                                //somewhere between
                                                gradient[2] = 0.5f * (volume[zpostIndex] - volume[zpreIndex]);

                                        gradients[gradIndex + 2] = (gradient[2] * 1.0f) / 2.0f; //transform gradients into [0...1]

                                        if(gradient[2] < min)
                                                min = gradient[2];
                                        if(gradient[2] > max)
                                                max = gradient[2];

                                        
                                }//end for i
                                
                                
                                
                        }//end for j
                        
                        
                        std::cout << "\r- Computing gradients (" <<(kThird * 100) / (m_iDepth-1) << "%) ...";
                        
                }//end for k

                std::cout << std::endl << "- Gradients computed - "<< std::endl;

                
                
        };//end computeGradients()

        float *getGradients()
        {
                return gradients;       
        };

protected:

private:
        std::vector<Voxel> m_vecVoxels;
        int m_iWidth,m_iHeight,m_iDepth;

        float *gradients;
};

---