3d/datafield.hh

Go to the documentation of this file.

/* -*- mia-c++  -*-
 *
 * This file is part of MIA - a toolbox for medical image analysis 
 * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
 *
 * MIA is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with MIA; if not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef __MIA_3DDATAFIELD_HH
#define __MIA_3DDATAFIELD_HH 1

#include <cstdio>
#include <vector>
#include <cmath>
#include <cassert>

#include <mia/3d/vector.hh>
#include <mia/3d/defines3d.hh>
#include <mia/3d/iterator.hh>
#include <mia/2d/datafield.hh>
#include <mia/core/msgstream.hh>
#include <mia/core/parameter.hh>
#include <mia/core/typedescr.hh>
#include <miaconfig.h>

NS_MIA_BEGIN


#define DECLARE_EXTERN_ITERATORS(TYPE)                                          \
        extern template class  EXPORT_3D range3d_iterator<std::vector<TYPE>::iterator>; \
        extern template class  EXPORT_3D range3d_iterator<std::vector<TYPE>::const_iterator>; \
        extern template class  EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::iterator>; \
        extern template class  EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::const_iterator>; \
        extern template class  EXPORT_3D range2d_iterator<std::vector<TYPE>::iterator>; \
        extern template class  EXPORT_3D range2d_iterator<std::vector<TYPE>::const_iterator>;


#ifdef __GNUC__
#pragma GCC diagnostic push
#ifndef __clang__
#pragma GCC diagnostic ignored "-Wattributes"
#endif
#endif

DECLARE_EXTERN_ITERATORS(double);
DECLARE_EXTERN_ITERATORS(float);
DECLARE_EXTERN_ITERATORS(uint32_t);
DECLARE_EXTERN_ITERATORS(int32_t);
DECLARE_EXTERN_ITERATORS(int16_t);
DECLARE_EXTERN_ITERATORS(uint16_t);
DECLARE_EXTERN_ITERATORS(int8_t);
DECLARE_EXTERN_ITERATORS(uint8_t);
DECLARE_EXTERN_ITERATORS(bool);
DECLARE_EXTERN_ITERATORS(int64_t);
DECLARE_EXTERN_ITERATORS(uint64_t);


DECLARE_EXTERN_ITERATORS(C3DFVector)
DECLARE_EXTERN_ITERATORS(C3DDVector)

#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif 

template <class T>
class  EXPORT_3D T3DDatafield {

        typedef  ::std::vector<typename __holder_type_dispatch<T>::type> data_array;

public:
        

        void make_single_ref() __attribute__((deprecated));

        bool holds_unique_data()const __attribute__((deprecated)){ 
                return true; 
        }
                        


        typedef typename data_array::iterator iterator;
        typedef typename data_array::const_iterator const_iterator;
        typedef typename data_array::const_reference const_reference;
        typedef typename data_array::reference reference;
        typedef typename data_array::const_pointer const_pointer;
        typedef typename data_array::pointer pointer;
        typedef typename data_array::value_type value_type;
        typedef typename data_array::size_type size_type;
        typedef typename data_array::difference_type difference_type;
        typedef typename atomic_data<value_type>::type atomic_type; 
        typedef range3d_iterator<iterator> range_iterator; 
        typedef range3d_iterator<const_iterator> const_range_iterator; 

        typedef range3d_iterator_with_boundary_flag<iterator> range_iterator_with_boundary_flag; 
        typedef range3d_iterator_with_boundary_flag<const_iterator> const_range_iterator_with_boundary_flag; 

        typedef C3DBounds dimsize_type;

        class EXPORT_3D Range {
                friend class T3DDatafield<T>;
                friend class ConstRange;
        public:
                
                typedef T3DDatafield<T>::range_iterator iterator;
                
                iterator begin();
                
                iterator end();
                
        private:
                Range(const C3DBounds& start, const C3DBounds& end, T3DDatafield<T>& field);

                iterator m_begin;
                iterator m_end;
        }; 
        
        class EXPORT_3D ConstRange {
        public:
                friend class T3DDatafield<T>;

                typedef T3DDatafield<T>::const_range_iterator iterator;
                
                iterator begin() const;
                
                iterator end() const;

        private:
                ConstRange(const C3DBounds& start, const C3DBounds& end, const T3DDatafield<T>& field);

                ConstRange(const Range& range); 
                
                iterator m_begin;
                iterator m_end;
        }; 

        
        T3DDatafield();

        explicit T3DDatafield(const C3DBounds& _Size);

        T3DDatafield(const C3DBounds& size, const T *data);


        T3DDatafield(const C3DBounds& size, const data_array& data);


        T3DDatafield(const T3DDatafield<T>& org);

        T3DDatafield(T3DDatafield<T>&& org);

        virtual ~T3DDatafield();

        template <typename Out>
        T3DVector<Out> get_gradient(const T3DVector<float >& p) const;

        template <typename Out>
        T3DVector<Out> get_gradient(size_t  x, size_t  y, size_t  z) const;

        template <typename Out>
        T3DVector<Out> get_gradient(int index) const;

        value_type get_interpol_val_at(const T3DVector<float >& p) const __attribute__((deprecated));

        value_type get_block_avrg(const C3DBounds& Start, const C3DBounds& BlockSize) const;

        T3DDatafield& operator = (const T3DDatafield& org);

        T3DDatafield& operator = (T3DDatafield&& org);

        const C3DBounds&  get_size() const
        {
                return m_size;
        }

        void clear();

        size_type size()const
        {
                return m_data.size();
        }

        void swap(T3DDatafield& other);

        value_type get_avg();

        value_type strip_avg();

        const_reference operator()(size_t  x, size_t  y, size_t  z) const
        {
                // Look if we are inside, and give reference, else give the zero
                if (x < m_size.x && y < m_size.y && z < m_size.z) {
                        return m_data[x+ m_size.x * (y  + m_size.y * z)];
                }
                return Zero;
        }


        const_reference operator()(const C3DBounds& l)const
        {
                return (*this)(l.x,l.y,l.z);
        }

        reference operator()(size_t  x, size_t  y, size_t  z)
        {
                // Look if we are inside, and give reference, else throw exception
                // since write access is wanted
                assert(x < m_size.x && y < m_size.y && z < m_size.z);
                return m_data[x + m_size.x *(y + m_size.y * z)];
        }



        reference operator()(const C3DBounds& l)
        {
                return (*this)(l.x,l.y,l.z);
        }

        void get_data_line_x(int y, int z, std::vector<T>& buffer)const;

        void get_data_line_y(int x, int z, std::vector<T>& buffer)const;

        void get_data_line_z(int x, int y, std::vector<T>& buffer)const;

        void put_data_line_x(int y, int z, const std::vector<T> &buffer);

        void put_data_line_y(int x, int z, const std::vector<T> &buffer);

        void put_data_line_z(int x, int y, const std::vector<T> &buffer);

        template <class TMask>
        void mask(const TMask& m);

        void read_xslice_flat(size_t x, std::vector<atomic_type>& buffer) const;

        void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
        
        void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
        
        void write_zslice_flat(size_t z, const std::vector<atomic_type>& buffer); 


        void write_yslice_flat(size_t y, const std::vector<atomic_type>& buffer); 

        void write_xslice_flat(size_t x, const std::vector<atomic_type>& buffer); 

        T2DDatafield<T> get_data_plane_xy(size_t  z)const;
        
        T2DDatafield<T> get_data_plane_yz(size_t  x)const;

        T2DDatafield<T> get_data_plane_xz(size_t  y)const;

        void put_data_plane_xy(size_t  z, const T2DDatafield<T>& p);

        void put_data_plane_yz(size_t  x, const T2DDatafield<T>& p);

        void put_data_plane_xz(size_t  y, const T2DDatafield<T>& p);

        const_iterator begin()const
        {
                return m_data.begin();
        }
        
        const_iterator begin_at(size_t x, size_t y, size_t z)const
        {
                return m_data.begin() + (z * m_size.y + y) * m_size.x + x;
        }


        const_iterator end()const
        {
                return m_data.end();
        }

        iterator begin()
        {
                return m_data.begin();
        }
        
        Range get_range(const C3DBounds& start, const C3DBounds& end);

        ConstRange get_range(const C3DBounds& start, const C3DBounds& end) const; 
        
        range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end); 

        range_iterator end_range(const C3DBounds& begin, const C3DBounds& end); 


        const_range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end)const; 

        const_range_iterator end_range(const C3DBounds& begin, const C3DBounds& end)const; 


        range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end); 

        range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end); 


        const_range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const; 

        const_range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const; 


        iterator begin_at(size_t x, size_t y, size_t z)
        {
                return m_data.begin() + (z * m_size.y + y) * m_size.x + x;
        }

        iterator end()
        {
                return m_data.end();
        }

        const_reference operator[](int i)const
        {
                return m_data[i];
        }

        reference operator[](int i)
        {
                return m_data[i];
        }


        size_t  get_plane_size_xy()const
        {
                return m_xy;
        };

private:
        C3DBounds  m_size;

        size_t  m_xy;

        data_array m_data;

        static const value_type Zero;
        
        static const size_t m_elements; 

};

typedef T3DDatafield<float>  C3DFDatafield;

typedef T3DDatafield<uint32_t> C3DUIDatafield;

typedef T3DDatafield<int32_t>  C3DSIDatafield;

typedef T3DDatafield<uint16_t> C3DUSDatafield;

typedef T3DDatafield<int16_t>  C3DSSDatafield;

typedef T3DDatafield<uint64_t> C3DULDatafield;

typedef T3DDatafield<int64_t>  C3DLDatafield;

typedef T3DDatafield<uint8_t>  C3DUBDatafield;

typedef T3DDatafield<int8_t>  C3DSBDatafield;


typedef T3DDatafield<bool>  C3DBitDatafield;

typedef  CTParameter<C3DBounds> C3DBoundsParameter;

typedef  CTParameter<C3DFVector> C3DFVectorParameter;

typedef  TTranslator<C3DFVector> C3DFVectorTranslator; 

DECLARE_TYPE_DESCR(C3DBounds); 
DECLARE_TYPE_DESCR(C3DFVector);

extern template class EXPORT_3D TAttribute<C3DFVector>; 

// some implementations

template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(size_t  x, size_t  y, size_t  z) const
{
        const int sizex = m_size.x;
        // Look if we are inside the used space
        if (x - 1 < m_size.x - 2 &&  y - 1 < m_size.y - 2 &&  z - 1 < m_size.z - 2) {

                // Lookup all neccessary Values
                const T *H  = &m_data[x + m_size.x * (y + m_size.y * z)];

                return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
                                       Out((H[sizex] - H[-sizex]) * 0.5),
                                       Out((H[m_xy] - H[-m_xy]) * 0.5));
        }

        return T3DVector<Out>();
}


template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(int hardcode) const
{
        const int sizex = m_size.x;
        // Lookup all neccessary Values
        const T *H  = &m_data[hardcode];

        return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
                               Out((H[sizex] - H[-sizex]) * 0.5),
                               Out((H[m_xy] - H[-m_xy]) * 0.5));
}


template <>
template <typename Out>
T3DVector<Out> T3DDatafield<bool>::get_gradient(int hardcode) const
{

        // Lookup all neccessary Values
        return T3DVector<Out> (Out((m_data[hardcode + 1] - m_data[hardcode -1]) * 0.5),
                               Out((m_data[hardcode + m_size.x] - m_data[hardcode -m_size.x]) * 0.5),
                               Out((m_data[hardcode + m_xy] - m_data[hardcode -m_xy]) * 0.5));
}

template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(const T3DVector<float >& p) const
{
        // This will become really funny
        const int sizex = m_size.x;
        // Calculate the int coordinates near the POI
        // and the distances
        size_t  x = size_t (p.x);
        float  dx = p.x - x;
        float  xm = 1 - dx;
        size_t  y = size_t (p.y);
        float  dy = p.y - y;
        float  ym = 1 - dy;
        size_t  z = size_t (p.z);
        float  dz = p.z - z;
        float  zm = 1 - dz;

        // Look if we are inside the used space
        if (x-1 < m_size.x-3 &&  y -1 < m_size.y-3 && z - 1 < m_size.z-3 ) {
                // Lookup all neccessary Values
                const T *H000  = &m_data[x + sizex * y + m_xy * z];

                const T* H_100 = &H000[-m_xy];
                const T* H_101 = &H_100[1];
                const T* H_110 = &H_100[sizex];
                const T* H_111 = &H_110[1];

                const T* H0_10 = &H000[-sizex];
                const T* H0_11 = &H0_10[1];

                const T* H00_1 = &H000[-1];
                const T* H001  = &H000[ 1];
                const T* H002  = &H000[ 2];


                const T* H010  = &H000[sizex];
                const T* H011  = &H010[ 1];
                const T* H012  = &H010[ 2];
                const T* H01_1 = &H010[-1];

                const T* H020 = &H010[sizex];
                const T* H021 = &H020[ 1];

                const T* H100 = &H000[m_xy];

                const T* H1_10 = &H100[sizex];
                const T* H1_11 = &H1_10[1];

                const T* H10_1 = &H100[-1];
                const T* H101  = &H100[ 1];
                const T* H102  = &H100[ 2];

                const T* H110  = &H100[sizex];
                const T* H111  = &H110[ 1];
                const T* H112  = &H110[ 2];
                const T* H11_1 = &H110[-1];

                const T* H120 = &H110[sizex];
                const T* H121 = &H120[ 1];

                const T* H200 = &H100[m_xy];
                const T* H201 = &H200[1];
                const T* H210 = &H200[sizex];
                const T* H211 = &H210[1];

                // use trilinear interpolation to calc the gradient
                return T3DVector<Out> (
                        Out((zm * (ym * (dx * (*H002 - *H000) + xm * (*H001 - *H00_1))+
                                   dy * (dx * (*H012 - *H010) + xm * (*H011 - *H01_1)))+
                             dz * (ym * (dx * (*H102 - *H100) + xm * (*H101 - *H10_1))+
                                   dy * (dx * (*H112 - *H110) + xm * (*H111 - *H11_1)))) * 0.5),

                        Out((zm * (ym * (xm * (*H010 - *H0_10) + dx * (*H011 - *H0_11))+
                                   dy * (xm * (*H020 - *H000)  + dx * (*H021 - *H001)))+
                             dz * (ym * (xm * (*H110 - *H1_10) + dx * (*H111 - *H1_11))+
                                   dy * (xm * (*H120 - *H100)  + dx * (*H121 - *H101)))) * 0.5),
                        Out((zm * (ym * (xm * (*H100 - *H_100) + dx * (*H101 - *H_101))+
                                   dy * (xm * (*H110 - *H_110) + dx * (*H111 - *H_111)))+
                             dz * (ym * (xm * (*H200 - *H000)  + dx * (*H201 - *H001))+
                                   dy * (xm * (*H210 - *H010)  + dx * (*H211 - *H011)))) * 0.5));
        }
        return T3DVector<Out>();

}

#ifdef __GNUC__
#pragma GCC diagnostic push
#ifndef __clang__
#pragma GCC diagnostic ignored "-Wattributes"
#endif
#endif

#define DECLARE_EXTERN(TYPE)                                            \
        extern template class  EXPORT_3D T3DDatafield<TYPE>; 


DECLARE_EXTERN(double);
DECLARE_EXTERN(float);
DECLARE_EXTERN(uint8_t);
DECLARE_EXTERN(uint16_t);
DECLARE_EXTERN(uint32_t);
DECLARE_EXTERN(uint64_t);
DECLARE_EXTERN(int8_t);
DECLARE_EXTERN(int16_t);
DECLARE_EXTERN(int32_t);
DECLARE_EXTERN(int64_t);

DECLARE_EXTERN(C3DFVector);
DECLARE_EXTERN(C3DDVector);

extern template class EXPORT_3D CTParameter<C3DBounds>;
extern template class EXPORT_3D CTParameter<C3DFVector>;
extern template class EXPORT_3D TTranslator<C3DFVector>; 
extern template class EXPORT_3D TAttribute<C3DFVector>; 


#undef DECLARE_EXTERN

#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif 

NS_MIA_END

#endif