Container types should support arithmetic operations

[keenancrane]

Container types such as VertexData<T> f( M ) effectively describe a map f: M -> T. If the type T supports arithmetic operations, it would therefore be natural to support the same arithmetic on the container type. E.g.,

VertexData<double> f(M), g(M), h(M);
double c;
h = f + c*g;

would be equivalent to

VertexData<T> f(M), g(M), h(M);
double c;
for( Vertex i : M.vertices() )
{
   h[i] = f[i] + c*g[i];
}

This syntax would greatly improve readability of code for, e.g., linear algebra and optimization, and bring geometry-central in better parity with packages that provide array-style manipulation.

The Eigen vector interoperability does not provide an attractive solution because:

• it requires the user to be familiar with Eigen,
• it still results in messy code, because of all the extra declarations and conversions,
• extra time is spent making copies, and
• most importantly, it can cause data to go out of synch, since methods like MeshData::toVector() make copies by value rather than by reference

This last point is particularly annoying for code that wants to intermingle container-wise and element-wise operations. E.g., a very common use case is: I want to treat a whole container like a vector for doing some kind of gradient descent or line search, but I want to use the nice Geometry methods (such as computing the area of a triangle, say) for expressing the objective function. Currently the only feasible usage pattern is:

• copy data from containers into Eigen vectors
• perform vector-wise operations on Eigen vectors
• copy Eigen values back into containers
• evaluate objective function using idiosyncratic geometry-central code
• repeat

or

• copy data from containers into Eigen vectors
• perform vector-wise operations on Eigen vectors
• copy Eigen values back into containers
• evaluate objective function by writing custom code that operates directly on Eigen vectors
• repeat

(But at this point, why bother using geometry-central? Just do everything in Eigen...)

[keenancrane]

Since MeshData uses an Eigen::Matrix as internal storage, containers can piggyback on Eigen to define arithmetic operations, e.g.,

template <typename E, typename T>
MeshData<E, T>  MeshData<E, T>::operator+( const MeshData<E, T>& u ) const
{
   assert( u.mesh == mesh );
   MeshData<E, T> sum( *mesh );
   sum.data = data + u.data;
   return sum;
}

Though I am not sure how well Eigen supports matrix arithmetic for types beyond standard atomic types (floats, doubles, etc.).

[ctlee]

Native arithmetic operations between MeshData types would be very helpful.

As an aside, it's possible to get a reference to the underlying Eigen objects using MeshData::raw(). This at least prevents the copying out and back.

Another point which may be worth considering is the preservation of an interface to the raw buffers. Eigen provides an Eigen::Map object which binds a raw buffer and allows manipulation with Eigen methods which can be very helpful.

For better or for worse, @cuzhucuncheng and I have some utility functions to map MeshData<gc::Vector3> to Eigen::Map<double, Dynamic, 3, RowMajor> which acts as a normal Eigen object such that we can benefit from Eigen solvers and optimizations (such as lazy evaluation). I think this is ok since Vector3 is POD. Also since the Map object shares the same memory, as MeshData<gc::Vector3>, changes are shared across typed objects. There is a danger where compressing a mesh could invalidate the Eigen::Map though.

[nmwsharp]

Totally agree!

In fact, in a previous version this was implemented, via a somewhat-ridiculous collection of macros to provide operators for each VertexData, FaceData, etc. When I reworked the MeshData<> containers to be templated on an element and avoid ridiculous macros, I just never got around to restoring all the various arithmetic functions :)

I think this is straightforward; just needs to be implemented. Delegating to Eigen like @keenancrane mentions should be the best strategy, and was part of the motivation for the recent change to use Eigen under the hood!

[nmwsharp]

@ctlee that Eigen::Map<> utility sounds very useful! If you'd like to share it we can definitely add it to the utilities somewhere.

[ctlee]

Happy to share the code. Some guidance on where/how to include would be helpful.

There are a few different utilities we use. Given an original data type T made up of k components of type D. Note that in 2 and 3, T is cast into k*D which can be dangerous if certain conditions (checked via static_assert) aren't met. Even if the conditions are met, there's no guarantee that aggregate type T is meaningful when decomposed into Ds. That's up to the user.

1. Normal map (Nx1 -> Nx1): MeshData<T> -> Map<Vector<T, Dynamic 1>>

This is deprecated by the move to Eigen objects in GC. Although it's still useful for binding data in an std::vector or other data types with a raw buffer.

2. Data to matrix (Nx1 -> Nxk): MeshData<T> -> Map<Eigen::Matrix<D, Dynamic, k>

This is what we're using for Vector 3. We check for the validity of the expansion using a few static assertions.

3. Flattening data (Nx1->kNx1): `MeshData -> Map<Eigen::Vector<D, Dynamic (kN), 1>>`

[nmwsharp]

MeshData<> arithmetic is now implemented in 9d04d76765a07e79f52ecfafd56aef79162cf1aa!

@ctlee that makes sense! If you'd like, we could just drop some utility converter functions in a utilities/eigen_interop_helpers.h or something like that. It'd also be nice to add a few examples to the docs, since Eigen template expressions can be scary. If you want to just share any examples here I can Markdown-ify them and get them online.

[nmwsharp]

By the way, one other comment here for posterity:

Unfortunately I ended up not implementing this with Eigen's builtin operations. The reason why is that Eigen's arithmetic templates don't really support operations between matrices of different types, so expressions like

VertexData<float> scales(*mesh, 2.);
VertexData<Vector3> vecs(*mesh, Vector3{1., 2., 3.});
VertexData<Vector3> scaledVecs = scales * vecs ;

wouldn't work out of the box.

Instead, the operations are implemented with a manual for-loop over the underlying data arrays. I once again ended up using a few simple macros, to keep code duplication down. It's definitely much nicer than the old version, at least!

[keenancrane]

Terrific. Another reason it may be necessary to implement these operations natively (i.e., rather than through Eigen) is that I'm not sure Eigen supports all the operators supported in C++. (Say for instance I want to perform a bitwise XOR using the caret operator (^), or have overloaded this operator in a custom type.)

[nmwsharp]

Eigen interop now also added in #42. Thanks!