add max operation to DataContainer

elsa/core/DataContainer.h

@@ -590,6 +590,45 @@ namespace elsa
         }
     }
 
+    /// Element-wise maximum value operation between two operands
+    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
+    auto cwiseMax(LHS const& lhs, RHS const& rhs)
+    {
+        constexpr bool isLHSComplex = isComplex<GetOperandDataType_t<LHS>>;
+        constexpr bool isRHSComplex = isComplex<GetOperandDataType_t<RHS>>;
+
+#ifdef ELSA_CUDA_VECTOR
+        auto cwiseMaxGPU = [](auto const& lhs, auto const& rhs, bool) {
+            return quickvec::cwiseMax(lhs, rhs);
+        };
+#endif
+        auto cwiseMax = [] {
+            if constexpr (isLHSComplex && isRHSComplex) {
+                return [](auto const& left, auto const& right) {
+                    return (left.array().abs().max(right.array().abs())).matrix();
+                };
+            } else if constexpr (isLHSComplex) {
+                return [](auto const& left, auto const& right) {
+                    return (left.array().abs().max(right.array())).matrix();
+                };
+            } else if constexpr (isRHSComplex) {
+                return [](auto const& left, auto const& right) {
+                    return (left.array().max(right.array().abs())).matrix();
+                };
+            } else {
+                return [](auto const& left, auto const& right) {
+                    return (left.array().max(right.array())).matrix();
+                };
+            }
+        }();
+
+#ifdef ELSA_CUDA_VECTOR
+        return Expression{Callables{cwiseMax, cwiseMaxGPU}, lhs, rhs};
+#else
+        return Expression{cwiseMax, lhs, rhs};
+#endif
+    }
+
     /// Element-wise absolute value operation
     template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
     auto cwiseAbs(Operand const& operand)
@@ -654,5 +693,4 @@ namespace elsa
         return Expression{log, operand};
 #endif
     }
-
 } // namespace elsa

elsa/core/ExpressionPredicates.h

@@ -73,6 +73,9 @@ namespace elsa
         using data_t = typename Expression<Callable, Operands...>::data_t;
     };
 
+    template <typename T>
+    using GetOperandDataType_t = typename GetOperandDataType<T>::data_t;
+
     /* Uses the data type used in the first or last operand depending on whether the first operand
      * is an anrithmetic type
      */

elsa/core/tests/test_DataContainer.cpp

@@ -378,6 +378,59 @@ TEST_CASE_TEMPLATE_DEFINE("DataContainer: Testing element-wise access", TestType
                         REQUIRE_UNARY(checkApproxEq(dc[i], randVec(i) / randVec2(i)));
             }
         }
+
+        WHEN("having two containers with real and complex data each")
+        {
+            auto [dcReals1, realsVec1] = generateRandomContainer<real_t>(desc, TestType::handler_t);
+            auto [dcComps1, compsVec1] =
+                generateRandomContainer<std::complex<real_t>>(desc, TestType::handler_t);
+
+            THEN("the element-wise maximum operation works as expected for two real DataContainers")
+            {
+                auto [dcReals2, realsVec2] =
+                    generateRandomContainer<real_t>(desc, TestType::handler_t);
+
+                DataContainer dcCWiseMax = cwiseMax(dcReals1, dcReals2);
+                for (index_t i = 0; i < dcCWiseMax.getSize(); ++i)
+                    REQUIRE_UNARY(
+                        checkApproxEq(dcCWiseMax[i], realsVec1.array().max(realsVec2.array())[i]));
+            }
+
+            THEN("the element-wise maximum operation works as expected for a real and a complex "
+                 "DataContainer")
+            {
+                auto [dcComps2, compsVec2] =
+                    generateRandomContainer<std::complex<real_t>>(desc, TestType::handler_t);
+
+                DataContainer dcCWiseMax = cwiseMax(dcReals1, dcComps2);
+                for (index_t i = 0; i < dcCWiseMax.getSize(); ++i)
+                    REQUIRE_UNARY(checkApproxEq(dcCWiseMax[i],
+                                                realsVec1.array().max(compsVec2.array().abs())[i]));
+            }
+
+            THEN("the element-wise maximum operation works as expected for a complex and a real "
+                 "DataContainer")
+            {
+                auto [dcComps2, compsVec2] =
+                    generateRandomContainer<std::complex<real_t>>(desc, TestType::handler_t);
+
+                DataContainer dcCWiseMax = cwiseMax(dcComps2, dcReals1);
+                for (index_t i = 0; i < dcCWiseMax.getSize(); ++i)
+                    REQUIRE_UNARY(checkApproxEq(dcCWiseMax[i],
+                                                compsVec2.array().abs().max(realsVec1.array())[i]));
+            }
+
+            THEN("the element-wise maximum operation works as expected for two DataContainers")
+            {
+                auto [dcComps2, compsVec2] =
+                    generateRandomContainer<std::complex<real_t>>(desc, TestType::handler_t);
+
+                DataContainer dcCWiseMax = cwiseMax(dcComps1, dcComps2);
+                for (index_t i = 0; i < dcCWiseMax.getSize(); ++i)
+                    REQUIRE_UNARY(checkApproxEq(
+                        dcCWiseMax[i], compsVec1.array().abs().max(compsVec2.array().abs())[i]));
+            }
+        }
     }
 }