Python bindings for chunked array functions

include/vigra/multi_array_chunked.hxx

@@ -139,6 +139,7 @@
 #include "metaprogramming.hxx"
 #include "threading.hxx"
 #include "compression.hxx"
+#include "error.hxx"
 
 #ifdef _WIN32
 # include "windows.h"
@@ -3483,6 +3484,25 @@ class ChunkIterator
     shape_type start_, stop_, chunk_shape_, array_point_;
 };
 
+/**
+    Compare if two ChunkedArray references point to the same object
+    and throw a precondition violation if the comparison evaluates to true.
+*/
+template <unsigned int N, class T1, class T2>
+struct CompareChunkedArrays
+{
+    CompareChunkedArrays(const ChunkedArray<N,T1> &, const ChunkedArray<N,T2> &) {}
+};
+
+template <unsigned int N, class T>
+struct CompareChunkedArrays<N,T,T>
+{
+    CompareChunkedArrays(const ChunkedArray<N,T> & source, const ChunkedArray<N,T> & dest)
+    {
+        vigra_precondition(&source != &dest, "&source must be unequal &destination");
+    }
+};
+
 //@}
 
 } // namespace vigra

include/vigra/multi_chunked.hxx

@@ -0,0 +1,146 @@
+/************************************************************************/
+/*                                                                      */
+/*     Copyright 2016 by Ullrich Koethe and Kevin Kiefer                */
+/*                                                                      */
+/*    This file is part of the VIGRA computer vision library.           */
+/*    The VIGRA Website is                                              */
+/*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
+/*    Please direct questions, bug reports, and contributions to        */
+/*        [email protected]    or                    */
+/*        [email protected]                               */
+/*                                                                      */
+/*    Permission is hereby granted, free of charge, to any person       */
+/*    obtaining a copy of this software and associated documentation    */
+/*    files (the "Software"), to deal in the Software without           */
+/*    restriction, including without limitation the rights to use,      */
+/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
+/*    sell copies of the Software, and to permit persons to whom the    */
+/*    Software is furnished to do so, subject to the following          */
+/*    conditions:                                                       */
+/*                                                                      */
+/*    The above copyright notice and this permission notice shall be    */
+/*    included in all copies or substantial portions of the             */
+/*    Software.                                                         */
+/*                                                                      */
+/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
+/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
+/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
+/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
+/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
+/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
+/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
+/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
+/*                                                                      */
+/************************************************************************/
+
+#ifndef VIGRA_CHUNKED_ARRAY_BLOCKWISE_CONVOLUTION_HXX
+#define VIGRA_CHUNKED_ARRAY_BLOCKWISE_CONVOLUTION_HXX
+
+#include "multi_array_chunked.hxx"
+#include "multi_blocking.hxx"
+#include "multi_blockwise.hxx"
+#include "multi_convolution.hxx"
+#include "threadpool.hxx"
+
+
+namespace vigra {
+
+namespace chunked_blockwise {
+
+/**
+    Helper function to create blockwise parallel filters for chunked arrays.
+*/
+template <unsigned int N, class T1, class T2, class FUNC, class C>
+void chunkedBlockwiseCaller(
+        const ChunkedArray<N,T1> & source,
+        ChunkedArray<N,T2> & dest,
+        FUNC & func,
+        const MultiBlocking<N,C> & blocking,
+        const typename MultiBlocking<N,C>::Shape & borderWidth,
+        const BlockwiseConvolutionOptions<N> & opt
+)
+{
+    typedef typename MultiBlocking<N,C>::BlockWithBorder BlockWithBorder;
+
+    auto begin = blocking.blockWithBorderBegin(borderWidth);
+    auto end = blocking.blockWithBorderEnd(borderWidth);
+
+    parallel_foreach(opt.getNumThreads(), begin, end,
+        [&](const int /*threadId*/, const BlockWithBorder bwb)
+        {
+            // copy input of the block into a new allocated array
+            MultiArray<N,T1> tmp(bwb.border().end() - bwb.border().begin());
+            source.checkoutSubarray(bwb.border().begin(), tmp);
+
+            // get the output as new allocated array
+            MultiArray<N,T2> tmpOut(bwb.core().end() - bwb.core().begin());
+
+            func(tmp, tmpOut, bwb.localCore().begin(), bwb.localCore().end());
+
+            // copy output into destination
+            dest.commitSubarray(bwb.core().begin(), tmpOut);
+        },
+        blocking.numBlocks()
+    );
+}
+
+} // END NAMESPACE chunked_blockwise
+
+/**
+    Overload the functions listed below for chunked arrays.
+    NOTE: Even if the MultiArrayView version may work in place
+    the ChunkedArray overload does not.
+*/
+#define VIGRA_CHUNKED_BLOCKWISE(FUNCTOR, FUNCTION, ORDER, USES_OUTER_SCALE)                 \
+template <unsigned int N, class T1, class T2>                                               \
+void FUNCTION(                                                                              \
+        const ChunkedArray<N,T1> & source,                                                  \
+        ChunkedArray<N,T2> & dest,                                                          \
+        BlockwiseConvolutionOptions<N> const & opt                                          \
+)                                                                                           \
+{                                                                                           \
+    typedef MultiBlocking<N, MultiArrayIndex> Blocking;                                     \
+    typedef typename Blocking::Shape Shape;                                                 \
+\
+    CompareChunkedArrays<N,T1,T2>(source, dest);                                            \
+\
+    const Shape border = blockwise::getBorder(opt, ORDER, USES_OUTER_SCALE);                \
+    const Blocking blocking(source.shape(), opt.template getBlockShapeN<N>());              \
+\
+    BlockwiseConvolutionOptions<N> subOpt(opt);                                             \
+    subOpt.subarray(Shape(0), Shape(0));                                                    \
+\
+    blockwise::FUNCTOR<N> func(subOpt);                                                     \
+    chunked_blockwise::chunkedBlockwiseCaller(source, dest, func, blocking, border, opt);   \
+}
+
+// Reuse the blockwise functors from \<vigra/multi_blockwise.hxx\>
+VIGRA_CHUNKED_BLOCKWISE(GaussianSmoothFunctor,                   gaussianSmoothMultiArray,                   0, false);
+VIGRA_CHUNKED_BLOCKWISE(GaussianGradientFunctor,                 gaussianGradientMultiArray,                 1, false);
+VIGRA_CHUNKED_BLOCKWISE(SymmetricGradientFunctor,                symmetricGradientMultiArray,                1, false);
+VIGRA_CHUNKED_BLOCKWISE(GaussianDivergenceFunctor,               gaussianDivergenceMultiArray,               1, false);
+VIGRA_CHUNKED_BLOCKWISE(HessianOfGaussianFunctor,                hessianOfGaussianMultiArray,                2, false);
+VIGRA_CHUNKED_BLOCKWISE(HessianOfGaussianEigenvaluesFunctor,     hessianOfGaussianEigenvaluesMultiArray,     2, false);
+VIGRA_CHUNKED_BLOCKWISE(HessianOfGaussianFirstEigenvalueFunctor, hessianOfGaussianFirstEigenvalueMultiArray, 2, false);
+VIGRA_CHUNKED_BLOCKWISE(HessianOfGaussianLastEigenvalueFunctor,  hessianOfGaussianLastEigenvalueMultiArray,  2, false);
+VIGRA_CHUNKED_BLOCKWISE(LaplacianOfGaussianFunctor,              laplacianOfGaussianMultiArray,              2, false);
+VIGRA_CHUNKED_BLOCKWISE(GaussianGradientMagnitudeFunctor,        gaussianGradientMagnitudeMultiArray,        1, false);
+VIGRA_CHUNKED_BLOCKWISE(StructureTensorFunctor,                  structureTensorMultiArray,                  1, true);
+
+#undef VIGRA_CHUNKED_BLOCKWISE
+
+
+// Alternative name for backward compatibility.
+template <unsigned int N, class T1, class T2>
+inline void gaussianGradientMagnitude(
+        ChunkedArray<N,T1> const & source,
+        ChunkedArray<N,T2> & dest,
+        BlockwiseConvolutionOptions<N> const & opt
+)
+{
+    gaussianGradientMagnitudeMultiArray(source, dest, opt);
+}
+
+} // END NAMESPACE vigra
+
+#endif // VIGRA_CHUNKED_ARRAY_BLOCKWISE_CONVOLUTION_HXX

include/vigra/python_utility.hxx

@@ -41,6 +41,7 @@
 #include <string>
 #include "vigra/error.hxx"
 #include "vigra/tinyvector.hxx"
+#include "array_vector.hxx"
 
 namespace vigra {
 

test/CMakeLists.txt

@@ -4,6 +4,7 @@ INCLUDE(VigraAddTest)
 
 ADD_SUBDIRECTORY(adjacency_list_graph)
 ADD_SUBDIRECTORY(blockwisealgorithms)
+ADD_SUBDIRECTORY(chunkedconvolution)
 ADD_SUBDIRECTORY(classifier)
 ADD_SUBDIRECTORY(colorspaces)
 ADD_SUBDIRECTORY(convolution)

test/chunkedconvolution/CMakeLists.txt

@@ -0,0 +1,23 @@
+VIGRA_CONFIGURE_THREADING()
+
+if(NOT ${VIGRA_CPP_VERSION})
+    message(FATAL_ERROR
+            "cmake error: VIGRA_CPP_VERSION not defined yet. "
+            "Call VIGRA_DETECT_CPP_VERSION() from the main CMakeLists file.")
+endif()
+
+# c++11 'auto' is used in the test as well as in multi_chunked.hxx
+string(COMPARE LESS ${VIGRA_CPP_VERSION} "201103" NO_CXX11)
+if(NO_CXX11)
+    MESSAGE(STATUS "** WARNING: You are compiling in c++98 mode.")
+    MESSAGE(STATUS "**          Chunked convolution tests will be skipped.")
+    MESSAGE(STATUS "**          Add -std=c++11 to CMAKE_CXX_FLAGS to enable this tests.")
+else()
+    if(NOT THREADING_FOUND)
+        MESSAGE(STATUS "** WARNING: No threading implementation found.")
+        MESSAGE(STATUS "**          Chunked convolution tests will be skipped.")
+    else()
+        VIGRA_ADD_TEST(test_chunked_convolution test_chunked_convolution.cxx
+            LIBRARIES vigraimpex ${THREADING_LIBRARIES})
+    endif()
+endif()

test/chunkedconvolution/test_chunked_convolution.cxx

@@ -0,0 +1,142 @@
+/************************************************************************/
+/*                                                                      */
+/*     Copyright 2016 by Ullrich Koethe and Kevin Kiefer                */
+/*                                                                      */
+/*    This file is part of the VIGRA computer vision library.           */
+/*    The VIGRA Website is                                              */
+/*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
+/*    Please direct questions, bug reports, and contributions to        */
+/*        [email protected]    or                    */
+/*        [email protected]                               */
+/*                                                                      */
+/*    Permission is hereby granted, free of charge, to any person       */
+/*    obtaining a copy of this software and associated documentation    */
+/*    files (the "Software"), to deal in the Software without           */
+/*    restriction, including without limitation the rights to use,      */
+/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
+/*    sell copies of the Software, and to permit persons to whom the    */
+/*    Software is furnished to do so, subject to the following          */
+/*    conditions:                                                       */
+/*                                                                      */
+/*    The above copyright notice and this permission notice shall be    */
+/*    included in all copies or substantial portions of the             */
+/*    Software.                                                         */
+/*                                                                      */
+/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
+/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
+/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
+/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
+/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
+/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
+/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
+/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
+/*                                                                      */
+/************************************************************************/
+
+#include <vigra/multi_array_chunked.hxx>
+#include <vigra/multi_blockwise.hxx>
+#include <vigra/multi_chunked.hxx>
+#include <vigra/multi_array.hxx>
+#include <vigra/tinyvector.hxx>
+#include <vigra/unittest.hxx>
+
+#include <iostream>
+#include "utils.hxx"
+
+using namespace std;
+using namespace vigra;
+
+#define MAKE_TEST(TEST_NAME, FUNCTION_NAME, USES_OUTER_SCALE)           \
+template <unsigned int N, class T1, class T2>                           \
+void TEST_NAME()                                                        \
+{                                                                       \
+    int max = 256;                                                      \
+    TinyVector<int,N> shape(512,512);                                   \
+    TinyVector<int,N> chunk_shape(32,32);                               \
+    TinyVector<int,N> block_shape(30,30);                               \
+\
+    MultiArray<N,T1> source_array(shape);                               \
+    MultiArray<N,T2> dest_array(shape);                                 \
+\
+    ChunkedArrayLazy<N,T1> source_lazy(shape, chunk_shape);             \
+    ChunkedArrayLazy<N,T2> dest_lazy(shape, chunk_shape);               \
+\
+    ChunkedArrayCompressed<N,T1> source_comp(shape, chunk_shape);       \
+    ChunkedArrayCompressed<N,T2> dest_comp(shape, chunk_shape);         \
+\
+    fillRandom(source_array.begin(), source_array.end(), max);          \
+    source_lazy.commitSubarray(TinyVector<int,N>(0,0), source_array);   \
+    source_comp.commitSubarray(TinyVector<int,N>(0,0), source_array);   \
+\
+    BlockwiseConvolutionOptions<N> opt;                                 \
+    opt.stdDev(2.0);                                                    \
+    opt.blockShape(block_shape);                                        \
+    if (opt.getNumThreads() > 4)                                        \
+        opt.numThreads(4);                                              \
+    if (USES_OUTER_SCALE)                                               \
+        opt.outerScale(1.0, 1.0);                                       \
+\
+    FUNCTION_NAME(source_array, dest_array, opt);                       \
+    FUNCTION_NAME(source_lazy, dest_lazy, opt);                         \
+    FUNCTION_NAME(source_comp, dest_comp, opt);                         \
+\
+    compare(dest_lazy.begin(), dest_lazy.end(), dest_array.begin());    \
+    compare(dest_comp.begin(), dest_comp.end(), dest_array.begin());    \
+}
+
+MAKE_TEST(gaussianSmoothTestImpl,                   gaussianSmoothMultiArray,                   false)
+MAKE_TEST(gaussianGradientMagnitudeTestImpl,        gaussianGradientMagnitudeMultiArray,        false)
+MAKE_TEST(laplacianOfGaussianTestImpl,              laplacianOfGaussianMultiArray,              false)
+MAKE_TEST(hessianOfGaussianFirstEigenvalueTestImpl, hessianOfGaussianFirstEigenvalueMultiArray, false)
+MAKE_TEST(hessianOfGaussianLastEigenvalueTestImpl,  hessianOfGaussianLastEigenvalueMultiArray,  false)
+MAKE_TEST(gaussianGradientTestImpl,                 gaussianGradientMultiArray,                 false)
+MAKE_TEST(symmetricGradientTestImpl,                symmetricGradientMultiArray,                false)
+MAKE_TEST(hessianOfGaussianEigenvaluesTestImpl,     hessianOfGaussianEigenvaluesMultiArray,     false)
+MAKE_TEST(gaussianDivergenceTestImpl,               gaussianDivergenceMultiArray,               false)
+MAKE_TEST(hessianOfGaussianTestImpl,                hessianOfGaussianMultiArray,                false)
+MAKE_TEST(structureTensorTestImpl,                  structureTensorMultiArray,                  true)
+
+struct ChunkedConvolutionTestSuite : public test_suite
+{
+    ChunkedConvolutionTestSuite()
+    : test_suite("chunked blockwise convolution test")
+    {
+        typedef float T1;
+        typedef float T2;
+        constexpr int N = 2;
+
+        auto gaussianSmoothTest                     = gaussianSmoothTestImpl<N,T1,T2>;
+        auto gaussianGradientMagnitudeTest          = gaussianGradientMagnitudeTestImpl<N,T1,T2>;
+        auto laplacianOfGaussianTest                = laplacianOfGaussianTestImpl<N,T1,T2>;
+        auto hessianOfGaussianFirstEigenvalueTest   = hessianOfGaussianFirstEigenvalueTestImpl<N,T1,T2>;
+        auto hessianOfGaussianLastEigenvalueTest    = hessianOfGaussianLastEigenvalueTestImpl<N,T1,T2>;
+        auto gaussianGradientTest                   = gaussianGradientTestImpl<N,T1,TinyVector<T2,N> >;
+        auto symmetricGradientTest                  = symmetricGradientTestImpl<N,T1,TinyVector<T2,N> >;
+        auto hessianOfGaussianEigenvaluesTest       = hessianOfGaussianEigenvaluesTestImpl<N,T1,TinyVector<T2,N> >;
+        auto gaussianDivergenceTest                 = gaussianDivergenceTestImpl<N,TinyVector<T1,N>,T2>;
+        auto hessianOfGaussianTest                  = hessianOfGaussianTestImpl<N,T1,TinyVector<T2,int(N*(N+1)/2)> >;
+        auto structureTensorTest                    = structureTensorTestImpl<N,T1,TinyVector<T2,int(N*(N+1)/2)> >;
+
+        add(testCase(gaussianSmoothTest));
+        add(testCase(gaussianGradientMagnitudeTest));
+        add(testCase(laplacianOfGaussianTest));
+        add(testCase(hessianOfGaussianFirstEigenvalueTest));
+        add(testCase(hessianOfGaussianLastEigenvalueTest));
+        add(testCase(gaussianGradientTest));
+        add(testCase(symmetricGradientTest));
+        add(testCase(hessianOfGaussianEigenvaluesTest));
+        add(testCase(gaussianDivergenceTest));
+        add(testCase(hessianOfGaussianTest));
+        add(testCase(structureTensorTest));
+    }
+};
+
+int main(int argc, char** argv)
+{
+    ChunkedConvolutionTestSuite test;
+    int failed = test.run(testsToBeExecuted(argc, argv));
+
+    cout << test.report() << endl;
+
+    return failed != 0;
+}