OpenCMT/uLib
1
0
Fork 0
mirror of https://github.com/OpenCMT/uLib.git synced 2025-12-06 07:21:31 +01:00
Code Issues Packages Projects Releases Wiki Activity

[uLib Geometry]

Browse Source 
adds ImageSpace ImageMap and ImageData
stats ProgrammableAccess to data
This commit is contained in:
Andrea Rigoni
2014-11-06 15:05:54 +00:00
parent d5531caa0c
commit 4b4e0b2959
23 changed files with 962 additions and 378 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Math/CMakeLists.txt
View File

@@ -7,6 +7,7 @@ set(HEADERS
Transform.h
StructuredData.h
# StructuredGrid.h
ImageSpace.h
ImageMap.h
ImageData.h
VoxImage.h
@@ -32,6 +33,7 @@ set(SOURCES
VoxRaytracer.cpp
StructuredData.cpp
# StructuredGrid.cpp
ImageSpace.cpp
ImageMap.cpp
ImageData.cpp
VoxImage.cpp

52
src/Math/DataSet.h
View File

@@ -35,37 +35,61 @@
#include "Math/StructuredData.h"
#include <Core/ProgrammableAccessor.h>
#include <vtkInformation.h>
namespace uLib {
class DataSet_Base {
public:
DataSet_Base() {}
virtual ~DataSet_Base() {}
virtual const void * GetPtr(Id_t id) const = 0;
virtual void * GetPtr(Id_t id) = 0;
virtual size_t Size() const = 0;
class AbstractArray {
public:
virtual void * GetDataPointer(Id_t id) const = 0;
virtual void SetSize(const size_t size) = 0;
virtual const size_t GetSize() const = 0;
virtual ~AbstractArray() {}
};
class DataSetAttributes {
public:
template < typename GeT, typename SeT >
inline void SetScalars(GeT get, SeT set) {
m_Scalars.SetAccessFunctions(get,set);
}
private:
ProgrammableAccessor<double> m_Scalars;
};
template <typename T>
class VectorData : public DataSet_Base, public Vector<T> {
template < typename T >
class DataVector : public AbstractArray {
public:
void * GetDataPointer(Id_t id) const { return (void *)&m_Data.at(id); }
inline void SetSize(const size_t size) { m_Data.resize(size); }
inline const size_t GetSize() const { return m_Data.size(); }
protected:
const void *GetPtr(Id_t id) const { return &this->at(id); }
void * GetPtr(Id_t id) { return &this->at(id); }
uLibRefMacro(Data,Vector<T>)
uLibRefMacro(Attributes,DataSetAttributes)
private:
DataSetAttributes m_Attributes;
Vector<T> m_Data;
};
//class DataSet {
//public:
// virtual Vector3d GetPoint(const Id_t) const = 0;
// virtual Vector<Id_t> GetCell(const Id_t) const = 0;
//private:
//};

57
src/Math/Dense.h
View File

@@ -53,7 +53,9 @@
#include <stdlib.h>
#include <Eigen/Core>
#include <Eigen/Dense>
#include "Core/Types.h"
#include "Core/Serializable.h"
@@ -106,6 +108,61 @@ EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<double>, cd)
#undef EIGEN_MAKE_TYPEDEFS_ALL_SIZES
#undef EIGEN_MAKE_TYPEDEFS
//#define EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix) \
///** \ingroup arraytypedefs */ \
//typedef Eigen::Array<Type, Size, Size> Array##SizeSuffix##SizeSuffix##TypeSuffix; \
///** \ingroup arraytypedefs */ \
//typedef Eigen::Array<Type, Size, 1> Array##SizeSuffix##TypeSuffix;
//#define EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, Size) \
///** \ingroup arraytypedefs */ \
//typedef Eigen::Array<Type, Size, Dynamic> Array##Size##X##TypeSuffix; \
///** \ingroup arraytypedefs */ \
//typedef Eigen::Array<Type, Dynamic, Size> Array##X##Size##TypeSuffix;
//#define EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
//EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 2, 2) \
//EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 3, 3) \
//EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 4, 4) \
//EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, Dynamic, X) \
//EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 2) \
//EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 3) \
//EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
//EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(int, i)
//EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(float, f)
//EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(double, d)
//EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(std::complex<float>, cf)
//EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(std::complex<double>, cd)
//#undef EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES
//#undef EIGEN_MAKE_ARRAY_TYPEDEFS
//#undef EIGEN_MAKE_ARRAY_TYPEDEFS_LARGE
//#define EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, SizeSuffix) \
//using Eigen::Matrix##SizeSuffix##TypeSuffix; \
//using Eigen::Vector##SizeSuffix##TypeSuffix; \
//using Eigen::RowVector##SizeSuffix##TypeSuffix;
//#define EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(TypeSuffix) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 2) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 3) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 4) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, X) \
//#define EIGEN_USING_ARRAY_TYPEDEFS \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(i) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(f) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(d) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(cf) \
//EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(cd)
} // uLib

181
src/Math/ImageData.cpp
View File

@@ -0,0 +1,181 @@
#include "Math/ImageData.h"
using namespace uLib;
void ImageData::SetSize(const Vector3f v)
{
ImageSpace::SetSize( v.array() / this->GetDims().array().cast<float>() );
}
Vector3f ImageData::GetSize() const
{
return ImageSpace::GetSize().array() * this->GetDims().array().cast<float>();
}
bool ImageData::IsInsideBounds(const Vector4f pt) const
{
Vector4f ptl = ImageSpace::GetLocalPoint(pt);
int result = 0;
for ( int i=0; i<3 ;++i) {
result += ptl(i) > (float)this->GetDims()(i);
result += ptl(i) < 0;
}
return result == 0;
}
void ImageData::ExportToVtk (const char *file, bool density_type)
{
FILE * vtk_file = fopen(file,"wb");
assert(vtk_file);
ImageData *voxels = this;
float norm;
if (density_type) {
norm = 40000;
} else norm = 1.E6;
int nx = voxels->GetDims()(0);
int ny = voxels->GetDims()(1);
int nz = voxels->GetDims()(2);
fprintf(vtk_file,
"# vtk DataFile Version 2.0\n"
"Image Builder vtk output\n"
"ASCII\n"
"DATASET STRUCTURED_POINTS\n"
"DIMENSIONS %d %d %d\n"
"SPACING %f %f %f\n"
"ORIGIN %f %f %f\n"
"POINT_DATA %d\n"
"SCALARS volume_scalars float 1\n"
"LOOKUP_TABLE default\n",
nx , ny , nz ,
voxels->GetSpacing()(0),
voxels->GetSpacing()(1),
voxels->GetSpacing()(2),
voxels->GetPosition()(0), // FIX FOR ORIGIN //
voxels->GetPosition()(1),
voxels->GetPosition()(2),
nx * ny * nz
);
Vector3i index(0,0,0);
for (int zv = 0; zv < nz; ++zv) {
for (int yv = 0; yv < ny; ++yv) {
for (int xv = 0; xv < nx; ++xv) {
index << xv,yv,zv;
// write voxel density in mrad^2/cm //
float density = fabs(voxels->GetValue(index)) * norm;
fprintf(vtk_file, "%f ", density);
}
}
}
fclose(vtk_file);
printf("%s vtk file saved\n",file);
}
int ImageData::ImportFromVtk(const char *file)
{
FILE * vtk_file = fopen(file, "r");
if (!vtk_file) return false;
char word[20];
int dx, dy, dz, n_tot;
float sx, sy, sz, ox, oy, oz;
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "DIMENSIONS"));
fscanf(vtk_file, "%d %d %d", &dx, &dy, &dz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "SPACING"));
fscanf(vtk_file, "%f %f %f", &sx, &sy, &sz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "ORIGIN"));
fscanf(vtk_file, "%f %f %f", &ox, &oy, &oz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "POINT_DATA"));
fscanf(vtk_file, "%d", &n_tot);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "default"));
this->SetDims(Vector3i(dx,dy,dz));
this->SetSpacing(Vector3f(sx,sy,sz));
this->SetPosition(Vector3f(ox,oy,oz));
for (int k = 0; k < dz; ++k) {
for (int j = 0; j < dy; ++j) {
for (int i = 0; i < dx; ++i) {
Vector3i idx(i, j, k);
float tmp_val;
fscanf(vtk_file, "%f", &tmp_val);
this->SetValue(idx,fabs(tmp_val)*1E-6);
}
}
}
fclose(vtk_file);
return true;
}
void ImageData::ExportToVtkXml(const char *file, bool density_type)
{
// Not implemented yet //
FILE * vtk_file = fopen(file,"wb");
assert(vtk_file);
ImageData *voxels = this;
float norm;
if (density_type) {
norm = 40000;
} else norm = 1.E6;
int nx = voxels->GetDims()(0);
int ny = voxels->GetDims()(1);
int nz = voxels->GetDims()(2);
fprintf(vtk_file,
"<?xml version=\"1.0\"?>\n"
"<VTKFile type=\"ImageData\" version=\"0.1\" "
"byte_order=\"LittleEndian\">\n"
"<ImageData WholeExtent=\"0 %d 0 %d 0 %d\" "
"Origin=\"%f %f %f\" Spacing=\"%f %f %f\">\n"
"<Piece Extent=\"0 %d 0 %d 0 %d\">\n"
"<CellData Scalars=\"lambda\">\n"
"<DataArray format=\"ascii\" Name=\"lambda\" type=\"Float32\" >\n",
nx, ny, nz,
voxels->GetPosition()(0), // FIX FOR ORIGIN //
voxels->GetPosition()(1),
voxels->GetPosition()(2),
voxels->GetSpacing()(0),
voxels->GetSpacing()(1),
voxels->GetSpacing()(2),
nx, ny, nz
);
Vector3i index(0,0,0);
for (int zv = 0; zv < nz; ++zv) {
for (int yv = 0; yv < ny; ++yv) {
for (int xv = 0; xv < nx; ++xv) {
index << xv,yv,zv;
// write voxel density in mrad^2/cm //
float density = fabs(voxels->GetValue(index)) * norm;
fprintf(vtk_file, "%f ", density);
}
}
}
fprintf(vtk_file,
"\n </DataArray>\n </CellData>\n");
fprintf(vtk_file,
"</Piece>\n </ImageData>\n </VTKFile>\n");
fclose(vtk_file);
printf("%s vtk file saved\n",file);
}

48
src/Math/ImageData.h
View File

@@ -37,18 +37,64 @@ namespace uLib {
class ImageData : public ImageSpace, public ImageMap {
ImageData() {}
public:
ImageData() : ImageMap(Vector3i::Zero()) {}
ImageData(const Vector3i &size) : ImageMap(size) {}
void SetSize(const Vector3f v);
Vector3f GetSize() const;
bool IsInsideBounds(const Vector4f pt) const;
virtual void * GetDataPointer(const Id_t id) const = 0;
virtual float GetValue(const Vector3i id) const = 0;
virtual float GetValue(const int id) const = 0;
virtual void SetValue(const Vector3i id, float value) = 0;
virtual void SetValue(const int id, float value) = 0;
virtual void SetDims(const Vector3i &size) = 0;
void ExportToVtk(const char *file, bool density_type = 0);
void ExportToVtkXml(const char *file, bool density_type = 0);
int ImportFromVtk(const char *file);
};
class DataVectorImage : public ImageData {
public:
inline void * GetDataPointer(const Id_t id) const {
return m_Data->GetDataPointer(id);
}
float GetValue(const Vector3i id) const { return m_Scalars.Get<float>(GetDataPointer(Map(id))); }
float GetValue(const int id) const { return m_Scalars.Get<float>(GetDataPointer(id)); }
void SetValue(const Vector3i id, float value) { m_Scalars.Set<float>(GetDataPointer(Map(id)),value); }
void SetValue(const int id, float value) { m_Scalars.Set<float>(GetDataPointer(id),value); }
void SetDims(const Vector3i &size) {
ImageMap::SetDims(size);
m_Data->SetSize(size.prod());
}
uLibRefMacro(Data,SmartPointer<AbstractArray>)
uLibRefMacro(Scalars,ProgrammableAccessor<double>)
private:
SmartPointer<AbstractArray> m_Data;
ProgrammableAccessor<double> m_Scalars;
};
} // uLib
#endif // IMAGEDATA_H

14
src/Math/ImageMap.h
View File

@@ -31,8 +31,19 @@
#include "Math/Dense.h"
#include "Math/DataSet.h"
//#include "boost/range.hpp"
//#include "boost/range/iterator.hpp"
namespace uLib {
class ImageMap {
public:
typedef enum {
@@ -69,6 +80,9 @@ public:
Vector3i UnMap(int index) const;
private:
Order m_DataOrder;
Vector3i m_Dims;

27
src/Math/ImageSpace.cpp
View File

@@ -28,19 +28,10 @@
namespace uLib {
ImageSpace::ImageSpace(const Vector3i size) :
StructuredData(size)
{}
//ImageSpace::ImageSpace(const Vector3i size) :
// StructuredData(size)
//{}
//void ImageSpace::SetSize(const Vector3f v)
//{
// ContainerBox::SetSize( v.array() / this->GetDims().array().cast<float>() );
//}
//Vector3f ImageSpace::GetSize() const
//{
// return ContainerBox::GetSize().array() * this->GetDims().array().cast<float>();
//}
void ImageSpace::SetSpacing(const Vector3f spacing)
{
@@ -52,18 +43,6 @@ Vector3f ImageSpace::GetSpacing() const
return ContainerBox::GetSize();
}
//bool ImageSpace::IsInsideBounds(const Vector4f pt) const
//{
// Vector4f ptl = this->GetLocalPoint(pt);
// int result = 0;
// for ( int i=0; i<3 ;++i) {
// result += ptl(i) > (float)this->GetDims()(i);
// result += ptl(i) < 0;
// }
// return result == 0;
//}
Vector3i ImageSpace::Find(const Vector4f pt) const
{
Vector3i out;

9
src/Math/ImageSpace.h
View File

@@ -35,25 +35,18 @@ namespace uLib {
class ImageSpace : public ContainerBox {
typedef ContainerBox BaseClass;
public:
ImageSpace(const Vector3i size);
using ContainerBox::SetOrigin;
// void SetSize(const Vector3f v);
// Vector3f GetSize() const;
void SetSpacing(const Vector3f spacing);
Vector3f GetSpacing() const;
// bool IsInsideBounds(const Vector4f pt) const;
Vector3i Find(const Vector4f pt) const;
void PrintSelf(std::ostream &o);
private:
protected:
using ContainerBox::GetSize;
using ContainerBox::SetSize;
};

4
src/Math/Line.h
View File

@@ -50,6 +50,10 @@ public:
Vector & direction() { return this->m_Data[1]; }
const Vector & origin() const { return this->m_Data[0]; }
const Vector & direction() const { return this->m_Data[1]; }
_Scalar & origin(Id_t id) { return this->m_Data[0](id); }
_Scalar & direction(Id_t id) { return this->m_Data[1](id); }
const _Scalar & origin(Id_t id) const { return this->m_Data[0](id); }
const _Scalar & direction(Id_t id) const { return this->m_Data[1](id); }
};
typedef Line<float,2> Line2f;

84
src/Math/StructuredData.cpp
View File

@@ -27,54 +27,54 @@
#include "StructuredData.h"
using namespace uLib;
//using namespace uLib;
StructuredData::StructuredData(const Vector3i &size) :
m_Dims(size)
{
SetDataOrder();
}
//StructuredData::StructuredData(const Vector3i &size) :
// m_Dims(size)
//{
// SetDataOrder();
//}
void StructuredData::SetDims(const Vector3i &size)
{
this->m_Dims = size;
SetDataOrder();
}
//void StructuredData::SetDims(const Vector3i &size)
//{
// this->m_Dims = size;
// SetDataOrder();
//}
void StructuredData::SetDataOrder(StructuredData::Order order)
{
int i = order & 0x3;
int j = (order >> 2) & 0x3;
int k = (order >> 4) & 0x3;
this->m_Increments[i] = 1;
this->m_Increments[j] = m_Dims[i];
this->m_Increments[k] = m_Dims[i] * m_Dims[j];
this->m_DataOrder = order;
}
//void StructuredData::SetDataOrder(StructuredData::Order order)
//{
// int i = order & 0x3;
// int j = (order >> 2) & 0x3;
// int k = (order >> 4) & 0x3;
// this->m_Increments[i] = 1;
// this->m_Increments[j] = m_Dims[i];
// this->m_Increments[k] = m_Dims[i] * m_Dims[j];
// this->m_DataOrder = order;
//}
bool StructuredData::IsInsideGrid(const Vector3i &v) const
{
int vok = 1;
vok *= (v(0) >= 0 && v(0) < m_Dims[0]);
vok *= (v(1) >= 0 && v(1) < m_Dims[1]);
vok *= (v(2) >= 0 && v(2) < m_Dims[2]);
return vok;
}
//bool StructuredData::IsInsideGrid(const Vector3i &v) const
//{
// int vok = 1;
// vok *= (v(0) >= 0 && v(0) < m_Dims[0]);
// vok *= (v(1) >= 0 && v(1) < m_Dims[1]);
// vok *= (v(2) >= 0 && v(2) < m_Dims[2]);
// return vok;
//}
Vector3i StructuredData::UnMap(int index) const
{
Vector3i v( 0,0,0 );
Vector3i iv = m_Increments;
int id = 0;
for(int k=0; k<3; ++k) {
int inc = iv.maxCoeff(&id);
v(id) = index / inc;
index -= v(id) * inc;
iv(id) = 0;
}
return v;
}
//Vector3i StructuredData::UnMap(int index) const
//{
// Vector3i v( 0,0,0 );
// Vector3i iv = m_Increments;
// int id = 0;
// for(int k=0; k<3; ++k) {
// int inc = iv.maxCoeff(&id);
// v(id) = index / inc;
// index -= v(id) * inc;
// iv(id) = 0;
// }
// return v;
//}

86
src/Math/StructuredData.h
View File

@@ -28,74 +28,74 @@
#ifndef STRUCTUREDDATA_H
#define STRUCTUREDDATA_H
#include "Core/Macros.h"
#include "Core/Object.h"
#include "Math/Dense.h"
#include "Math/DataSet.h"
//#include "Core/Macros.h"
//#include "Core/Object.h"
//#include "Math/Dense.h"
//#include "Math/DataSet.h"
namespace uLib {
//namespace uLib {
class StructuredData {
public:
enum _Order
{
CustomOrder = 0,
XYZ = 0 | 1 << 2 | 2 << 4,
XZY = 0 | 2 << 2 | 1 << 4,
YXZ = 1 | 0 << 2 | 2 << 4,
YZX = 2 | 0 << 2 | 1 << 4,
ZXY = 1 | 2 << 2 | 0 << 4,
ZYX = 2 | 1 << 2 | 0 << 4
};
//class StructuredData {
//public:
// enum _Order
// {
// CustomOrder = 0,
// XYZ = 0 | 1 << 2 | 2 << 4,
// XZY = 0 | 2 << 2 | 1 << 4,
// YXZ = 1 | 0 << 2 | 2 << 4,
// YZX = 2 | 0 << 2 | 1 << 4,
// ZXY = 1 | 2 << 2 | 0 << 4,
// ZYX = 2 | 1 << 2 | 0 << 4
// };
typedef enum _Order Order;
// typedef enum _Order Order;
StructuredData(const Vector3i &size);
// StructuredData(const Vector3i &size);
StructuredData(const StructuredData &copy) :
m_DataOrder(copy.m_DataOrder),
m_Dims(copy.m_Dims),
m_Increments(copy.m_Increments)
{}
// StructuredData(const StructuredData &copy) :
// m_DataOrder(copy.m_DataOrder),
// m_Dims(copy.m_Dims),
// m_Increments(copy.m_Increments)
// {}
uLibGetMacro(Dims,Vector3i)
// uLibGetMacro(Dims,Vector3i)
void SetDims(const Vector3i &size);
// void SetDims(const Vector3i &size);
uLibGetSetMacro(Increments,Vector3i)
// uLibGetSetMacro(Increments,Vector3i)
void SetDataOrder(Order order = YXZ);
// void SetDataOrder(Order order = YXZ);
uLibGetMacro(DataOrder,Order)
// uLibGetMacro(DataOrder,Order)
bool IsInsideGrid(const Vector3i &v) const;
// bool IsInsideGrid(const Vector3i &v) const;
inline int Map(Vector3i index) const;
// inline int Map(Vector3i index) const;
Vector3i UnMap(int index) const;
// Vector3i UnMap(int index) const;
private:
Order m_DataOrder;
Vector3i m_Dims;
Vector3i m_Increments; //TODO: make this line matrix //
};
//private:
// Order m_DataOrder;
// Vector3i m_Dims;
// Vector3i m_Increments; //TODO: make this line matrix //
//};
// --- INLINES -------------------------------------------------------------- //
//// --- INLINES -------------------------------------------------------------- //
inline int StructuredData::Map(Vector3i index) const
{
return (m_Increments.transpose() * index);
}
//inline int StructuredData::Map(Vector3i index) const
//{
// return (m_Increments.transpose() * index);
//}
}
//}
#endif // STRUCTUREDDATA_H

275
src/Math/VoxImage.cpp
View File

@@ -27,163 +27,162 @@
#include <stdlib.h>
#include <stdio.h>
#include "VoxImage.h"
namespace uLib {
void Abstract::VoxImage::ExportToVtk (const char *file, bool density_type)
{
FILE * vtk_file = fopen(file,"wb");
assert(vtk_file);
Abstract::VoxImage *voxels = this;
float norm;
if (density_type) {
norm = 40000;
} else norm = 1.E6;
int nx = voxels->GetDims()(0);
int ny = voxels->GetDims()(1);
int nz = voxels->GetDims()(2);
//void Abstract::VoxImage::ExportToVtk (const char *file, bool density_type)
//{
// FILE * vtk_file = fopen(file,"wb");
// assert(vtk_file);
// ImageData *voxels = this;
// float norm;
// if (density_type) {
// norm = 40000;
// } else norm = 1.E6;
// int nx = voxels->GetDims()(0);
// int ny = voxels->GetDims()(1);
// int nz = voxels->GetDims()(2);
fprintf(vtk_file,
"# vtk DataFile Version 2.0\n"
"Image Builder vtk output\n"
"ASCII\n"
"DATASET STRUCTURED_POINTS\n"
"DIMENSIONS %d %d %d\n"
"SPACING %f %f %f\n"
"ORIGIN %f %f %f\n"
"POINT_DATA %d\n"
"SCALARS volume_scalars float 1\n"
"LOOKUP_TABLE default\n",
nx , ny , nz ,
voxels->GetSpacing()(0),
voxels->GetSpacing()(1),
voxels->GetSpacing()(2),
voxels->GetPosition()(0), // FIX FOR ORIGIN //
voxels->GetPosition()(1),
voxels->GetPosition()(2),
nx * ny * nz
);
// fprintf(vtk_file,
// "# vtk DataFile Version 2.0\n"
// "Image Builder vtk output\n"
// "ASCII\n"
// "DATASET STRUCTURED_POINTS\n"
// "DIMENSIONS %d %d %d\n"
// "SPACING %f %f %f\n"
// "ORIGIN %f %f %f\n"
// "POINT_DATA %d\n"
// "SCALARS volume_scalars float 1\n"
// "LOOKUP_TABLE default\n",
// nx , ny , nz ,
// voxels->GetSpacing()(0),
// voxels->GetSpacing()(1),
// voxels->GetSpacing()(2),
// voxels->GetPosition()(0), // FIX FOR ORIGIN //
// voxels->GetPosition()(1),
// voxels->GetPosition()(2),
// nx * ny * nz
// );
Vector3i index(0,0,0);
for (int zv = 0; zv < nz; ++zv) {
for (int yv = 0; yv < ny; ++yv) {
for (int xv = 0; xv < nx; ++xv) {
index << xv,yv,zv;
// write voxel density in mrad^2/cm //
float density = fabs(voxels->GetValue(index)) * norm;
fprintf(vtk_file, "%f ", density);
}
}
}
// Vector3i index(0,0,0);
// for (int zv = 0; zv < nz; ++zv) {
// for (int yv = 0; yv < ny; ++yv) {
// for (int xv = 0; xv < nx; ++xv) {
// index << xv,yv,zv;
// // write voxel density in mrad^2/cm //
// float density = fabs(voxels->GetValue(index)) * norm;
// fprintf(vtk_file, "%f ", density);
// }
// }
// }
fclose(vtk_file);
printf("%s vtk file saved\n",file);
}
// fclose(vtk_file);
// printf("%s vtk file saved\n",file);
//}
int Abstract::VoxImage::ImportFromVtk(const char *file)
{
FILE * vtk_file = fopen(file, "r");
if (!vtk_file) return false;
//int Abstract::VoxImage::ImportFromVtk(const char *file)
//{
// FILE * vtk_file = fopen(file, "r");
// if (!vtk_file) return false;
char word[20];
int dx, dy, dz, n_tot;
float sx, sy, sz, ox, oy, oz;
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "DIMENSIONS"));
fscanf(vtk_file, "%d %d %d", &dx, &dy, &dz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "SPACING"));
fscanf(vtk_file, "%f %f %f", &sx, &sy, &sz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "ORIGIN"));
fscanf(vtk_file, "%f %f %f", &ox, &oy, &oz);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "POINT_DATA"));
fscanf(vtk_file, "%d", &n_tot);
do {
fscanf(vtk_file, "%s", word);
} while (strcmp(word, "default"));
// char word[20];
// int dx, dy, dz, n_tot;
// float sx, sy, sz, ox, oy, oz;
// do {
// fscanf(vtk_file, "%s", word);
// } while (strcmp(word, "DIMENSIONS"));
// fscanf(vtk_file, "%d %d %d", &dx, &dy, &dz);
// do {
// fscanf(vtk_file, "%s", word);
// } while (strcmp(word, "SPACING"));
// fscanf(vtk_file, "%f %f %f", &sx, &sy, &sz);
// do {
// fscanf(vtk_file, "%s", word);
// } while (strcmp(word, "ORIGIN"));
// fscanf(vtk_file, "%f %f %f", &ox, &oy, &oz);
// do {
// fscanf(vtk_file, "%s", word);
// } while (strcmp(word, "POINT_DATA"));
// fscanf(vtk_file, "%d", &n_tot);
// do {
// fscanf(vtk_file, "%s", word);
// } while (strcmp(word, "default"));
this->SetDims(Vector3i(dx,dy,dz));
this->SetSpacing(Vector3f(sx,sy,sz));
this->SetPosition(Vector3f(ox,oy,oz));
// this->SetDims(Vector3i(dx,dy,dz));
// this->SetSpacing(Vector3f(sx,sy,sz));
// this->SetPosition(Vector3f(ox,oy,oz));
for (int k = 0; k < dz; ++k) {
for (int j = 0; j < dy; ++j) {
for (int i = 0; i < dx; ++i) {
Vector3i idx(i, j, k);
float tmp_val;
fscanf(vtk_file, "%f", &tmp_val);
this->SetValue(idx,fabs(tmp_val)*1E-6);
}
}
}
fclose(vtk_file);
return true;
}
// for (int k = 0; k < dz; ++k) {
// for (int j = 0; j < dy; ++j) {
// for (int i = 0; i < dx; ++i) {
// Vector3i idx(i, j, k);
// float tmp_val;
// fscanf(vtk_file, "%f", &tmp_val);
// this->SetValue(idx,fabs(tmp_val)*1E-6);
// }
// }
// }
// fclose(vtk_file);
// return true;
//}
void Abstract::VoxImage::ExportToVtkXml(const char *file, bool density_type)
{
// Not implemented yet //
FILE * vtk_file = fopen(file,"wb");
assert(vtk_file);
Abstract::VoxImage *voxels = this;
float norm;
if (density_type) {
norm = 40000;
} else norm = 1.E6;
int nx = voxels->GetDims()(0);
int ny = voxels->GetDims()(1);
int nz = voxels->GetDims()(2);
//void Abstract::VoxImage::ExportToVtkXml(const char *file, bool density_type)
//{
// // Not implemented yet //
// FILE * vtk_file = fopen(file,"wb");
// assert(vtk_file);
// ImageData *voxels = this;
// float norm;
// if (density_type) {
// norm = 40000;
// } else norm = 1.E6;
// int nx = voxels->GetDims()(0);
// int ny = voxels->GetDims()(1);
// int nz = voxels->GetDims()(2);
fprintf(vtk_file,
"<?xml version=\"1.0\"?>\n"
// fprintf(vtk_file,
// "<?xml version=\"1.0\"?>\n"
"<VTKFile type=\"ImageData\" version=\"0.1\" "
"byte_order=\"LittleEndian\">\n"
// "<VTKFile type=\"ImageData\" version=\"0.1\" "
// "byte_order=\"LittleEndian\">\n"
"<ImageData WholeExtent=\"0 %d 0 %d 0 %d\" "
"Origin=\"%f %f %f\" Spacing=\"%f %f %f\">\n"
// "<ImageData WholeExtent=\"0 %d 0 %d 0 %d\" "
// "Origin=\"%f %f %f\" Spacing=\"%f %f %f\">\n"
"<Piece Extent=\"0 %d 0 %d 0 %d\">\n"
// "<Piece Extent=\"0 %d 0 %d 0 %d\">\n"
"<CellData Scalars=\"lambda\">\n"
"<DataArray format=\"ascii\" Name=\"lambda\" type=\"Float32\" >\n",
nx, ny, nz,
voxels->GetPosition()(0), // FIX FOR ORIGIN //
voxels->GetPosition()(1),
voxels->GetPosition()(2),
voxels->GetSpacing()(0),
voxels->GetSpacing()(1),
voxels->GetSpacing()(2),
nx, ny, nz
);
Vector3i index(0,0,0);
for (int zv = 0; zv < nz; ++zv) {
for (int yv = 0; yv < ny; ++yv) {
for (int xv = 0; xv < nx; ++xv) {
index << xv,yv,zv;
// write voxel density in mrad^2/cm //
float density = fabs(voxels->GetValue(index)) * norm;
fprintf(vtk_file, "%f ", density);
}
}
}
fprintf(vtk_file,
"\n </DataArray>\n </CellData>\n");
// "<CellData Scalars=\"lambda\">\n"
// "<DataArray format=\"ascii\" Name=\"lambda\" type=\"Float32\" >\n",
// nx, ny, nz,
// voxels->GetPosition()(0), // FIX FOR ORIGIN //
// voxels->GetPosition()(1),
// voxels->GetPosition()(2),
// voxels->GetSpacing()(0),
// voxels->GetSpacing()(1),
// voxels->GetSpacing()(2),
// nx, ny, nz
// );
// Vector3i index(0,0,0);
// for (int zv = 0; zv < nz; ++zv) {
// for (int yv = 0; yv < ny; ++yv) {
// for (int xv = 0; xv < nx; ++xv) {
// index << xv,yv,zv;
// // write voxel density in mrad^2/cm //
// float density = fabs(voxels->GetValue(index)) * norm;
// fprintf(vtk_file, "%f ", density);
// }
// }
// }
// fprintf(vtk_file,
// "\n </DataArray>\n </CellData>\n");
fprintf(vtk_file,
"</Piece>\n </ImageData>\n </VTKFile>\n");
fclose(vtk_file);
printf("%s vtk file saved\n",file);
}
// fprintf(vtk_file,
// "</Piece>\n </ImageData>\n </VTKFile>\n");
// fclose(vtk_file);
// printf("%s vtk file saved\n",file);
//}

52
src/Math/VoxImage.h
View File

@@ -45,28 +45,28 @@ namespace uLib {
// ABSTRACT VOX IMAGE //////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
namespace Abstract {
class VoxImage : public uLib::ImageData {
public:
typedef uLib::ImageData BaseClass;
//namespace Abstract {
//class VoxImage : public uLib::ImageData {
//public:
// typedef uLib::ImageData BaseClass;
virtual float GetValue(const Vector3i id) const = 0;
virtual float GetValue(const int id) const = 0;
virtual void SetValue(const Vector3i id, float value) = 0;
virtual void SetValue(const int id, float value) = 0;
// virtual float GetValue(const Vector3i id) const = 0;
// virtual float GetValue(const int id) const = 0;
// virtual void SetValue(const Vector3i id, float value) = 0;
// virtual void SetValue(const int id, float value) = 0;
virtual void SetDims(const Vector3i &size) = 0;
// virtual void SetDims(const Vector3i &size) = 0;
void ExportToVtk(const char *file, bool density_type = 0);
void ExportToVtkXml(const char *file, bool density_type = 0);
int ImportFromVtk(const char *file);
// void ExportToVtk(const char *file, bool density_type = 0);
// void ExportToVtkXml(const char *file, bool density_type = 0);
// int ImportFromVtk(const char *file);
protected:
//protected:
virtual ~VoxImage() {}
VoxImage(const Vector3i &size) : BaseClass(size) {}
};
}
// virtual ~VoxImage() {}
// VoxImage(const Vector3i &size) : BaseClass(size) {}
//};
//}
////////////////////////////////////////////////////////////////////////////////
// VOXEL ////////////////////////////////////////////////////////////////////
@@ -92,9 +92,9 @@ std::basic_ostream<T> & operator << (std::basic_ostream<T> &o, const Voxel &v) {
template< typename T >
class VoxImage : public Abstract::VoxImage {
class VoxImage : public ImageData {
public:
typedef Abstract::VoxImage BaseClass;
typedef ImageData BaseClass;
VoxImage();
@@ -120,14 +120,16 @@ public:
inline T& operator[](const Vector3i &id) { return m_Data[Map(id)]; }
// this implements Abstract interface //
inline Scalarf GetValue(const Vector3i id) const { return this->At(id)(); }
inline Scalarf GetValue(const int id) const { return this->At(id)(); }
inline void SetValue(const Vector3i id, Scalarf value) { this->operator [](id)() = value; }
inline void SetValue(const int id, float value) { this->operator [](id)() = value; }
inline void * GetDataPointer(const Id_t id) const { return (void *)&this->At(id); }
inline Scalarf GetValue(const Vector3i id) const { return this->At(id).Value; }
inline Scalarf GetValue(const int id) const { return this->At(id).Value; }
inline void SetValue(const Vector3i id, Scalarf value) { this->operator [](id).Value = value; }
inline void SetValue(const int id, float value) { this->operator [](id).Value = value; }
// fix needed
inline void SetDims(const Vector3i &size) {
this->m_Data.resize(size.prod());
StructuredData::SetDims(size);
ImageMap::SetDims(size);
}
inline VoxImage<T> clipImage(const Vector3i begin, const Vector3i end) const;
@@ -213,13 +215,11 @@ template<typename T>
VoxImage<T>::VoxImage() :
m_Data(0),
BaseClass(Vector3i(0,0,0)) {}
//{ Interface::IsA <T,Interface::Voxel>(); /* structural check for T */ }
template<typename T>
VoxImage<T>::VoxImage(const Vector3i &size) :
m_Data(size.prod()),
BaseClass(size) {}
//{ Interface::IsA <T,Interface::Voxel>(); /* structural check for T */ }
template <typename T>

4
src/Math/VoxImageFilter.h
View File

@@ -55,7 +55,7 @@ class VoxImageFilter {
public:
virtual void Run() = 0;
virtual void SetImage(Abstract::VoxImage *image) = 0;
virtual void SetImage(ImageData *image) = 0;
protected:
virtual ~VoxImageFilter() {}
@@ -88,7 +88,7 @@ public:
uLibGetMacro(Image,VoxImage<VoxelT> *)
void SetImage(Abstract::VoxImage *image);
void SetImage(ImageData *image);
protected:

10
src/Math/VoxImageFilter.hpp
View File

@@ -39,8 +39,8 @@ namespace uLib {
// KERNEL //////////////////////////////////////////////////////////////////////
template < typename T >
class Kernel : public StructuredData {
typedef StructuredData BaseClass;
class Kernel : public ImageMap {
typedef ImageMap BaseClass;
public:
Kernel(const Vector3i &size);
@@ -63,7 +63,7 @@ Kernel<T>::Kernel(const Vector3i &size) :
BaseClass(size),
m_Data(size.prod())
{
Interface::IsA<T,Interface::Voxel>();
// Interface::IsA<T,Interface::Voxel>();
}
template < typename T >
@@ -166,7 +166,7 @@ _TPL_
void VoxImageFilter<_TPLT_>::SetKernelNumericXZY(const Vector<float> &numeric)
{
// set data order //
StructuredData::Order order = m_KernelData.GetDataOrder();
ImageMap::Order order = m_KernelData.GetDataOrder();
//m_KernelData.SetDataOrder(StructuredData::XZY);
Vector3i id;
int index = 0;
@@ -258,7 +258,7 @@ void VoxImageFilter<_TPLT_>::SetKernelWeightFunction(ShapeT shape)
_TPL_
void VoxImageFilter<_TPLT_>::SetImage(Abstract::VoxImage *image)
void VoxImageFilter<_TPLT_>::SetImage(ImageData *image)
{
this->m_Image = reinterpret_cast<VoxImage<VoxelT> *> (image);
this->SetKernelOffset();

1
src/Math/testing/CMakeLists.txt
View File

@@ -3,6 +3,7 @@ set(TESTS
MathVectorTest
GeometryTest
ContainerBoxTest
ImageDataTest
VoxImageTest
VoxRaytracerTest
StructuredDataTest

78
src/Math/testing/ImageDataTest.cpp Normal file
View File

@@ -0,0 +1,78 @@
/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
All rights reserved
Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
------------------------------------------------------------------
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3.0 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library.
//////////////////////////////////////////////////////////////////////////////*/
#include "testing-prototype.h"
#include "Math/ImageData.h"
using namespace uLib;
namespace {
struct MyVoxel {
MyVoxel() : value(0), count(0) {}
float value;
int count;
};
struct VoxelMean : public MyVoxel {
VoxelMean() {}
void SetValue(const float value) { this->value += value; ++count; }
float GetValue() const { return value/count; }
};
}
int main() {
DataVector<MyVoxel> v;
DataVectorImage img;
img.Data() = v;
img.Scalars().SetAccessFunctions(&MyVoxel::value);
img.SetDims(Vector3i(3,3,3));
for (int x=0; x<img.GetDims().prod(); ++x){
img.SetValue(x,x);
std::cout << img.UnMap(x).transpose() << " -> " << img.GetValue(x) << "\n";
}
DataVector<VoxelMean> vm;
img.Data() = vm;
img.SetDims(Vector3i(3,3,3));
img.Scalars().SetAccessFunctions(&VoxelMean::GetValue,&VoxelMean::SetValue);
for (int x=0; x<img.GetDims().prod(); ++x){
img.SetValue(x,x);
img.SetValue(x,x+2);
img.SetValue(x,x-1);
std::cout << img.UnMap(x).transpose() << " -> " << img.GetValue(x) << "\n";
}
}