44 #ifndef _INCLUDED_Field3D_SparseField_H_
45 #define _INCLUDED_Field3D_SparseField_H_
51 #include <boost/thread/mutex.hpp>
52 #include <boost/lexical_cast.hpp>
57 #define BLOCK_ORDER 4 // 2^BLOCK_ORDER is the block size along each axis
69 template <
class Field_T>
71 template <
class Field_T>
85 template <
typename Data_T>
93 typedef boost::intrusive_ptr<LinearSparseFieldInterp>
Ptr;
102 return "LinearSparseFieldInterp";
117 FIELD3D_VEC3_T<double> p(vsP - FIELD3D_VEC3_T<double>(0.5));
120 V3i c1(static_cast<int>(
floor(p.x)),
121 static_cast<int>(
floor(p.y)),
122 static_cast<int>(
floor(p.z)));
126 FIELD3D_VEC3_T<double> f1(
static_cast<FIELD3D_VEC3_T<double>
>(c2) - p);
128 FIELD3D_VEC3_T<double> f2(
static_cast<FIELD3D_VEC3_T<double>
>(1.0) - f1);
141 int i = c1.x, j = c1.y, k = c1.z, vi, vj, vk, bi, bj, bk;
148 if (vi < blockSize - 1 && vj < blockSize - 1 && vk < blockSize - 1) {
151 const int blockId = field.
blockId(bi, bj, bk);
158 const Data_T *
const p = field.
blockData(bi, bj, bk);
159 const Data_T *
const c111 =
160 p + vi + vj * blockSize + vk * blockSize * blockSize;
161 const Data_T *
const c121 = c111 + blockSize * (c2.y - c1.y);
163 c112 = c111 + blockSize * blockSize * (c2.z - c1.z);
164 const Data_T *
const c122 = c112 + blockSize * (c2.y - c1.y);
165 int xInc = c2.x - c1.x;
166 Data_T value = static_cast<Data_T>
167 (f1.x * (f1.y * (f1.z * *c111 +
169 f2.y * (f1.z * *c121 +
171 f2.x * (f1.y * (f1.z * *(c111 + xInc) +
172 f2.z * *(c112 + xInc)) +
173 f2.y * (f1.z * *(c121 + xInc) +
174 f2.z * *(c122 + xInc))));
185 return static_cast<Data_T>
186 (f1.x * (f1.y * (f1.z * field.
fastValue(c1.x, c1.y, c1.z) +
187 f2.z * field.
fastValue(c1.x, c1.y, c2.z)) +
188 f2.y * (f1.z * field.
fastValue(c1.x, c2.y, c1.z) +
189 f2.z * field.
fastValue(c1.x, c2.y, c2.z))) +
190 f2.x * (f1.y * (f1.z * field.
fastValue(c2.x, c1.y, c1.z) +
191 f2.z * field.
fastValue(c2.x, c1.y, c2.z)) +
192 f2.y * (f1.z * field.
fastValue(c2.x, c2.y, c1.z) +
193 f2.z * field.
fastValue(c2.x, c2.y, c2.z))));
226 template <
typename Data_T>
249 inline Data_T&
value(
int i,
int j,
int k,
int blockOrder)
251 {
return data[(k << blockOrder << blockOrder) + (j << blockOrder) + i]; }
255 inline const Data_T&
value(
int i,
int j,
int k,
int blockOrder)
const
256 {
return data[(k << blockOrder << blockOrder) + (j << blockOrder) + i]; }
267 data =
new Data_T[n];
349 template <
class Data_T>
357 typedef boost::intrusive_ptr<SparseField>
Ptr;
358 typedef std::vector<Ptr>
Vec;
370 return "SparseField";
438 template <
typename Functor_T>
442 int blockId(
int blockI,
int blockJ,
int blockK)
const;
447 void getBlockCoord(
int i,
int j,
int k,
int &bi,
int &bj,
int &bk)
const;
452 void getVoxelInBlock(
int i,
int j,
int k,
int &vi,
int &vj,
int &vk)
const;
487 virtual Data_T
value(
int i,
int j,
int k)
const;
488 virtual long long int memSize()
const;
496 virtual Data_T&
lvalue(
int i,
int j,
int k);
502 Data_T
fastValue(
int i,
int j,
int k)
const;
508 Data_T*
blockData(
int bi,
int bj,
int bk)
const;
528 class const_iterator;
531 const_iterator
cbegin()
const;
535 const_iterator
cend()
const;
538 const_iterator
cend(
const Box3i &subset)
const;
558 class block_iterator;
570 void addReference(
const std::string &filename,
const std::string &layerPath,
571 int valuesPerBlock,
int numVoxels,
int occupiedBlocks);
655 template <
typename Data_T>
675 template <
typename Data_T>
683 Box3i dovsBounds = dvsBounds;
693 f->
getBlockCoord(dovsBounds.min.x, dovsBounds.min.y, dovsBounds.min.z,
694 dbsBounds.min.x, dbsBounds.min.y, dbsBounds.min.z);
695 f->
getBlockCoord(dovsBounds.max.x, dovsBounds.max.y, dovsBounds.max.z,
696 dbsBounds.max.x, dbsBounds.max.y, dbsBounds.max.z);
709 template <
typename Data_T>
721 const V3i &validSize,
const V3i &blockSize)
724 Data_T first = block.
data[0];
727 size_t len = blockSize.x * blockSize.y * blockSize.z;
728 if (validSize == blockSize) {
730 for (
size_t i = 0; i < len; i++) {
731 if (block.
data[i] != first) {
739 for (
size_t i = 0; i < len; i++, x++) {
740 if (x >= blockSize.x) {
743 if (y >= blockSize.y) {
748 if (x >= validSize.x || y >= validSize.y || z >= validSize.z) {
751 if (block.
data[i] != first) {
759 retEmptyValue = first;
769 template <
typename Data_T>
770 inline bool isAnyLess(
const Data_T &left,
const Data_T &right)
772 return (std::abs(left) < right);
780 return (std::abs(left.x) < right.x ||
781 std::abs(left.y) < right.y ||
782 std::abs(left.z) < right.z );
790 return (std::abs(left.x) < right.x ||
791 std::abs(left.y) < right.y ||
792 std::abs(left.z) < right.z );
800 return (std::abs(left.x) < right.x ||
801 std::abs(left.y) < right.y ||
802 std::abs(left.z) < right.z );
810 template <
typename Data_T>
826 const V3i &validSize,
const V3i &blockSize)
829 Data_T first = block.
data[0];
831 bool allGreater =
true;
832 size_t len = blockSize.x * blockSize.y * blockSize.z;
834 if (validSize == blockSize) {
836 for (
size_t i = 0; i < len; i++) {
845 for (
size_t i = 0; i < len; i++, x++) {
846 if (x >= blockSize.x) {
849 if (y >= blockSize.y) {
854 if (x >= validSize.x || y >= validSize.y || z >= validSize.z) {
865 retEmptyValue = first;
884 template <
class Data_T>
888 #if defined(WIN32) || __MAC_OS_X_VERSION_MIN_REQUIRED >= 1090
889 typedef std::forward_iterator_tag iterator_category;
891 typedef ptrdiff_t difference_type;
892 typedef ptrdiff_t distance_type;
893 typedef Data_T *pointer;
894 typedef Data_T& reference;
901 : x(currentPos.x), y(currentPos.y), z(currentPos.z),
902 m_p(NULL), m_blockIsActivated(false),
904 m_blockId(-1), m_window(window), m_field(&field)
906 m_manager = m_field->m_fileManager;
907 setupNextBlock(x, y, z);
910 if (m_manager && m_blockId >= 0 &&
911 m_blockId < static_cast<int>(m_field->m_numBlocks)) {
912 if (m_field->m_blocks[m_blockId].isAllocated)
913 m_manager->decBlockRef<Data_T>(m_field->m_fileId, m_blockId);
918 bool resetPtr =
false;
920 if (x == m_window.max.x) {
921 if (y == m_window.max.y) {
935 ++m_blockStepsTicker;
937 if (!m_isEmptyBlock && (!m_manager || m_blockIsActivated))
944 m_blockStepsTicker = 0;
945 setupNextBlock(x, y, z);
949 template <
class Iter_T>
950 inline bool operator == (
const Iter_T &rhs)
const
952 return x == rhs.x && y == rhs.y && z == rhs.z;
954 template <
class Iter_T>
955 inline bool operator != (
const Iter_T &rhs)
const
957 return x != rhs.x || y != rhs.y || z != rhs.z;
961 if (!m_isEmptyBlock && m_manager && !m_blockIsActivated) {
962 m_manager->activateBlock<Data_T>(m_field->m_fileId, m_blockId);
963 m_blockIsActivated =
true;
964 const Block &block = m_field->m_blocks[m_blockId];
966 m_field->getVoxelInBlock(x, y, z, vi, vj, vk);
971 inline const Data_T* operator -> ()
const
973 if (!m_isEmptyBlock && m_manager && !m_blockIsActivated) {
975 manager->
activateBlock<Data_T>(m_field->m_fileId, m_blockId);
976 m_blockIsActivated =
true;
977 const Block &block = m_field->m_blocks[m_blockId];
979 m_field->getVoxelInBlock(x, y, z, vi, vj, vk);
1000 m_field->applyDataWindowOffset(i, j, k);
1001 m_field->getBlockCoord(i, j, k, m_blockI, m_blockJ, m_blockK);
1002 int oldBlockId = m_blockId;
1003 m_blockId = m_field->blockId(m_blockI, m_blockJ, m_blockK);
1004 if (m_manager && oldBlockId != m_blockId &&
1006 oldBlockId < static_cast<int>(m_field->m_numBlocks) &&
1007 m_field->m_blocks[oldBlockId].isAllocated) {
1008 m_manager->decBlockRef<Data_T>(m_field->m_fileId, oldBlockId);
1010 if (m_blockId >= m_field->m_blockXYSize * m_field->m_blockRes.z) {
1011 m_isEmptyBlock =
true;
1015 const Block &block = m_field->m_blocks[m_blockId];
1017 m_field->getVoxelInBlock(i, j, k, vi, vj, vk);
1018 m_blockStepsTicker = vi;
1020 if (m_manager && oldBlockId != m_blockId && m_blockId >= 0) {
1021 m_manager->incBlockRef<Data_T>(m_field->m_fileId, m_blockId);
1030 m_isEmptyBlock =
false;
1033 m_isEmptyBlock =
true;
1035 if (m_field->m_fileManager) {
1036 m_blockIsActivated =
false;
1066 template <
class Data_T>
1070 #if defined(WIN32) || __MAC_OS_X_VERSION_MIN_REQUIRED >= 1090
1071 typedef std::forward_iterator_tag iterator_category;
1073 typedef ptrdiff_t difference_type;
1074 typedef ptrdiff_t distance_type;
1075 typedef Data_T *pointer;
1076 typedef Data_T& reference;
1081 const Box3i &window,
1083 : x(currentPos.x), y(currentPos.y), z(currentPos.z),
1085 m_blockId(-1), m_window(window), m_field(&field)
1087 setupNextBlock(x, y, z);
1091 bool resetPtr =
false;
1093 if (x == m_window.max.x) {
1094 if (y == m_window.max.y) {
1108 ++m_blockStepsTicker;
1110 if (!m_isEmptyBlock)
1117 m_blockStepsTicker = 0;
1118 setupNextBlock(x, y, z);
1124 return x == rhs.
x && y == rhs.
y && z == rhs.
z;
1128 return x != rhs.
x || y != rhs.
y || z != rhs.
z;
1132 if (m_field->m_fileManager) {
1133 assert(
false &&
"Dereferencing iterator on a dynamic-read sparse field");
1139 if (m_isEmptyBlock) {
1141 m_field->lvalue(x, y, z);
1143 setupNextBlock(x, y, z);
1147 inline Data_T* operator -> ()
1149 if (m_field->m_fileManager) {
1150 assert(
false &&
"Dereferencing iterator on a dynamic-read sparse field");
1156 if (m_isEmptyBlock) {
1158 m_field->lvalue(x, y, z);
1160 setupNextBlock(x, y, z);
1171 m_field->applyDataWindowOffset(i, j, k);
1172 m_field->getBlockCoord(i, j, k, m_blockI, m_blockJ, m_blockK);
1173 m_blockId = m_field->blockId(m_blockI, m_blockJ, m_blockK);
1174 if (m_blockId >= m_field->m_blockXYSize * m_field->m_blockRes.z) {
1175 m_isEmptyBlock =
true;
1178 Block &block = m_field->m_blocks[m_blockId];
1180 m_field->getVoxelInBlock(i, j, k, vi, vj, vk);
1181 m_blockStepsTicker = vi;
1184 m_isEmptyBlock =
false;
1187 m_isEmptyBlock =
true;
1212 template <
class Data_T>
1220 const V3i ¤tPos)
1221 : x(currentPos.x), y(currentPos.y), z(currentPos.z),
1222 m_window(window), m_field(field)
1224 recomputeBlockBoundingBox();
1229 if (x == m_window.max.x) {
1230 if (y == m_window.max.y) {
1241 recomputeBlockBoundingBox();
1247 return x == rhs.
x && y == rhs.
y && z == rhs.
z;
1252 return x != rhs.
x || y != rhs.
y || z != rhs.
z;
1257 return m_currentBlockWindow;
1272 m_currentBlockWindow = box;
1286 template <
class Data_T>
1298 template <
class Data_T>
1301 m_blockOrder(o.m_blockOrder),
1303 m_fileManager(o.m_fileManager)
1310 template <
class Data_T>
1313 if (m_fileManager) {
1317 m_fileManager->removeFieldFromCache<Data_T>(m_fileId);
1326 template <
class Data_T>
1331 this->base::operator=(o);
1339 template <
class Data_T>
1349 m_fileManager->reference<Data_T>(o.
m_fileId);
1355 setupReferenceBlocks();
1364 m_blocks =
new Block[m_numBlocks];
1365 for (
size_t i = 0; i < m_numBlocks; ++i) {
1369 1 << m_blockOrder << m_blockOrder << m_blockOrder);
1372 m_fileManager = NULL;
1378 template <
class Data_T>
1380 const std::string &layerPath,
1386 m_fileId = m_fileManager->getNextId<Data_T>(filename, layerPath);
1389 m_fileManager->reference<Data_T>(m_fileId);
1391 reference->numVoxels = numVoxels;
1392 reference->occupiedBlocks = occupiedBlocks;
1393 reference->setNumBlocks(m_numBlocks);
1398 template <
class Data_T>
1403 for (
size_t i = 0; i < m_numBlocks; ++i) {
1406 m_blocks[i].clear();
1412 template <
class Data_T>
1415 if (!m_fileManager || m_fileId < 0)
return;
1418 m_fileManager->reference<Data_T>(m_fileId);
1420 #if F3D_NO_BLOCKS_ARRAY
1422 reference->
blocks = m_blocks;
1423 int nextBlockIdx = 0;
1424 for (
size_t i = 0; i < m_numBlocks; ++i, ++fb) {
1425 if (m_blocks[i].isAllocated) {
1435 reference->
blocks.begin();
1436 int nextBlockIdx = 0;
1437 for (
size_t i = 0; i < m_numBlocks; ++i, ++fb, ++bp) {
1438 if (m_blocks[i].isAllocated) {
1451 template <
class Data_T>
1456 Block *p = m_blocks, *end = m_blocks + m_numBlocks;
1465 template <
class Data_T>
1468 m_blockOrder = order;
1474 template <
class Data_T>
1477 return m_blockOrder;
1482 template <
class Data_T>
1485 return 1 << m_blockOrder;
1490 template <
class Data_T>
1494 applyDataWindowOffset(i, j, k);
1495 getBlockCoord(i, j, k, bi, bj, bk);
1496 return blockIsAllocated(bi, bj, bk);
1501 template <
class Data_T>
1504 const Block &block = m_blocks[blockId(bi, bj, bk)];
1510 template <
class Data_T>
1513 return m_blocks[blockId(bi, bj, bk)].emptyValue;
1518 template <
class Data_T>
1522 Block &block = m_blocks[blockId(bi, bj, bk)];
1524 deallocBlock(block, val);
1532 template <
class Data_T>
1535 return bi >= 0 && bj >= 0 && bk >= 0 &&
1536 bi < m_blockRes.x && bj < m_blockRes.y && bk < m_blockRes.z;
1541 template <
class Data_T>
1549 template <
class Data_T>
1550 template <
typename Functor_T>
1554 int numDeallocs = 0;
1562 V3i blockAllocSize(blockSize());
1564 int bx = 0, by = 0, bz = 0;
1565 for (
size_t i = 0; i < m_numBlocks; ++i, ++bx) {
1566 if (bx >= m_blockRes.x) {
1569 if (by >= m_blockRes.y) {
1574 validSize = blockAllocSize;
1575 if (bx == m_blockRes.x-1) {
1576 validSize.x = dataRes.x - bx * blockAllocSize.x;
1578 if (by == m_blockRes.y-1) {
1579 validSize.y = dataRes.y - by * blockAllocSize.y;
1581 if (bz == m_blockRes.z-1) {
1582 validSize.z = dataRes.z - bz * blockAllocSize.z;
1585 if (m_blocks[i].isAllocated) {
1586 if (func.check(m_blocks[i], emptyValue, validSize, blockAllocSize)) {
1587 deallocBlock(m_blocks[i], emptyValue);
1597 template <
class Data_T>
1600 return fastValue(i, j, k);
1605 template <
class Data_T>
1608 return fastLValue(i, j, k);
1613 template <
class Data_T>
1616 assert (i >= base::m_dataWindow.
min.x);
1617 assert (i <= base::m_dataWindow.
max.x);
1618 assert (j >= base::m_dataWindow.
min.y);
1619 assert (j <= base::m_dataWindow.
max.y);
1620 assert (k >= base::m_dataWindow.
min.z);
1621 assert (k <= base::m_dataWindow.
max.z);
1623 applyDataWindowOffset(i, j, k);
1626 getBlockCoord(i, j, k, bi, bj, bk);
1629 getVoxelInBlock(i, j, k, vi, vj, vk);
1631 int id = blockId(bi, bj, bk);
1632 const Block &block = m_blocks[id];
1635 if (m_fileManager) {
1636 m_fileManager->incBlockRef<Data_T>(m_fileId, id);
1637 m_fileManager->activateBlock<Data_T>(m_fileId, id);
1638 Data_T tmpValue = block.
value(vi, vj, vk, m_blockOrder);
1639 m_fileManager->decBlockRef<Data_T>(m_fileId, id);
1642 return block.
value(vi, vj, vk, m_blockOrder);
1652 template <
class Data_T>
1655 assert (i >= base::m_dataWindow.
min.x);
1656 assert (i <= base::m_dataWindow.
max.x);
1657 assert (j >= base::m_dataWindow.
min.y);
1658 assert (j <= base::m_dataWindow.
max.y);
1659 assert (k >= base::m_dataWindow.
min.z);
1660 assert (k <= base::m_dataWindow.
max.z);
1662 if (m_fileManager) {
1663 assert(
false &&
"Called fastLValue() on a dynamic-read sparse field");
1670 applyDataWindowOffset(i, j, k);
1673 getBlockCoord(i, j, k, bi, bj, bk);
1676 getVoxelInBlock(i, j, k, vi, vj, vk);
1678 int id = blockId(bi, bj, bk);
1679 Block &block = m_blocks[id];
1682 return block.
value(vi, vj, vk, m_blockOrder);
1685 size_t blockSize = 1 << m_blockOrder << m_blockOrder << m_blockOrder;
1687 return block.
value(vi, vj, vk, m_blockOrder);
1693 template <
class Data_T>
1696 int id = blockId(bi, bj, bk);
1697 const Block &block = m_blocks[id];
1707 template <
class Data_T>
1710 long long int blockSize = m_numBlocks *
sizeof(
Block);
1711 long long int dataSize = 0;
1713 for (
size_t i = 0; i < m_numBlocks; ++i) {
1714 if (m_blocks[i].data) {
1715 dataSize += (1 << m_blockOrder << m_blockOrder << m_blockOrder) *
1720 return sizeof(*this) + dataSize + blockSize;
1725 template <
class Data_T>
1729 const size_t blockSize = (1 << m_blockOrder << m_blockOrder << m_blockOrder);
1731 for (
size_t i = 0; i < m_numBlocks; ++i) {
1732 if (m_blocks[i].isAllocated) {
1742 template <
class Data_T>
1754 template <
class Data_T>
1758 if (subset.isEmpty())
1759 return cend(subset);
1765 template <
class Data_T>
1770 V3i(base::m_dataWindow.
min.x,
1771 base::m_dataWindow.min.y,
1772 base::m_dataWindow.max.z + 1),
1778 template <
class Data_T>
1785 subset.max.z + 1), m_blockOrder);
1790 template <
class Data_T>
1796 return iterator(*
this, base::m_dataWindow,
1797 base::m_dataWindow.
min, m_blockOrder); }
1801 template <
class Data_T>
1805 if (subset.isEmpty())
1807 return iterator(*
this, subset, subset.min, m_blockOrder);
1812 template <
class Data_T>
1816 return iterator(*
this, base::m_dataWindow,
1817 V3i(base::m_dataWindow.
min.x,
1818 base::m_dataWindow.min.y,
1819 base::m_dataWindow.max.z + 1), m_blockOrder);
1824 template <
class Data_T>
1829 V3i(subset.min.x, subset.min.y, subset.max.z + 1),
1835 template <
class Data_T>
1847 template <
class Data_T>
1852 V3i(0, 0, m_blockRes.z));
1857 template <
class Data_T>
1861 V3f res(base::m_dataWindow.size() +
V3i(1));
1862 V3f blockRes(res / (1 << m_blockOrder));
1863 blockRes.x =
ceil(blockRes.x);
1864 blockRes.y =
ceil(blockRes.y);
1865 blockRes.z =
ceil(blockRes.z);
1866 V3i intBlockRes(static_cast<int>(blockRes.x),
1867 static_cast<int>(blockRes.y),
1868 static_cast<int>(blockRes.z));
1869 m_blockRes = intBlockRes;
1870 m_blockXYSize = m_blockRes.x * m_blockRes.y;
1874 m_numBlocks = intBlockRes.x * intBlockRes.y * intBlockRes.z;
1875 m_blocks =
new Block[m_numBlocks];
1880 template <
class Data_T>
1883 return blockK * m_blockXYSize + blockJ * m_blockRes.x + blockI;
1889 template <
class Data_T>
1891 int &bi,
int &bj,
int &bk)
const
1896 bi = i >> m_blockOrder;
1897 bj = j >> m_blockOrder;
1898 bk = k >> m_blockOrder;
1904 template <
class Data_T>
1906 int &vi,
int &vj,
int &vk)
const
1911 vi = i & ((1 << m_blockOrder) - 1);
1912 vj = j & ((1 << m_blockOrder) - 1);
1913 vk = k & ((1 << m_blockOrder) - 1);
1918 template <
class Data_T>
1921 m_fileManager->incBlockRef<Data_T>(m_fileId, blockId);
1926 template <
class Data_T>
1929 m_fileManager->activateBlock<Data_T>(m_fileId, blockId);
1934 template <
class Data_T>
1937 m_fileManager->decBlockRef<Data_T>(m_fileId, blockId);
1942 template <
class Data_T>
1950 template <
class Data_T>
1954 const size_t blockSide = (1 << m_blockOrder);
1958 const V3i start(bCoord * blockSide + base::m_dataWindow.
min);
1961 const Box3i unclipped(start, end);
1962 bounds =
clipBounds(unclipped, base::m_dataWindow);
1964 return bounds == unclipped;
1969 template <
class Data_T>
1984 #endif // Include guard