mask.h

Go to the documentation of this file.

/*
 *  mask.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef MASK_H
#define MASK_H
 
// C++ includes:
#include <memory>
 
// Includes from libnestutil:
#include "dictionary.h"
#include "numerics.h"
 
// Includes from nestkernel:
#include "exceptions.h"
#include "nest_names.h"
#include "nest_types.h"
 
// Includes from spatial:
#include "position.h"
 
namespace nest
{
class AbstractMask;
using MaskPTR = std::shared_ptr< AbstractMask >;
 
 
class AbstractMask
{
public:
  virtual ~AbstractMask()
  {
  }
 
  virtual bool inside( const std::vector< double >& ) const = 0;
 
  virtual Dictionary
  get_dict() const
  {
    throw KernelException( "Can not convert mask to dict" );
  }
 
  virtual AbstractMask* intersect_mask( const AbstractMask& other ) const = 0;
 
  virtual AbstractMask* union_mask( const AbstractMask& other ) const = 0;
 
  virtual AbstractMask* minus_mask( const AbstractMask& other ) const = 0;
};
 
template < int D >
class Mask : public AbstractMask
{
public:
  using AbstractMask::inside;
 
  ~Mask() override
  {
  }
 
  virtual bool inside( const Position< D >& ) const = 0;
 
  bool inside( const std::vector< double >& pt ) const override;
 
  virtual bool inside( const Box< D >& ) const = 0;
 
  virtual bool outside( const Box< D >& b ) const;
 
  virtual Box< D > get_bbox() const = 0;
 
  virtual Mask* clone() const = 0;
 
  AbstractMask* intersect_mask( const AbstractMask& other ) const override;
  AbstractMask* union_mask( const AbstractMask& other ) const override;
  AbstractMask* minus_mask( const AbstractMask& other ) const override;
};
 
template < int D >
class AllMask : public Mask< D >
{
public:
  ~AllMask()
  {
  }
 
  using Mask< D >::inside;
 
  bool
  inside( const Position< D >& ) const
  {
    return true;
  }
 
  bool
  inside( const Box< D >& ) const
  {
    return true;
  }
 
  bool
  outside( const Box< D >& ) const
  {
    return false;
  }
 
  Box< D >
  get_bbox() const
  {
    const double inf = std::numeric_limits< double >::infinity();
    return Box< D >( Position< D >( -inf, -inf ), Position< D >( inf, inf ) );
  }
 
  Mask< D >*
  clone() const
  {
    return new AllMask();
  }
};
 
template < int D >
class BoxMask : public Mask< D >
{
public:
  BoxMask( const Dictionary& );
 
  BoxMask( const Position< D >& lower_left,
    const Position< D >& upper_right,
    const double azimuth_angle = 0.0,
    const double polar_angle = 0.0 );
 
  ~BoxMask() override
  {
  }
 
  using Mask< D >::inside;
 
  bool inside( const Position< D >& p ) const override;
 
  bool inside( const Box< D >& b ) const override;
 
  bool outside( const Box< D >& b ) const override;
 
  Box< D > get_bbox() const override;
 
  Dictionary get_dict() const override;
 
  Mask< D >* clone() const override;
 
  static std::string get_name();
 
protected:
  void calculate_min_max_values_();
 
  Position< D > lower_left_;
  Position< D > upper_right_;
 
  Position< D > min_values_;
  Position< D > max_values_;
 
  double azimuth_angle_;
  double polar_angle_;
  double azimuth_cos_;
  double azimuth_sin_;
  double polar_cos_;
  double polar_sin_;
 
  Position< D > cntr_;
  Position< D > eps_;
  double cntr_x_az_cos_;
  double cntr_x_az_sin_;
  double cntr_y_az_cos_;
  double cntr_y_az_sin_;
  double cntr_z_pol_cos_;
  double cntr_z_pol_sin_;
  double cntr_x_az_cos_pol_cos_;
  double cntr_x_az_cos_pol_sin_;
  double cntr_y_az_sin_pol_cos_;
  double cntr_y_az_sin_pol_sin_;
  double az_cos_pol_cos_;
  double az_cos_pol_sin_;
  double az_sin_pol_cos_;
  double az_sin_pol_sin_;
 
  bool is_rotated_;
};
 
template < int D >
class BallMask : public Mask< D >
{
public:
  BallMask( Position< D > center, double radius )
    : center_( center )
    , radius_( radius )
  {
  }
 
  BallMask( const Dictionary& );
 
  ~BallMask() override
  {
  }
 
  using Mask< D >::inside;
 
  bool inside( const Position< D >& p ) const override;
 
  bool inside( const Box< D >& ) const override;
 
  bool outside( const Box< D >& b ) const override;
 
  Box< D > get_bbox() const override;
 
  Dictionary get_dict() const override;
 
  Mask< D >* clone() const override;
 
  static std::string get_name();
 
protected:
  Position< D > center_;
  double radius_;
};
 
template < int D >
class EllipseMask : public Mask< D >
{
public:
  EllipseMask( Position< D > center,
    double major_axis,
    double minor_axis,
    double polar_axis,
    double azimuth_angle,
    double polar_angle )
    : center_( center )
    , major_axis_( major_axis )
    , minor_axis_( minor_axis )
    , polar_axis_( polar_axis )
    , azimuth_angle_( azimuth_angle )
    , polar_angle_( polar_angle )
    , x_scale_( 4.0 / ( major_axis_ * major_axis_ ) )
    , y_scale_( 4.0 / ( minor_axis_ * minor_axis_ ) )
    , z_scale_( 4.0 / ( polar_axis_ * polar_axis_ ) )
    , azimuth_cos_( std::cos( azimuth_angle_ * numerics::pi / 180. ) )
    , azimuth_sin_( std::sin( azimuth_angle_ * numerics::pi / 180. ) )
    , polar_cos_( std::cos( polar_angle_ * numerics::pi / 180. ) )
    , polar_sin_( std::sin( polar_angle_ * numerics::pi / 180. ) )
  {
    if ( major_axis_ <= 0 or minor_axis_ <= 0 or polar_axis_ <= 0 )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "All axis > 0 required." );
    }
    if ( major_axis_ < minor_axis_ )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "major_axis greater than minor_axis required." );
    }
    if ( D == 2 and not( polar_angle_ == 0.0 ) )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "polar_angle not defined in 2D." );
    }
 
    create_bbox_();
  }
 
  EllipseMask( const Dictionary& );
 
  ~EllipseMask() override
  {
  }
 
  using Mask< D >::inside;
 
  bool inside( const Position< D >& p ) const override;
 
  bool inside( const Box< D >& ) const override;
 
  bool outside( const Box< D >& b ) const override;
 
  Box< D > get_bbox() const override;
 
  Dictionary get_dict() const override;
 
  Mask< D >* clone() const override;
 
  static std::string get_name();
 
private:
  void create_bbox_();
 
  Position< D > center_;
  double major_axis_;
  double minor_axis_;
  double polar_axis_;
  double azimuth_angle_;
  double polar_angle_;
 
  double x_scale_;
  double y_scale_;
  double z_scale_;
 
  double azimuth_cos_;
  double azimuth_sin_;
  double polar_cos_;
  double polar_sin_;
 
  Box< D > bbox_;
};
 
template < int D >
class IntersectionMask : public Mask< D >
{
public:
  using Mask< D >::inside;
 
  IntersectionMask( const Mask< D >& m1, const Mask< D >& m2 )
    : mask1_( m1.clone() )
    , mask2_( m2.clone() )
  {
  }
 
  IntersectionMask( const IntersectionMask& m )
    : Mask< D >( m )
    , mask1_( m.mask1_->clone() )
    , mask2_( m.mask2_->clone() )
  {
  }
 
  ~IntersectionMask()
  {
    delete mask1_;
    delete mask2_;
  }
 
  bool inside( const Position< D >& p ) const;
 
  bool inside( const Box< D >& b ) const;
 
  bool outside( const Box< D >& b ) const;
 
  Box< D > get_bbox() const;
 
  Mask< D >* clone() const;
 
protected:
  Mask< D >*mask1_, *mask2_;
};
 
template < int D >
class UnionMask : public Mask< D >
{
public:
  using Mask< D >::inside;
 
  UnionMask( const Mask< D >& m1, const Mask< D >& m2 )
    : mask1_( m1.clone() )
    , mask2_( m2.clone() )
  {
  }
 
  UnionMask( const UnionMask& m )
    : Mask< D >( m )
    , mask1_( m.mask1_->clone() )
    , mask2_( m.mask2_->clone() )
  {
  }
 
  ~UnionMask()
  {
    delete mask1_;
    delete mask2_;
  }
 
  bool inside( const Position< D >& p ) const;
 
  bool inside( const Box< D >& b ) const;
 
  bool outside( const Box< D >& b ) const;
 
  Box< D > get_bbox() const;
 
  Mask< D >* clone() const;
 
protected:
  Mask< D >*mask1_, *mask2_;
};
 
template < int D >
class DifferenceMask : public Mask< D >
{
public:
  using Mask< D >::inside;
 
  DifferenceMask( const Mask< D >& m1, const Mask< D >& m2 )
    : mask1_( m1.clone() )
    , mask2_( m2.clone() )
  {
  }
 
  DifferenceMask( const DifferenceMask& m )
    : Mask< D >( m )
    , mask1_( m.mask1_->clone() )
    , mask2_( m.mask2_->clone() )
  {
  }
 
  ~DifferenceMask()
  {
    delete mask1_;
    delete mask2_;
  }
 
  bool inside( const Position< D >& p ) const;
 
  bool inside( const Box< D >& b ) const;
 
  bool outside( const Box< D >& b ) const;
 
  Box< D > get_bbox() const;
 
  Mask< D >* clone() const;
 
protected:
  Mask< D >*mask1_, *mask2_;
};
 
 
template < int D >
class ConverseMask : public Mask< D >
{
public:
  using Mask< D >::inside;
 
  ConverseMask( const Mask< D >& m )
    : m_( m.clone() )
  {
  }
 
  ConverseMask( const ConverseMask& m )
    : Mask< D >( m )
    , m_( m.m_->clone() )
  {
  }
 
  ~ConverseMask()
  {
    delete m_;
  }
 
  bool inside( const Position< D >& p ) const;
 
  bool inside( const Box< D >& b ) const;
 
  bool outside( const Box< D >& b ) const;
 
  Box< D > get_bbox() const;
 
  Mask< D >* clone() const;
 
protected:
  Mask< D >* m_;
};
 
 
template < int D >
class AnchoredMask : public Mask< D >
{
public:
  using Mask< D >::inside;
 
  AnchoredMask( const Mask< D >& m, Position< D > anchor )
    : m_( m.clone() )
    , anchor_( anchor )
  {
  }
 
  AnchoredMask( const AnchoredMask& m )
    : Mask< D >( m )
    , m_( m.m_->clone() )
    , anchor_( m.anchor_ )
  {
  }
 
  ~AnchoredMask()
  {
    delete m_;
  }
 
  bool inside( const Position< D >& p ) const;
 
  bool inside( const Box< D >& b ) const;
 
  bool outside( const Box< D >& b ) const;
 
  Box< D > get_bbox() const;
 
  Dictionary get_dict() const;
 
  Mask< D >* clone() const;
 
protected:
  Mask< D >* m_;
  Position< D > anchor_;
};
 
template <>
inline std::string
BoxMask< 2 >::get_name()
{
  return names::rectangular;
}
 
template <>
inline std::string
BoxMask< 3 >::get_name()
{
  return names::box;
}
 
template < int D >
BoxMask< D >::BoxMask( const Dictionary& d )
{
  lower_left_ = d.get< std::vector< double > >( names::lower_left );
  upper_right_ = d.get< std::vector< double > >( names::upper_right );
 
  if ( not( lower_left_ < upper_right_ ) )
  {
    throw BadProperty(
      "BoxMask<D>: "
      "Upper right must be strictly to the right and above lower left." );
  }
 
  if ( d.known( names::azimuth_angle ) )
  {
    azimuth_angle_ = d.get< double >( names::azimuth_angle );
  }
  else
  {
    azimuth_angle_ = 0.0;
  }
 
  if ( d.known( names::polar_angle ) )
  {
    if ( D == 2 )
    {
      throw BadProperty(
        "BoxMask<D>: "
        "polar_angle not defined in 2D." );
    }
    polar_angle_ = d.get< double >( names::polar_angle );
  }
  else
  {
    polar_angle_ = 0.0;
  }
 
  azimuth_cos_ = std::cos( azimuth_angle_ * numerics::pi / 180. );
  azimuth_sin_ = std::sin( azimuth_angle_ * numerics::pi / 180. );
  polar_cos_ = std::cos( polar_angle_ * numerics::pi / 180. );
  polar_sin_ = std::sin( polar_angle_ * numerics::pi / 180. );
 
  cntr_ = ( upper_right_ + lower_left_ ) * 0.5;
  for ( int i = 0; i != D; ++i )
  {
    eps_[ i ] = 1e-12;
  }
 
  cntr_x_az_cos_ = cntr_[ 0 ] * azimuth_cos_;
  cntr_x_az_sin_ = cntr_[ 0 ] * azimuth_sin_;
  cntr_y_az_cos_ = cntr_[ 1 ] * azimuth_cos_;
  cntr_y_az_sin_ = cntr_[ 1 ] * azimuth_sin_;
  if ( D == 3 )
  {
    cntr_z_pol_cos_ = cntr_[ 2 ] * polar_cos_;
    cntr_z_pol_sin_ = cntr_[ 2 ] * polar_sin_;
    cntr_x_az_cos_pol_cos_ = cntr_x_az_cos_ * polar_cos_;
    cntr_x_az_cos_pol_sin_ = cntr_x_az_cos_ * polar_sin_;
    cntr_y_az_sin_pol_cos_ = cntr_y_az_sin_ * polar_cos_;
    cntr_y_az_sin_pol_sin_ = cntr_y_az_sin_ * polar_sin_;
    az_cos_pol_cos_ = azimuth_cos_ * polar_cos_;
    az_cos_pol_sin_ = azimuth_cos_ * polar_sin_;
    az_sin_pol_cos_ = azimuth_sin_ * polar_cos_;
    az_sin_pol_sin_ = azimuth_sin_ * polar_sin_;
  }
  else
  {
    cntr_z_pol_cos_ = 0.0;
    cntr_z_pol_sin_ = 0.0;
    cntr_x_az_cos_pol_cos_ = 0.0;
    cntr_x_az_cos_pol_sin_ = 0.0;
    cntr_y_az_sin_pol_cos_ = 0.0;
    cntr_y_az_sin_pol_sin_ = 0.0;
    az_cos_pol_cos_ = 0.0;
    az_cos_pol_sin_ = 0.0;
    az_sin_pol_cos_ = 0.0;
    az_sin_pol_sin_ = 0.0;
  }
 
  is_rotated_ = azimuth_angle_ != 0.0 or polar_angle_ != 0.0;
 
  calculate_min_max_values_();
}
 
template < int D >
inline BoxMask< D >::BoxMask( const Position< D >& lower_left,
  const Position< D >& upper_right,
  const double azimuth_angle,
  const double polar_angle )
  : lower_left_( lower_left )
  , upper_right_( upper_right )
  , azimuth_angle_( azimuth_angle )
  , polar_angle_( polar_angle )
  , azimuth_cos_( std::cos( azimuth_angle_ * numerics::pi / 180. ) )
  , azimuth_sin_( std::sin( azimuth_angle_ * numerics::pi / 180. ) )
  , polar_cos_( std::cos( polar_angle_ * numerics::pi / 180. ) )
  , polar_sin_( std::sin( polar_angle_ * numerics::pi / 180. ) )
  , cntr_( ( upper_right_ + lower_left_ ) * 0.5 )
  , cntr_x_az_cos_( cntr_[ 0 ] * azimuth_cos_ )
  , cntr_x_az_sin_( cntr_[ 0 ] * azimuth_sin_ )
  , cntr_y_az_cos_( cntr_[ 1 ] * azimuth_cos_ )
  , cntr_y_az_sin_( cntr_[ 1 ] * azimuth_sin_ )
{
  if ( D == 2 and not( polar_angle_ == 0.0 ) )
  {
    throw BadProperty(
      "BoxMask<D>: "
      "polar_angle not defined in 2D." );
  }
 
  for ( int i = 0; i != D; ++i )
  {
    eps_[ i ] = 1e-12;
  }
 
  if ( D == 3 )
  {
    cntr_z_pol_cos_ = cntr_[ 2 ] * polar_cos_;
    cntr_z_pol_sin_ = cntr_[ 2 ] * polar_sin_;
    cntr_x_az_cos_pol_cos_ = cntr_x_az_cos_ * polar_cos_;
    cntr_x_az_cos_pol_sin_ = cntr_x_az_cos_ * polar_sin_;
    cntr_y_az_sin_pol_cos_ = cntr_y_az_sin_ * polar_cos_;
    cntr_y_az_sin_pol_sin_ = cntr_y_az_sin_ * polar_sin_;
    az_cos_pol_cos_ = azimuth_cos_ * polar_cos_;
    az_cos_pol_sin_ = azimuth_cos_ * polar_sin_;
    az_sin_pol_cos_ = azimuth_sin_ * polar_cos_;
    az_sin_pol_sin_ = azimuth_sin_ * polar_sin_;
  }
  else
  {
    cntr_z_pol_cos_ = 0.0;
    cntr_z_pol_sin_ = 0.0;
    cntr_x_az_cos_pol_cos_ = 0.0;
    cntr_x_az_cos_pol_sin_ = 0.0;
    cntr_y_az_sin_pol_cos_ = 0.0;
    cntr_y_az_sin_pol_sin_ = 0.0;
    az_cos_pol_cos_ = 0.0;
    az_cos_pol_sin_ = 0.0;
    az_sin_pol_cos_ = 0.0;
    az_sin_pol_sin_ = 0.0;
  }
 
  is_rotated_ = azimuth_angle_ != 0.0 or polar_angle_ != 0.0;
 
  calculate_min_max_values_();
}
 
template <>
inline std::string
BallMask< 2 >::get_name()
{
  return names::circular;
}
 
template <>
inline std::string
BallMask< 3 >::get_name()
{
  return names::spherical;
}
 
template < int D >
BallMask< D >::BallMask( const Dictionary& d )
{
  radius_ = d.get< double >( names::radius );
  if ( radius_ <= 0 )
  {
    throw BadProperty( "BallMask<D>: radius > 0 required." );
  }
 
  if ( d.known( names::anchor ) )
  {
    center_ = d.get< std::vector< double > >( names::anchor );
  }
}
 
template <>
inline std::string
EllipseMask< 2 >::get_name()
{
  return names::elliptical;
}
 
template <>
inline std::string
EllipseMask< 3 >::get_name()
{
  return names::ellipsoidal;
}
 
template < int D >
EllipseMask< D >::EllipseMask( const Dictionary& d )
{
  major_axis_ = d.get< double >( names::major_axis );
  minor_axis_ = d.get< double >( names::minor_axis );
  if ( major_axis_ <= 0 or minor_axis_ <= 0 )
  {
    throw BadProperty(
      "EllipseMask<D>: "
      "All axis > 0 required." );
  }
  if ( major_axis_ < minor_axis_ )
  {
    throw BadProperty(
      "EllipseMask<D>: "
      "major_axis greater than minor_axis required." );
  }
 
  x_scale_ = 4.0 / ( major_axis_ * major_axis_ );
  y_scale_ = 4.0 / ( minor_axis_ * minor_axis_ );
 
  if ( d.known( names::polar_axis ) )
  {
    if ( D == 2 )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "polar_axis not defined in 2D." );
    }
    polar_axis_ = d.get< double >( names::polar_axis );
 
    if ( polar_axis_ <= 0 )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "All axis > 0 required." );
    }
 
    z_scale_ = 4.0 / ( polar_axis_ * polar_axis_ );
  }
  else
  {
    polar_axis_ = 0.0;
    z_scale_ = 0.0;
  }
 
  if ( d.known( names::anchor ) )
  {
    center_ = d.get< std::vector< double > >( names::anchor );
  }
 
  if ( d.known( names::azimuth_angle ) )
  {
    azimuth_angle_ = d.get< double >( names::azimuth_angle );
  }
  else
  {
    azimuth_angle_ = 0.0;
  }
 
  if ( d.known( names::polar_angle ) )
  {
    if ( D == 2 )
    {
      throw BadProperty(
        "EllipseMask<D>: "
        "polar_angle not defined in 2D." );
    }
    polar_angle_ = d.get< double >( names::polar_angle );
  }
  else
  {
    polar_angle_ = 0.0;
  }
 
  azimuth_cos_ = std::cos( azimuth_angle_ * numerics::pi / 180. );
  azimuth_sin_ = std::sin( azimuth_angle_ * numerics::pi / 180. );
  polar_cos_ = std::cos( polar_angle_ * numerics::pi / 180. );
  polar_sin_ = std::sin( polar_angle_ * numerics::pi / 180. );
 
  create_bbox_();
}
 
}  // namespace nest
 
#endif