diff --git a/Core/Algorithms/inc/Mask.h b/Core/Algorithms/inc/Mask.h
index 3f3b446fd619c8999102fe1a6e2a40119bff3038..ba3507fb23b6a475a5339d4a217561fe51fce71f 100644
--- a/Core/Algorithms/inc/Mask.h
+++ b/Core/Algorithms/inc/Mask.h
@@ -1,6 +1,5 @@
 #ifndef MASK_H_
 #define MASK_H_
-#include <cstddef>
 // ********************************************************************
 // * The BornAgain project                                            *
 // * Simulation of neutron and x-ray scattering at grazing incidence  *
@@ -15,6 +14,10 @@
 //! @author Scientific Computing Group at FRM II
 //! @date   Nov 15, 2012
 
+#include "MaskCoordinateFunction.h"
+
+#include <cstddef>
+
 //- -------------------------------------------------------------------
 //! @class Mask
 //! @brief Definition of base class for masking OutputData elements
@@ -23,6 +26,7 @@ class Mask
 {
 public:
     template <class TValue, class TContainer> friend class OutputDataIterator;
+    template <class TValue> friend class OutputData;
     Mask(Mask *p_submask=0);
     virtual ~Mask() { delete mp_submask; }
     virtual Mask *clone() const;
@@ -31,10 +35,13 @@ public:
     size_t getNextIndex(size_t total_index);
     void setMaxIndex(size_t max_index);
 
+    size_t getOwnIndex() const { return m_own_index; }
+    size_t getMaxIndex() const { return m_max_index; }
+
 protected:
     virtual bool isMasked(size_t total_index) const;
-    size_t m_mask_index;
-    size_t m_end_index;
+    size_t m_own_index;
+    size_t m_max_index;
     Mask *mp_submask;
 private:
     size_t nextSubIndex(size_t total_index);
@@ -58,4 +65,26 @@ private:
     size_t m_remainder;
 };
 
+//- -------------------------------------------------------------------
+//! @class MaskCoordinates
+//! @brief Mask based on the coordinates
+//- -------------------------------------------------------------------
+class MaskCoordinates : public Mask
+{
+public:
+    MaskCoordinates(size_t rank, int *dims, Mask *p_submask=0);
+    virtual ~MaskCoordinates();
+    virtual MaskCoordinates *clone() const;
+
+    void setMaskCoordinateFunction(MaskCoordinateFunction *p_mask_function);
+
+private:
+    virtual bool isMasked(size_t total_index) const;
+    void setCachedCoordinates(size_t index) const;
+    size_t m_rank;
+    int *m_dims;
+    mutable int *m_coordinates;
+    MaskCoordinateFunction *mp_mask_function;
+};
+
 #endif /* MASK_H_ */
diff --git a/Core/Algorithms/inc/MaskCoordinateFunction.h b/Core/Algorithms/inc/MaskCoordinateFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..de3251fc1501c3d18d33ad698ddd2b7d02bb2e9c
--- /dev/null
+++ b/Core/Algorithms/inc/MaskCoordinateFunction.h
@@ -0,0 +1,64 @@
+#ifndef MASKCOORDINATEFUNCTION_H_
+#define MASKCOORDINATEFUNCTION_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   MaskCoordinateFunction.h
+//! @brief  Definition of MaskCoordinateFunction classes
+//! @author Scientific Computing Group at FRM II
+//! @date   Nov 20, 2012
+
+#include <cstddef>
+
+class MaskCoordinateFunction
+{
+public:
+    MaskCoordinateFunction(size_t rank);
+    virtual MaskCoordinateFunction *clone() const;
+    virtual ~MaskCoordinateFunction() {}
+
+    bool isMasked(size_t rank, const int *coordinates) const;
+    void setInvertFlag(bool invert) { m_invert = invert; }
+protected:
+    virtual bool isInStandardMaskedArea(const int *coordinates) const;
+    size_t m_rank;
+    bool m_invert;  //!< if true, the complement is masked instead
+};
+
+class MaskCoordinateRectangleFunction : public MaskCoordinateFunction
+{
+public:
+    MaskCoordinateRectangleFunction(size_t rank, const int *minima, const int *maxima);
+    virtual MaskCoordinateRectangleFunction *clone() const;
+    virtual ~MaskCoordinateRectangleFunction();
+
+protected:
+    virtual bool isInStandardMaskedArea(const int *coordinates) const;
+
+private:
+    int *m_minima;
+    int *m_maxima;
+};
+
+class MaskCoordinateEllipseFunction : public MaskCoordinateFunction
+{
+public:
+    MaskCoordinateEllipseFunction(size_t rank, const int *center, const int *radii);
+    virtual MaskCoordinateEllipseFunction *clone() const;
+    virtual ~MaskCoordinateEllipseFunction();
+
+protected:
+    virtual bool isInStandardMaskedArea(const int *coordinates) const;
+
+private:
+    int *m_center;
+    int *m_radii;
+};
+
+#endif /* MASKCOORDINATEFUNCTION_H_ */
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index 35aaa3ae7b2070c77f95a6203bf1bba301dd7af6..f29880967514b97b7ebab7a5521c8302c49e6cb5 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -18,6 +18,9 @@ void DWBASimulation::init(const Experiment& experiment)
     for (size_t dim=0; dim<detector_dimension; ++dim) {
         m_dwba_intensity.addAxis(detector.getAxis(dim).clone());
     }
+    if (experiment.getOutputData()->getMask()) {
+        m_dwba_intensity.setMask(*experiment.getOutputData()->getMask());
+    }
     Beam beam = experiment.getBeam();
     m_ki = beam.getCentralK();
     kvector_t ki_real(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index 9af59f30cecfa4f8e13dc7570e7076647527e1df..93f631468b80652299f4519e5c1b031933ba8525 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -39,6 +39,8 @@ void GISASExperiment::runSimulation()
     }
 
     m_intensity_map.setAllTo(0.0);
+//    MaskCoordinates circle_mask;
+//    m_intensity_map.addMask(circle_mask);
     if(n_threads_total<0) {
         DWBASimulation *p_dwba_simulation = mp_sample->createDWBASimulation();
         if (!p_dwba_simulation) throw NullPointerException("GISASExperiment::runSimulation() -> No dwba simulation");
@@ -84,8 +86,8 @@ void GISASExperiment::runSimulation()
             delete simulations[i];
             delete threads[i];
         }
-
     }
+//    m_intensity_map.removeAllMasks();
     m_detector.applyDetectorResolution(&m_intensity_map);
 }
 
diff --git a/Core/Algorithms/src/Mask.cpp b/Core/Algorithms/src/Mask.cpp
index f3b6a2be61161104b0f0c74edd442acefc2046fe..b02eff9f8defe0c5cc8ba71e7df783960fb67db4 100644
--- a/Core/Algorithms/src/Mask.cpp
+++ b/Core/Algorithms/src/Mask.cpp
@@ -1,28 +1,34 @@
 #include "Mask.h"
 
+#include <iostream>
+#include <algorithm>
+
 Mask::Mask(Mask* p_submask)
-: m_mask_index(0), m_end_index(0), mp_submask(p_submask)
+: m_own_index(0), m_max_index(0), mp_submask(p_submask)
 {
 }
 
 Mask* Mask::clone() const
 {
-    Mask *p_new_submask = mp_submask ? mp_submask->clone() : 0;
+    Mask *p_new_submask = 0;
+    if (mp_submask) {
+        p_new_submask = mp_submask->clone();
+    }
     Mask *p_new = new Mask(p_new_submask);
-    p_new->m_mask_index = m_mask_index;
-    p_new->setMaxIndex(m_end_index);
+    p_new->m_own_index = m_own_index;
+    p_new->setMaxIndex(m_max_index);
     return p_new;
 }
 
 size_t Mask::getFirstValidIndex()
 {
-    m_mask_index = 0;
+    m_own_index = 0;
     size_t result = 0;
     if (mp_submask) result = mp_submask->getFirstValidIndex();
     while (isMasked(result)) {
         result = nextSubIndex(result);
-        if (result >= m_end_index) {
-            return m_end_index;
+        if (result >= m_max_index) {
+            return m_max_index;
         }
     }
     return result;
@@ -33,8 +39,8 @@ size_t Mask::getNextIndex(size_t total_index)
     size_t result = nextSubIndex(total_index);
     while (isMasked(result)) {
         result = nextSubIndex(result);
-        if (result >= m_end_index) {
-            return m_end_index;
+        if (result >= m_max_index) {
+            return m_max_index;
         }
     }
     return result;
@@ -42,7 +48,7 @@ size_t Mask::getNextIndex(size_t total_index)
 
 void Mask::setMaxIndex(size_t max_index)
 {
-    m_end_index = max_index;
+    m_max_index = max_index;
     if (mp_submask) mp_submask->setMaxIndex(max_index);
 }
 
@@ -54,28 +60,100 @@ bool Mask::isMasked(size_t total_index) const
 
 size_t Mask::nextSubIndex(size_t total_index)
 {
-    size_t result=0;
-    if (!mp_submask) {
-        result = total_index+1;
+    size_t result=total_index;
+    if (total_index < m_max_index) {
+        if (!mp_submask) {
+            ++result;
+        }
+        else {
+            result = mp_submask->getNextIndex(total_index);
+        }
+        ++m_own_index;
     }
     else {
-        result = mp_submask->getNextIndex(total_index);
+        return m_max_index;
     }
-    ++m_mask_index;
     return result;
 }
 
 MaskIndexModulus* MaskIndexModulus::clone() const
 {
-    Mask *p_new_submask = mp_submask ? mp_submask->clone() : 0;
+
+    Mask *p_new_submask = 0;
+    if (mp_submask) {
+        p_new_submask = mp_submask->clone();
+    }
     MaskIndexModulus *p_new = new MaskIndexModulus(m_modulus, m_remainder, p_new_submask);
-    p_new->m_mask_index = m_mask_index;
-    p_new->setMaxIndex(m_end_index);
+    p_new->m_own_index = m_own_index;
+    p_new->setMaxIndex(m_max_index);
     return p_new;
 }
 
 bool MaskIndexModulus::isMasked(size_t total_index) const
 {
     (void)total_index;
-    return m_mask_index % m_modulus != m_remainder;
+    return m_own_index % m_modulus != m_remainder;
+}
+
+MaskCoordinates::MaskCoordinates(size_t rank, int* dims, Mask* p_submask)
+: Mask(p_submask)
+, m_rank(rank)
+, m_dims(0)
+, m_coordinates(0)
+, mp_mask_function(0)
+{
+    m_dims = new int[m_rank];
+    m_coordinates = new int[m_rank];
+    std::copy(dims, dims + m_rank, m_dims);
+    size_t max_index = (m_rank==0) ? 0 : 1;
+    for (size_t i=0; i<m_rank; ++i) {
+        max_index *= m_dims[i];
+    }
+    setMaxIndex(max_index);
+}
+
+MaskCoordinates::~MaskCoordinates()
+{
+    delete m_dims;
+    delete m_coordinates;
+    if (mp_mask_function) delete mp_mask_function;
+}
+
+MaskCoordinates* MaskCoordinates::clone() const
+{
+    Mask *p_new_submask = 0;
+    if (mp_submask) {
+        p_new_submask = mp_submask->clone();
+    }
+    MaskCoordinates *p_new = new MaskCoordinates(m_rank, m_dims, p_new_submask);
+    p_new->m_own_index = m_own_index;
+    p_new->setMaxIndex(m_max_index);
+    if (mp_mask_function) {
+        p_new->setMaskCoordinateFunction(mp_mask_function->clone());
+    }
+    return p_new;
+}
+
+void MaskCoordinates::setMaskCoordinateFunction(MaskCoordinateFunction* p_mask_function)
+{
+    if(mp_mask_function && mp_mask_function != p_mask_function) {
+        delete mp_mask_function;
+    }
+    mp_mask_function = p_mask_function;
+}
+
+bool MaskCoordinates::isMasked(size_t total_index) const
+{
+    if (!mp_mask_function) return false;
+    setCachedCoordinates(total_index);
+    return mp_mask_function->isMasked(m_rank, m_coordinates);
+}
+
+void MaskCoordinates::setCachedCoordinates(size_t index) const
+{
+    int remainder = (int)index;
+    for (size_t i=m_rank; i>0; --i) {
+        m_coordinates[i-1] = remainder % m_dims[i-1];
+        remainder = remainder / m_dims[i-1];
+    }
 }
diff --git a/Core/Algorithms/src/MaskCoordinateFunction.cpp b/Core/Algorithms/src/MaskCoordinateFunction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f87d79eaf5c3d3a6b830660ea94254e121218ad5
--- /dev/null
+++ b/Core/Algorithms/src/MaskCoordinateFunction.cpp
@@ -0,0 +1,106 @@
+#include "MaskCoordinateFunction.h"
+#include "Exceptions.h"
+
+#include <algorithm>
+
+MaskCoordinateFunction::MaskCoordinateFunction(size_t rank)
+: m_rank(rank)
+, m_invert(false)
+{
+}
+
+MaskCoordinateFunction* MaskCoordinateFunction::clone() const
+{
+    MaskCoordinateFunction *p_clone = new MaskCoordinateFunction(m_rank);
+    p_clone->setInvertFlag(m_invert);
+    return p_clone;
+}
+
+bool MaskCoordinateFunction::isMasked(size_t rank, const int* coordinates) const
+{
+    if (rank != m_rank) {
+        throw LogicErrorException("Mask function must have same rank as data structure");
+    }
+    switch (m_invert)
+    {
+    case true:
+        return !isInStandardMaskedArea(coordinates);
+    case false:
+        return isInStandardMaskedArea(coordinates);
+    }
+    throw LogicErrorException("Invalid instruction point");
+}
+
+bool MaskCoordinateFunction::isInStandardMaskedArea(const int* coordinates) const
+{
+    (void)coordinates;
+    return false;
+}
+
+MaskCoordinateRectangleFunction::MaskCoordinateRectangleFunction(size_t rank, const int* minima, const int* maxima)
+: MaskCoordinateFunction(rank)
+, m_minima(0)
+, m_maxima(0)
+{
+    m_minima = new int[m_rank];
+    m_maxima = new int[m_rank];
+    std::copy(minima, minima + m_rank, m_minima);
+    std::copy(maxima, maxima + m_rank, m_maxima);
+}
+
+MaskCoordinateRectangleFunction* MaskCoordinateRectangleFunction::clone() const
+{
+    MaskCoordinateRectangleFunction *p_clone = new MaskCoordinateRectangleFunction(m_rank, m_minima, m_maxima);
+    p_clone->setInvertFlag(m_invert);
+    return p_clone;
+}
+
+MaskCoordinateRectangleFunction::~MaskCoordinateRectangleFunction()
+{
+    delete m_minima;
+    delete m_maxima;
+}
+
+bool MaskCoordinateRectangleFunction::isInStandardMaskedArea(const int* coordinates) const
+{
+    for (size_t i=0; i<m_rank; ++i) {
+        if (m_minima[i] > coordinates[i]) return false;
+        if (m_maxima[i] < coordinates[i]) return false;
+    }
+    return true;
+}
+
+MaskCoordinateEllipseFunction::MaskCoordinateEllipseFunction(size_t rank, const int* center, const int* radii)
+: MaskCoordinateFunction(rank)
+, m_center(0)
+, m_radii(0)
+{
+    m_center = new int[m_rank];
+    m_radii = new int[m_rank];
+    std::copy(center, center + m_rank, m_center);
+    std::copy(radii, radii + m_rank, m_radii);
+}
+
+MaskCoordinateEllipseFunction* MaskCoordinateEllipseFunction::clone() const
+{
+    MaskCoordinateEllipseFunction *p_clone = new MaskCoordinateEllipseFunction(m_rank, m_center, m_radii);
+    p_clone->setInvertFlag(m_invert);
+    return p_clone;
+}
+
+MaskCoordinateEllipseFunction::~MaskCoordinateEllipseFunction()
+{
+    delete m_center;
+    delete m_radii;
+}
+
+bool MaskCoordinateEllipseFunction::isInStandardMaskedArea(const int* coordinates) const
+{
+    double weighted_squares = 0.0;
+    for (size_t i=0; i<m_rank; ++i) {
+        double dist = (double)coordinates[i] - (double)m_center[i];
+        double d_radius = (double)m_radii[i];
+        weighted_squares += dist*dist/d_radius/d_radius;
+    }
+    return weighted_squares <= 1.0;
+}
diff --git a/Core/Core.pro b/Core/Core.pro
index bada46dcbdd7ca9083c0542f544306c44b03fe99..5edadb77df9d97296086bf6844715527bc987a63 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -35,6 +35,7 @@ SOURCES += \
     Algorithms/src/LayerDWBASimulation.cpp \
     Algorithms/src/LocalMonodisperseApproximationStrategy.cpp \
     Algorithms/src/Mask.cpp \
+    Algorithms/src/MaskCoordinateFunction.cpp \
     Algorithms/src/MultiLayerDWBASimulation.cpp \
     Algorithms/src/MultiLayerRoughnessDWBASimulation.cpp \
     Algorithms/src/OpticalFresnel.cpp \
@@ -145,6 +146,7 @@ HEADERS += \
     Algorithms/inc/LayerDWBASimulation.h \
     Algorithms/inc/LocalMonodisperseApproximationStrategy.h \
     Algorithms/inc/Mask.h \
+    Algorithms/inc/MaskCoordinateFunction.h \
     Algorithms/inc/MultiLayerDWBASimulation.h \
     Algorithms/inc/MultiLayerRoughnessDWBASimulation.h \
     Algorithms/inc/OpticalFresnel.h \
diff --git a/Core/Tools/inc/OutputData.h b/Core/Tools/inc/OutputData.h
index 3c37d81ab7af90e205cf77cd2635039e29ea6a2f..ce3ba7b7ea8f1ff19c31248536dddb203b81eea3 100644
--- a/Core/Tools/inc/OutputData.h
+++ b/Core/Tools/inc/OutputData.h
@@ -74,7 +74,7 @@ public:
     T totalSum() const;
 
     // ---------------------------------
-    // external iterators
+    // external iterators (with their possible masking)
     // ---------------------------------
 
     friend class OutputDataIterator<T, OutputData<T> >;
@@ -87,16 +87,28 @@ public:
     typedef OutputDataIterator<const T, const OutputData<T> > const_iterator;
 
     //! return a read/write iterator that points to the first element
-    iterator begin() { return iterator(this); }
+    iterator begin();
 
     //! return a read-only iterator that points to the first element
-    const_iterator begin() const { return const_iterator(this); }
+    const_iterator begin() const;
 
     //! return a read/write iterator that points to the one past last element
-    const iterator end() { return iterator(this, getAllocatedSize()); }
+    iterator end() { return iterator(this, getAllocatedSize()); }
 
     //! return a read-only iterator that points to the one past last element
-    const const_iterator end() const  { return const_iterator(this, getAllocatedSize()); }
+    const_iterator end() const  { return const_iterator(this, getAllocatedSize()); }
+
+    //! get mask that will be used by iterators
+    Mask *getMask() const { return mp_mask; }
+
+    //! set mask (or a stack of masks)
+    void setMask(const Mask &mask);
+
+    //! add mask that will be used by iterators
+    void addMask(const Mask &mask);
+
+    //! remove all masks
+    void removeAllMasks();
 
     // ---------------------------------
     // coordinate and index functions
@@ -171,6 +183,7 @@ private:
 
     std::vector<NamedVectorBase*> m_value_axes;
     LLData<T> *mp_ll_data;
+    Mask *mp_mask;
 };
 
 
@@ -196,6 +209,7 @@ OutputData<double> *getModulusPart(const OutputData<complex_t> &source);
 
 template <class T> OutputData<T>::OutputData()
 : mp_ll_data(0)
+, mp_mask(0)
 {
     allocate();
 }
@@ -214,6 +228,10 @@ template <class T> OutputData<T>* OutputData<T>::clone() const
     }
     (*p_result->mp_ll_data) = *mp_ll_data;
 
+    if (mp_mask) {
+        p_result->mp_mask = mp_mask->clone();
+    }
+
 	return p_result;
 }
 
@@ -225,6 +243,9 @@ template <class T> void OutputData<T>::copyFrom(const OutputData<T> &other)
         addAxis(other.getAxis(i)->clone());
     }
     (*mp_ll_data) = *other.mp_ll_data;
+    if (other.getMask()) {
+        mp_mask = other.getMask()->clone();
+    }
 }
 
 
@@ -273,7 +294,7 @@ inline std::vector<size_t> OutputData<T>::getAllSizes() const
     return result;
 }
 
-template<class T>
+template <class T>
 inline std::vector<T> OutputData<T>::getRawDataVector() const
 {
     std::vector<T> result;
@@ -291,6 +312,50 @@ template <class T> void OutputData<T>::fillRawDataArray(T *destination) const
     return;
 }
 
+template <class T> typename OutputData<T>::iterator OutputData<T>::begin()
+{
+    typename OutputData<T>::iterator result(this);
+    if (mp_mask) {
+        result.setMask(*mp_mask);
+    }
+    return result;
+}
+
+template <class T> typename OutputData<T>::const_iterator OutputData<T>::begin() const
+{
+    typename OutputData<T>::const_iterator result(this);
+    if (mp_mask) {
+        result.setMask(*mp_mask);
+    }
+    return result;
+}
+
+template <class T> void OutputData<T>::setMask(const Mask &mask)
+{
+    if (mp_mask != &mask) {
+        delete mp_mask;
+        mp_mask = mask.clone();
+        mp_mask->setMaxIndex(getAllocatedSize());
+    }
+}
+
+template <class T> void OutputData<T>::addMask(const Mask &mask)
+{
+    if (mask.mp_submask) {
+        throw RuntimeErrorException("One can only add single masks to OutputDataIterator at a time");
+    }
+    Mask *p_old_mask = getMask();
+    mp_mask = mask.clone();
+    mp_mask->mp_submask = p_old_mask;
+    mp_mask->setMaxIndex(getAllocatedSize());
+}
+
+template<class T> void OutputData<T>::removeAllMasks()
+{
+    delete mp_mask;
+    mp_mask = 0;
+}
+
 template<class T> std::vector<int> OutputData<T>::toCoordinates(size_t index) const
 {
     size_t remainder = index;
@@ -368,6 +433,8 @@ template <class T> void OutputData<T>::clear()
     m_value_axes.clear();
     delete mp_ll_data;
     mp_ll_data = 0;
+    delete mp_mask;
+    mp_mask = 0;
 }
 
 template <class T> void OutputData<T>::setAllTo(const T &value)
diff --git a/Core/Tools/inc/OutputDataIterator.h b/Core/Tools/inc/OutputDataIterator.h
index 99799c94808b8318e71db7872b1a1e320174d3bc..a8c87591e28a6b1b8b32dacd5c46c244ce4ac80e 100644
--- a/Core/Tools/inc/OutputDataIterator.h
+++ b/Core/Tools/inc/OutputDataIterator.h
@@ -30,17 +30,23 @@ public:
     template<class TValue2, class TContainer2> OutputDataIterator(
             const OutputDataIterator<TValue2, TContainer2> &other);
 
+    //! non-templated copy construction
+    OutputDataIterator(const OutputDataIterator<TValue, TContainer> &other);
+
     //! templated copy assignment
     template<class TValue2, class TContainer2> OutputDataIterator<TValue, TContainer> &operator=(
             const OutputDataIterator<TValue2, TContainer2> &right);
 
+    //! non-templated copy asssignment
+    OutputDataIterator<TValue, TContainer> &operator=(const OutputDataIterator<TValue, TContainer> &right);
+
     virtual ~OutputDataIterator();
 
     //! prefix increment
     virtual OutputDataIterator<TValue, TContainer> &operator++();
 
     //! postfix increment
-    virtual OutputDataIterator<const TValue, const TContainer> operator++(int);
+    virtual OutputDataIterator<TValue, TContainer> operator++(int);
 
     //! retrieve current element
     virtual TValue &operator*() const;
@@ -52,20 +58,17 @@ public:
     const size_t getIndex() const { return m_current_index; }
 
     //! get container pointer
-    TContainer *const getContainer() const { return mp_output_data; }
+    TContainer *getContainer() const { return mp_output_data; }
 
     //! get mask
-    Mask *const getMask() const { return mp_mask; }
+    Mask *getMask() const { return mp_mask; }
+
+    //! set mask (or a stack of masks)
+    void setMask(const Mask &mask);
 
     //! add mask (also resets index to first available element)
     void addMask(const Mask &mask);
 
-    //! comparison
-    template <class TValue2, class TContainer2> bool operator==(
-            const OutputDataIterator<TValue2, TContainer2> &other);
-    template <class TValue2, class TContainer2> bool operator!=(
-            const OutputDataIterator<TValue2, TContainer2> &other) { return !(*this == other); }
-
 protected:
     virtual void swapContents(OutputDataIterator<TValue, TContainer> &other);
     size_t m_current_index;
@@ -73,6 +76,19 @@ protected:
     Mask *mp_mask;
 };
 
+//! comparison
+template <class TValue1, class TContainer1, class TValue2, class TContainer2> bool operator==(
+        const OutputDataIterator<TValue1, TContainer1> &left,
+        const OutputDataIterator<TValue2, TContainer2> &right) {
+    return left.getContainer()==right.getContainer() && left.getIndex()==right.getIndex();
+}
+
+template <class TValue1, class TContainer1, class TValue2, class TContainer2> bool operator!=(
+        const OutputDataIterator<TValue1, TContainer1> &left,
+        const OutputDataIterator<TValue2, TContainer2> &right) {
+    return !(left == right);
+}
+
 template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer>::OutputDataIterator(
         TContainer *p_output_data, size_t start_at_index)
 : m_current_index(start_at_index)
@@ -87,15 +103,40 @@ OutputDataIterator<TValue, TContainer>::OutputDataIterator(const OutputDataItera
 : m_current_index(0)
 , mp_output_data(0)
 , mp_mask(0)
+{
+    mp_output_data = static_cast<TContainer *>(other.getContainer());
+    m_current_index = other.getIndex();
+    if (other.getMask()) {
+        mp_mask = other.getMask()->clone();
+    }
+}
+
+template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer>::OutputDataIterator(
+        const OutputDataIterator<TValue, TContainer> &other)
+: m_current_index(0)
+, mp_output_data(0)
+, mp_mask(0)
 {
     mp_output_data = other.getContainer();
     m_current_index = other.getIndex();
-    mp_mask = other.getMask() ? other.getMask()->clone() : 0;
+    if (other.getMask()) {
+        mp_mask = other.getMask()->clone();
+    }
 }
 
 template<class TValue, class TContainer>
 template<class TValue2, class TContainer2>
-OutputDataIterator<TValue, TContainer> &OutputDataIterator<TValue, TContainer>::operator=(const OutputDataIterator<TValue2, TContainer2> &right)
+OutputDataIterator<TValue, TContainer> &OutputDataIterator<TValue, TContainer>::operator=(
+        const OutputDataIterator<TValue2, TContainer2> &right)
+{
+    OutputDataIterator<TValue, TContainer> copy(right);
+    swapContents(copy);
+    return *this;
+}
+
+template<class TValue, class TContainer>
+OutputDataIterator<TValue, TContainer> &OutputDataIterator<TValue, TContainer>::operator=(
+        const OutputDataIterator<TValue, TContainer>& right)
 {
     OutputDataIterator<TValue, TContainer> copy(right);
     swapContents(copy);
@@ -104,7 +145,7 @@ OutputDataIterator<TValue, TContainer> &OutputDataIterator<TValue, TContainer>::
 
 template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer>::~OutputDataIterator()
 {
-    delete mp_mask;
+    if (mp_mask) delete mp_mask;
 }
 
 template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer> &OutputDataIterator<TValue, TContainer>::operator++()
@@ -113,16 +154,18 @@ template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer>
         m_current_index = mp_mask->getNextIndex(m_current_index);
     }
     else {
-        ++m_current_index;
+        if (m_current_index<mp_output_data->getAllocatedSize()) {
+            ++m_current_index;
+        }
     }
     return *this;
 }
 
-template<class TValue, class TContainer> OutputDataIterator<const TValue, const TContainer> OutputDataIterator<TValue, TContainer>::operator++(int dummy)
+template<class TValue, class TContainer> OutputDataIterator<TValue, TContainer> OutputDataIterator<TValue, TContainer>::operator++(int dummy)
 {
     (void)dummy;
-    OutputDataIterator<const TValue, const TContainer> result = *this;
-    ++(*this);
+    OutputDataIterator<TValue, TContainer> result(*this);
+    this->operator++();
     return result;
 }
 
@@ -136,10 +179,20 @@ template<class TValue, class TContainer> TValue* OutputDataIterator<TValue, TCon
     return &((*mp_output_data)[m_current_index]);
 }
 
+template<class TValue, class TContainer> void OutputDataIterator<TValue, TContainer>::setMask(const Mask &mask)
+{
+    if (mp_mask != &mask) {
+        delete mp_mask;
+        mp_mask = mask.clone();
+        mp_mask->setMaxIndex(mp_output_data->getAllocatedSize());
+    }
+    m_current_index = mp_mask->getFirstValidIndex();
+}
+
 template<class TValue, class TContainer> void OutputDataIterator<TValue, TContainer>::addMask(const Mask &mask)
 {
     if (mask.mp_submask) {
-        throw RuntimeErrorException("One can only add single masks to OUtputDataIterator at a time");
+        throw RuntimeErrorException("One can only add single masks to OutputDataIterator at a time");
     }
     Mask *p_old_mask = getMask();
     mp_mask = mask.clone();
@@ -148,13 +201,6 @@ template<class TValue, class TContainer> void OutputDataIterator<TValue, TContai
     m_current_index = mp_mask->getFirstValidIndex();
 }
 
-template<class TValue, class TContainer>
-template<class TValue2, class TContainer2> bool OutputDataIterator<TValue, TContainer>::operator==(
-        const OutputDataIterator<TValue2, TContainer2>& other)
-{
-    return mp_output_data==other.getContainer() && m_current_index==other.getIndex();
-}
-
 template<class TValue, class TContainer> void OutputDataIterator<TValue, TContainer>::swapContents(OutputDataIterator<TValue, TContainer>& other)
 {
     std::swap(this->m_current_index, other.m_current_index);
diff --git a/UnitTests/TestCore/MaskTest.h b/UnitTests/TestCore/MaskTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..4199c32569f74153aa5d96b9741dd844d7793caf
--- /dev/null
+++ b/UnitTests/TestCore/MaskTest.h
@@ -0,0 +1,298 @@
+#ifndef MASKTEST_H_
+#define MASKTEST_H_
+
+#include "Mask.h"
+#include "OutputData.h"
+
+#include "gtest/gtest.h"
+
+class MaskTest : public ::testing::Test
+{
+protected:
+    MaskTest();
+    virtual ~MaskTest();
+
+    Mask *p_default_mask;
+    MaskIndexModulus *p_modulo_mask;
+    MaskCoordinates *p_rectangular_mask;
+    MaskCoordinates *p_ellipse_mask;
+};
+
+MaskTest::MaskTest()
+{
+    p_default_mask = new Mask();
+    p_default_mask->setMaxIndex(10);
+
+    p_modulo_mask = new MaskIndexModulus(3, 1);
+    p_modulo_mask->setMaxIndex(10);
+
+    int dims[] = { 5, 5 };
+    p_rectangular_mask = new MaskCoordinates(2, dims);
+    int minima[] = { 1, 1 };
+    int maxima[] = { 2, 2 };
+    MaskCoordinateRectangleFunction *p_rectangle = new MaskCoordinateRectangleFunction(2, minima, maxima);
+    p_rectangular_mask->setMaskCoordinateFunction(p_rectangle);
+
+    p_ellipse_mask = new MaskCoordinates(2, dims);
+    int center[] = { 2, 2 };
+    int radii[] = { 1, 1 };
+    MaskCoordinateEllipseFunction *p_ellipse = new MaskCoordinateEllipseFunction(2, center, radii);
+    p_ellipse_mask->setMaskCoordinateFunction(p_ellipse);
+}
+
+MaskTest::~MaskTest()
+{
+    delete p_default_mask;
+    delete p_modulo_mask;
+    delete p_rectangular_mask;
+    delete p_ellipse_mask;
+}
+
+TEST_F(MaskTest, DefaultFirstIndex)
+{
+    EXPECT_EQ((size_t)0, p_default_mask->getFirstValidIndex());
+}
+
+TEST_F(MaskTest, DefaultNextUninitialized)
+{
+    EXPECT_EQ((size_t)1, p_default_mask->getNextIndex(0));
+    EXPECT_EQ((size_t)10, p_default_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_default_mask->getOwnIndex());
+    EXPECT_EQ((size_t)2, p_default_mask->getNextIndex(1));
+    EXPECT_EQ((size_t)10, p_default_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_default_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, DefaultNextEnd)
+{
+    EXPECT_EQ((size_t)9, p_default_mask->getNextIndex(8));
+    EXPECT_EQ((size_t)10, p_default_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_default_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_default_mask->getNextIndex(9));
+    EXPECT_EQ((size_t)10, p_default_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_default_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_default_mask->getNextIndex(12));
+    EXPECT_EQ((size_t)10, p_default_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_default_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, DefaultClone)
+{
+    Mask *p_mask  = new Mask();
+    p_mask->setMaxIndex(10);
+    EXPECT_EQ((size_t)0, p_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)1, p_mask->getNextIndex(0));
+    Mask *p_clone = p_mask->clone();
+    EXPECT_EQ((size_t)2, p_clone->getNextIndex(1));
+    EXPECT_EQ((size_t)2, p_clone->getOwnIndex());
+    delete p_mask;
+    EXPECT_EQ((size_t)3, p_clone->getNextIndex(2));
+    EXPECT_EQ((size_t)3, p_clone->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_clone->getMaxIndex());
+    delete p_clone;
+}
+
+TEST_F(MaskTest, ModuloFirstIndex)
+{
+    EXPECT_EQ((size_t)1, p_modulo_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+}
+
+TEST_F(MaskTest, ModuloNextUninitialized)
+{
+    EXPECT_EQ((size_t)1, p_modulo_mask->getNextIndex(0));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_modulo_mask->getOwnIndex());
+    EXPECT_EQ((size_t)4, p_modulo_mask->getNextIndex(1));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)4, p_modulo_mask->getOwnIndex());
+    EXPECT_EQ((size_t)7, p_modulo_mask->getNextIndex(4));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)7, p_modulo_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_modulo_mask->getNextIndex(7));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)10, p_modulo_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, ModuloNextEnd)
+{
+    EXPECT_EQ((size_t)9, p_modulo_mask->getNextIndex(8));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_modulo_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_modulo_mask->getNextIndex(9));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_modulo_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_modulo_mask->getNextIndex(12));
+    EXPECT_EQ((size_t)10, p_modulo_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_modulo_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, ModuloClone)
+{
+    Mask *p_mask  = new MaskIndexModulus(3,1);
+    p_mask->setMaxIndex(10);
+    EXPECT_EQ((size_t)1, p_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)4, p_mask->getNextIndex(1));
+    Mask *p_clone = p_mask->clone();
+    EXPECT_EQ((size_t)7, p_clone->getNextIndex(4));
+    EXPECT_EQ((size_t)7, p_clone->getOwnIndex());
+    delete p_mask;
+    EXPECT_EQ((size_t)10, p_clone->getNextIndex(7));
+    EXPECT_EQ((size_t)10, p_clone->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_clone->getMaxIndex());
+    delete p_clone;
+}
+
+TEST_F(MaskTest, RectangularFirstIndex)
+{
+    EXPECT_EQ((size_t)0, p_rectangular_mask->getFirstValidIndex());
+}
+
+TEST_F(MaskTest, RectangularNextUninitialized)
+{
+    EXPECT_EQ((size_t)1, p_rectangular_mask->getNextIndex(0));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_rectangular_mask->getOwnIndex());
+    EXPECT_EQ((size_t)2, p_rectangular_mask->getNextIndex(1));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_rectangular_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, RectangularMaskedElements)
+{
+    EXPECT_EQ((size_t)5, p_rectangular_mask->getNextIndex(4));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_rectangular_mask->getOwnIndex());
+    EXPECT_EQ((size_t)8, p_rectangular_mask->getNextIndex(5));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)4, p_rectangular_mask->getOwnIndex());
+    EXPECT_EQ((size_t)13, p_rectangular_mask->getNextIndex(12));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)5, p_rectangular_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, RectangularNextEnd)
+{
+    EXPECT_EQ((size_t)24, p_rectangular_mask->getNextIndex(23));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_rectangular_mask->getOwnIndex());
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getNextIndex(24));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_rectangular_mask->getOwnIndex());
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getNextIndex(25));
+    EXPECT_EQ((size_t)25, p_rectangular_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_rectangular_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, RectangularClone)
+{
+    Mask *p_mask  = p_rectangular_mask->clone();
+    EXPECT_EQ((size_t)0, p_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)1, p_mask->getNextIndex(0));
+    Mask *p_clone = p_mask->clone();
+    EXPECT_EQ((size_t)2, p_clone->getNextIndex(1));
+    EXPECT_EQ((size_t)2, p_clone->getOwnIndex());
+    delete p_mask;
+    EXPECT_EQ((size_t)3, p_clone->getNextIndex(2));
+    EXPECT_EQ((size_t)3, p_clone->getOwnIndex());
+    EXPECT_EQ((size_t)25, p_clone->getMaxIndex());
+    delete p_clone;
+}
+
+TEST_F(MaskTest, EllipseFirstIndex)
+{
+    EXPECT_EQ((size_t)0, p_ellipse_mask->getFirstValidIndex());
+}
+
+TEST_F(MaskTest, EllipseNextUninitialized)
+{
+    EXPECT_EQ((size_t)1, p_ellipse_mask->getNextIndex(0));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_ellipse_mask->getOwnIndex());
+    EXPECT_EQ((size_t)2, p_ellipse_mask->getNextIndex(1));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_ellipse_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, EllipseMaskedElements)
+{
+    EXPECT_EQ((size_t)6, p_ellipse_mask->getNextIndex(5));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_ellipse_mask->getOwnIndex());
+    EXPECT_EQ((size_t)8, p_ellipse_mask->getNextIndex(6));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)3, p_ellipse_mask->getOwnIndex());
+    EXPECT_EQ((size_t)14, p_ellipse_mask->getNextIndex(12));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)5, p_ellipse_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, EllipseNextEnd)
+{
+    EXPECT_EQ((size_t)24, p_ellipse_mask->getNextIndex(23));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)1, p_ellipse_mask->getOwnIndex());
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getNextIndex(24));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_ellipse_mask->getOwnIndex());
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getNextIndex(25));
+    EXPECT_EQ((size_t)25, p_ellipse_mask->getMaxIndex());
+    EXPECT_EQ((size_t)2, p_ellipse_mask->getOwnIndex());
+}
+
+TEST_F(MaskTest, EllipseClone)
+{
+    Mask *p_mask  = p_ellipse_mask->clone();
+    p_mask->setMaxIndex(10);
+    EXPECT_EQ((size_t)0, p_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)1, p_mask->getNextIndex(0));
+    Mask *p_clone = p_mask->clone();
+    EXPECT_EQ((size_t)2, p_clone->getNextIndex(1));
+    EXPECT_EQ((size_t)2, p_clone->getOwnIndex());
+    delete p_mask;
+    EXPECT_EQ((size_t)3, p_clone->getNextIndex(2));
+    EXPECT_EQ((size_t)3, p_clone->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_clone->getMaxIndex());
+    delete p_clone;
+}
+
+TEST_F(MaskTest, CompoundMask)
+{
+    Mask *p_compound_mask = new MaskIndexModulus(3, 0, p_default_mask->clone());
+    p_compound_mask->setMaxIndex(10);
+    EXPECT_EQ((size_t)0, p_compound_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)3, p_compound_mask->getNextIndex(0));
+    EXPECT_EQ((size_t)6, p_compound_mask->getNextIndex(3));
+    EXPECT_EQ((size_t)9, p_compound_mask->getNextIndex(6));
+    EXPECT_EQ((size_t)10, p_compound_mask->getNextIndex(9));
+    EXPECT_EQ((size_t)10, p_compound_mask->getOwnIndex());
+    EXPECT_EQ((size_t)10, p_compound_mask->getMaxIndex());
+    delete p_compound_mask;
+}
+
+TEST_F(MaskTest, CompoundMaskClone)
+{
+    Mask *p_compound_mask = new MaskIndexModulus(3, 1, p_default_mask->clone());
+    p_compound_mask->setMaxIndex(20);
+    EXPECT_EQ((size_t)1, p_compound_mask->getFirstValidIndex());
+    EXPECT_EQ((size_t)4, p_compound_mask->getNextIndex(1));
+    EXPECT_EQ((size_t)7, p_compound_mask->getNextIndex(4));
+    EXPECT_EQ((size_t)7, p_compound_mask->getOwnIndex());
+    EXPECT_EQ((size_t)20, p_compound_mask->getMaxIndex());
+    Mask *p_compound_clone = p_compound_mask->clone();
+    delete p_compound_mask;
+    EXPECT_EQ((size_t)1, p_compound_clone->getFirstValidIndex());
+    EXPECT_EQ((size_t)4, p_compound_clone->getNextIndex(1));
+    EXPECT_EQ((size_t)7, p_compound_clone->getNextIndex(4));
+    EXPECT_EQ((size_t)7, p_compound_clone->getOwnIndex());
+    EXPECT_EQ((size_t)20, p_compound_clone->getMaxIndex());
+    size_t index = p_compound_clone->getNextIndex(7);
+    for (size_t i=0; i<3; ++i) {
+        index = p_compound_clone->getNextIndex(index);
+    }
+    EXPECT_EQ((size_t)19, p_compound_clone->getOwnIndex());
+    EXPECT_EQ((size_t)20, p_compound_clone->getMaxIndex());
+    delete p_compound_clone;
+}
+
+#endif /* MASKTEST_H_ */
diff --git a/UnitTests/TestCore/NamedVectorTest.h b/UnitTests/TestCore/NamedVectorTest.h
index 16e6ea5523ae3aa4066ed6a9fb1ba0d11af82481..5e04730a0f951e36529b54ea9b0d77f2ca557c9a 100644
--- a/UnitTests/TestCore/NamedVectorTest.h
+++ b/UnitTests/TestCore/NamedVectorTest.h
@@ -3,6 +3,8 @@
 
 #include "NamedVector.h"
 
+#include "gtest/gtest.h"
+
 class NamedVectorTest : public ::testing::Test
 {
 protected:
@@ -19,4 +21,25 @@ NamedVectorTest::NamedVectorTest()
 {
 }
 
+TEST_F(NamedVectorTest, DefaultIsEmpty)
+{
+    EXPECT_EQ((size_t)0, floatAngleVector.getSize());
+}
+
+TEST_F(NamedVectorTest, AddElementsToEmpty)
+{
+    floatAngleVector.push_back(1.0f);
+    ASSERT_EQ((size_t)1, floatAngleVector.getSize());
+    EXPECT_FLOAT_EQ(1.0f, floatAngleVector[0]);
+}
+
+TEST_F(NamedVectorTest, ExtendedConstructor)
+{
+    ASSERT_EQ((size_t)100, doubleLengthVector.getSize());
+    EXPECT_DOUBLE_EQ(0.0, doubleLengthVector[0]);
+    EXPECT_DOUBLE_EQ(0.1, doubleLengthVector[1]);
+    EXPECT_DOUBLE_EQ(6.5, doubleLengthVector[65]);
+    EXPECT_DOUBLE_EQ(9.9, doubleLengthVector[99]);
+}
+
 #endif // NAMEDVECTORTEST_H
diff --git a/UnitTests/TestCore/OutputDataIteratorTest.h b/UnitTests/TestCore/OutputDataIteratorTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed594499c859dc206400103055dc4f8553798176
--- /dev/null
+++ b/UnitTests/TestCore/OutputDataIteratorTest.h
@@ -0,0 +1,137 @@
+#ifndef OUTPUTDATAITERATORTEST_H_
+#define OUTPUTDATAITERATORTEST_H_
+
+#include "OutputDataIterator.h"
+
+#include "gtest/gtest.h"
+
+class OutputDataIteratorTest : public ::testing::Test
+{
+protected:
+    OutputDataIteratorTest();
+    virtual ~OutputDataIteratorTest();
+
+    OutputData<double> *p_data;
+    Mask *p_default_mask;
+    Mask *p_modulo_mask;
+};
+
+OutputDataIteratorTest::OutputDataIteratorTest()
+{
+    p_data = new OutputData<double>();
+    int *dims = new int[2];
+    dims[0] = 3;
+    dims[1] = 5;
+    p_data->setAxisSizes(2, dims);
+    double value = 0.0;
+    OutputData<double>::iterator it = p_data->begin();
+    while (it != p_data->end()) {
+        *it = value;
+        value += 1.0;
+        ++it;
+    }
+
+    p_default_mask = new Mask();
+
+    p_modulo_mask = new MaskIndexModulus(3, 1);
+
+}
+
+OutputDataIteratorTest::~OutputDataIteratorTest()
+{
+    delete p_data;
+}
+
+TEST_F(OutputDataIteratorTest, IterateWithoutMask)
+{
+    OutputData<double>::iterator it = p_data->begin();
+    EXPECT_EQ(0.0, *it);
+    for (size_t i=0; i<14; ++i) {
+        ++it;
+    }
+    EXPECT_DOUBLE_EQ(14.0, *it);
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+}
+
+TEST_F(OutputDataIteratorTest, ConstIterateWithoutMask)
+{
+    OutputData<double>::const_iterator it = p_data->begin();
+    EXPECT_EQ(0.0, *it);
+    for (size_t i=0; i<14; ++i) {
+        ++it;
+    }
+    EXPECT_DOUBLE_EQ(14.0, *it);
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+}
+
+TEST_F(OutputDataIteratorTest, IterateDefaultMask)
+{
+    p_data->addMask(*p_default_mask);
+    OutputData<double>::iterator it = p_data->begin();
+    EXPECT_EQ(0.0, *it);
+    for (size_t i=0; i<14; ++i) {
+        ++it;
+    }
+    EXPECT_DOUBLE_EQ(14.0, *it);
+    it++;
+    EXPECT_EQ(it, p_data->end());
+    it++;
+    EXPECT_EQ(it, p_data->end());
+}
+
+TEST_F(OutputDataIteratorTest, ConstIterateDefaultMask)
+{
+    p_data->addMask(*p_default_mask);
+    OutputData<double>::const_iterator it = p_data->begin();
+    EXPECT_EQ(0.0, *it);
+    for (size_t i=0; i<14; ++i) {
+        ++it;
+    }
+    EXPECT_DOUBLE_EQ(14.0, *it);
+    it++;
+    EXPECT_EQ(it, p_data->end());
+    it++;
+    EXPECT_EQ(it, p_data->end());
+}
+
+TEST_F(OutputDataIteratorTest, IterateCompoundMask)
+{
+    p_data->addMask(*p_default_mask);
+    OutputData<double>::iterator it = p_data->begin();
+    it.addMask(*p_modulo_mask);
+    EXPECT_EQ(1.0, *it);
+    for (size_t i=0; i<4; ++i) {
+        ++it;
+    }
+    EXPECT_DOUBLE_EQ(13.0, *it);
+    it++;
+    EXPECT_EQ(it, p_data->end());
+    it++;
+    EXPECT_EQ(it, p_data->end());
+}
+
+TEST_F(OutputDataIteratorTest, ConstIterateCompoundMask)
+{
+    p_data->addMask(*p_default_mask);
+    OutputData<double>::const_iterator it = p_data->begin();
+    it.addMask(*p_modulo_mask);
+    EXPECT_EQ(1.0, *it);
+    double value;
+    for (size_t i=0; i<4; ++i) {
+        value = *it++;
+    }
+    EXPECT_DOUBLE_EQ(10.0, value);
+    EXPECT_DOUBLE_EQ(13.0, *it);
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+    ++it;
+    EXPECT_EQ(it, p_data->end());
+}
+
+#endif /* OUTPUTDATAITERATORTEST_H_ */
diff --git a/UnitTests/TestCore/OutputDataTest.h b/UnitTests/TestCore/OutputDataTest.h
index a42966899aed1c86c0afbfd562c2676ea7f1606a..50ea14dfba9968f6417ad83f1fc0abf97f1812f6 100644
--- a/UnitTests/TestCore/OutputDataTest.h
+++ b/UnitTests/TestCore/OutputDataTest.h
@@ -3,6 +3,8 @@
 
 #include "OutputData.h"
 
+#include "gtest/gtest.h"
+
 class OutputDataTest : public ::testing::Test
 {
 protected:
@@ -37,4 +39,25 @@ OutputDataTest::~OutputDataTest()
 {
 }
 
+TEST_F(OutputDataTest, SingleElementAfterConstruction)
+{
+    EXPECT_EQ((size_t)1, int_data_0d.getAllocatedSize());
+}
+
+TEST_F(OutputDataTest, SizeAfterAddingAxes)
+{
+    EXPECT_EQ((size_t)1, int_data_0d.getAllocatedSize());
+    EXPECT_EQ((size_t)20, fl_data_1d.getAllocatedSize());
+    EXPECT_EQ((size_t)2000, db_data_3d.getAllocatedSize());
+}
+
+TEST_F(OutputDataTest, DataInitialization)
+{
+    std::vector<int> coordinates;
+    coordinates.push_back(11);
+    coordinates.push_back(4);
+    coordinates.push_back(3);
+    EXPECT_DOUBLE_EQ((double)1143, db_data_3d[db_data_3d.toIndex(coordinates)]);
+}
+
 #endif // OUTPUTDATATEST_H
diff --git a/UnitTests/TestCore/TestCore.pro b/UnitTests/TestCore/TestCore.pro
index e0da4d750f859638bd6ef0b15f3037c7a1aeccb6..30ccbf799c6e7730a9574f296940ae92ab4ab978 100644
--- a/UnitTests/TestCore/TestCore.pro
+++ b/UnitTests/TestCore/TestCore.pro
@@ -12,8 +12,10 @@ include($$PWD/../../shared.pri)
 SOURCES += main.cpp
 
 HEADERS += \
+	MaskTest.h \
     NamedVectorTest.h \
-    OutputDataTest.h
+    OutputDataTest.h \
+    OutputDataIteratorTest.h \
 
 OBJECTS_DIR = obj
 
diff --git a/UnitTests/TestCore/main.cpp b/UnitTests/TestCore/main.cpp
index af04d971889948c36401d3bb10d19d90e2906f87..e27e9996333b3d35efa19c54ed0e9072553db5ae 100644
--- a/UnitTests/TestCore/main.cpp
+++ b/UnitTests/TestCore/main.cpp
@@ -1,48 +1,9 @@
 #include "gtest/gtest.h"
+
+#include "MaskTest.h"
 #include "NamedVectorTest.h"
 #include "OutputDataTest.h"
-
-TEST_F(NamedVectorTest, DefaultIsEmpty)
-{
-    EXPECT_EQ((size_t)0, floatAngleVector.getSize());
-}
-
-TEST_F(NamedVectorTest, AddElementsToEmpty)
-{
-    floatAngleVector.push_back(1.0f);
-    ASSERT_EQ((size_t)1, floatAngleVector.getSize());
-    EXPECT_FLOAT_EQ(1.0f, floatAngleVector[0]);
-}
-
-TEST_F(NamedVectorTest, ExtendedConstructor)
-{
-    ASSERT_EQ((size_t)100, doubleLengthVector.getSize());
-    EXPECT_DOUBLE_EQ(0.0, doubleLengthVector[0]);
-    EXPECT_DOUBLE_EQ(0.1, doubleLengthVector[1]);
-    EXPECT_DOUBLE_EQ(6.5, doubleLengthVector[65]);
-    EXPECT_DOUBLE_EQ(9.9, doubleLengthVector[99]);
-}
-
-TEST_F(OutputDataTest, SingleElementAfterConstruction)
-{
-    EXPECT_EQ((size_t)1, int_data_0d.getAllocatedSize());
-}
-
-TEST_F(OutputDataTest, SizeAfterAddingAxes)
-{
-    EXPECT_EQ((size_t)1, int_data_0d.getAllocatedSize());
-    EXPECT_EQ((size_t)20, fl_data_1d.getAllocatedSize());
-    EXPECT_EQ((size_t)2000, db_data_3d.getAllocatedSize());
-}
-
-TEST_F(OutputDataTest, DataInitialization)
-{
-    std::vector<int> coordinates;
-    coordinates.push_back(11);
-    coordinates.push_back(4);
-    coordinates.push_back(3);
-    EXPECT_DOUBLE_EQ((double)1143, db_data_3d[db_data_3d.toIndex(coordinates)]);
-}
+#include "OutputDataIteratorTest.h"
 
 int main(int argc, char** argv)
 {