cmtkUniformVolume.h

Go to the documentation of this file.

/*
//
//  Copyright 1997-2010 Torsten Rohlfing
//
//  Copyright 2004-2011 SRI International
//
//  This file is part of the Computational Morphometry Toolkit.
//
//  http://www.nitrc.org/projects/cmtk/
//
//  The Computational Morphometry Toolkit is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  The Computational Morphometry Toolkit is distributed in the hope that it
//  will be useful, but WITHOUT ANY WARRANTY; without even the implied
//  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with the Computational Morphometry Toolkit.  If not, see
//  <http://www.gnu.org/licenses/>.
//
//  $Revision: 2752 $
//
//  $LastChangedDate: 2011-01-17 11:33:31 -0800 (Mon, 17 Jan 2011) $
//
//  $LastChangedBy: torstenrohlfing $
//
*/

#ifndef __cmtkUniformVolume_h_included_
#define __cmtkUniformVolume_h_included_

#include <cmtkconfig.h>

#include <Base/cmtkVolume.h>

#include <System/cmtkSmartPtr.h>
#include <System/cmtkSmartConstPtr.h>

#include <algorithm>

namespace
cmtk
{


class VolumeGridToGridLookup;

class UniformVolume : 
  public Volume 
{
public:
  typedef UniformVolume Self;

  typedef Volume Superclass;

  typedef SmartPointer<Self> SmartPtr;

  typedef SmartConstPointer<Self> SmartConstPtr;

  typedef Superclass:: CoordinateRegionType CoordinateRegionType;

  typedef Superclass::CoordinateVectorType CoordinateVectorType;

  CoordinateVectorType m_Delta; 

  const CoordinateVectorType& Deltas() const
  {
    return this->m_Delta;
  }

  CoordinateVectorType& Deltas()
  {
    return this->m_Delta;
  }

  virtual ~UniformVolume() {}

  UniformVolume( const UniformVolume& other, const Types::Coordinate resolution = 0, const bool allowUpsampling = false );

  UniformVolume( const DataGrid::IndexType& dims, const Self::CoordinateVectorType& size, TypedArray::SmartPtr& data = TypedArray::SmartPtr::Null );

  UniformVolume( const DataGrid::IndexType& dims, const Types::Coordinate deltaX, const Types::Coordinate deltaY, const Types::Coordinate deltaZ, TypedArray::SmartPtr& data = TypedArray::SmartPtr::Null );

  bool GridMatches( const Self& other ) const
  {
    return Superclass::GridMatches( other ) && ((this->m_Delta-other.m_Delta).MaxAbsValue() < 1e-5) && ((this->m_Offset-other.m_Offset).MaxAbsValue() < 1e-5);
  }

  AffineXform::MatrixType m_IndexToPhysicalMatrix;

  virtual void ChangeCoordinateSpace( const std::string& newSpace );
  
  std::string GetOrientationFromDirections() const;

  virtual void CreateDefaultIndexToPhysicalMatrix();

  virtual AffineXform::MatrixType GetImageToPhysicalMatrix() const;

  Self::SmartPtr Clone( const bool copyData )
  {
    return Self::SmartPtr( this->CloneVirtual( copyData ) );
  }

  Self::SmartPtr Clone() const
  {
    return Self::SmartPtr( this->CloneVirtual() );
  }

  Self::SmartPtr CloneGrid() const
  {
    return Self::SmartPtr( this->CloneGridVirtual() );
  }

  virtual Types::Coordinate GetDelta( const int axis ) const
  {
    return this->m_Delta[axis];
  }

  virtual Types::Coordinate GetMinDelta () const 
  {
    return this->m_Delta.MinValue();
  }
  
  virtual Types::Coordinate GetMaxDelta () const 
  {
    return this->m_Delta.MaxValue();
  }

  virtual TypedArray::SmartPtr Resample ( const UniformVolume& ) const;

  const UniformVolume::SmartPtr GetReoriented ( const char* newOrientation = AnatomicalOrientation::ORIENTATION_STANDARD ) const;

  virtual UniformVolume* GetDownsampledAndAveraged( const int downsample, const bool approxIsotropic = false ) const;
 
  virtual UniformVolume* GetDownsampledAndAveraged( const int (&downsample)[3] ) const;

  UniformVolume* GetInterleavedSubVolume( const int axis, const int factor, const int idx ) const;

  UniformVolume* GetInterleavedPaddedSubVolume( const int axis, const int factor, const int idx ) const;

  virtual void Mirror ( const int axis = AXIS_X );

  virtual ScalarImage* GetOrthoSlice( const int axis, const unsigned int plane ) const;

  Self::SmartPtr ExtractSlice( const int axis, const int plane ) const
  {
    return Self::SmartPtr( this->ExtractSliceVirtual( axis, plane ) );
  }

  virtual ScalarImage* GetOrthoSliceInterp( const int axis, const Types::Coordinate location ) const;

  virtual ScalarImage* GetNearestOrthoSlice( const int axis, const Types::Coordinate location ) const;

  virtual const Self::CoordinateVectorType GetGradientAt( const int i, const int j, const int k );

  virtual void GetHessianAt( Matrix3x3<Types::DataItem>& H, const int i, const int j, const int k );

  template<class TData> inline bool ProbeData( TData& result, const TData* dataPtr, const Self::CoordinateVectorType& location ) const;

  inline bool ProbeNoXform ( ProbeInfo&, const Self::CoordinateVectorType& ) const;

  inline bool FindVoxel( const Self::CoordinateVectorType& location, int *const idx, Types::Coordinate *const from, Types::Coordinate *const to ) const;
  
  inline bool FindVoxel( const Self::CoordinateVectorType& location, int *const idx ) const;

  inline void GetVoxelIndexNoBounds( const Self::CoordinateVectorType& location, int *const idx ) const;

  inline bool FindVoxelByIndex( const Self::CoordinateVectorType& fracIndex, int *const idx, Types::Coordinate *const frac ) const;

  const UniformVolume::RegionType GetGridRange( const Self::CoordinateVectorType& fromVOI, const Self::CoordinateVectorType& toVOI ) const;
 
  virtual Types::Coordinate GetPlaneCoord( const int axis, const int plane ) const 
  {
    return this->m_Offset[axis] + plane * this->m_Delta[axis];
  }

  virtual int GetCoordIndex( const int axis, const Types::Coordinate location ) const 
  {
    return std::max<int>( 0, std::min<int>( (int) ((location-this->m_Offset[axis]) / this->m_Delta[axis]), this->m_Dims[axis]-1 ) );
  }
  
  virtual int GetClosestCoordIndex( const int axis, const Types::Coordinate location ) const 
  {
    const int idx = (int)MathUtil::Round((location-this->m_Offset[axis]) / this->m_Delta[axis]);
    return std::max<int>( 0, std::min<int>( idx, this->m_Dims[axis]-1 ) );
  }

  virtual bool GetClosestGridPointIndex( const Self::CoordinateVectorType v, int *const xyz ) const 
  {
    for ( int dim = 0; dim < 3; ++dim )
      {
      xyz[dim] = (int) MathUtil::Round((v[dim]-this->m_Offset[dim]) / this->m_Delta[dim]);
      if ( (xyz[dim] < 0) || ( xyz[dim] > this->m_Dims[dim]-1) )
        return false;
      }
    return true;
  }

  virtual int GetTruncCoordIndex( const int axis, const Types::Coordinate location ) const 
  {
    const int idx = static_cast<int>((location-this->m_Offset[axis]) / this->m_Delta[axis]);
    return std::max<int>( 0, std::min<int>( idx, this->m_Dims[axis]-1 ) );
  }

  virtual bool GetTruncGridPointIndex( const Self::CoordinateVectorType v , 
                                       int *const xyz  ) const 
  {
    for ( int dim = 0; dim < 3; ++dim )
      {
      xyz[dim] = static_cast<int>((v[dim]-this->m_Offset[dim]) / this->m_Delta[dim]);
      if ( (xyz[dim] < 0) || ( xyz[dim] > this->m_Dims[dim]-1) )
        return false;
      }
    return true;
  }

  virtual const Self::CoordinateVectorType GetGridLocation( const int x, const int y, const int z ) const 
  {
    const Types::Coordinate loc[3] = { this->m_Offset[0] + x * this->m_Delta[0], this->m_Offset[1] + y * this->m_Delta[1], this->m_Offset[2] + z * this->m_Delta[2] };
    return Self::CoordinateVectorType( loc );
  }
  
  virtual const Self::CoordinateVectorType IndexToPhysical( const Self::CoordinateVectorType& idxV ) const 
  {
    const Self::CoordinateVectorType v( idxV );
    return v * this->m_IndexToPhysicalMatrix;
  }
  
  virtual const Self::CoordinateVectorType GetGridLocation( const size_t idx ) const 
  {
    const Types::Coordinate loc[3] = { this->m_Offset[0] +  (idx % this->nextJ) * this->m_Delta[0], 
                                       this->m_Offset[1] +  (idx % this->nextK) / this->nextJ * this->m_Delta[1], 
                                       this->m_Offset[2] +  (idx / this->nextK) * this->m_Delta[2] };
    return Self::CoordinateVectorType( loc );
  }


  void SetHighResCropRegion( const Self::CoordinateRegionType& crop );

  const Self::CoordinateRegionType GetHighResCropRegion() const;

  virtual void SetCropRegion( const Self::RegionType& region )
  {
    this->m_HighResCropRegion = Self::CoordinateRegionType::SmartPtr( NULL );
    Superclass::SetCropRegion( region );
  }

  Self::CoordinateVectorType GetCenterCropRegion() const 
  {
    const Self::CoordinateRegionType region = this->GetHighResCropRegion();
    return 0.5 * ( region.From() + region.To() );
  }
  
  Self::SmartPtr GetCroppedVolume() const;

  template<class TAccumulator>
  ScalarImage* ComputeProjection( const int axis ) const;
  
  virtual Self::CoordinateVectorType GetCenterOfMass() const
  {
    Self::CoordinateVectorType com = this->Superclass::GetCenterOfMassGrid();
    for ( int dim = 0; dim < 3; ++dim )
      (com[dim] *= this->m_Delta[dim]) += this->m_Offset[dim];
    return com;
  }
  
  virtual Self::CoordinateVectorType GetCenterOfMass( Self::CoordinateVectorType& firstOrderMoment ) const
  {
    Self::CoordinateVectorType com = this->Superclass::GetCenterOfMassGrid( firstOrderMoment );
    for ( int dim = 0; dim < 3; ++dim )
      {
      (com[dim] *= this->m_Delta[dim]) += this->m_Offset[dim];
      firstOrderMoment[dim] *= this->m_Delta[dim];
      }
    return com;
  }

  void GetPrincipalAxes( Matrix3x3<Types::Coordinate>& directions, Self::CoordinateVectorType& centerOfMass ) const;

protected:
  virtual Self* CloneVirtual( const bool copyData );
  virtual Self* CloneVirtual() const;

  virtual Self* CloneGridVirtual() const;

private:
  Self::CoordinateRegionType::SmartPtr m_HighResCropRegion;

  friend class WarpXform;
  
  friend class VolumeAxesHash;

  friend class VolumeGridToGridLookup;

  class ResampleThreadInfo;

  friend class ResampleThreadInfo;

  class ResampleTaskInfo :
    public ThreadParameters<const Self> 
  {
  public:
    Types::DataItem *ResampledData;
    const VolumeGridToGridLookup *GridLookup;
    const Self* OtherVolume;
    const TypedArray* FromData;
  };

  static void ResampleThreadPoolExecuteGrey( void *const arg, const size_t taskIdx, const size_t taskCnt, const size_t threadIdx, const size_t threadCnt );

  static void ResampleThreadPoolExecuteLabels( void *const arg, const size_t taskIdx, const size_t taskCnt, const size_t threadIdx, const size_t threadCnt );

  virtual Self* ExtractSliceVirtual( const int axis, const int plane ) const;
};


} // namespace cmtk

#include "cmtkUniformVolume.txx"

#endif // #ifndef __cmtkUniformVolume_h_included_