Commit 615e3dd9 authored by Jørgen Lind's avatar Jørgen Lind
Browse files

Add VolumeIndexer and CoordinateTransformer

and add functionality for adding noice to data pages. Also add a test
for the InMemory IOManager that uses the nocie functionality
parent ee22e835
......@@ -29,6 +29,7 @@ set(SOURCE_FILES
VDS/FSE/fse_decompress.cpp
VDS/Rle.cpp
VDS/VolumeDataRequestProcessor.cpp
VDS/VolumeIndexer.cpp
)
set (PRIVATE_HEADER_FILES
......@@ -70,6 +71,7 @@ set (PRIVATE_HEADER_FILES
VDS/Rle.h
VDS/VolumeDataRequestProcessor.h
VDS/ThreadPool.h
VDS/SimplexNoiceKernel.h
)
set (EXPORTED_HEADER_FILES
......@@ -88,6 +90,8 @@ set (EXPORTED_HEADER_FILES
OpenVDS/VolumeDataLayout.h
OpenVDS/VolumeDataLayoutDescriptor.h
OpenVDS/VolumeSampler.h
OpenVDS/CoordinateTransformer.h
OpenVDS/VolumeIndexer.h
)
add_library(openvds_objects OBJECT
......
#ifndef COORDINATETRANSFORMER_H
#define COORDINATETRANSFORMER_H
#include <OpenVDS/Vector.h>
#include <OpenVDS/VolumeData.h>
#include <cstring>
#include <cmath>
namespace OpenVDS
{
struct M4
{
DoubleVector4 data[4];
};
static void fastInvert(M4 &m)
{
DoubleVector3 cT(m.data[3].X, m.data[3].Y, m.data[3].Z);
double fA0 = m.data[1].Y*m.data[2].Z-m.data[2].Y*m.data[1].Z;
double fA1 = m.data[0].Y*m.data[2].Z-m.data[2].Y*m.data[0].Z;
double fA2 = m.data[0].Y*m.data[1].Z-m.data[1].Y*m.data[0].Z;
double fA3 = m.data[1].X*m.data[2].Z-m.data[2].X*m.data[1].Z;
double fA4 = m.data[0].X*m.data[2].Z-m.data[2].X*m.data[0].Z;
double fA5 = m.data[0].X*m.data[1].Z-m.data[1].X*m.data[0].Z;
double fA6 = m.data[1].X*m.data[2].Y-m.data[2].X*m.data[1].Y;
double fA7 = m.data[0].X*m.data[2].Y-m.data[2].X*m.data[0].Y;
double fA8 = m.data[0].X*m.data[1].Y-m.data[1].X*m.data[0].Y;
double fDet = m.data[0].X*fA0 - m.data[0].Y*fA3 + m.data[0].Z*fA6;
const double rEpsilon = 1e-60;
if (fabs(fDet) <= rEpsilon)
{
memset(&m, 0, sizeof(m));
return;
}
double fInvDet = ((double)1.0) / fDet;
m.data[0].X = fA0*fInvDet;
m.data[0].Y = -fA1*fInvDet;
m.data[0].Z = fA2*fInvDet;
m.data[1].X = -fA3*fInvDet;
m.data[1].Y = fA4*fInvDet;
m.data[1].Z = -fA5*fInvDet;
m.data[2].X = fA6*fInvDet;
m.data[2].Y = -fA7*fInvDet;
m.data[2].Z = fA8*fInvDet;
m.data[3].X = -cT.X*m.data[0].X-cT.Y*m.data[1].X-cT.Z*m.data[2].X;
m.data[3].Y = -cT.X*m.data[0].Y-cT.Y*m.data[1].Y-cT.Z*m.data[2].Y;
m.data[3].Z = -cT.X*m.data[0].Z-cT.Y*m.data[1].Z-cT.Z*m.data[2].Z;
}
struct IJKGridDefinition
{
IJKGridDefinition()
: origin(0,0,0)
, iUnitStep(1,0,0)
, jUnitStep(0,1,0)
, kUnitStep(0,0,1)
{}
IJKGridDefinition(DoubleVector3 origin, DoubleVector3 iUnitStep, DoubleVector3 jUnitStep, DoubleVector3 kUnitStep)
: origin(origin)
, iUnitStep(iUnitStep)
, jUnitStep(jUnitStep)
, kUnitStep(kUnitStep)
{}
DoubleVector3 origin;
DoubleVector3 iUnitStep;
DoubleVector3 jUnitStep;
DoubleVector3 kUnitStep;
};
struct VDSIJKGridDefinition : public IJKGridDefinition
{
VDSIJKGridDefinition()
: IJKGridDefinition()
, dimensionMap(2,1,0)
{}
VDSIJKGridDefinition(const IJKGridDefinition &ijkGridDefinition, IntVector3 dimensionMap)
: IJKGridDefinition(ijkGridDefinition)
, dimensionMap(dimensionMap)
{}
IntVector3 dimensionMap;
};
template<int N>
struct VDSCoordinateTransformerBase
{
float ijkToWorldTransform[3][3]; ///< Coordinate transfrom matrix
float worldToIJKTransform[3][3]; ///< Coordinate transform matrix
float ijkToWorldTranslation[3]; ///< Coordinate translation vector
float worldToIJKTranslation[3]; ///< Coordinate translation vector
int ijkDimensionMap[3]; ///< Map to determine which dimension map to I, J, and K
///
/// Converts a voxel index to world coordinates using the indexer's IJK grid definition and IJK dimension map
/// @param iVoxelIndex the voxel index to convert
/// @return the world coordinates
///
FloatVector<3> VoxelIndexToWorldCoordinates(const IntVector<N> &iVoxelIndex) const
{
//remap volume index to ijk
int ijk[3];
ijk[0] = iVoxelIndex[ijkDimensionMap[0]];
ijk[1] = iVoxelIndex[ijkDimensionMap[1]];
if (ijkDimensionMap[2] >= 0)
ijk[2] = iVoxelIndex[ijkDimensionMap[2]];
else
ijk[2] = 0;
FloatVector<3> rWorldCoords;
rWorldCoords[0] = ijkToWorldTranslation[0] +
ijk[0] * ijkToWorldTransform[0][0] +
ijk[1] * ijkToWorldTransform[1][0] +
ijk[2] * ijkToWorldTransform[2][0];
rWorldCoords[1] = ijkToWorldTranslation[1] +
ijk[0] * ijkToWorldTransform[0][1] +
ijk[1] * ijkToWorldTransform[1][1] +
ijk[2] * ijkToWorldTransform[2][1];
rWorldCoords[2] = ijkToWorldTranslation[2] +
ijk[0] * ijkToWorldTransform[0][2] +
ijk[1] * ijkToWorldTransform[1][2] +
ijk[2] * ijkToWorldTransform[2][2];
return rWorldCoords;
}
///
/// Converts a float voxel index to world coordinates using the indexer's IJK grid definition and IJK dimension map
/// @param rVoxelIndex the float voxel index to convert
/// @return the world coordinates
///
FloatVector<3> VoxelIndexToWorldCoordinates(const FloatVector<N> &rVoxelIndex) const
{
//remap volume index to ijk
float ijk[3];
ijk[0] = rVoxelIndex[ijkDimensionMap[0]];
ijk[1] = rVoxelIndex[ijkDimensionMap[1]];
if (ijkDimensionMap[2] >= 0)
ijk[2] = rVoxelIndex[ijkDimensionMap[2]];
else
ijk[2] = 0;
FloatVector<3> rWorldCoords;
rWorldCoords[0] = ijkToWorldTranslation[0] +
ijk[0] * ijkToWorldTransform[0][0] +
ijk[1] * ijkToWorldTransform[1][0] +
ijk[2] * ijkToWorldTransform[2][0];
rWorldCoords[1] = ijkToWorldTranslation[1] +
ijk[0] * ijkToWorldTransform[0][1] +
ijk[1] * ijkToWorldTransform[1][1] +
ijk[2] * ijkToWorldTransform[2][1];
rWorldCoords[2] = ijkToWorldTranslation[2] +
ijk[0] * ijkToWorldTransform[0][2] +
ijk[1] * ijkToWorldTransform[1][2] +
ijk[2] * ijkToWorldTransform[2][2];
return rWorldCoords;
}
///
/// Converts world coordinates to a voxel index, rounding to the nearest integer
/// @param rWorldCoords the world coordinates to convert
/// @param anVoxelMin
/// @return the voxel index
///
IntVector<N> WorldCoordinatesToVoxelIndex(const FloatVector<3> &rWorldCoords, const int (&anVoxelMin)[Dimensionality_Max]) const
{
int ijk[3];
ijk[0] = (int)((worldToIJKTranslation[0] +
rWorldCoords[0] * worldToIJKTransform[0][0] +
rWorldCoords[1] * worldToIJKTransform[1][0] +
rWorldCoords[2] * worldToIJKTransform[2][0]) + 0.5f);
ijk[1] = (int)((worldToIJKTranslation[1] +
rWorldCoords[0] * worldToIJKTransform[0][1] +
rWorldCoords[1] * worldToIJKTransform[1][1] +
rWorldCoords[2] * worldToIJKTransform[2][1]) + 0.5f);
ijk[2] = (int)((worldToIJKTranslation[2] +
rWorldCoords[0] * worldToIJKTransform[0][2] +
rWorldCoords[1] * worldToIJKTransform[1][2] +
rWorldCoords[2] * worldToIJKTransform[2][2]) + 0.5f);
//remap ijk to volume index
IntVector<N> iVoxelIndex;
for (int i = 0; i < N; i++)
iVoxelIndex[i] = anVoxelMin[i];
iVoxelIndex[ijkDimensionMap[0]] = ijk[0];
iVoxelIndex[ijkDimensionMap[1]] = ijk[1];
if (ijkDimensionMap[2] >= 0)
iVoxelIndex[ijkDimensionMap[2]] = ijk[2];
return iVoxelIndex;
}
///
/// Converts world coordinates to a float voxel index without rounding
/// @param rWorldCoords the world coordinates to convert
/// @param anVoxelMin
/// @return the float voxel index
///
FloatVector<N> WorldCoordinatesToVoxelIndexFloat(const FloatVector<3> &rWorldCoords, const int (&anVoxelMin)[Dimensionality_Max]) const
{
float ijk[3];
ijk[0] = worldToIJKTranslation[0] +
rWorldCoords[0] * worldToIJKTransform[0][0] +
rWorldCoords[1] * worldToIJKTransform[1][0] +
rWorldCoords[2] * worldToIJKTransform[2][0];
ijk[1] = worldToIJKTranslation[1] +
rWorldCoords[0] * worldToIJKTransform[0][1] +
rWorldCoords[1] * worldToIJKTransform[1][1] +
rWorldCoords[2] * worldToIJKTransform[2][1];
ijk[2] = worldToIJKTranslation[2] +
rWorldCoords[0] * worldToIJKTransform[0][2] +
rWorldCoords[1] * worldToIJKTransform[1][2] +
rWorldCoords[2] * worldToIJKTransform[2][2];
//remap ijk to volume index
FloatVector<N> rVoxelIndex;
for (int i = 0; i < N; i++)
rVoxelIndex[i] = (float)anVoxelMin[i];
rVoxelIndex[ijkDimensionMap[0]] = ijk[0];
rVoxelIndex[ijkDimensionMap[1]] = ijk[1];
if (ijkDimensionMap[2] >= 0)
rVoxelIndex[ijkDimensionMap[2]] = ijk[2];
return rVoxelIndex;
}
///////////////////////////// Constructors /////////////////////////////
///
/// Constructor
/// @param cIJKGridAndDImensionMap the IJK grid definition and IJK dimension map
///
VDSCoordinateTransformerBase(const VDSIJKGridDefinition& vdsIjkGridDefinition)
{
SetGridDefinition(vdsIjkGridDefinition);
}
///
/// Constructor
///
VDSCoordinateTransformerBase()
{
VDSIJKGridDefinition unitDefinition(IJKGridDefinition(DoubleVector3(0.0, 0.0, 0.0),
DoubleVector3(1.0, 0.0, 0.0),
DoubleVector3(0.0, 1.0, 0.0),
DoubleVector3(0.0, 0.0, 1.0)),
IntVector3(0, 1, 2));
SetGridDefinition(unitDefinition);
}
private:
///
/// Sets the IJK grid definition and IJK dimension map for use in world coordinate conversions
/// @see VoxelIndexToWorldCoordinates
/// @see WorldCoordinatesToVoxelIndex
/// @see WorldCoordinatesToVoxelIndexFloat
/// @see LocalIndexToWorldCoordinates
/// @see WorldCoordinatesToLocalIndex
/// @param vdsIjkGridDefinition the IJK grid definition and IJK dimension map
///
void SetGridDefinition(const VDSIJKGridDefinition& vdsIjkGridDefinition)
{
ijkDimensionMap[0] = vdsIjkGridDefinition.dimensionMap.X;
ijkDimensionMap[1] = vdsIjkGridDefinition.dimensionMap.Y;
ijkDimensionMap[2] = vdsIjkGridDefinition.dimensionMap.Z;
ijkToWorldTransform[0][0] = (float)vdsIjkGridDefinition.iUnitStep.X;
ijkToWorldTransform[0][1] = (float)vdsIjkGridDefinition.iUnitStep.Y;
ijkToWorldTransform[0][2] = (float)vdsIjkGridDefinition.iUnitStep.Z;
ijkToWorldTransform[1][0] = (float)vdsIjkGridDefinition.jUnitStep.X;
ijkToWorldTransform[1][1] = (float)vdsIjkGridDefinition.jUnitStep.Y;
ijkToWorldTransform[1][2] = (float)vdsIjkGridDefinition.jUnitStep.Z;
ijkToWorldTransform[2][0] = (float)vdsIjkGridDefinition.kUnitStep.X;
ijkToWorldTransform[2][1] = (float)vdsIjkGridDefinition.kUnitStep.Y;
ijkToWorldTransform[2][2] = (float)vdsIjkGridDefinition.kUnitStep.Z;
ijkToWorldTranslation[0] = (float)vdsIjkGridDefinition.origin.X;
ijkToWorldTranslation[1] = (float)vdsIjkGridDefinition.origin.Y;
ijkToWorldTranslation[2] = (float)vdsIjkGridDefinition.origin.Z;
M4 matrix = {
{vdsIjkGridDefinition.iUnitStep.X, vdsIjkGridDefinition.iUnitStep.Y, vdsIjkGridDefinition.iUnitStep.Z, 0},
{vdsIjkGridDefinition.jUnitStep.X, vdsIjkGridDefinition.jUnitStep.Y, vdsIjkGridDefinition.jUnitStep.Z, 0},
{vdsIjkGridDefinition.kUnitStep.X, vdsIjkGridDefinition.kUnitStep.Y, vdsIjkGridDefinition.kUnitStep.Z, 0},
{vdsIjkGridDefinition.origin.X, vdsIjkGridDefinition.origin.Y, vdsIjkGridDefinition.origin.Z, 1}
};
fastInvert(matrix);
worldToIJKTransform[0][0] = (float)matrix.data[0].X;
worldToIJKTransform[0][1] = (float)matrix.data[0].Y;
worldToIJKTransform[0][2] = (float)matrix.data[0].Z;
worldToIJKTransform[1][0] = (float)matrix.data[1].X;
worldToIJKTransform[1][1] = (float)matrix.data[1].Y;
worldToIJKTransform[1][2] = (float)matrix.data[1].Z;
worldToIJKTransform[2][0] = (float)matrix.data[2].X;
worldToIJKTransform[2][1] = (float)matrix.data[2].Y;
worldToIJKTransform[2][2] = (float)matrix.data[2].Z;
worldToIJKTranslation[0] = (float)matrix.data[3].X;
worldToIJKTranslation[1] = (float)matrix.data[3].Y;
worldToIJKTranslation[2] = (float)matrix.data[3].Z;
}
};
typedef VDSCoordinateTransformerBase<2> VDSCoordinateTransformer2D;
typedef VDSCoordinateTransformerBase<3> VDSCoordinateTransformer3D;
typedef VDSCoordinateTransformerBase<4> VDSCoordinateTransformer4D;
typedef VDSCoordinateTransformerBase<5> VDSCoordinateTransformer5D;
typedef VDSCoordinateTransformerBase<6> VDSCoordinateTransformer6D;
}
#endif //COORDINATETRANSFORMER_H
......@@ -21,8 +21,25 @@
namespace OpenVDS
{
template<typename T, size_t N>
struct Vector
{
T data[N];
Vector()
: data{}
{}
template<typename ...Args>
Vector(Args... args)
: data{args...}
{}
inline T &operator[] (size_t n) { return data[n]; }
inline const T &operator[] (size_t n) const { return data[n]; }
};
template<typename T>
struct Vector2
struct Vector<T,2>
{
typedef T element_type;
enum { element_count = 2 };
......@@ -36,15 +53,15 @@ struct Vector2
T data[2];
};
Vector2() : X(), Y() {}
Vector2(T X, T Y) : X(X), Y(Y) {}
Vector() : X(), Y() {}
Vector(T X, T Y) : X(X), Y(Y) {}
inline T &operator[] (size_t n) { return data[n]; }
inline const T &operator[] (size_t n) const { return data[n]; }
};
template<typename T>
struct Vector3
struct Vector<T,3>
{
typedef T element_type;
enum { element_count = 3 };
......@@ -58,15 +75,15 @@ struct Vector3
T data[3];
};
Vector3() : X(), Y(), Z() {}
Vector3(T X, T Y, T Z) : X(X), Y(Y), Z(Z) {}
Vector() : X(), Y(), Z() {}
Vector(T X, T Y, T Z) : X(X), Y(Y), Z(Z) {}
inline T &operator[] (size_t n) { return data[n]; }
inline const T &operator[] (size_t n) const { return data[n]; }
};
template<typename TYPE>
struct Vector4
struct Vector<TYPE, 4>
{
typedef TYPE element_type;
enum { element_count = 4 };
......@@ -80,8 +97,8 @@ struct Vector4
TYPE data[4];
};
Vector4() : X(), Y(), Z(), T() {}
Vector4(TYPE X, TYPE Y, TYPE Z, TYPE T) : X(X), Y(Y), Z(Z), T(T) {}
Vector() : X(), Y(), Z(), T() {}
Vector(TYPE X, TYPE Y, TYPE Z, TYPE T) : X(X), Y(Y), Z(Z), T(T) {}
inline TYPE &operator[] (size_t n) { return data[n]; }
inline const TYPE &operator[] (size_t n) const { return data[n]; }
......@@ -117,17 +134,23 @@ bool operator!=(const Vector<T, N> &a, const Vector<T, N> &b)
return true;
}
*/
using IntVector2 = Vector2<int>;
using IntVector3 = Vector3<int>;
using IntVector4 = Vector4<int>;
using FloatVector2 = Vector2<float>;
using FloatVector3 = Vector3<float>;
using FloatVector4 = Vector4<float>;
using DoubleVector2 = Vector2<double>;
using DoubleVector3 = Vector3<double>;
using DoubleVector4 = Vector4<double>;
using IntVector2 = Vector<int,2>;
using IntVector3 = Vector<int,3>;
using IntVector4 = Vector<int,4>;
using FloatVector2 = Vector<float, 2>;
using FloatVector3 = Vector<float, 3>;
using FloatVector4 = Vector<float, 4>;
using DoubleVector2 = Vector<double, 2>;
using DoubleVector3 = Vector<double, 3>;
using DoubleVector4 = Vector<double, 4>;
template<size_t N>
using IntVector = Vector<int, N>;
template<size_t N>
using FloatVector = Vector<float, N>;
template<size_t N>
using DoubleVector = Vector<double, N>;
}
This diff is collapsed.
/****************************************************************************
** Copyright 2019 The Open Group
** Copyright 2019 Bluware, Inc.
**
** Licensed under the Apache License, Version 2.0 (the "License");
** you may not use this file except in compliance with the License.
** You may obtain a copy of the License at
**
** http://www.apache.org/licenses/LICENSE-2.0
**
** Unless required by applicable law or agreed to in writing, software
** distributed under the License is distributed on an "AS IS" BASIS,
** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
** See the License for the specific language governing permissions and
** limitations under the License.
****************************************************************************/
#include <OpenVDS/Vector.h>
#include <OpenVDS/ValueConversion.h>
#include <OpenVDS/VolumeIndexer.h>
#define NOISE_OCTAVES (5)
#define NOISE_SCALE (0.06f * 0.1f)
namespace OpenVDS
{
template<int dimensionality, int octaves>
float SimplexNoise(float *position, unsigned int randomSeed)
{
int integerPosition[4];
int index[4];
float distance[4];
float sortedDistance[4];
int dimension;
float skew = (sqrtf((float)dimensionality + 1) - 1.f) / (float)dimensionality;
float unskew = (dimensionality + 1 - sqrtf((float)dimensionality + 1)) / ((float)dimensionality * (dimensionality + 1));
float noise = 0.f;
float scale = 1.f;
float amplitude = 1.f;
float amplitudeSum = 0.f;
for(int octave = 0; octave < octaves; octave++)
{
float positionSum = position[0] * scale;
for(dimension = 1; dimension < dimensionality; dimension++)
{
positionSum += position[dimension] * scale;
}
integerPosition[0] = (int)floorf(position[0] * scale + positionSum * skew);
int integerPositionSum = integerPosition[0];
for(dimension = 1; dimension < dimensionality; dimension++)
{
integerPosition[dimension] = (int)floorf(position[dimension] * scale + positionSum * skew);
integerPositionSum += integerPosition[dimension];
}
for(dimension = 0; dimension < dimensionality; dimension++)
{
distance[dimension] = position[dimension] * scale - (integerPosition[dimension] - integerPositionSum * unskew);
sortedDistance[dimension] = distance[dimension];
index[dimension] = dimension;
}
for(dimension = 0; dimension < dimensionality - 1; dimension++)
{
if(sortedDistance[dimension] < sortedDistance[dimension + 1])
{
float rTemp = sortedDistance[dimension];
int iTemp = index[dimension];
sortedDistance[dimension] = sortedDistance[dimension + 1];
index[dimension] = index[dimension + 1];
sortedDistance[dimension + 1] = rTemp;
index[dimension + 1] = iTemp;
if(dimension > 0) dimension-=2;
}
}
for(int vertex = 0; true; vertex++)
{
unsigned int
random = randomSeed,
u0 = 79555561 + octave * 322147;
for(dimension = 0; dimension < dimensionality; dimension++)
{
random += integerPosition[dimension] * random;
u0 += random;
random ^= u0 + 8953453 * random;
}
float weight = 0.6f - distance[0] * distance[0]; // This vertex contributes relative to the distance of the point from this vertex
for(dimension = 1; dimension < dimensionality; dimension++)
{
weight -= distance[dimension] * distance[dimension];
}
if(weight > 0.f)
{
weight *= weight;
weight *= weight;
// We use a gradient vector which is the midpoint of some edge of an N-dimensional cube centered on origo
// To get this, we chose a vertex of the N-dimensional cube, and then chose one of it's coordinates to be 0
// This gives us two ways of making every edge, but that is completely unproblematic
int iNull = random % dimensionality;
random /= dimensionality;
// We loop through each component to form the dot product between the gradient vector and the distance vector
float dot = 0.f;
for(dimension = 0; dimension < dimensionality; dimension++)
{
if(dimension != iNull)
{
dot += distance[dimension] * ((int)(random & 2) - 1);
}
random >>= 1;
}
noise += dot * weight * amplitude;
}