Libs/Optimize/Domain/ImageDomain.h

Namespaces

Name
shapeworks User usage reporting (telemetry)

Classes

	Name
class	shapeworks::ImageDomain

Source code

```cpp

pragma once

include

include "ParticleRegionDomain.h"

// we have to undef foreach here because both Qt and OpenVDB define foreach

undef foreach

ifndef Q_MOC_RUN

include

endif

namespace shapeworks { template class ImageDomain : public ParticleRegionDomain { public: using Pointer = std::shared_ptr;

using ImageType = itk::Image;

using PointType = ParticleRegionDomain::PointType;

void SetImage(ImageType* I, double narrow_band) { this->m_FixedDomain = false; // this->Modified();

openvdb::initialize();  // It is safe to initialize multiple times.

// Set a large background value, so that we quickly catch particles outside or on the edge the narrow band.
// (Downside: its more difficult to display the correct location of the point of failure.)
m_VDBImage = openvdb::FloatGrid::create(1e8);
m_VDBImage->setGridClass(openvdb::GRID_LEVEL_SET);
auto vdbAccessor = m_VDBImage->getAccessor();

// Save properties of the Image needed for the optimizer
m_Size = I->GetRequestedRegion().GetSize();
m_Spacing = I->GetSpacing();
m_Origin = I->GetOrigin();
m_Index = I->GetRequestedRegion().GetIndex();

// Transformation from index space to world space
openvdb::math::Mat4f mat;
mat.setIdentity();
mat.postScale(openvdb::Vec3f(m_Spacing[0], m_Spacing[1], m_Spacing[2]));
mat.postTranslate(openvdb::Vec3f(m_Origin[0], m_Origin[1], m_Origin[2]));
const auto xform = openvdb::math::Transform::createLinearTransform(mat);
m_VDBImage->setTransform(xform);

itk::ImageRegionIterator<ImageType> it(I, I->GetRequestedRegion());
it.GoToBegin();

while (!it.IsAtEnd()) {
  const auto idx = it.GetIndex();
  const auto pixel = it.Get();
  if (abs(pixel) > narrow_band) {
    ++it;
    continue;
  }
  const auto coord = openvdb::Coord(idx[0], idx[1], idx[2]);
  vdbAccessor.setValue(coord, pixel);
  ++it;
}

typename ImageType::PointType l0;
I->TransformIndexToPhysicalPoint(m_Index, l0);
for (unsigned int i = 0; i < DIMENSION; i++) m_Index[i] += m_Size[i] - 1;

typename ImageType::PointType u0;
I->TransformIndexToPhysicalPoint(m_Index, u0);

// Cast points to higher precision if needed.  Parent class uses doubles
// because they are compared directly with points in the particle system,
// which are always double precision.
PointType l;
PointType u;

for (unsigned int i = 0; i < DIMENSION; i++) {
  l[i] = static_cast<double>(l0[i]);
  u[i] = static_cast<double>(u0[i]);
}

this->SetLowerBound(l);
this->SetUpperBound(u);

// Precompute and save values that are used in parts of the optimizer
this->SetupImageForCrossingPointUpdate(I);  // this->UpdateZeroCrossingPoint(I);
this->UpdateSurfaceArea(I);

}

double GetSurfaceArea() const override { return m_SurfaceArea; }

inline PointType GetOrigin() const { return m_Origin; }

inline typename ImageType::SizeType GetSize() const { return m_Size; }

inline typename ImageType::SpacingType GetSpacing() const { return m_Spacing; }

inline typename ImageType::RegionType::IndexType GetIndex() const { return m_Index; }

inline PointType GetValidLocationNear(PointType p) const override { // todo why is this function ignoring the argument? Also see Optimize::AddSinglePoint return m_ZeroCrossingPoint; }

inline T Sample(const PointType& p) const { if (this->IsInsideBuffer(p)) { const auto coord = this->ToVDBCoord(p); return openvdb::tools::BoxSampler::sample(m_VDBImage->tree(), coord); } else { std::ostringstream message; message << "Domain " << m_DomainID << ": " << m_DomainName << " : Distance transform queried for a Point, " << p << ", outside the given image domain. Consider increasing the padding in grooming or the narrow band " "optimization parameter"; throw std::runtime_error(message.str()); } }

inline double GetMaxDiameter() const override { double bestRadius = 0; double maxdim = 0; for (unsigned int i = 0; i < ImageType::ImageDimension; i++) { if (GetSize()[i] > maxdim) { maxdim = GetSize()[i]; bestRadius = maxdim * GetSpacing()[i]; } } return bestRadius; }

void DeleteImages() override { m_VDBImage = 0; }

// Updates zero crossing points. Raster scans candidate zero crossing points, and finds one that does not violate any // constraints. void UpdateZeroCrossingPoint() override { for (size_t i = 0; i < m_possible_zero_crossings.size(); i++) { this->m_ZeroCrossingPoint = m_possible_zero_crossings[i]; if (!this->GetConstraints()->isAnyViolated(this->m_ZeroCrossingPoint)) { // std::cout << "Chosen initial point " << this->m_ZeroCrossingPoint << std::endl; break; } }

if (this->GetConstraints()->isAnyViolated(this->m_ZeroCrossingPoint)) {
  std::cerr << "A particle initialization violates at least one constraint. Make sure at least one point satisfies "
               "all constraints"
            << std::endl;
}

}

protected: openvdb::FloatGrid::Ptr GetVDBImage() const { return m_VDBImage; }

ImageDomain() {} virtual ~ImageDomain() {};

void PrintSelf(std::ostream& os, itk::Indent indent) const { ParticleRegionDomain::PrintSelf(os, indent); os << indent << "VDB Active Voxels = " << m_VDBImage->activeVoxelCount() << std::endl; }

inline openvdb::math::Transform::Ptr transform() const { return this->m_VDBImage->transformPtr(); }

// Converts a coordinate from an ITK Image point in world space to the corresponding // coordinate in OpenVDB Index space. Raises an exception if the narrow band is not // sufficiently large to sample the point. inline openvdb::Vec3R ToVDBCoord(const PointType& p) const { const auto worldCoord = openvdb::Vec3R(p[0], p[1], p[2]); const auto idxCoord = this->transform()->worldToIndex(worldCoord);

// Make sure the coordinate is part of the narrow band
if (m_VDBImage->tree().isValueOff(
        openvdb::Coord::round(idxCoord))) {  // `isValueOff` requires an integer coordinate
  // If multiple threads crash here at the same time, the error message displayed is just "terminate called
  // recursively", which isn't helpful. So we std::cerr the error to make sure its printed to the console.
  std::cerr << "Sampled point outside the narrow band: " << p << std::endl;

  std::ostringstream message;
  message << "Attempt to sample at a point outside the narrow band: " << p
          << ". Consider increasing the narrow band";
  throw std::runtime_error(message.str());
}

return idxCoord;

}

private: openvdb::FloatGrid::Ptr m_VDBImage; typename ImageType::SizeType m_Size; typename ImageType::SpacingType m_Spacing; PointType m_Origin; PointType m_ZeroCrossingPoint; typename ImageType::RegionType::IndexType m_Index; // Index defining the corner of the region double m_SurfaceArea; std::vector m_possible_zero_crossings;

// Computes possible zero crossing points. Later on, one can find the ones that do not violate constraints. void SetupImageForCrossingPointUpdate(ImageType* I) { typename itk::ZeroCrossingImageFilter::Pointer zc = itk::ZeroCrossingImageFilter::New(); zc->SetInput(I); zc->Update(); typename itk::ImageRegionConstIteratorWithIndex zcIt(zc->GetOutput(), zc->GetOutput()->GetRequestedRegion());

for (zcIt.GoToReverseBegin(); !zcIt.IsAtReverseEnd(); --zcIt) {
  if (zcIt.Get() == 1.0) {
    PointType pos;
    I->TransformIndexToPhysicalPoint(zcIt.GetIndex(), pos);
    this->m_ZeroCrossingPoint = pos;
    m_possible_zero_crossings.push_back(pos);
  }
}

}

void UpdateSurfaceArea(ImageType* I) { Image image(I); m_SurfaceArea = image.toMesh(0).getSurfaceArea(); } };

} // end namespace shapeworks ```

Updated on 2026-03-31 at 16:02:11 +0000