libinsighttoolkit4-dev 4.5.0-3 / usr / include / ITK-4.5 / itkBioCellularAggregate.hxx

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  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.txt
 *
 *  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.
 *
 *=========================================================================*/
#ifndef __itkBioCellularAggregate_hxx
#define __itkBioCellularAggregate_hxx

#include "itkBioCellularAggregate.h"

namespace itk
{
namespace bio
{
template< unsigned int NSpaceDimension >
CellularAggregate< NSpaceDimension >
::CellularAggregate()
{
  typename BioCellType::ColorType color;
  color.SetRed(1.0);
  color.SetGreen(1.0);
  color.SetBlue(1.0);
  BioCellType::SetDefaultColor(color);

  m_Mesh = MeshType::New();

  m_Mesh->SetCellsAllocationMethod(MeshType::CellsAllocatedDynamicallyCellByCell);
  m_Mesh->SetPoints( PointsContainer::New() );
  m_Mesh->SetPointData( PointDataContainer::New() );
  m_Mesh->SetCells( VoronoiRegionsContainer::New() );

  m_Iteration = 0;
  m_ClosestPointComputationInterval = 5;

  m_FrictionForce = 1.0f;
}

template< unsigned int NSpaceDimension >
CellularAggregate< NSpaceDimension >
::~CellularAggregate()
{
  this->KillAll();
}

template< unsigned int NSpaceDimension >
unsigned int
CellularAggregate< NSpaceDimension >
::GetNumberOfCells(void) const
{
  return m_Mesh->GetPointData()->Size();
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::SetGrowthRadiusLimit(double value)
{
  BioCellType::SetGrowthRadiusLimit(value);
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::SetGrowthRadiusIncrement(double value)
{
  BioCellType::SetGrowthRadiusIncrement(value);
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::PrintSelf(std::ostream & os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << "Cellular aggregate " << m_Mesh << std::endl;
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::Remove(CellBase *cellbase)
{
  BioCellType *cell = dynamic_cast< BioCellType * >( cellbase );

  if ( !cell )
    {
    itk::ExceptionObject exception;
    exception.SetDescription("Trying to remove a null pointer to cell");
    exception.SetLocation("CellularAggregate::Remove(BioCellType*)");
    throw exception;
    }

  IdentifierType id = cell->GetSelfIdentifier();

  typename MeshType::CellAutoPointer region;
  bool regionExist = m_Mesh->GetCell(id, region);
  if ( regionExist )
    {
    VoronoiRegionType *realRegion =
      dynamic_cast< VoronoiRegionType * >( region.GetPointer() );

    if ( !realRegion )
      {
      itkExceptionMacro("CellularAggregate::Remove() couldn't dynamic_cast region ");
      }
    else
      {
      //
      // Notify all the neighbors that this cell is going away
      //
      typename VoronoiRegionType::PointIdIterator neighbor = realRegion->PointIdsBegin();
      typename VoronoiRegionType::PointIdIterator end      = realRegion->PointIdsEnd();
      while ( neighbor != end )
        {
        const IdentifierType neighborId = ( *neighbor );
        typename MeshType::CellAutoPointer cellPointer;
        bool neighborVoronoiExist = m_Mesh->GetCell(neighborId, cellPointer);
        if ( neighborVoronoiExist )
          {
          VoronoiRegionType *vregion =
            dynamic_cast< VoronoiRegionType * >( cellPointer.GetPointer() );
          if ( !vregion )
            {
            std::cerr << "CellularAggregate::Add() Failed to find a region"  << std::endl;
            }
          else
            {
            vregion->RemovePointId(id);
            }
          }
        neighbor++;
        }

      // We can now remove the entry from the list of cells in the Mesh
      m_Mesh->GetCells()->DeleteIndex(id);

      // and then release the memory used by the VoronoiRegion
      delete realRegion;
      }
    }
  else
    {
    itkExceptionMacro(" Region " << id << " doesn't exist ");
    }

  m_Mesh->GetPoints()->DeleteIndex(id);
  m_Mesh->GetPointData()->DeleteIndex(id);

  // Finally we can delete the BioCell;
  delete cell;
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::GetVoronoi(IdentifierType cellId, VoronoiRegionAutoPointer & voronoiPointer) const
{
  typename MeshType::CellAutoPointer cellPointer;

  bool voronoiExists = m_Mesh->GetCell(cellId, cellPointer);
  if ( !voronoiExists )
    {
    itk::ExceptionObject exception;
    exception.SetDescription("voronoi region does not exist in the container");
    exception.SetLocation("GetVoronoi( IdentifierType");
    throw exception;
    }

  VoronoiRegionType *region =
    dynamic_cast< VoronoiRegionType * >( cellPointer.GetPointer() );

  voronoiPointer.TakeNoOwnership(region);
  if ( cellPointer.IsOwner() )
    {
    voronoiPointer.TakeOwnership();
    cellPointer.ReleaseOwnership();
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::SetEgg(BioCellType *cell, const PointType & position)
{
  VectorType perturbation = position.GetVectorFromOrigin();

  this->Add(cell, perturbation);
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::Add(CellBase *cell)
{
  VectorType perturbation;

  perturbation.Fill(0.0);
  this->Add(cell, perturbation);
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::Add(CellBase *cellA, CellBase *cellB, double perturbationLength)
{
  // Create a perturbation for separating the daugther cells
  typename BioCellType::VectorType perturbationVector;
  for ( unsigned int d = 0; d < NSpaceDimension; d++ )
    {
    perturbationVector[d] = vnl_sample_uniform(-1.0f, 1.0f);
    }

  const double norm = perturbationVector.GetNorm();
  if ( vnl_math_abs(norm) > 1e-10 )
    {
    perturbationVector *= perturbationLength / norm;
    }
  else
    {
    // this event should rarely happen... very rarely
    perturbationVector[0] = perturbationLength;
    }

  this->Add(cellA,  perturbationVector);
  this->Add(cellB, -perturbationVector);

  const IdentifierType cellAId = cellA->GetSelfIdentifier();
  const IdentifierType cellBId = cellB->GetSelfIdentifier();

  VoronoiRegionAutoPointer cellAptr;
  VoronoiRegionAutoPointer cellBptr;

  this->GetVoronoi(cellAId, cellAptr);
  this->GetVoronoi(cellBId, cellBptr);

  cellAptr->AddPointId(cellBId);
  cellBptr->AddPointId(cellAId);
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::Add(CellBase *cellBase, const VectorType & perturbation)
{
  BioCellType *  cell = dynamic_cast< BioCellType * >( cellBase );
  IdentifierType newcellId       = cell->GetSelfIdentifier();
  IdentifierType newcellparentId = cell->GetParentIdentifier();

  PointType       position;
  CellAutoPointer selfVoronoi;

  // If the cell does not have a parent
  // from which receive a position
  if ( !newcellparentId )
    {
    position.Fill(0.0);
    selfVoronoi.TakeOwnership(new VoronoiRegionType);
    }
  else
    {
    bool parentPositionExist = m_Mesh->GetPoint(newcellparentId, &position);
    if ( !parentPositionExist )
      {
      itk::ExceptionObject exception;
      exception.SetDescription("Parent cell does not exist in the container");
      exception.SetLocation("CellularAggregate::Add( cell * )");
      throw exception;
      }

    VoronoiRegionAutoPointer parentVoronoi;
    this->GetVoronoi(newcellparentId, parentVoronoi);
    parentVoronoi->MakeCopy(selfVoronoi);
    }

  position += perturbation;

  m_Mesh->SetCell(newcellId, selfVoronoi);
  m_Mesh->SetPoint(newcellId, position);
  m_Mesh->SetPointData(newcellId, cell);

  cell->SetCellularAggregate(this);

  // Add this new cell as neighbor to cells in its neighborhood

  typename VoronoiRegionType::PointIdIterator neighbor = selfVoronoi->PointIdsBegin();
  typename VoronoiRegionType::PointIdIterator end      = selfVoronoi->PointIdsEnd();

  while ( neighbor != end )
    {
    const IdentifierType neighborId = ( *neighbor );
    typename MeshType::CellAutoPointer cellPointer;
    bool neighborVoronoiExist = m_Mesh->GetCell(neighborId, cellPointer);
    if ( neighborVoronoiExist )
      {
      VoronoiRegionType *region =
        dynamic_cast< VoronoiRegionType * >( cellPointer.GetPointer() );

      if ( !region )
        {
        std::cerr << "CellularAggregate::Add() Failed to find a region"  << std::endl;
        }
      else
        {
        region->AddPointId(newcellId);
        }
      }
    neighbor++;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::AdvanceTimeStep(void)
{
  if ( m_Iteration % m_ClosestPointComputationInterval == 0 )
    {
    this->ComputeClosestPoints();
    }

  this->ComputeForces();
  this->UpdatePositions();

  CellsIterator begin = m_Mesh->GetPointData()->Begin();
  CellsIterator end   = m_Mesh->GetPointData()->End();

  CellsIterator cell  = begin;

  while ( cell != end )
    {
    BioCellType *theCell = cell.Value();
    theCell->AdvanceTimeStep();
    ++cell;
    if ( theCell->MarkedForRemoval() )
      {
      this->Remove(theCell);
      }
    }

  this->InvokeEvent( IterationEvent() );

  m_Iteration++;
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::KillAll(void)
{
  if ( !m_Mesh  )
    {
    return;
    }

  CellsIterator cell = m_Mesh->GetPointData()->Begin();
  CellsIterator end  = m_Mesh->GetPointData()->End();

  while ( cell != end )
    {
    delete ( cell.Value() );
    ++cell;
    }

  BioCellType::ResetCounter();
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::ClearForces(void)
{
  CellsIterator cell = m_Mesh->GetPointData()->Begin();
  CellsIterator end  = m_Mesh->GetPointData()->End();

  while ( cell != end )
    {
    cell.Value()->ClearForce();
    ++cell;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::UpdatePositions(void)
{
  CellsConstIterator cellIt = m_Mesh->GetPointData()->Begin();
  CellsConstIterator end    = m_Mesh->GetPointData()->End();

  PointType position;

  position.Fill(0);

  while ( cellIt != end )
    {
    BioCellType *  cell = cellIt.Value();
    IdentifierType cellId = cell->GetSelfIdentifier();
    m_Mesh->GetPoint(cellId, &position);
    const VectorType force = cell->GetForce();
    if ( force.GetNorm() > m_FrictionForce )
      {
      position += force / 50.0;
      }
    m_Mesh->SetPoint(cellId, position);
    cellIt++;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::ComputeForces(void)
{
  // Clear all the force accumulators
  this->ClearForces();

  CellsConstIterator cell1It     = m_Mesh->GetPointData()->Begin();
  CellsConstIterator end         = m_Mesh->GetPointData()->End();

  // compute forces
  while ( cell1It != end )
    {
    const IdentifierType cell1Id = cell1It.Index();

    BioCellType *cell1      =  cell1It.Value();

    PointType position1;
    position1.Fill(0);
    m_Mesh->GetPoint(cell1Id, &position1);

    const double rA          = cell1->GetRadius();

    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1Id, voronoiRegion);

    typename VoronoiRegionType::PointIdIterator neighbor = voronoiRegion->PointIdsBegin();
    typename VoronoiRegionType::PointIdIterator vend      = voronoiRegion->PointIdsEnd();

    while ( neighbor != vend )
      {
      const IdentifierType cell2Id = ( *neighbor );

      BioCellType *cell2 = 0;
      PointType    position2;

      if ( !m_Mesh->GetPoint(cell2Id, &position2) )
        {
        neighbor++;  // if the neigbor has been removed, skip it
        continue;
        }
      m_Mesh->GetPointData(cell2Id, &cell2);

      const double rB      = cell2->GetRadius();

      typename BioCellType::VectorType relativePosition = position1 - position2;

      const double distance = relativePosition.GetNorm();

      if ( distance < ( rA + rB ) / 2.0 )
        {
        const double factor = 2.0 * BioCellType::GetGrowthRadiusLimit() / distance;
        typename BioCellType::VectorType force = relativePosition * factor;
        cell1->AddForce(force);
        cell2->AddForce(-force);
        }
      else if ( distance < rA + rB )
        {
        typename BioCellType::VectorType force = relativePosition;
        cell1->AddForce(force);
        cell2->AddForce(-force);
        }

      neighbor++;
      }

    cell1It++;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::ComputeClosestPoints(void)
{
  PointsConstIterator beginPoints   = m_Mesh->GetPoints()->Begin();
  PointsConstIterator endPoints     = m_Mesh->GetPoints()->End();

  PointsConstIterator point1It = beginPoints;

  while ( point1It != endPoints )
    {
    PointType position1 = point1It.Value();

    PointsConstIterator point2It   = beginPoints;

    BioCellType *cell1 = 0;

    IdentifierType cell1Id = point1It.Index();
    m_Mesh->GetPointData(cell1Id, &cell1);
    const double radius  = cell1->GetRadius();
    const double limitDistance = radius * 4.0;

    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1Id, voronoiRegion);

    voronoiRegion->ClearPoints();

    while ( point2It != endPoints )
      {
      if ( point2It == point1It )
        {
        point2It++;
        continue;
        }

      PointType position2 = point2It.Value();

      typename BioCellType::VectorType relativePosition = position1 - position2;

      const double distance = relativePosition.GetNorm();
      if ( distance < limitDistance )
        {
        voronoiRegion->AddPointId( point2It.Index() );
        }
      point2It++;
      }

    point1It++;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::DumpContent(std::ostream & os) const
{
  CellsConstIterator beginCell   = m_Mesh->GetPointData()->Begin();
  CellsConstIterator endCell     = m_Mesh->GetPointData()->End();

  CellsConstIterator cell1It = beginCell;

  os << "Cell Identifiers " << std::endl;
  while ( cell1It != endCell )
    {
    os << cell1It.Index() << " == ";
    os << cell1It.Value()->GetSelfIdentifier() << std::endl;
    cell1It++;
    }

  os << std::endl << "Points " << std::endl;
  PointsConstIterator pointIt     = m_Mesh->GetPoints()->Begin();
  PointsConstIterator endPoint    = m_Mesh->GetPoints()->End();

  while ( pointIt != endPoint )
    {
    os << pointIt.Index() << "  ";
    os << pointIt.Value() << std::endl;
    pointIt++;
    }

  os << std::endl << "Neighborhoods " << std::endl;
  cell1It = beginCell;
  while ( cell1It != endCell )
    {
    VoronoiRegionAutoPointer voronoiRegion;
    this->GetVoronoi(cell1It.Index(), voronoiRegion);
    typename VoronoiRegionType::PointIdIterator neighbor = voronoiRegion->PointIdsBegin();
    typename VoronoiRegionType::PointIdIterator end      = voronoiRegion->PointIdsEnd();

    while ( neighbor != end )
      {
      os << ( *neighbor ) << "   ";
      neighbor++;
      }
    os << std::endl;

    cell1It++;
    }
}

template< unsigned int NSpaceDimension >
void
CellularAggregate< NSpaceDimension >
::AddSubstrate(SubstrateType *substrate)
{
  SubstratePointer smartPointer(substrate);

  m_Substrates.push_back(smartPointer);
}

template< unsigned int NSpaceDimension >
typename CellularAggregate< NSpaceDimension >::SubstratesVector &
CellularAggregate< NSpaceDimension >
::GetSubstrates(void)
{
  return m_Substrates;
}

template< unsigned int NSpaceDimension >
typename CellularAggregate< NSpaceDimension >::SubstrateValueType
CellularAggregate< NSpaceDimension >
::GetSubstrateValue(IdentifierType cellId, unsigned int substrateId) const
{
  PointType cellPosition;
  bool      cellPositionExists = m_Mesh->GetPoint(cellId, &cellPosition);

  if ( !cellPositionExists )
    {
    std::cerr << " Cell position doesn't exist for cell Id = ";
    std::cerr << cellId << std::endl;
    return itk::NumericTraits< SubstrateValueType >::Zero;
    }

  SubstratePointer substrate = m_Substrates[substrateId];

  typename SubstrateType::IndexType index;
  typedef typename SubstrateType::IndexType::IndexValueType IndexValueType;

  substrate->TransformPhysicalPointToIndex(cellPosition, index);

  SubstrateValueType value = 0;
  if ( substrate->GetBufferedRegion().IsInside(index) )
    {
    value = substrate->GetPixel(index);
    }

  return value;
}
} // end namespace bio
} // end namespace itk

#endif