libinsighttoolkit3-dev 3.20.1-1 / usr / include / InsightToolkit / Algorithms / itkSphereMeshSource.txx

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkSphereMeshSource.txx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) Insight Software 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 __itkSphereMeshSource_txx
#define __itkSphereMeshSource_txx

#include "itkSphereMeshSource.h"

namespace itk
{

/**
 *
 */
template<class TOutputMesh>
SphereMeshSource<TOutputMesh>
::SphereMeshSource()
{
  /**
   * Create the output
   */
  typename TOutputMesh::Pointer output = TOutputMesh::New();
  this->ProcessObject::SetNumberOfRequiredOutputs(1);
  this->ProcessObject::SetNthOutput(0, output.GetPointer());
  m_Squareness1 = 1.0;
  m_Squareness2 = 1.0;
  m_Center.Fill(0);
  m_Scale.Fill(1);
  m_ResolutionX = 4;
  m_ResolutionY = 4;
}

/*
 *
 */
template<class TOutputMesh>
void
SphereMeshSource<TOutputMesh>
::GenerateData()
{
  unsigned long i, j, jn, p, numpts;
  double ustep, vstep, ubeg, vbeg, u, v; 
  int signu, signv; 

  // calculate the number os cells and points
  numpts = m_ResolutionX*m_ResolutionY + 2; 

  // calculate the steps using resolution
  ustep = vnl_math::pi / (m_ResolutionX+1); 
  vstep = 2.0*vnl_math::pi / m_ResolutionY; 
  ubeg = (-vnl_math::pi/2.0) + ustep; 
  vbeg = -vnl_math::pi; 

  ///////////////////////////////////////////////////////////////////////////
  // nodes allocation

  // the temporary container of nodes' connectness
  unsigned long tripoints[3] = {0,1,2};
  
  // memory allocation for nodes
  typename OutputMeshType::Pointer outputMesh = this->GetOutput();  

  outputMesh->GetPoints()->Reserve(numpts);

  outputMesh->SetCellsAllocationMethod( OutputMeshType::CellsAllocatedDynamicallyCellByCell );

  PointsContainerPointer  myPoints = outputMesh->GetPoints();
  typename PointsContainer::Iterator   point = myPoints->Begin();

  OPointType p1;

  // calculate all regular nodes
  while( point != myPoints->End() ) 
    { 
    for (u=ubeg, i=0; i < m_ResolutionX; u += ustep, i++) 
      { 
      for (v=vbeg, j=0; j < m_ResolutionY; v += vstep, j++) 
        { 
        if (vcl_cos(u) > 0) 
          {
          signu = 1;
          } 
        else 
          {
          signu = -1;
          }
        if (vcl_cos(v) > 0) 
          {
          signv = 1;
          } 
        else 
          {
          signv = -1;
          }
        
        p1[0] = m_Scale[0]*signu*(vcl_pow((float)(vcl_fabs(vcl_cos(u))), (float) m_Squareness1))*signv* 
          (vcl_pow((float)(vcl_fabs(vcl_cos(v))), (float) m_Squareness2)) + m_Center[0]; 
        
        if (vcl_sin(v) > 0) 
          {
          signv = 1;
          } 
        else 
          {
          signv = -1;
          }
        
        p1[1] = m_Scale[1]*signu*(vcl_pow((float)(vcl_fabs(vcl_cos(u))), (float) m_Squareness1))*signv* 
          (vcl_pow((float)(vcl_fabs(vcl_sin(v))), (float) m_Squareness2)) + m_Center[1]; 
        
        if (vcl_sin(u) > 0) 
          {
          signu = 1;
          } 
        else 
          {
          signu = -1;
          }
        
        p1[2] = m_Scale[2]*signu*(vcl_pow((float)(vcl_fabs(vcl_sin(u))), (float) m_Squareness1)) + 
          m_Center[2];
        
        point.Value() = p1;
        ++point;
        } 
      }   

    // calculate the south pole node
    p1[0] = (m_Scale[0]*(vcl_pow((float)(vcl_fabs(vcl_cos(-vnl_math::pi/2))),1.0f))* 
             (vcl_pow((float)(vcl_fabs(vcl_cos(0.0))),1.0f)) + m_Center[0]); 
    p1[1] = (m_Scale[1]*(vcl_pow((float)(vcl_fabs(vcl_cos(-vnl_math::pi/2))),1.0f))* 
             (vcl_pow((float)(vcl_fabs(vcl_sin(0.0))),1.0f)) + m_Center[1]); 
    p1[2] = (m_Scale[2]*-1*(vcl_pow((float)(vcl_fabs(vcl_sin(-vnl_math::pi/2))),1.0f)) 
             + m_Center[2]);
    point.Value() = p1;
    ++point;
    
    // calculate the north pole node
    p1[0] = (m_Scale[0]*(vcl_pow((float)(vcl_fabs(vcl_cos(vnl_math::pi/2))),1.0f))* 
             (vcl_pow(vcl_fabs(vcl_cos(0.0)),1.0)) + m_Center[0]); 
    p1[1] = (m_Scale[1]*(vcl_pow((float)(vcl_fabs(vcl_cos(vnl_math::pi/2))),1.0f))* 
             (vcl_pow(vcl_fabs(vcl_sin(0.0)),1.0)) + m_Center[1]); 
    p1[2] = (m_Scale[2]*(vcl_pow((float)(vcl_fabs(vcl_sin(vnl_math::pi/2))),1.0f)) 
             + m_Center[2]);
    point.Value() = p1;
    ++point;
    }
  
  ///////////////////////////////////////////////////////////////////////////
  // cells allocation
  p = 0;

  // store all regular cells
  CellAutoPointer testCell;
  for(unsigned int ii=0; ii+1 < m_ResolutionX; ii++) 
    {
    for (unsigned int jj=0; jj<m_ResolutionY; jj++) 
      {
      jn = (jj+1)%m_ResolutionY; 
      tripoints[0] = ii*m_ResolutionY+jj; 
      tripoints[1] = tripoints[0]-jj+jn; 
      tripoints[2] = tripoints[0]+m_ResolutionY; 
      testCell.TakeOwnership( new TriCellType );
      testCell->SetPointIds(tripoints);
      outputMesh->SetCell(p, testCell );
      outputMesh->SetCellData(p, (OPixelType)3.0);
      p++;
      testCell.TakeOwnership( new TriCellType );
      tripoints[0] = tripoints[1]; 
      tripoints[1] = tripoints[0]+m_ResolutionY; 
      testCell->SetPointIds(tripoints);
      outputMesh->SetCell(p, testCell );
      outputMesh->SetCellData(p, (OPixelType)3.0);
      p++;
      }
    }
  
  // store cells containing the south pole nodes
  for (unsigned int jj=0; jj<m_ResolutionY; jj++) 
    {
    jn = (jj+1)%m_ResolutionY; 
    tripoints[0] = numpts-2; 
    tripoints[1] = jn; 
    tripoints[2] = jj; 
    testCell.TakeOwnership( new TriCellType );
    testCell->SetPointIds(tripoints);
    outputMesh->SetCell(p, testCell );
    outputMesh->SetCellData(p, (OPixelType)1.0);
    p++;
    }
  
  // store cells containing the north pole nodes
  for (unsigned int jj=0; jj<m_ResolutionY; jj++) 
    {
    jn = (jj+1)%m_ResolutionY; 
    tripoints[2] = (m_ResolutionX-1)*m_ResolutionY+jj; 
    tripoints[1] = numpts-1; 
    tripoints[0] = tripoints[2]-jj+jn; 
    testCell.TakeOwnership( new TriCellType );
    testCell->SetPointIds(tripoints);
    outputMesh->SetCell(p, testCell );
    outputMesh->SetCellData(p, (OPixelType)2.0);
    p++;
    }
}

template<class TOutputMesh>
void
SphereMeshSource<TOutputMesh>
::PrintSelf( std::ostream& os, Indent indent ) const
{
  Superclass::PrintSelf(os,indent);

  os << indent << "Center: " << m_Center << std::endl;
  os << indent << "Scale: " << m_Scale << std::endl;
  os << indent << "ResolutionX: " << m_ResolutionX << std::endl;
  os << indent << "ResolutionX: " << m_ResolutionY << std::endl;
  os << indent << "Squareness1: " << m_Squareness1 << std::endl;
  os << indent << "Squareness2: " << m_Squareness2 << std::endl;
}

} // end namespace itk

#endif