Image.h

Go to the documentation of this file.

/*

Copyright (c) 2017, Dr Franck P. Vidal (franck.p.vidal@fpvidal.net),
http://www.fpvidal.net/
All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.

3. Neither the name of the Bangor University nor the names of its contributors
may be used to endorse or promote products derived from this software without
specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/


#ifndef __Image_h
#define __Image_h


//******************************************************************************
//  Include
//******************************************************************************
#ifndef NOMINMAX
#define NOMINMAX
#endif

#include <string>
#include <map>

#ifdef HAS_GDCM
#include <gdcmFile.h>
#endif

#ifdef HAS_TIFF
#include <tiff.h>
#endif

#ifndef __Types_h
#include "gVirtualXRay/Types.h"
#endif


//******************************************************************************
//  namespace
//******************************************************************************
namespace gVirtualXRay {


//******************************************************************************
//  Class/type declaration
//******************************************************************************
template<typename T> class Image;
template<typename T> class FFT;
template<typename T> class Sinogram;
class PolygonMesh;

// This is needed for MacOS
#ifndef TIFF_UINT16_T
#define TIFF_UINT16_T unsigned short
#endif

#ifndef TIFF_UINT32_T
#define TIFF_UINT32_T unsigned int
#endif


//******************************************************************************
//  Function declaration
//******************************************************************************
template<typename T> Image<T> operator*(const T& aValue,
                                        const Image<T>& anImage);


template<typename T> Image<T> operator/(const T& aValue,
                                        const Image<T>& anImage);


template<typename T> Image<T> operator+(const T& aValue,
                                        const Image<T>& anImage);


template<typename T> Image<T> operator-(const T& aValue,
                                        const Image<T>& anImage);


template<typename T> Image<T> log(const Image<T>& anImage);


template<typename T> Image<T> abs(const Image<T>& anImage);


std::string getPixelType(const std::string& aFileName);


std::string getPixelType(const char* aFileName);


template<typename T> Image<T> gauss2D(unsigned int aSize, double aSigmaValue);


//==============================================================================
//==============================================================================
template<typename T> class Image
//------------------------------------------------------------------------------
{
//******************************************************************************
public:
    //--------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    Image();


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image(const char* aFileName);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image(const std::string& aFileName);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image(const T* apData,
                 unsigned int aWidth,
                 unsigned int aHeight,
                 unsigned int aNumberOfSlices = 1,
                 const VEC3& aVoxelSize = VEC3(1, 1, 1));


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image(unsigned int aWidth,
                 unsigned int aHeight,
                 unsigned int aNumberOfSlices = 1,
                 T aDefaultValue = 0,
                 const VEC3& aVoxelSize = VEC3(1, 1, 1));


    //--------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    ~Image();


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void makeBlank(const T& aDefaultColour = 0);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image getROI(unsigned int i,
            unsigned int j,
            unsigned int k,
            unsigned int aWidth,
            unsigned int aHeight,
            unsigned int aDepth) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void setPixel(unsigned int i, unsigned int j, unsigned int k, T aValue);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T& getPixel(unsigned int i, unsigned int j, unsigned int k = 0);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    const T& getPixel(unsigned int i, unsigned int j, unsigned int k = 0) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T& getPixelValue(unsigned int i, unsigned int j, unsigned int k = 0);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    const T& getPixelValue(unsigned int i,
                           unsigned int j,
                           unsigned int k = 0) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    VEC3 getPixelPosition(unsigned int i,
                          unsigned int j,
                          unsigned int k = 0) const;


    T& operator[](unsigned anIndex);


    const T& operator[](unsigned anIndex) const;


    T& operator()(unsigned int i, unsigned int j, unsigned int k);


    const T& operator()(unsigned int i, unsigned int j, unsigned int k) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator=(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator=(const T* apImage);


    //--------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    void destroy();


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T* getRawData();


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    const T* getRawData() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    unsigned int getWidth() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    unsigned int getHeight() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    unsigned int getDepth() const;


    void setSpacing(const VEC3& aPixelSize);
    void setSpacing(const double& aPixelWidth,
            const double& aPixelHeight = 0.0,
            const double& aPixelDepth = 0.0);


    const VEC3& getSpacing() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T getMinValue() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T getMaxValue() const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T getMinValue(unsigned int aSliceID) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    T getMaxValue(unsigned int aSliceID) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image shiftScaleFilter(double aShiftValue,
            double aScaleValue,
            bool aRunOnGPUFlag = true) const;


    Image<unsigned char> threshold(const T& aLowerThreshold, const T& aUpperThreshold,
            unsigned char anInValue,
            unsigned char anOutValue,
            bool aRunOnGPUFlag = true) const;

    Image<T> convolution1D(const std::vector<float>& aKernel) const;
    Image<T> convolution2D(const Image<float>& aKernel) const;
    Image<T> convolution2D(const Image<double>& aKernel) const;

    Image<unsigned char> opening(unsigned int aRadius) const;
    Image<unsigned char> closing(unsigned int aRadius) const;
    Image<float> getIsoTropic() const;
    Image<float> antiAliasBinaryImageFilter(double aMaximumRMSChange, unsigned int aNumberOfIterations, double aVariance) const;


    Image constantPadFilter(const T& aPadValue) const;


    Image<unsigned char> threshold(const T& aThreshold,
        unsigned char aValueIn = 1,
        unsigned char aValueOut = 0) const;
    
    
    Image<unsigned char> threshold(const T& aLowerThreshold,
        const T& aUpperThreshold,
        unsigned char aValueIn = 1,
        unsigned char aValueOut = 0) const;


    //--------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    Image normalised() const;


    //--------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    Image normalized() const;


    void load(const char* aFileName);


    void load(const std::string& aFileName);


    void loadUsingITK(const char* aFileName);


    void loadUsingITK(const std::string& aFileName);


    void loadSeries(const char* aPattern,
                    int aFirstSliceIndex,
                    int aLastSliceIndex,
                    int anIncrementIndex = 1);


    void loadSeries(const std::string& aPattern,
                    int aFirstSliceIndex,
                    int aLastSliceIndex,
                    int anIncrementIndex = 1);


    void loadDicomSeries(const char* aDirectory);


    void loadDicomSeries(const std::string& aDirectory);


    void loadASCII(const char* aFileName);


    void loadASCII(const std::string& aFileName);


    void loadASCII(const char* aFileName, unsigned int aWidth, unsigned int aHeight, unsigned int aDepth, const VEC3& aVoxelSize);


    void loadASCII(const std::string& aFileName, unsigned int aWidth, unsigned int aHeight, unsigned int aDepth, const VEC3& aVoxelSize);


    void loadMHD(const char* aFileName);


    void loadMHD(const std::string& aFileName);


    void loadMHA(const char* aFileName);


    void loadMHA(const std::string& aFileName);


    unsigned int sizeOf(const char* anElementType) const;

    unsigned int sizeOf(const std::string& anElementType) const;

    void loadRAW(const char* aFileName,
                 unsigned int aWidth,
                 unsigned int aHeight,
                 unsigned int aNumberOfSlices = 1,
                 bool aChangeEndianessFlag = false,
                 const char* anElementType = 0);


    void loadRAW(const std::string& aFileName,
                 unsigned int aWidth,
                 unsigned int aHeight,
                 unsigned int aNumberOfSlices = 1,
                 bool aChangeEndianessFlag = false,
                 const std::string& anElementType = "");


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void save(const char* aFileName,
              bool anInvertLUTFlag = false,
              const char* aComment = "",
          bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void save(const std::string& aFileName,
              bool anInvertLUTFlag = false,
              const std::string& aComment = "",
              bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void savePGM(const char* aFileName) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void savePGM(const std::string& aFileName) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveRaw(const char* aFileName, bool useDeflateCompressionIfPossible) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveRaw(const std::string& aFileName, bool useDeflateCompressionIfPossible) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveRAW(const char* aFileName, bool useDeflateCompressionIfPossible) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveRAW(const std::string& aFileName, bool useDeflateCompressionIfPossible) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveASCII(const char* aFileName) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveASCII(const std::string& aFileName) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveMHD(const char* aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveMHD(const std::string& aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveMHA(const char* aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveMHA(const std::string& aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveDCM(const char* aFileName,
                 bool anInvertLUTFlag = false,
                 const char* aComment = "") const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveDCM(const std::string& aFileName,
                 bool anInvertLUTFlag = false,
                 const std::string& aComment = "") const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveTIFF(const char* aFileName,
                  bool anInvertLUTFlag = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveTIFF(const std::string& aFileName,
                  bool anInvertLUTFlag = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveUsingITK(const char* aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void saveUsingITK(const std::string& aFileName,
        bool useDeflateCompressionIfPossible = false) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    bool operator==(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    bool operator!=(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator+(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator-(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator*(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator/(const Image<T>& anImage) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator+=(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator-=(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator*=(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator/=(const Image<T>& anImage);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator+(T aValue) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator-(T aValue) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator*(T aValue) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator/(T aValue) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator+=(T aValue);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator-=(T aValue);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator*=(T aValue);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image<T>& operator/=(T aValue);


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    Image operator!() const;


    Image<T> abs() const;


    Image<T> log() const;


    Image<T> flipVertically() const;


    Image<T> flipHorizontally() const;


    double getSum() const;


    double getAverage() const;


    double getVariance() const;


    double getStandardDeviation() const;


    double computeSAE(const Image<T>& anImage) const;


    double computeMAE(const Image<T>& anImage) const;


    double computeSSE(const Image<T>& anImage) const;


    double computeMSE(const Image<T>& anImage) const;


    double computeRMSE(const Image<T>& anImage) const;


    double computeNCC(const Image<T>& anImage) const;


    gVirtualXRay::FFT<T> getFFT() const;


    Image<T> transpose() const;


    Image<T> resize(unsigned int aNewWidth,
                    unsigned int aNewHeightunsigned,
                    unsigned int aNewNumberOfSlices) const;


    Image<T> rotate(float anAngleInDegrees) const;


    Image<T> setCanvasSize(unsigned int aNewWidth,
                           unsigned int aNewHeight,
                           unsigned int aNewNumberOfSlices) const;


    gVirtualXRay::Sinogram<T> radonTransform(double aFirstAngle,
                                             double anAngleStep,
                                             double aLastAngle) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    PolygonMesh marchingCubes(const T& aFirstThresholdValue,
        const T& aSecondThresholdValue = 0,
        bool aTwoThresholdFlag = false,
        unsigned char aVerboseLevel = 0) const;


    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    void convertVoxelsToCubes(std::map<T, PolygonMesh>& aPolygonMeshSet,                     
        T aThreshold = -1000,
        unsigned char aVerboseLevel = 0) const;


    Image<T> convertHU2mu(double anEnergy) const;


    void allocateHostMemory();


    // Load the image as an OpenGL texture
    void transferHostToDevice();


    void transferDeviceToHost();


    // Destroy the OpenGL texture of the image if any
    void destroyOpenGLTexture();


    unsigned int getTextureId() const;


//******************************************************************************
protected:
    VEC3 VertexInterp(const T& aThresholdValue,
            const VEC3& p0,
            const VEC3& p1,
            const T& v0,
            const T& v1) const;


    unsigned int m_width;


    unsigned int m_height;


    unsigned int m_number_of_slices;


    VEC3 m_spacing;


    T* m_p_image;


    unsigned int m_texture_id;


    static unsigned int g_mesh_type;


    static unsigned int g_vbo_type;


//******************************************************************************
private:
    void loadMetaHeader(const char* aFileName);


    void loadMetaHeader(const std::string& aFileName);


    void loadRAW(std::ifstream& anInputStream,
                 unsigned int aWidth,
                 unsigned int aHeight,
                 unsigned int aNumberOfSlices,
                 bool aChangeEndianessFlag,
                 const char* aFileName,
                 const char* anElementType/*,
                 bool aCompressionFlag*/);


    void saveMetaHeader(const char* aFileName, bool useDeflateCompressionIfPossible) const;


    void saveMetaHeader(const std::string& aFileName, bool useDeflateCompressionIfPossible) const;


#ifdef HAS_GDCM
    void addTag(double aValue,
                unsigned int anAddress1,
                unsigned int anAddress2,
                gdcm::File& aFile) const;


    void addTag(const char* aValue,
                unsigned int anAddress1,
                unsigned int anAddress2,
                gdcm::File& aFile) const;


    void addTag(const std::string& aValue,
                unsigned int anAddress1,
                unsigned int anAddress2,
                gdcm::File& aFile) const;
#endif

};


} // namespace gVirtualXRay


//******************************************************************************
#include "Image.inl"


#endif // __Image_h