diff --git a/include/vigra/multi_array_chunked.hxx b/include/vigra/multi_array_chunked.hxx index dd8c7f666..18ab21cda 100644 --- a/include/vigra/multi_array_chunked.hxx +++ b/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 +struct CompareChunkedArrays +{ + CompareChunkedArrays(const ChunkedArray &, const ChunkedArray &) {} +}; + +template +struct CompareChunkedArrays +{ + CompareChunkedArrays(const ChunkedArray & source, const ChunkedArray & dest) + { + vigra_precondition(&source != &dest, "&source must be unequal &destination"); + } +}; + //@} } // namespace vigra diff --git a/include/vigra/multi_chunked.hxx b/include/vigra/multi_chunked.hxx new file mode 100644 index 000000000..25d3f8064 --- /dev/null +++ b/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 */ +/* ullrich.koethe@iwr.uni-heidelberg.de or */ +/* vigra@informatik.uni-hamburg.de */ +/* */ +/* 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 +void chunkedBlockwiseCaller( + const ChunkedArray & source, + ChunkedArray & dest, + FUNC & func, + const MultiBlocking & blocking, + const typename MultiBlocking::Shape & borderWidth, + const BlockwiseConvolutionOptions & opt +) +{ + typedef typename MultiBlocking::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 tmp(bwb.border().end() - bwb.border().begin()); + source.checkoutSubarray(bwb.border().begin(), tmp); + + // get the output as new allocated array + MultiArray 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 \ +void FUNCTION( \ + const ChunkedArray & source, \ + ChunkedArray & dest, \ + BlockwiseConvolutionOptions const & opt \ +) \ +{ \ + typedef MultiBlocking Blocking; \ + typedef typename Blocking::Shape Shape; \ +\ + CompareChunkedArrays(source, dest); \ +\ + const Shape border = blockwise::getBorder(opt, ORDER, USES_OUTER_SCALE); \ + const Blocking blocking(source.shape(), opt.template getBlockShapeN()); \ +\ + BlockwiseConvolutionOptions subOpt(opt); \ + subOpt.subarray(Shape(0), Shape(0)); \ +\ + blockwise::FUNCTOR func(subOpt); \ + chunked_blockwise::chunkedBlockwiseCaller(source, dest, func, blocking, border, opt); \ +} + +// Reuse the blockwise functors from \ +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 +inline void gaussianGradientMagnitude( + ChunkedArray const & source, + ChunkedArray & dest, + BlockwiseConvolutionOptions const & opt +) +{ + gaussianGradientMagnitudeMultiArray(source, dest, opt); +} + +} // END NAMESPACE vigra + +#endif // VIGRA_CHUNKED_ARRAY_BLOCKWISE_CONVOLUTION_HXX diff --git a/include/vigra/python_utility.hxx b/include/vigra/python_utility.hxx index 50908f4d8..0fb26fc6c 100644 --- a/include/vigra/python_utility.hxx +++ b/include/vigra/python_utility.hxx @@ -41,6 +41,7 @@ #include #include "vigra/error.hxx" #include "vigra/tinyvector.hxx" +#include "array_vector.hxx" namespace vigra { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1d1d6e0cc..0662a24a6 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/chunkedconvolution/CMakeLists.txt b/test/chunkedconvolution/CMakeLists.txt new file mode 100644 index 000000000..3e527bf79 --- /dev/null +++ b/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() diff --git a/test/chunkedconvolution/test_chunked_convolution.cxx b/test/chunkedconvolution/test_chunked_convolution.cxx new file mode 100644 index 000000000..7b3a68616 --- /dev/null +++ b/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 */ +/* ullrich.koethe@iwr.uni-heidelberg.de or */ +/* vigra@informatik.uni-hamburg.de */ +/* */ +/* 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 +#include +#include +#include +#include +#include + +#include +#include "utils.hxx" + +using namespace std; +using namespace vigra; + +#define MAKE_TEST(TEST_NAME, FUNCTION_NAME, USES_OUTER_SCALE) \ +template \ +void TEST_NAME() \ +{ \ + int max = 256; \ + TinyVector shape(512,512); \ + TinyVector chunk_shape(32,32); \ + TinyVector block_shape(30,30); \ +\ + MultiArray source_array(shape); \ + MultiArray dest_array(shape); \ +\ + ChunkedArrayLazy source_lazy(shape, chunk_shape); \ + ChunkedArrayLazy dest_lazy(shape, chunk_shape); \ +\ + ChunkedArrayCompressed source_comp(shape, chunk_shape); \ + ChunkedArrayCompressed dest_comp(shape, chunk_shape); \ +\ + fillRandom(source_array.begin(), source_array.end(), max); \ + source_lazy.commitSubarray(TinyVector(0,0), source_array); \ + source_comp.commitSubarray(TinyVector(0,0), source_array); \ +\ + BlockwiseConvolutionOptions 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; + auto gaussianGradientMagnitudeTest = gaussianGradientMagnitudeTestImpl; + auto laplacianOfGaussianTest = laplacianOfGaussianTestImpl; + auto hessianOfGaussianFirstEigenvalueTest = hessianOfGaussianFirstEigenvalueTestImpl; + auto hessianOfGaussianLastEigenvalueTest = hessianOfGaussianLastEigenvalueTestImpl; + auto gaussianGradientTest = gaussianGradientTestImpl >; + auto symmetricGradientTest = symmetricGradientTestImpl >; + auto hessianOfGaussianEigenvaluesTest = hessianOfGaussianEigenvaluesTestImpl >; + auto gaussianDivergenceTest = gaussianDivergenceTestImpl,T2>; + auto hessianOfGaussianTest = hessianOfGaussianTestImpl >; + auto structureTensorTest = structureTensorTestImpl >; + + 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; +} diff --git a/test/chunkedconvolution/utils.hxx b/test/chunkedconvolution/utils.hxx new file mode 100644 index 000000000..efce5cfe8 --- /dev/null +++ b/test/chunkedconvolution/utils.hxx @@ -0,0 +1,51 @@ +/************************************************************************/ +/* */ +/* 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 */ +/* ullrich.koethe@iwr.uni-heidelberg.de or */ +/* vigra@informatik.uni-hamburg.de */ +/* */ +/* 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 +#include + +template +void fillRandom(Iter begin, Iter end, int max) +{ + for (; begin != end; ++begin) + *begin = rand() % max; +} + +template +void compare(Iter1 b1, Iter1 e1, Iter2 b2) +{ + for (; b1 != e1; ++b1, ++b2) + shouldEqual(*b1, *b2); +} diff --git a/vigranumpy/lib/__init__.py b/vigranumpy/lib/__init__.py index 4873928d1..737587cb2 100644 --- a/vigranumpy/lib/__init__.py +++ b/vigranumpy/lib/__init__.py @@ -2616,18 +2616,37 @@ def convolutionOptions(blockShape, sigma=None,innerScale=None, outerScale=None, blockwise.convolutionOptions = convolutionOptions blockwise.convOpts = convolutionOptions + + def labelOptions(blockShape, neighborhood="direct", backgroundValue=None, numThreads=cpu_count()): + options = blockwise.BlockwiseLabelOptions() + options.blockShape = blockShape + options.neighborhood = neighborhood + options.numThreads = numThreads + + if backgroundValue is not None: + options.backgroundValue = backgroundValue + + return options + + labelOptions.__module__ = 'vigra.blockwise' + blockwise.labelOptions = labelOptions + blockwise.labelOpts = labelOptions + + def gaussianSmooth(image,options,out=None): out = blockwise._gaussianSmooth(image,options,out) return out gaussianSmooth.__module__ = 'vigra.blockwise' blockwise.gaussianSmooth = gaussianSmooth + def gaussianGradient(image,options,out=None): out = blockwise._gaussianGradient(image,options,out) return out gaussianGradient.__module__ = 'vigra.blockwise' blockwise.gaussianGradient = gaussianGradient + def gaussianGradientMagnitude(image,options,out=None): out = blockwise._gaussianGradientMagnitude(image,options,out) return out @@ -2635,18 +2654,34 @@ def gaussianGradientMagnitude(image,options,out=None): blockwise.gaussianGradientMagnitude = gaussianGradientMagnitude + def gaussianDivergence(image, options, out=None): + out = blockwise._gaussianDivergence(image, options, out) + return out + gaussianDivergence.__module__ = 'vigra.blockwise' + blockwise.gaussianDivergence = gaussianDivergence + + + def hessianOfGaussian(image, options, out=None): + out = blockwise._hessianOfGaussian(image, options, out) + return out + hessianOfGaussian.__module__ = 'vigra.blockwise' + blockwise.hessianOfGaussian = hessianOfGaussian + + def hessianOfGaussianEigenvalues(image,options,out=None): out = blockwise._hessianOfGaussianEigenvalues(image,options,out) return out hessianOfGaussianEigenvalues.__module__ = 'vigra.blockwise' blockwise.hessianOfGaussianEigenvalues = hessianOfGaussianEigenvalues + def hessianOfGaussianFirstEigenvalue(image,options,out=None): out = blockwise._hessianOfGaussianFirstEigenvalue(image,options,out) return out hessianOfGaussianFirstEigenvalue.__module__ = 'vigra.blockwise' blockwise.hessianOfGaussianFirstEigenvalue = hessianOfGaussianFirstEigenvalue + def hessianOfGaussianLastEigenvalue(image,options,out=None): out = blockwise._hessianOfGaussianLastEigenvalue(image,options,out) return out @@ -2654,10 +2689,46 @@ def hessianOfGaussianLastEigenvalue(image,options,out=None): blockwise.hessianOfGaussianLastEigenvalue = hessianOfGaussianLastEigenvalue + def laplacianOfGaussian(image,options,out=None): + out = blockwise._laplacianOfGaussian(image,options,out) + return out + laplacianOfGaussian.__module__ = 'vigra.blockwise' + blockwise.laplacianOfGaussian = laplacianOfGaussian + + + def symmetricGradient(image, options,out=None): + out = blockwise._symmetricGradient(image, options, out) + return out + symmetricGradient.__module__ = 'vigra.blockwise' + blockwise.symmetricGradient = symmetricGradient + + + def structureTensor(image, options, out = None): + out = blockwise._structureTensor(image, options, out) + return out + structureTensor.__module__ = 'vigra.blockwise' + blockwise.structureTensor = structureTensor + + + def unionFindWatersheds(data, options, labels=None): + labels = blockwise._unionFindWatersheds(data, options, labels) + return labels + unionFindWatersheds.__module__ = 'vigra.blockwise' + blockwise.unionFindWatersheds = unionFindWatersheds + + + def labelArray(data, options, labels=None): + labels = blockwise._labelArray(data, options, labels) + return labels + labelArray.__module__ = 'vigra.blockwise' + blockwise.labelArray = labelArray + + _genBlockwiseFunctions() del _genBlockwiseFunctions + def loadBSDGt(filename): import scipy.io as sio matContents = sio.loadmat(filename) diff --git a/vigranumpy/src/core/CMakeLists.txt b/vigranumpy/src/core/CMakeLists.txt index e003402f7..f096f4e80 100644 --- a/vigranumpy/src/core/CMakeLists.txt +++ b/vigranumpy/src/core/CMakeLists.txt @@ -96,8 +96,9 @@ VIGRA_ADD_NUMPY_MODULE(utilities SOURCES VIGRA_ADD_NUMPY_MODULE(blockwise SOURCES + chunked.cxx blockwise.cxx LIBRARIES - ${VIGRANUMPY_THREAD_LIBRARIES} + ${VIGRANUMPY_IMPEX_LIBRARIES} ${VIGRANUMPY_THREAD_LIBRARIES} VIGRANUMPY) diff --git a/vigranumpy/src/core/blockwise.cxx b/vigranumpy/src/core/blockwise.cxx index 4e20c529b..b3590e648 100644 --- a/vigranumpy/src/core/blockwise.cxx +++ b/vigranumpy/src/core/blockwise.cxx @@ -1,6 +1,6 @@ /************************************************************************/ /* */ -/* Copyright 2011 by Ullrich Koethe */ +/* Copyright 2011 by Ullrich Koethe and Kevin Kiefer */ /* */ /* This file is part of the VIGRA computer vision library. */ /* The VIGRA Website is */ @@ -38,237 +38,355 @@ #include #include - #include #include +#include +#include +#include + namespace python = boost::python; +namespace vigra{ -namespace vigra{ +template +NumpyAnyArray pyGaussianSmooth( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape()); + gaussianSmoothMultiArray(source, dest, opt); + return dest; +} - template - NumpyAnyArray pyBlockwiseGaussianSmoothMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - gaussianSmoothMultiArray(source, dest, opt); - return dest; - } - - template - NumpyAnyArray pyBlockwiseGaussianGradientMagnitudeMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - gaussianGradientMagnitudeMultiArray(source, dest, opt); - return dest; - } - - template - NumpyAnyArray pyBlockwiseGaussianGradientMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - gaussianGradientMultiArray(source, dest, opt); - return dest; - } - - template - NumpyAnyArray pyBlockwiseHessianOfGaussianEigenvaluesMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - hessianOfGaussianEigenvaluesMultiArray(source, dest, opt); - return dest; - } - - template - NumpyAnyArray pyBlockwiseHessianOfGaussianFirstEigenvalueMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - hessianOfGaussianFirstEigenvalueMultiArray(source, dest, opt); - return dest; - } - - template - NumpyAnyArray pyBlockwiseHessianOfGaussianLastEigenvalueMultiArray( - const NumpyArray & source, - const BlockwiseConvolutionOptions & opt, - NumpyArray dest - ){ - dest.reshapeIfEmpty(source.taggedShape()); - hessianOfGaussianLastEigenvalueMultiArray(source, dest, opt); - return dest; - } - - - - - template - void defineBlockwiseFilters(){ - //typedef BlockwiseConvolutionOptions Opt; - - python::def("_gaussianSmooth",registerConverters(&pyBlockwiseGaussianSmoothMultiArray), - ( - python::arg("source"), - python::arg("options"), - python::arg("out") = python::object() - ) - ); +template +NumpyAnyArray pyHessianOfGaussianFirstEigenvalue( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape()); + hessianOfGaussianFirstEigenvalueMultiArray(source, dest, opt); + return dest; +} - python::def("_gaussianGradientMagnitude",registerConverters(&pyBlockwiseGaussianGradientMagnitudeMultiArray), - ( - python::arg("source"), - python::arg("options"), - python::arg("out") = python::object() - ) - ); +template +NumpyAnyArray pyGaussianGradientMagnitude( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape()); + gaussianGradientMagnitudeMultiArray(source, dest, opt); + return dest; +} - python::def("_gaussianGradient",registerConverters(&pyBlockwiseGaussianGradientMultiArray >), - ( - python::arg("source"), - python::arg("options"), - python::arg("out") = python::object() - ) - ); +template +NumpyAnyArray pyHessianOfGaussianLastEigenvalue( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape()); + hessianOfGaussianLastEigenvalueMultiArray(source, dest, opt); + return dest; +} - python::def("_hessianOfGaussianEigenvalues",registerConverters(&pyBlockwiseHessianOfGaussianEigenvaluesMultiArray >), - ( - python::arg("source"), - python::arg("options"), - python::arg("out") = python::object() - ) - ); - python::def("_hessianOfGaussianFirstEigenvalue",registerConverters(&pyBlockwiseHessianOfGaussianFirstEigenvalueMultiArray), - ( - python::arg("source"), - python::arg("options"), - python::arg("out") = python::object() - ) - ); - python::def("_hessianOfGaussianLastEigenvalue",registerConverters(&pyBlockwiseHessianOfGaussianLastEigenvalueMultiArray), +template +NumpyAnyArray pyLaplacianOfGaussian( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape()); + laplacianOfGaussianMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pyGaussianDivergence( + const NumpyArray > & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(0)); + gaussianDivergenceMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pyGaussianGradient( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray > dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(DIM)); + gaussianGradientMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pyHessianOfGaussian( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray > dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(DIM*(DIM+1)/2)); + hessianOfGaussianMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pyHessianOfGaussianEigenvalues( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray > dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(DIM)); + hessianOfGaussianEigenvaluesMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pyStructureTensor( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray > dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(DIM*(DIM+1)/2)); + structureTensorMultiArray(source, dest, opt); + return dest; +} + +template +NumpyAnyArray pySymmetricGradient( + const NumpyArray & source, + const BlockwiseConvolutionOptions & opt, + NumpyArray > dest +){ + dest.reshapeIfEmpty(source.taggedShape().setChannelCount(DIM)); + symmetricGradientMultiArray(source, dest, opt); + return dest; +} + +#define VIGRA_NUMPY_FILTER_BINDINGS(NAME_STR, FUNCTOR) \ +python::def(NAME_STR, registerConverters(&FUNCTOR), \ + ( \ + python::arg("source"), \ + python::arg("options"), \ + python::arg("out") = python::object() \ + ) \ +); + +template +void defineBlockwiseFilters(){ + VIGRA_NUMPY_FILTER_BINDINGS("_gaussianSmooth", pyGaussianSmooth) + VIGRA_NUMPY_FILTER_BINDINGS("_gaussianGradient", pyGaussianGradient) + VIGRA_NUMPY_FILTER_BINDINGS("_gaussianGradientMagnitude", pyGaussianGradientMagnitude) + VIGRA_NUMPY_FILTER_BINDINGS("_gaussianDivergence", pyGaussianDivergence) + VIGRA_NUMPY_FILTER_BINDINGS("_hessianOfGaussian", pyHessianOfGaussian) + VIGRA_NUMPY_FILTER_BINDINGS("_hessianOfGaussianEigenvalues", pyHessianOfGaussianEigenvalues) + VIGRA_NUMPY_FILTER_BINDINGS("_hessianOfGaussianFirstEigenvalue",pyHessianOfGaussianFirstEigenvalue) + VIGRA_NUMPY_FILTER_BINDINGS("_hessianOfGaussianLastEigenvalue", pyHessianOfGaussianLastEigenvalue) + VIGRA_NUMPY_FILTER_BINDINGS("_laplacianOfGaussian", pyLaplacianOfGaussian) + VIGRA_NUMPY_FILTER_BINDINGS("_symmetricGradient", pySymmetricGradient) + VIGRA_NUMPY_FILTER_BINDINGS("_structureTensor", pyStructureTensor) +} + +#undef VIGRA_NUMPY_FILTER_BINDINGS + + + +template +python::tuple pyUnionFindWatersheds( + const NumpyArray & data, + const BlockwiseLabelOptions & opt, + NumpyArray labels +){ + labels.reshapeIfEmpty(data.taggedShape()); + auto res = unionFindWatershedsBlockwise(data, labels, opt); + return python::make_tuple(labels, res); +} + +template +void defineUnionFindWatershedsImpl() +{ + python::def("_unionFindWatersheds", registerConverters(&pyUnionFindWatersheds), + ( + python::arg("data"), + python::arg("options"), + python::arg("labels") = python::object() + ) + ); +} + + + +template +python::tuple pyLabelArray( + const NumpyArray & data, + const BlockwiseLabelOptions & opt, + NumpyArray labels +){ + labels.reshapeIfEmpty(data.taggedShape()); + auto res = labelMultiArrayBlockwise(data, labels, opt); + return python::make_tuple(labels, res); +} + +template +void defineLabelArrayImpl() +{ + python::def("_labelArray", registerConverters(&pyLabelArray), + ( + python::arg("data"), + python::arg("options"), + python::arg("labels") = python::object() + ) + ); +} + + + +template +NumpyAnyArray intersectingBlocks( + const MB & mb, + const typename MB::Shape begin, + const typename MB::Shape end, + NumpyArray<1, UInt32> out +){ + std::vector outVec = mb.intersectingBlocks(begin,end); + out.reshapeIfEmpty(typename NumpyArray<1,UInt32>::difference_type(outVec.size())); + std::copy(outVec.begin(),outVec.end(), out.begin()); + return out; +} + +template +python::tuple getBlock( + const MB & mb, + const UInt32 blockIndex +){ + const auto iter = mb.blockBegin(); + const auto & block = iter[blockIndex]; + auto tl = block.begin(); + auto br = block.end(); + return python::make_tuple(tl,br); +} + + +template +python::tuple getBlock2( + const MB & mb, + const typename MB::BlockDesc desc +){ + const auto block = mb.blockDescToBlock(desc); + auto tl = block.begin(); + auto br = block.end(); + return python::make_tuple(tl,br); +} + +template +typename BLOCK::Vector +blockBegin(const BLOCK & b){ + return b.begin(); +} +template +typename BLOCK::Vector +blockEnd(const BLOCK & b){ + return b.end(); +} + +template +typename BLOCK::Vector +blockShape(const BLOCK & b){ + return b.size(); +} + + +template +void defineMultiBlocking(const std::string & clsName){ + + typedef MultiBlocking Blocking; + typedef typename Blocking::Shape Shape; + typedef typename Blocking::Block Block; + + python::class_(clsName.c_str(), python::init()) + .def("intersectingBlocks",registerConverters(&intersectingBlocks), ( - python::arg("source"), - python::arg("options"), + python::arg("begin"), + python::arg("end"), python::arg("out") = python::object() ) - ); - } - - template - NumpyAnyArray intersectingBlocks( - const MB & mb, - const typename MB::Shape begin, - const typename MB::Shape end, - NumpyArray<1, UInt32> out - ){ - std::vector outVec = mb.intersectingBlocks(begin,end); - out.reshapeIfEmpty(typename NumpyArray<1,UInt32>::difference_type(outVec.size())); - std::copy(outVec.begin(),outVec.end(), out.begin()); - return out; - } - - template - python::tuple getBlock( - const MB & mb, - const UInt32 blockIndex - ){ - const auto iter = mb.blockBegin(); - const auto & block = iter[blockIndex]; - auto tl = block.begin(); - auto br = block.end(); - return python::make_tuple(tl,br); - } - - - template - python::tuple getBlock2( - const MB & mb, - const typename MB::BlockDesc desc - ){ - const auto block = mb.blockDescToBlock(desc); - auto tl = block.begin(); - auto br = block.end(); - return python::make_tuple(tl,br); - } - - template - typename BLOCK::Vector - blockBegin(const BLOCK & b){ - return b.begin(); - } - template - typename BLOCK::Vector - blockEnd(const BLOCK & b){ - return b.end(); - } - - template - typename BLOCK::Vector - blockShape(const BLOCK & b){ - return b.size(); - } - - - template - void defineMultiBlocking(const std::string & clsName){ - - typedef MultiBlocking Blocking; - typedef typename Blocking::Shape Shape; - typedef typename Blocking::Block Block; - - python::class_(clsName.c_str(), python::init()) - .def("intersectingBlocks",registerConverters(&intersectingBlocks), - ( - python::arg("begin"), - python::arg("end"), - python::arg("out") = python::object() - ) - ) - .def("__len__", &Blocking::numBlocks) - .def("__getitem__", &getBlock) - .def("__getitem__", &getBlock2) - ; + ) + .def("__len__", &Blocking::numBlocks) + .def("__getitem__", &getBlock) + .def("__getitem__", &getBlock2) + ; + + const std::string blockName = clsName + std::string("Block"); + + python::class_(blockName.c_str()) + .add_property("begin",&blockBegin) + .add_property("end", &blockEnd) + .add_property("shape",&blockShape) + ; +} + + - const std::string blockName = clsName + std::string("Block"); +template +void defineBlockwiseConvolutionOptions(const std::string & clsName){ + + typedef BlockwiseConvolutionOptions Opt; + python::class_(clsName.c_str(), python::init<>()) + .add_property("stdDev", &Opt::getStdDev, &Opt::setStdDev) + //.add_property("scale", &Opt::getScale, &Opt::setScale) + .add_property("innerScale", &Opt::getInnerScale, &Opt::setInnerScale) + .add_property("outerScale", &Opt::getOuterScale, &Opt::setOuterScale) + .add_property("blockShape", &Opt::readBlockShape, &Opt::setBlockShape) + .add_property("numThreads", &Opt::getNumThreads, &Opt::setNumThreads) + ; +} + + +void pySetNeighborhood(BlockwiseLabelOptions& self, std::string str) +{ + if (str == "direct") + self.neighborhood(NeighborhoodType::DirectNeighborhood); + else if (str == "indirect") + self.neighborhood(NeighborhoodType::IndirectNeighborhood); + else + vigra_precondition(false, "Neighborhood must be either 'direct' or 'indirect'."); +} - python::class_(blockName.c_str()) - .add_property("begin",&blockBegin) - .add_property("end", &blockEnd) - .add_property("shape",&blockShape) - ; - } +std::string pyGetNeighborhood(const BlockwiseLabelOptions& self) +{ + if (self.getNeighborhood() == NeighborhoodType::DirectNeighborhood) + return "direct"; + return "indirect"; +} - template - void defineBlockwiseConvolutionOptions(const std::string & clsName){ +void defineBlockwiseLabelOptions() +{ + typedef BlockwiseLabelOptions Opt; + + python::class_("BlockwiseLabelOptions", python::init<>()) + .add_property("blockShape", &Opt::readBlockShape, &Opt::setBlockShape) + .add_property("numThreads", &Opt::getNumThreads, &Opt::setNumThreads) + .add_property("backgroundValue", &Opt::getBackgroundValue, + python::make_function(&Opt::ignoreBackgroundValue, python::return_internal_reference<>())) + .add_property("neighborhood", &pyGetNeighborhood, &pySetNeighborhood) + .def("hasBackgroundValue", &Opt::hasBackgroundValue) + ; +} - typedef BlockwiseConvolutionOptions Opt; - python::class_(clsName.c_str(), python::init<>()) - .add_property("stdDev", &Opt::getStdDev, &Opt::setStdDev) - //.add_property("scale", &Opt::getScale, &Opt::setScale) - .add_property("innerScale", &Opt::getInnerScale, &Opt::setInnerScale) - .add_property("outerScale", &Opt::getOuterScale, &Opt::setOuterScale) - .add_property("blockShape", &Opt::readBlockShape, &Opt::setBlockShape) - .add_property("numThreads", &Opt::getNumThreads, &Opt::setNumThreads) - ; - } +//import from chunked.cxx +void defineChunkedFunctions(); } @@ -293,8 +411,28 @@ BOOST_PYTHON_MODULE_INIT(blockwise) defineBlockwiseConvolutionOptions<2>("BlockwiseConvolutionOptions2D"); defineBlockwiseConvolutionOptions<3>("BlockwiseConvolutionOptions3D"); defineBlockwiseConvolutionOptions<4>("BlockwiseConvolutionOptions4D"); - defineBlockwiseConvolutionOptions<5>("BlockwiseConvolutionOptions4D"); + defineBlockwiseConvolutionOptions<5>("BlockwiseConvolutionOptions5D"); + + defineBlockwiseLabelOptions(); + + defineBlockwiseFilters<2, npy_float32, npy_float32>(); + defineBlockwiseFilters<3, npy_float32, npy_float32>(); + + defineUnionFindWatershedsImpl<2, npy_uint8, npy_uint32>(); + defineUnionFindWatershedsImpl<3, npy_uint8, npy_uint32>(); + defineUnionFindWatershedsImpl<2, npy_uint32, npy_uint32>(); + defineUnionFindWatershedsImpl<3, npy_uint32, npy_uint32>(); + defineUnionFindWatershedsImpl<2, npy_float32, npy_uint32>(); + defineUnionFindWatershedsImpl<3, npy_float32, npy_uint32>(); + + + defineLabelArrayImpl<2, npy_uint8, npy_uint32>(); + defineLabelArrayImpl<3, npy_uint8, npy_uint32>(); + defineLabelArrayImpl<2, npy_uint32, npy_uint32>(); + defineLabelArrayImpl<3, npy_uint32, npy_uint32>(); + defineLabelArrayImpl<2, npy_float32, npy_uint32>(); + defineLabelArrayImpl<3, npy_float32, npy_uint32>(); + - defineBlockwiseFilters<2, float>(); - defineBlockwiseFilters<3, float>(); + defineChunkedFunctions(); } diff --git a/vigranumpy/src/core/chunked.cxx b/vigranumpy/src/core/chunked.cxx new file mode 100644 index 000000000..a0e8c5f5a --- /dev/null +++ b/vigranumpy/src/core/chunked.cxx @@ -0,0 +1,206 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace python = boost::python; + + +namespace vigra { + + +namespace detail { + + template + python::object makeOwningHolder(T* p) + { + typedef python::detail::make_owning_holder owning; + auto nonzero = python_ptr::new_nonzero_reference; + + python_ptr ptr(python::to_python_indirect()(p), nonzero); + python::handle<> h(ptr.release()); + return python::object(h); + } + + template + python::object makeChunkedArray(const ChunkedArray & source) + { + std::string backend = source.backend(); + ChunkedArrayOptions opt; + opt.cacheMax(source.cacheMaxSize()); + ChunkedArray * out = nullptr; + + if (backend == "ChunkedArrayFull") { + out = new ChunkedArrayFull(source.shape(), opt); + } + else if (backend == "ChunkedArrayLazy") { + out = new ChunkedArrayLazy(source.shape(), source.chunkShape(), opt); + } + else if (backend == "ChunkedArrayTmpFile") { + out = new ChunkedArrayTmpFile(source.shape(), source.chunkShape(), opt); + } + else if (backend.find("ChunkedArrayCompressed") != std::string::npos) { + const ChunkedArrayCompressed & sourceComp = dynamic_cast &>(source); + opt.compression(sourceComp.compression_method_); + out = new ChunkedArrayCompressed(sourceComp.shape(), sourceComp.chunkShape(), opt); + } + else if (backend.find("ChunkedArrayHDF5") != std::string::npos) { + vigra_precondition(false, "The 'out' parameter is mandatory for ChunkedArrayHDF5"); + } + else { + vigra_fail("Unable to derive backend from 'input' parameter"); + } + + return makeOwningHolder(out); + } + +} // END NAMESPACE detail + + + +#define PY_VIGRA_FILTERS(FUNCTOR, FUNCTION) \ +template \ +python::object FUNCTOR( \ + const ChunkedArray & source, \ + const BlockwiseConvolutionOptions & opt, \ + python::object out \ +){ \ + if (out == python::object()) out = detail::makeChunkedArray(source); \ + ChunkedArray & dest = python::extract&>(out)(); \ + FUNCTION(source, dest, opt); \ +\ + return out; \ +} + +PY_VIGRA_FILTERS(pyGaussianSmooth, gaussianSmoothMultiArray) +PY_VIGRA_FILTERS(pyGaussianGradient, gaussianGradientMultiArray) +PY_VIGRA_FILTERS(pyGaussianGradientMagnitude, gaussianGradientMagnitudeMultiArray) +PY_VIGRA_FILTERS(pyGaussianDivergence, gaussianDivergenceMultiArray) +PY_VIGRA_FILTERS(pyHessianOfGaussian, hessianOfGaussianMultiArray) +PY_VIGRA_FILTERS(pyHessianOfGaussianEigenvalues, hessianOfGaussianEigenvaluesMultiArray) +PY_VIGRA_FILTERS(pyHessianOfGaussianFirstEigenvalue,hessianOfGaussianFirstEigenvalueMultiArray) +PY_VIGRA_FILTERS(pyHessianOfGaussianLastEigenvalue, hessianOfGaussianLastEigenvalueMultiArray) +PY_VIGRA_FILTERS(pyLaplacianOfGaussian, laplacianOfGaussianMultiArray) +PY_VIGRA_FILTERS(pySymmetricGradient, symmetricGradientMultiArray) +PY_VIGRA_FILTERS(pyStructureTensor, structureTensorMultiArray) + +#undef PY_VIGRA_CONVOLUTION + + +#define PY_VIGRA_BINDINGS(PY_NAME, FUNCTOR, T_IN, T_OUT) \ +python::def(PY_NAME, &FUNCTOR, \ + ( \ + python::arg("source"), \ + python::arg("options"), \ + python::arg("out") = python::object() \ + ) \ +); + +template +void defineChunkedFiltersImpl() +{ +// typedef TinyVector in_type; +// typedef TinyVector out_type; +// typedef TinyVector out_type_2; + + PY_VIGRA_BINDINGS("_gaussianSmooth", pyGaussianSmooth, T_IN, T_OUT) +// PY_VIGRA_BINDINGS("_gaussianGradient", pyGaussianGradient, T_IN, out_type) + PY_VIGRA_BINDINGS("_gaussianGradientMagnitude", pyGaussianGradientMagnitude, T_IN, T_OUT) +// PY_VIGRA_BINDINGS("_gaussianDivergence", pyGaussianDivergence, in_type, T_OUT) +// PY_VIGRA_BINDINGS("_hessianOfGaussian", pyHessianOfGaussian, T_IN, out_type_2) +// PY_VIGRA_BINDINGS("_hessianOfGaussianEigenvalues", pyHessianOfGaussianEigenvalues, T_IN, out_type) + PY_VIGRA_BINDINGS("_hessianOfGaussianFirstEigenvalue", pyHessianOfGaussianFirstEigenvalue, T_IN, T_OUT) + PY_VIGRA_BINDINGS("_hessianOfGaussianLastEigenvalue", pyHessianOfGaussianLastEigenvalue, T_IN, T_OUT) + PY_VIGRA_BINDINGS("_laplacianOfGaussian", pyLaplacianOfGaussian, T_IN, T_OUT) +// PY_VIGRA_BINDINGS("_symmetricGradient", pySymmetricGradient, T_IN, out_type) +// PY_VIGRA_BINDINGS("_structureTensor", pyStructureTensor, T_IN, out_type_2) +} + +#undef PY_VIGRA_BINDINGS + + + +template +python::tuple pyUnionFindWatersheds( + const ChunkedArray & source, + const BlockwiseLabelOptions & opt, + python::object out +){ + if (out == python::object()) out = detail::makeChunkedArray(source); + ChunkedArray & dest = python::extract&>(out)(); + auto res = unionFindWatershedsBlockwise(source, dest, opt); + + return python::make_tuple(out, res); +} + +template +void defineChunkedWatershedsImpl() +{ + python::def("_unionFindWatersheds", &pyUnionFindWatersheds, + ( + python::arg("data"), + python::arg("options"), + python::arg("labels") = python::object() + ) + ); +} + + +template +python::tuple pyLabelArray( + const ChunkedArray & source, + const BlockwiseLabelOptions & opt, + python::object out +){ + if (out == python::object()) out = detail::makeChunkedArray(source); + ChunkedArray & dest = python::extract&>(out)(); + auto res = labelMultiArrayBlockwise(source, dest, opt); + + return python::make_tuple(out, res); +} + +template +void defineChunkedLabelImpl() +{ + python::def("_labelArray", &pyLabelArray, + ( + python::arg("data"), + python::arg("options"), + python::arg("labels") = python::object() + ) + ); +} + + + +void defineChunkedFunctions() +{ + defineChunkedFiltersImpl<2, npy_float32, npy_float32>(); + defineChunkedFiltersImpl<3, npy_float32, npy_float32>(); + + defineChunkedWatershedsImpl<2, npy_uint8, npy_uint32>(); + defineChunkedWatershedsImpl<3, npy_uint8, npy_uint32>(); + defineChunkedWatershedsImpl<2, npy_uint32, npy_uint32>(); + defineChunkedWatershedsImpl<3, npy_uint32, npy_uint32>(); + defineChunkedWatershedsImpl<2, npy_float32, npy_uint32>(); + defineChunkedWatershedsImpl<3, npy_float32, npy_uint32>(); + + + defineChunkedLabelImpl<2, npy_uint8, npy_uint32>(); + defineChunkedLabelImpl<3, npy_uint8, npy_uint32>(); + defineChunkedLabelImpl<2, npy_uint32, npy_uint32>(); + defineChunkedLabelImpl<3, npy_uint32, npy_uint32>(); + defineChunkedLabelImpl<2, npy_float32, npy_uint32>(); + defineChunkedLabelImpl<3, npy_float32, npy_uint32>(); +} + + +} // END NAMESPACE VIGRA diff --git a/vigranumpy/src/core/multi_array_chunked.cxx b/vigranumpy/src/core/multi_array_chunked.cxx index ac6d3056e..102292737 100644 --- a/vigranumpy/src/core/multi_array_chunked.cxx +++ b/vigranumpy/src/core/multi_array_chunked.cxx @@ -242,6 +242,17 @@ int numpyScalarTypeNumber(python::object obj) return typeNum; } +int numpyScalarTypeNumber(PyObject* obj) +{ + PyArray_Descr* dtype; + if(!PyArray_DescrConverter(obj, &dtype)) + return NPY_NOTYPE; + int typeNum = dtype->type_num; + Py_DECREF(dtype); + return typeNum; +} + + template PyObject * ptr_to_python(Array * array, python::object axistags) @@ -287,6 +298,45 @@ construct_ChunkedArrayFull(TinyVector const & shape, python::object dtype, double fill_value, python::object axistags) { + /* + if (PyTuple_Check(dtype.ptr()) && (PyTuple_Size(dtype.ptr()) == 2)) + { + PyObject* type = PyTuple_GetItem(dtype.ptr(), 0); + PyObject* pySize = PyTuple_GetItem(dtype.ptr(), 1); + long size = PyInt_AsLong(pySize); + + if (size == N) + { + switch(numpyScalarTypeNumber(type)) + { + case: NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + case: NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayFull(): unsupported dtype"); + } + } + else if (size == N*(N+1)/2) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayFullImpl >(shape, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayFull(): unsupported dtype"); + } + } + + vigra_precondition(false, "ChunkedArrayFull(): unsupported data shape"); + } + */ switch(numpyScalarTypeNumber(dtype)) { case NPY_UINT8: @@ -319,6 +369,46 @@ construct_ChunkedArrayLazy(TinyVector const & shape, double fill_value, python::object axistags) { + + /* + if (PyTuple_Check(dtype.ptr()) && (PyTuple_Size(dtype.ptr()) == 2)) + { + PyObject* type = PyTuple_GetItem(dtype.ptr(), 0); + PyObject* pySize = PyTuple_GetItem(dtype.ptr(), 1); + long size = PyInt_AsLong(pySize); + + if (size == N) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayLazy(): unsupported dtype"); + } + } + else if (size == N*(N+1)/2) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayLazyImpl >(shape, chunk_shape, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayLazy(): unsupported dtype"); + } + } + + vigra_precondition(false, "ChunkedArrayLazy(): unsupported data shape"); + } + */ switch(numpyScalarTypeNumber(dtype)) { case NPY_UINT8: @@ -355,6 +445,51 @@ construct_ChunkedArrayCompressed(TinyVector const & shape, double fill_value, python::object axistags) { + /* + if (PyTuple_Check(dtype.ptr()) && (PyTuple_Size(dtype.ptr()) == 2)) + { + PyObject* type = PyTuple_GetItem(dtype.ptr(), 0); + PyObject* pySize = PyTuple_GetItem(dtype.ptr(), 1); + long size = PyInt_AsLong(pySize); + + if (size == N) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayCompressed(): unsupported dtype"); + } + } + else if (size == N*(N+1)/2) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayCompressedImpl >( + shape, method, chunk_shape, cache_max, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayCompressed(): unsupported dtype"); + } + } + + vigra_precondition(false, "ChunkedArrayCompressed(): unsupported data shape"); + } + */ switch(numpyScalarTypeNumber(dtype)) { case NPY_UINT8: @@ -394,6 +529,51 @@ construct_ChunkedArrayTmpFile(TinyVector const & shape, double fill_value, python::object axistags) { + /* + if (PyTuple_Check(dtype.ptr()) && (PyTuple_Size(dtype.ptr()) == 2)) + { + PyObject* type = PyTuple_GetItem(dtype.ptr(), 0); + PyObject* pySize = PyTuple_GetItem(dtype.ptr(), 1); + long size = PyInt_AsLong(pySize); + + if (size == N) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayTmpFile(): unsupported dtype"); + } + } + else if (size == N*(N+1)/2) + { + switch(numpyScalarTypeNumber(type)) + { + case NPY_UINT8: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + case NPY_UINT32: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + case NPY_FLOAT32: + return ptr_to_python(construct_ChunkedArrayTmpFileImpl >( + shape, chunk_shape, cache_max, path, fill_value), axistags); + default: + vigra_precondition(false, "ChunkedArrayTmpFile(): unsupported dtype"); + } + } + + vigra_precondition(false, "ChunkedArrayTmpFile(): unsupported data shape"); + } + */ switch(numpyScalarTypeNumber(dtype)) { case NPY_UINT8: @@ -886,7 +1066,12 @@ void defineChunkedArray() defineChunkedArrayImpl<3, npy_float32>(); defineChunkedArrayImpl<4, npy_float32>(); defineChunkedArrayImpl<5, npy_float32>(); - +/* + defineChunkedArrayImpl<2, TinyVector >(); + defineChunkedArrayImpl<2, TinyVector >(); + defineChunkedArrayImpl<3, TinyVector >(); + defineChunkedArrayImpl<3, TinyVector >(); +*/ defineChunkedArrayFactories<2>(); defineChunkedArrayFactories<3>(); defineChunkedArrayFactories<4>();