cmtkVector.h

Go to the documentation of this file.

/*
//
//  Copyright 1997-2010 Torsten Rohlfing
//
//  Copyright 2004-2010 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 __cmtkVector_h_included_
#define __cmtkVector_h_included_

#include <cmtkconfig.h>

#include <stdio.h>
#include <stdlib.h>
#include <cassert>
#include <math.h>
#include <string.h>

#include <algorithm>

#include <Base/cmtkMathUtil.h>
#include <Base/cmtkTypes.h>
#include <System/cmtkSmartPtr.h>

namespace
cmtk
{


template<class T>
class Vector
{
public:
  size_t Dim;

  T *Elements;

  typedef Vector<T> Self;

  typedef SmartPointer<Self> SmartPtr;


  Vector ( const size_t dim = 0, const T value = 0 ) 
  {
    Dim = dim;
    if ( Dim ) 
      {
      Elements = Memory::AllocateArray<T>( Dim );
      FreeElements = true;
      if ( value==0 )
        memset( Elements, 0, Dim * sizeof(T) );
      else
        for ( size_t i=0; i<Dim; ++i )
          Elements[i]=value;
      } 
    else
      {
      Elements = NULL;
      FreeElements = false;
      }
  }

  Vector ( const size_t dim, T *const elements, const bool freeElements = true ) 
  {
    Dim = dim;
    Elements = elements;
    FreeElements = freeElements;
  }
  
  Vector ( const Vector& other, const size_t len = 0, const size_t from = 0 ) 
  {
    if ( len )
      Dim = std::min( len, other.Dim - from );
    else
      Dim = other.Dim - from;
    
    Elements = Memory::AllocateArray<T>( Dim );
    FreeElements = true;
    memcpy( Elements, other.Elements + from, Dim * sizeof(T) );
  }

  Vector* Clone( const size_t len = 0, const size_t from = 0 ) const
  { 
    return new Vector( *this, len, from ); 
  }
  
  ~Vector () 
  {
    if ( Elements && FreeElements ) 
      {
      Memory::DeleteArray( this->Elements );
      }
  }
  
  Vector& SetDim ( const size_t dim, const bool zero = true ) 
  {
    if ( Dim != dim ) 
      {
      if ( Elements ) 
        {
        Memory::DeleteArray( this->Elements );
        }

      Dim = dim;
      
      if ( Dim ) 
        {
        Elements = Memory::AllocateArray<T>( Dim );
        } 
      else
        Elements = NULL;
      }
    
    if ( zero && Dim ) 
      {
      memset( Elements, 0, Dim * sizeof(T) );
      }
    
    return *this;
  }
  
  Vector& AdjustDimension( const size_t dim, const bool zero = true ) 
  {
    // If old and new size are the same, there is nothing to do.
    if ( Dim != dim ) 
      {
      T* newElements = Memory::AllocateArray<T>( dim );
      // copy common elements
      memcpy( newElements, this->Elements, sizeof(T) * std::min( dim, Dim ) );

      // reset new elements if so desired
      if ( zero && (dim > Dim) )
        {
        memset( newElements + Dim, 0, sizeof(T) * (dim-Dim) );
        }

      // new set new array.
      this->Dim = dim;
      if ( this->FreeElements )
        {
        Memory::DeleteArray( this->Elements );
        }
      this->Elements = newElements;
      this->FreeElements = true;
      } 
    
    return *this;
  }
  
  Vector& operator = ( const Vector& other ) 
  {
    if ( Dim != other.Dim ) {
    if (Elements) 
      {
      Memory::DeleteArray( this->Elements );
      Elements = NULL;
      }
    
    Dim = other.Dim;
    }
    
    if ( Elements == NULL ) 
      {
      Elements = Memory::AllocateArray<T>( Dim );
      }
    
    memcpy( Elements, other.Elements, Dim * sizeof(T) );
    return *this; 
  }

  void CopyToOffset( const Vector& other, const size_t offs = 0, size_t len = 0 )
  {
    if ( ! len ) 
      len = std::min( this->Dim - offs, other.Dim );
    for ( size_t idx=0; idx<len; ++idx )
      Elements[offs+idx] = other.Elements[idx];
  }

  int operator== ( const Vector& other ) const 
  {
    if ( Dim != other.Dim )
      return 0;
    
    for ( size_t i=0; i<Dim; ++i )
      if ( Elements[i] != other.Elements[i] )
        return 0;
    
    return 1;
  }
  
  int operator!= ( const Vector& other ) const 
  {
    return !(*this == other );
  }
  
  T EuclidNorm () const 
  { 
    T Result = 0;
    
#ifndef __SUNPRO_CC
#pragma omp parallel for if (Dim>1e4) reduction(+:Result)
#endif
    for ( size_t i=0; i<this->Dim; ++i ) 
      {
      const T e = Elements[i];
      Result+=e*e;
      }
    
    return sqrt(Result);
  }

  T MaxNorm () const 
  { 
    T Result = 0;
    
    for ( size_t i=0; i<Dim; ++i ) 
      {
      Result = std::max<T>( Result, fabs( Elements[i] ) );
      }
    
    return Result;
  }

  void Clear() 
  { 
    memset( Elements, 0, Dim * sizeof( *Elements ) ); 
  }

  void SetAll( const T value )
  {
#ifndef __SUNPRO_CC
#pragma omp parallel for if (Dim>1e5)
#endif
    for ( size_t i=0; i < this->Dim; ++i ) 
      this->Elements[i] = value;
  }

  T& operator [] ( const size_t index ) 
  {
    return this->Elements[index];
  }

  const T& operator [] ( const size_t index ) const 
  {
    return this->Elements[index];
  }

  Vector<T>& operator+= ( const Vector<T>& delta ) 
  {
    assert( Dim == delta.Dim );

#ifndef __SUNPRO_CC
#pragma omp parallel for if (Dim>1e4)
#endif
    for ( size_t i=0; i<this->Dim; ++i )
      Elements[i] += delta.Elements[i];
    
    return *this;
  }

  Vector<T>& operator-= ( const Vector<T>& delta ) 
  {
    assert( Dim == delta.Dim );
    
#ifndef __SUNPRO_CC
#pragma omp parallel for if (Dim>1e4)
#endif
    for ( size_t i=0; i < this->Dim; ++i )
      Elements[i] -= delta.Elements[i];
    
    return *this;
  }

  Vector<T>& operator*= ( const T a ) 
  {
#ifndef __SUNPRO_CC
#pragma omp parallel for if (Dim>1e4)
#endif
    for ( size_t i=0; i<this->Dim; ++i )
      this->Elements[i] *= a;
    
    return *this;
  }

  void Print ( FILE *const fp = NULL, const char* format = " %f" ) const 
  {
    if ( fp ) 
      {
      for ( size_t idx=0; idx < Dim; ++idx )
        fprintf( fp, format, (float) Elements[idx] );
      fputs( "\n", fp );
      } 
    else
      {
      for ( size_t idx=0; idx < Dim; ++idx )
        printf( format, (float) Elements[idx] );
      puts( "" );
      }
  }
  
  void Sort( const size_t from = 0, const size_t len = 0 ) 
  {
    T *ptr = Elements+from;
    if ( len )
      qsort( ptr, len, sizeof( T ), Vector<T>::Compare );
    else
      qsort( ptr, Dim-from, sizeof( T ), Vector<T>::Compare );
  }
  
private:
  bool FreeElements;

  static int Compare( const void* a, const void* b ) 
  {
    const T *Ta = (const T *) a;
    const T *Tb = (const T *) b;
    return (*Ta > *Tb) - (*Ta < *Tb);
  }
};

typedef Vector<Types::Coordinate> CoordinateVector;

typedef Vector<float> FloatVector;


} // namespace cmtk

#include "cmtkVector.txx"

#endif // #ifndef __cmtkVector_h_included_