Skip to content

Libs/Optimize/ParticleSystem.h

Namespaces

Name
shapeworks
User usage reporting (telemetry)

Classes

Name
class shapeworks::ParticleSystem
A facade class managing interactions with a particle system.

Source code

```cpp

pragma once

include

include

include

include "Libs/Optimize/Container/GenericContainer.h"

include "Libs/Optimize/Domain/ParticleDomain.h"

include "Libs/Optimize/Neighborhood/ParticleNeighborhood.h"

include "Observer.h"

include "ParticleEvents.h"

include "itkCommand.h"

include "itkDataObject.h"

include "itkEventObject.h"

include "itkObjectFactory.h"

include "itkPoint.h"

include "itkWeakPointer.h"

include "vnl/vnl_inverse.h"

include "vnl/vnl_matrix_fixed.h"

include "vnl/vnl_vector_fixed.h"

namespace shapeworks { class ParticleSystem : public itk::DataObject { public: static constexpr int VDimension = 3;

typedef ParticleSystem Self; typedef DataObject Superclass; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; typedef itk::WeakPointer ConstWeakPointer;

itkNewMacro(Self);

itkTypeMacro(ParticleSystem, itk::DataObject);

itkStaticConstMacro(Dimension, unsigned int, VDimension);

using DomainType = shapeworks::ParticleDomain;

typedef itk::Point PointType;

typedef GenericContainer PointContainerType;

// typedef Transform TransformType; typedef vnl_matrix_fixed TransformType; typedef vnl_vector_fixed VectorType; typedef vnl_matrix VnlMatrixType;

void RegisterObserver(Observer *attr);

void SynchronizePositions() { for (unsigned int d = 0; d < this->GetNumberOfDomains(); d++) { for (unsigned int p = 0; p < this->GetNumberOfParticles(d); p++) { this->SetPosition(this->GetPosition(p, d), p, d); } } }

unsigned long int GetNumberOfParticles(unsigned int d = 0) const { return m_Positions[d]->GetSize(); }

const PointType &AddPosition(const PointType &, unsigned int d = 0); const PointType &SetPosition(const PointType &, unsigned long int k, unsigned int d = 0);

void RemovePosition(unsigned long int k, unsigned int d = 0);

PointType &GetPosition(unsigned long int k, unsigned int d = 0) { return m_Positions[d]->operator; } const PointType &GetPosition(unsigned long int k, unsigned int d = 0) const { return m_Positions[d]->operator; } PointType GetTransformedPosition(unsigned long int k, unsigned int d = 0) const { return this->TransformPoint(m_Positions[d]->operator, m_Transforms[d] * m_PrefixTransforms[d]); } PointType GetPrefixTransformedPosition(unsigned long int k, unsigned int d = 0) const { return this->TransformPoint(m_Positions[d]->operator, m_PrefixTransforms[d]); }

void SplitAllParticles(double epsilon); void SplitParticle(double epsilon, unsigned int idx, unsigned int d = 0); void AdvancedAllParticleSplitting(double epsilon, unsigned int domains_per_shape, unsigned int dom_to_process); // Debug function void PrintParticleSystem();

std::shared_ptr GetNeighborhood(unsigned int k) const { return m_Neighborhoods[k]; }

using PointVectorType = std::vector;

PointVectorType FindNeighborhoodPoints(const PointType &p, int idx, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(p, idx, r); } inline PointVectorType FindNeighborhoodPoints(const PointType &p, int idx, std::vector &w, std::vector &distances, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(p, idx, w, distances, r); } inline PointVectorType FindNeighborhoodPoints(const PointType &p, int idx, std::vector &w, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(p, idx, w, r); } inline PointVectorType FindNeighborhoodPoints(unsigned int idx, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(GetPosition(idx, d), idx, r); } inline PointVectorType FindNeighborhoodPoints(unsigned int idx, std::vector &w, std::vector &distances, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(GetPosition(idx, d), idx, w, distances, r); } inline PointVectorType FindNeighborhoodPoints(unsigned int idx, std::vector &w, double r, unsigned int d = 0) const { return m_Neighborhoods[d]->find_neighborhood_points(GetPosition(idx, d), idx, w, r); }

// inline int FindNeighborhoodPoints(const PointType &p, double r, PointVectorType &vec, unsigned int d = 0) const // { return m_Neighborhoods[d]->FindNeighborhoodPoints(p, r, vec); }

// PointVectorType FindTransformedNeighborhoodPoints(const PointType &p, double r, unsigned int d = 0) const // { // PointVectorType ans = m_Neighborhoods[d] // ->FindNeighborhoodPoints(this->TransformPoint(p, InverseTransform[d]), r); // for (unsigned int i = 0; i < ans.size(); i++) // { // ans.Point[i] = this->TransformPoint(ans.Point[i], m_Transform[d]); // } // return ans; // }

void AddDomain(DomainType::Pointer input);

std::vector::const_iterator GetDomainsBegin() const { return m_Domains.begin(); }

std::vector::const_iterator GetDomainsEnd() const { return m_Domains.end(); }

DomainType *GetDomain(unsigned int i) { return m_Domains[i].get(); }

DomainType *GetDomain() { return m_Domains[0].get(); }

const DomainType *GetDomain(unsigned int i) const { return m_Domains[i].get(); }

const DomainType *GetDomain() const { return m_Domains[0].get(); }

unsigned int GetNumberOfDomains() const { return m_Domains.size(); }

void SetTransform(unsigned int i, const TransformType &); void SetTransform(const TransformType &p) { this->SetTransform(0, p); } void SetPrefixTransform(unsigned int i, const TransformType &); void SetPrefixTransform(const TransformType &p) { this->SetPrefixTransform(0, p); }

std::vector::const_iterator GetTransformsBegin() const { return m_Transforms.begin(); }

std::vector::const_iterator GetTransformsEnd() const { return m_Transforms.end(); }

const TransformType &GetTransform(unsigned int i) const { return m_Transforms[i]; }

const TransformType &GetTransform() const { return m_Transforms[0]; }

TransformType GetTransform(unsigned int i) { return m_Transforms[i]; }

TransformType GetTransform() { return m_Transforms[0]; }

const TransformType &GetPrefixTransform(unsigned int i) const { return m_PrefixTransforms[i]; }

const TransformType &GetPrefixTransform() const { return m_PrefixTransforms[0]; }

TransformType GetPrefixTransform(unsigned int i) { return m_PrefixTransforms[i]; }

TransformType GetPrefixTransform() { return m_PrefixTransforms[0]; }

std::vector::const_iterator GetInverseTransformsBegin() const { return m_InverseTransforms.begin(); }

std::vector::const_iterator GetInverseTransformsEnd() const { return m_InverseTransforms.end(); }

const TransformType &GetInverseTransform(unsigned int i) const { return m_InverseTransforms[i]; }

const TransformType &GetInverseTransform() const { return m_InverseTransforms[0]; }

const TransformType &GetInversePrefixTransform(unsigned int i) const { return m_InversePrefixTransforms[i]; }

const TransformType &GetInversePrefixTransform() const { return m_InversePrefixTransforms[0]; }

const std::vector &GetPositions() const { return m_Positions; } const PointContainerType::Pointer &GetPositions(unsigned int d) const { return m_Positions[d]; }

void AddPositionList(const std::vector &, unsigned int d = 0);

PointType TransformPoint(const PointType &, const TransformType &) const;

VectorType TransformVector(const VectorType &, const TransformType &) const;

VnlMatrixType TransformNormalDerivative(const VnlMatrixType &, const TransformType &) const;

inline TransformType InvertTransform(const TransformType &T) const { // Note, vnl_inverse is optimized for small matrices 1x1 - 4x4 return vnl_inverse(T); }

void FlagDomain(unsigned int i) { // ensure large enough while (i >= this->m_DomainFlags.size()) { m_DomainFlags.push_back(false); }

// set the flag
m_DomainFlags[i] = true;

} void UnflagDomain(unsigned int i) { m_DomainFlags[i] = false; } bool GetDomainFlag(unsigned int i) const { if (i >= m_DomainFlags.size()) { // not set return false; } return m_DomainFlags[i]; } const std::vector &GetDomainFlags() const { return m_DomainFlags; } void SetDomainFlags() { for (unsigned int i = 0; i < m_DomainFlags.size(); i++) { m_DomainFlags[i] = true; } } void ResetDomainFlags() { for (unsigned int i = 0; i < m_DomainFlags.size(); i++) { m_DomainFlags[i] = false; } }

std::vector GetFixedShapes(bool &has_fixed) const { unsigned int num_shapes = GetNumberOfDomains() / m_DomainsPerShape; std::vector is_fixed(num_shapes, false); has_fixed = false; for (unsigned int shape = 0; shape < num_shapes; shape++) { bool shape_is_fixed = true; for (unsigned int d = 0; d < m_DomainsPerShape; d++) { if (!GetDomainFlag(shape * m_DomainsPerShape + d)) { shape_is_fixed = false; break; } } if (shape_is_fixed && m_DomainsPerShape > 0) { is_fixed[shape] = true; has_fixed = true; } } return is_fixed; }

void SetFixedParticleFlag(unsigned int d, unsigned int i) { m_FixedParticleFlags[d][i] = true; } void ResetFixedParticleFlag(unsigned int d, unsigned int i) { m_FixedParticleFlags[d][i] = false; } bool GetFixedParticleFlag(unsigned int d, unsigned int i) const { return m_FixedParticleFlags[d][i]; } void ResetFixedParticleFlags() { for (unsigned d = 0; d < m_FixedParticleFlags.size(); d++) { for (unsigned int i = 0; i < m_FixedParticleFlags[d].size(); i++) m_FixedParticleFlags[d][i] = false; } }

void SetDomainsPerShape(unsigned int num) { m_DomainsPerShape = num; m_FixedParticleFlags.resize(m_DomainsPerShape); } unsigned int GetDomainsPerShape() const { return m_DomainsPerShape; }

void SetNumberOfDomains(unsigned int);

// Returns the maximum distance between nearest neighbors in domain dom double ComputeMaxDistNearestNeighbors(size_t dom);

void SetFieldAttributes(const std::vector &field_attributes) { m_FieldAttributes = field_attributes; }

const std::vector &GetFieldAttributes() const { return m_FieldAttributes; }

protected: ParticleSystem(); void PrintSelf(std::ostream &os, itk::Indent indent) const; virtual ~ParticleSystem(){};

TransformType &GetInverseTransform(unsigned int i) { return m_InverseTransforms[i]; }

TransformType &GetInverseTransform() { return m_InverseTransforms[0]; }

TransformType &GetInversePrefixTransform(unsigned int i) { return m_InversePrefixTransforms[i]; }

TransformType &GetInversePrefixTransform() { return m_InversePrefixTransforms[0]; }

private: ParticleSystem(const Self &); // purposely not implemented void operator=(const Self &); // purposely not implemented

std::vector m_Positions;

std::vector m_Domains;

unsigned int m_DomainsPerShape;

std::vector> m_Neighborhoods;

std::vector m_Transforms;

std::vector m_InverseTransforms;

std::vector m_PrefixTransforms;

std::vector m_InversePrefixTransforms;

std::vector m_IndexCounters;

std::vector m_DomainFlags;

std::vector> m_FixedParticleFlags;

std::vector m_FieldAttributes;

std::mt19937 m_rand{42}; };

} // end namespace shapeworks ```


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