/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkBioCell.txx,v $
  Language:  C++
  Date:      $Date: 2006/03/19 04:36:53 $
  Version:   $Revision: 1.7 $

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkBioCell_txx
#define __itkBioCell_txx

#include "itkBioCell.h"
#include "vnl/vnl_math.h"
#include "vnl/vnl_sample.h"
#include <new>


namespace itk {

namespace bio {


  
/**
 *    Constructor Lonely Cell
 */ 
template<unsigned int NSpaceDimension>
Cell<NSpaceDimension>
::Cell()
{

  m_Force.Fill( 0.0f );

  // Genome pointers are set to NULL in the superclass.
}




/**
 *    Destructor   
 */ 
template<unsigned int NSpaceDimension>
Cell<NSpaceDimension>
::~Cell()
{
  // Genomes are released in the destructor of the superclass.
}



/**
 *    Cell Division
 */ 

template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::Mitosis(void) 
{
  // Create the two daughters.
  Cell * siblingA = new Cell;
  Cell * siblingB = new Cell;

  // Broad compensation for Volume distribution among daugthers.
  // The type of root should depend on the Dimension...
  siblingA->m_Radius   = m_Radius / vcl_sqrt(2.0f );
  siblingB->m_Radius   = m_Radius / vcl_sqrt(2.0f );

  // Update Teleomeres
  siblingA->m_Generation = m_Generation + 1; 
  siblingB->m_Generation = m_Generation + 1;

  // Prepare to separate them by a specific distance.
  // This helps to avoid infinite interaction forces 
  // just after cellular division.
  const double perturbationLength = m_Radius * 0.75;

  // Register the parent
  siblingA->m_ParentIdentifier = m_SelfIdentifier;
  siblingB->m_ParentIdentifier = m_SelfIdentifier;

  // Pass the genome to each daughter cell
  siblingA->m_Genome = m_Genome;
  siblingB->m_Genome = m_GenomeCopy;

  // Mark that the genome pointer is not owned by this cell anymore.
  m_Genome     = NULL;
  m_GenomeCopy = NULL;

  // Register both daughter cells with the CellularAggregate.
  CellularAggregateBase * aggregate = this->GetCellularAggregate();
  aggregate->Add( siblingA, siblingB, perturbationLength );

  // Mark this cell for being removed from the Aggregate and deleted.
  this->MarkForRemoval();

}


/**
 *    Create a New Egg Cell
 *    this method behave like a factory, it is 
 *    intended to be overloaded in any class 
 *    deriving from Cell.
 */ 
template<unsigned int NSpaceDimension>
Cell<NSpaceDimension> *
Cell<NSpaceDimension>
::CreateEgg(void) 
{
  Superclass::Initialize();
  Cell * cell = new Cell;
  cell->m_ParentIdentifier = 0;
  cell->m_SelfIdentifier = 1;
  cell->m_Generation = 0;

  cell->m_Genome = new GenomeType;

  cell->ComputeGeneNetwork();
  cell->SecreteProducts();

  return cell;
}





/**
 *    Clear the cumulator for applied forces
 */ 
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::ClearForce(void) 
{
  m_Force.Fill( 0.0f );
  m_Pressure  = 0.0f;
}


/**
 *    Return the cumulated force
 */ 
template<unsigned int NSpaceDimension>
const typename Cell<NSpaceDimension>::VectorType &
Cell<NSpaceDimension>
::GetForce(void) const
{
  return m_Force;
}



/**
 *    Return a pointer to the Cellular Aggregate
 */ 
template<unsigned int NSpaceDimension>
CellularAggregateBase *
Cell<NSpaceDimension>
::GetCellularAggregate(void) 
{
  return m_Aggregate;
}



/**
 *    Return a const pointer to the Cellular Aggregate
 */ 
template<unsigned int NSpaceDimension>
const CellularAggregateBase *
Cell<NSpaceDimension>
::GetCellularAggregate(void) const
{
  return m_Aggregate;
}






/**
 *   Set Cellular Aggregate
 */ 
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::SetCellularAggregate( CellularAggregateBase * cells ) 
{
  m_Aggregate = cells;
}




/**
 *    Add a force to the cumulator
 */ 
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::AddForce( const VectorType & force )
{
  if( m_ChemoAttractantLevel > ChemoAttractantLowThreshold &&
      m_ChemoAttractantLevel < ChemoAttractantHighThreshold   )
    {
    double factor = 1.0 / vcl_pow(m_Radius, (double)(NSpaceDimension) );
    m_Force    += force;
    m_Pressure += force.GetNorm() * factor;
    }
  else
    {
    // no force so it is fixed in place....
    }
}


/**
 *    Programmed Cell Death 
 *    This is the cellular equivalent of suicide.
 */ 
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::Apoptosis(void) 
{
  // This call will release the Genomes
  this->Superclass::Apoptosis(); 

  CellularAggregateBase * aggregate = GetCellularAggregate();

  // "this" cell will be destroyed here
  aggregate->Remove( this ); 

}




/**
 *    Execute a time step in the life of the cell.
 *    This is one step in the cell cycle.
 *
 *    Nutrients are acquired
 *    Energy is acquired
 *    If conditions allow it, the cell will grow
 *    The position will be updated according to
 *    applied forces
 */ 
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::AdvanceTimeStep(void) 
{
  // get input from the environment
  this->ReceptorsReading(); 

  // update the level of expression of all the
  // genes in the gene network
  this->ComputeGeneNetwork();

  // this methods produce the effects of gene
  // activation and protein synthesis. It is
  // mostly used for secreting proteins already
  // synthetized in the ComputeGeneNetwork method.
  this->SecreteProducts();

  // If this happens, it is an
  // emergency situation: Do it first.
  if( this->CheckPointApoptosis() )
    {
    m_CycleState = Apop;
    }


  switch( m_CycleState )
  {
  case M: // Mitosis
    m_CycleState = Gap1;
    break;
  case Gap1: // Gap 1 : growing
    {
    if( this->CheckPointDNAReplication() )
      {
      m_CycleState = S;
      }

    break;
    }
  case S:
    m_CycleState = Gap2;
    break;
  case Gap2:
    if( this->CheckPointMitosis() )
      {
      m_CycleState = M;
      }
    break;
  case Gap0:
    // The cell is in cell cycle arrest
    m_CycleState = Gap0;
    break;
  case Apop:
    m_CycleState = Apop;
    break;
  }



  // Atomaton : Execute action
  switch( m_CycleState )
  {
  case M:  // Mitosis
    // This is a terminal action. The implementation of the cell 
    // is destroyed after division. Our abstraction assumes that 
    // the cell disapears and two new cell are created.
    this->Mitosis();
    break;
  case Gap1:
    // Eat and grow
    this->NutrientsIntake();
    this->EnergyIntake();
    this->Grow();
    break;
  case Gap0:
    this->NutrientsIntake();
    this->EnergyIntake();
    break;
  case S:
    this->DNAReplication();
    break; 
  case Gap2:
    break;
  case Apop:
    this->Apoptosis();
    break;
  }
}

/**
 *    Reading substrate using receptors
 */
template<unsigned int NSpaceDimension>
void
Cell<NSpaceDimension>
::ReceptorsReading(void) 
{
  m_Genome->SetExpressionLevel( Pressurin, m_Pressure );
  
  float substrate0 = m_Aggregate->GetSubstrateValue( m_SelfIdentifier, 0 );

  m_ChemoAttractantLevel = substrate0;

}




}  // end namespace bio

}  // end namespace itk

#endif

