diff --git a/source/simple/geom/iterator.hpp b/source/simple/geom/iterator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c3a01b1db0252d563e8075274e05f24ec16d658f
--- /dev/null
+++ b/source/simple/geom/iterator.hpp
@@ -0,0 +1,348 @@
+// key ideas:
+// offsets are the distances from the beginnings of the current row, page, book (just like a pointer is the distance from beginning of memory)
+// step is the distance to get from the end of a row to begin of next row (pitch - flat_size)
+// (pitch here is an image processing jargon, it the image size including padding)
+// flat_size is the total number of elements in given row/page/book, just like offset
+// flat_size can be obtained with (end - begin) of a given range, the relevant calculation(width*height*depth*...) performed to when forming the end pointer, unless it might be otherwise formed, through iteration for example
+// done_mask (it == end) indicates which dimensions are at their end and need to be stepped forward,
+// for example in 2D case
+// (false, false) - we are safely inside the range
+// (true, false) - we reached the end of the row
+// (true, true) - we reached the end of the whole range
+// (false, true) - not possible, invariant broken
+// in boolean context it reduces by conjunction
+// iterator::equality is an optimization of this logic, exploiting the patterns of iteration
+// begin != end is similar and thus the following pattern
+// while (begin != end) ++begin;
+// will only take you till the end of the row,
+// afterwords the step will need to be added (next, advance) according to done_mask
+// in iterator this is also optimized to a plain bool, so the dimensional information is lost
+
+// TODO: compare to Boost.GIL which stores a pointer and pitch, various advantages and disadvantages, reimplement all the relevant examples from there!!!
+
+#ifndef SIMPLE_GEOM_ITERATOR_HPP
+#define SIMPLE_GEOM_ITERATOR_HPP
+
+#include "vector.hpp"
+
+namespace simple::geom
+{
+
+ //TODO: in image processing step/pitch can be in bytes not in sizeof(value_type)s, for memory alignment purposes, which means that offsets will need to be in bytes as well, but maybe that's out of scope here, in general step/pitch is necessary for sub-ranges/regions
+
+ // a more sophisticated and somewhat optimized version
+ // the underlying iterator FlatIt serves as the last offset
+ // so the logical structure is
+ // 2D - vector(row_offset, flat_iterator)
+ // 3D - vector(row_offset, page_offset, flat_iterator)
+ // 4D - vector(row_offset, page_offset, book_offset, flat_iterator)
+ // etc
+ // advantages:
+ // performs better in practice,
+ // disadvantages:
+ // complex implementation
+ // very specific interface
+ template <typename FlatIt, size_t Dimensions = 2>
+ class iterator
+ {
+
+ public:
+
+ using offset_type = vector<std::uintptr_t, Dimensions-1>;
+ using scalar_difference_type = typename std::iterator_traits<FlatIt>::difference_type;
+ using difference_type = vector<scalar_difference_type, Dimensions>;
+ using value_type = typename std::iterator_traits<FlatIt>::value_type;
+ using reference = typename std::iterator_traits<FlatIt>::reference;
+ using pointer = typename std::iterator_traits<FlatIt>::pointer;
+
+ // NOTE: if FlatIt is random access, we are partially random access - adding and subtracting offsets is cheap, but offsets are multidimensional, so the required equivalence with repeated increments doesn't work, that is you can decrement the offset while it's not equal (by conjunction) to default constructed offset, incrementing the pointer, and it will be a valid operation, but not the same as adding the offset directly, the (conjunctive) inequality check will hit a wall on the end of first dimension, while the offset addition can get you across
+ // in more general terms something like binary search wouldn't work, cause halfing the offset is meaningless, when pitch is not zero
+ // TODO: we are actually input if FlatIt is input, output if output, forward if forward, bidirectional if anything higher, and otherwise not an iterator at all... need the meta-programming here and the proper free functions used on FlatIt member instead of operators
+ using iterator_category = std::bidirectional_iterator_tag;
+ //TODO: iterator_concept
+
+ constexpr iterator(FlatIt flat, offset_type offsets)
+ : offsets(offsets), flat(flat)
+ {}
+ constexpr iterator(FlatIt flat)
+ : offsets{}, flat(flat)
+ {}
+
+ struct equality
+ {
+ size_t first_false;
+ operator bool() { return first_false == Dimensions; }
+ };
+
+ [[nodiscard]] constexpr
+ equality operator==(const iterator& other) const
+ {
+ for(size_t i = 0; i < offsets.dimensions; ++i)
+ if(offsets[i] != other.offsets[i])
+ return {i};
+ if(flat != other.flat)
+ return {Dimensions - 1};
+ return {Dimensions};
+ }
+
+ [[nodiscard]] constexpr
+ bool operator!=(const iterator& other) const
+ {
+ return offsets.x() != other.offsets.x();
+ }
+
+ // TODO: relational operators, defining a lexicographic order (can just compare FlatIt), as opposed to partial order of vectors
+
+ [[nodiscard]] constexpr
+ decltype(auto) operator*() const { return *flat; }
+
+ constexpr FlatIt operator->() const { return flat; }
+
+ constexpr iterator& operator++()
+ {
+ ++offsets;
+ ++flat;
+ return *this;
+ }
+
+ [[nodiscard]] constexpr iterator operator++(int)
+ {
+ iterator old = *this;
+ ++(*this);
+ return old;
+ }
+
+ constexpr iterator& operator--()
+ {
+ --offsets;
+ --flat;
+ return *this;
+ }
+
+ [[nodiscard]] constexpr iterator operator--(int)
+ {
+ iterator old = *this;
+ --(*this);
+ return old;
+ }
+
+ [[nodiscard]] constexpr
+ difference_type operator-(const iterator& other) const
+ {
+ return combine_binary_op<scalar_difference_type>(other, std::minus<>{});
+ }
+
+ [[nodiscard]] constexpr
+ friend iterator operator+(iterator it, const difference_type& diff)
+ {
+ it.split_binary_op(diff, std::plus<>());
+ return it;
+ }
+
+ [[nodiscard]] constexpr
+ friend iterator operator+(const difference_type& diff, iterator it)
+ { return it + diff; }
+
+ [[nodiscard]] constexpr
+ friend iterator operator-(iterator it, const difference_type& diff)
+ {
+ it.split_binary_op(diff, std::minus<>());
+ return it;
+ }
+
+ [[nodiscard]] constexpr
+ friend iterator operator-(const difference_type& diff, iterator it)
+ { return it - diff; }
+
+ // TODO: in place arithmetic ops
+
+ constexpr void next(const offset_type& step, size_t dimension)
+ {
+ // dimension == 0 means we are in the beginning or somewhere in the middle, and just need to increment
+ assert(dimension > 0);
+
+ // dimension == Dimensions means we are at the end of the range, nowhere to go
+ assert(dimension < Dimensions);
+
+ for(size_t i = 0; i < dimension; ++i)
+ offsets[i] = 0;
+
+ for(size_t i = dimension; i < offset_type::dimensions; ++i)
+ offsets[i] += step[dimension - 1];
+
+ flat += step[dimension - 1];
+ }
+
+ constexpr void next(const offset_type& step, equality done)
+ {
+ next(step, done.first_false);
+ }
+
+ constexpr void prev(offset_type step, size_t dimension); // TODO
+ constexpr void prev(offset_type step, equality done); // TODO
+
+ constexpr void advance(vector<scalar_difference_type, offset_type::dimensions> step, size_t dimentins); // TODO
+ constexpr void advance(vector<scalar_difference_type, offset_type::dimensions> step, equality done); // TODO
+
+ private:
+ offset_type offsets;
+ FlatIt flat;
+
+ template <typename Result, typename BinaryOp>
+ [[nodiscard]] constexpr
+ vector<Result,Dimensions> combine_binary_op(const iterator& other, BinaryOp&& bop) const
+ {
+ return vector<Result, offset_type::dimensions>(bop(offsets,other.offsets))
+ .template concat(bop(flat,other.flat));
+ }
+
+ template <typename Combined, typename BinaryOp>
+ constexpr void split_binary_op(Combined value, BinaryOp&& bop)
+ {
+ offsets = bop(offsets, value.template first<offset_type::dimensions>());
+ flat = bop(flat, get<offset_type::dimensions>(value));
+ }
+
+ };
+
+ // a more straight forward version
+ // all offset as present, along with the underlying iterator that serves as a reference to the beginning of the range
+ // so the logical structure is
+ // 2D - vector(row_offset, page_offset), begin_iterator
+ // 3D - vector(row_offset, page_offset, book_offset), begin_iterator
+ // 4D - vector(row_offset, page_offset, book_offset, shelf_offset), begin_iterator
+ // etc
+ // advantages:
+ // simple implementation (in theory may be easier for compilers to optimize, and any future geom::vector optimizations will adhere)
+ // inevitably carries more information, which may be of use in some contexts
+ // generic interface
+ // disadvantages:
+ // slow in practice
+ // inevitably carries more information, which makes it more bulky
+ // TODO: compatibility with other version geom::iterator
+ template <typename BeginIt, size_t Dimensions = 2>
+ class vectirator
+ {
+
+ public:
+
+ using offset_type = vector<std::uintptr_t, Dimensions>;
+ using scalar_difference_type = typename std::iterator_traits<BeginIt>::difference_type;
+ using difference_type = vector<scalar_difference_type, Dimensions>;
+ using value_type = typename std::iterator_traits<BeginIt>::value_type;
+ using reference = typename std::iterator_traits<BeginIt>::reference;
+ using pointer = typename std::iterator_traits<BeginIt>::pointer;
+
+ // TODO: constrain BeginIt to not be input or output iterator
+ // TODO: use generic iterator function instead of the operators in the implementation
+ // NOTE: category for this is always bidirectional or partially random access(see geom::iterator), however if BeginIt is forward iterator dereference will be expensive
+ using iterator_category = std::bidirectional_iterator_tag;
+ //TODO: iterator_concept
+
+ constexpr vectirator(BeginIt origin, offset_type offsets)
+ : offsets(offsets), origin(origin)
+ {}
+ constexpr vectirator(BeginIt origin)
+ : offsets{}, origin(origin)
+ {}
+
+ [[nodiscard]] constexpr
+ auto operator==(const vectirator& other) const
+ {
+ return offsets == other.offsets;
+ }
+
+ [[nodiscard]] constexpr
+ auto operator!=(const vectirator& other) const
+ {
+ return to_conjunction(offsets != other.offsets);
+ }
+
+ // TODO: relational operators (as geom::iterator)
+
+ constexpr BeginIt operator->() const { return origin + *(offsets.end()-1); }
+
+ [[nodiscard]] constexpr
+ decltype(auto) operator*() const { return *(operator->()); }
+
+ constexpr vectirator& operator++()
+ {
+ ++offsets;
+ return *this;
+ }
+
+ [[nodiscard]] constexpr vectirator operator++(int)
+ {
+ vectirator old = *this;
+ ++(*this);
+ return old;
+ }
+
+ constexpr vectirator& operator--()
+ {
+ --offsets;
+ return *this;
+ }
+
+ [[nodiscard]] constexpr vectirator operator--(int)
+ {
+ vectirator old = *this;
+ --(*this);
+ return old;
+ }
+
+ [[nodiscard]] constexpr
+ difference_type operator-(const vectirator& other) const
+ {
+ return offsets - other.offsets;
+ }
+
+ [[nodiscard]] constexpr friend
+ vectirator operator+(vectirator it, const difference_type& diff)
+ {
+ it.offsets += diff;
+ return it;
+ }
+
+ [[nodiscard]] constexpr friend
+ vectirator operator+(const difference_type& diff, vectirator it)
+ { return it + diff; }
+
+ [[nodiscard]] constexpr friend
+ vectirator operator-(vectirator it, const difference_type& diff)
+ {
+ it.offsets -= diff;
+ return it;
+ }
+
+ [[nodiscard]] constexpr friend
+ vectirator operator-(const difference_type& diff, vectirator it)
+ { return it - diff; }
+
+ constexpr vectirator& operator+=(const difference_type& diff)
+ { offsets += diff; return *this; }
+ constexpr vectirator& operator-=(const difference_type& diff)
+ { offsets -= diff; return *this; }
+
+ constexpr void next(const offset_type& step, vector<bool, Dimensions> done_mask)
+ {
+ auto begin = std::begin(done_mask);
+ auto end = std::end(done_mask);
+ auto first_false = std::find(begin, end, false) - begin;
+
+ offsets += step[first_false];
+ offsets *= ~done_mask;
+ }
+
+ constexpr void prev(offset_type step, vector<bool, Dimensions> done_mask); // TODO
+
+ constexpr void advance(const offset_type& step, vector<bool, Dimensions> done_mask); // TODO
+
+ private:
+ offset_type offsets;
+ BeginIt origin;
+ };
+
+} // namespace simple::geom
+
+#endif /* end of include guard */
diff --git a/source/simple/geom/pointer.hpp b/source/simple/geom/pointer.hpp
deleted file mode 100644
index cf923ac3d4371809cdee69ef52ccd9229d6416a0..0000000000000000000000000000000000000000
--- a/source/simple/geom/pointer.hpp
+++ /dev/null
@@ -1,28 +0,0 @@
-// pointer into 2D span is (ptr, offset) and dimesntion_step common to all pointers of given range
-// ptr is a pointer to exact memory location
-// offset is the offset from the beginning of the current row
-// dimesntion_step is the step to get from the end of row to begin of next row (pitch - width)
-//
-// begin != end returns a vector<bool,2> (should it be conjunction or disjunction??), indicating in which directions the end has been reached
-// (begin.offset != end.offset, begin.ptr != end.ptr)
-// (false, false) - we are safely inside the range
-// (true, false) - we reached the end of the row
-// (true, true) - we reached the end of the range
-// (false, true) - not possible, invariant broken
-// more complicated with more dimensions, when end is reached not all dimentins maybe true, only the first offset and the ptr
-// ooorrrr, might be much simpler if say y oggset represents the number of pixel in current plane instead of y offset, so it'll go till width*height, and also increment with each step, hmmm... then z would be width*height*depth... it could escalate, however it's the same problem as with size_t or ptrdiff_t, the image still needs to fit in memory so these should not overflow... but definitely a potential limitation in more generic context, where image might not be in memory actually... with this we can get mote consistent and meanigful !=, and ++, however strictly speaking we will be doing more operations (simd though... simd can't do 3 integers and a pointer genious -_- ... but but I can dreaaam) aaand other comparison operators would also wooork, so cool
-//
-// end - begin returns a vector<ptrdiff_t, 2>
-// (end.ptr - begin.ptr, end.offset - begin.offset)
-// what is the meaning of this?
-// da heck, why did i order of ptr and offset, now I have to think which order is better -_-
-//
-// ++ increments both ptr and offset
-// how to jump to next row? maybe with a reset function? or an advance function?
-// even increment not so simple in more dimensions so maye just use a separate step and +=, instead of increment
-//
-// basic missing concept is a vector with different type of coordinates with reduced set of operations and interoperability with normal numeric vectors
-// might be worth it's own template, since the difference with generic vector is vast, though there must be a lowest common denominator
-// TODO: implement, generalize to N dimensions
-// compare to Boost.GIL which stores a pointer and pitch, various advantages and disatvantages, reimplement all the examples there!!!
-// mixing does not make sense, order maybe does, but ptr is always out of order
diff --git a/unit_tests/iterator.cpp b/unit_tests/iterator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..62a0acc39bdc92e005fb0667682aa96ad34d61ac
--- /dev/null
+++ b/unit_tests/iterator.cpp
@@ -0,0 +1,101 @@
+#include "simple/geom/iterator.hpp"
+#include "simple/geom/algorithm.hpp"
+
+#include <cassert>
+#include <memory>
+
+using namespace simple::geom;
+
+template <unsigned width, unsigned height, unsigned step>
+void BasicReadWrite()
+{
+ constexpr vector size{width,height};
+ constexpr unsigned pitch = size.x() + step;
+ auto image_ptr = std::make_unique<std::array<unsigned, pitch * size.y()>>();
+ auto& image = *image_ptr;
+
+ loop(size, [&, position = unsigned()] (auto i) mutable
+ {
+ image[i.y() * pitch + i.x()] = position;
+ ++position;
+ });
+
+ {
+ iterator begin{image.begin(), 0u};
+ iterator end{image.end()-step, size.x()};
+
+ unsigned i = 0;
+ while(true)
+ {
+ auto done = begin == end;
+ if(done) break;
+ if(done.first_false != 0)
+ begin.next(step,done);
+
+ assert(*begin == i);
+ ++begin;
+ ++i;
+ }
+ assert(i == size.x() * size.y());
+ }
+
+ {
+ vectirator begin{image.begin(), {0u, 0u}};
+ vectirator end{image.begin(), {size.x(), size.y()*pitch - step}};
+
+ unsigned i = 0;
+ while(true)
+ {
+ auto done = begin == end;
+ if(done) break;
+ if(to_disjunction(done))
+ begin.next({1u,step}, done);
+
+ assert(*begin == i);
+ ++begin;
+ ++i;
+ }
+ assert(i == size.x() * size.y());
+ }
+
+ {
+ iterator begin{image.begin(), 0u};
+ iterator end{image.end()-step, size.x()};
+
+ unsigned i = 0;
+ while(true)
+ {
+ auto done = begin == end;
+ if(done) break;
+ if(done.first_false != 0)
+ begin.next(step,done);
+
+ while(begin != end)
+ {
+ assert(*begin == i);
+ ++begin;
+ ++i;
+ };
+ }
+
+ assert(i == size.x() * size.y());
+ }
+
+}
+
+int main()
+{
+ BasicReadWrite<1,2,3>();
+ BasicReadWrite<1,1,1>();
+ BasicReadWrite<1,2,0>();
+ BasicReadWrite<1,1,0>();
+ // TODO: figure out what's wrong with these and if needs fixing
+ // BasicReadWrite<1,0,3>();
+ // BasicReadWrite<1,0,0>();
+ BasicReadWrite<0,0,0>();
+ BasicReadWrite<1000,1000,200>();
+ BasicReadWrite<2000,2000,4000>();
+ BasicReadWrite<1000,1000,0>();
+ BasicReadWrite<2000,2000,0>();
+ return 0;
+}