diff --git a/include/vigra/multi_localminmax.hxx b/include/vigra/multi_localminmax.hxx index eddbab0ae..7a733374c 100644 --- a/include/vigra/multi_localminmax.hxx +++ b/include/vigra/multi_localminmax.hxx @@ -44,60 +44,248 @@ #include "multi_gridgraph.hxx" #include "multi_labeling.hxx" #include "metaprogramming.hxx" +#include "openmp_def.h" namespace vigra { namespace boost_graph { - // Attempt without LValue propmaps, using only the free functions - // to access ReadablePropertyMap (input) and WritablePropertyMap (label) +#ifdef OPENMP + +//TODO: Get rid of break statement inside the parallel region? + +// Attempt without LValue propmaps, using only the free functions +// to access ReadablePropertyMap (input) and WritablePropertyMap (label) +template +unsigned int +localMinMaxGraph(Graph const &g, + T1Map const &src, + T2Map &dest, + typename property_traits::value_type marker, + typename property_traits::value_type threshold, + Compare const &compare, + bool allowAtBorder = true) +{ + typedef typename graph_traits::vertex_iterator graph_scanner; + typedef typename graph_traits::adjacency_iterator neighbor_iterator; + + typedef typename property_traits::value_type T1; + + graph_scanner node, srcend; + neighbor_iterator arc, nbend; + + unsigned int count = 0; + + #pragma omp parallel shared(count) + #pragma omp single + { + tie(node, srcend) = vertices(g); + + for (; node != srcend; ++node) + { + const T1 current = get(src, *node); + + if (!compare(current, threshold)) + continue; + + if(!allowAtBorder && node.atBorder()) + continue; + + tie(arc, nbend) = adjacent_vertices(*node, g); + + #pragma omp task firstprivate(node, current, arc, nbend) + { + for (;arc != nbend; ++arc) + if (!compare(current, get(src, *arc))) + break; + + if (arc == nbend) + { + put(dest, *node, marker); + #pragma omp atomic + ++count; + } + } + } + } + return count; +} + +#else //#ifndef OPENMP + +// Attempt without LValue propmaps, using only the free functions +// to access ReadablePropertyMap (input) and WritablePropertyMap (label) template unsigned int localMinMaxGraph(Graph const &g, + T1Map const &src, + T2Map &dest, + typename property_traits::value_type marker, + typename property_traits::value_type threshold, + Compare const &compare, + bool allowAtBorder = true) +{ + typedef typename graph_traits::vertex_iterator graph_scanner; + typedef typename graph_traits::adjacency_iterator neighbor_iterator; + + typedef typename property_traits::value_type T1; + + graph_scanner node, srcend; + neighbor_iterator arc, nbend; + + unsigned int count = 0; + tie(node, srcend) = vertices(g); + for (; node != srcend; ++node) + { + const T1 current = get(src, *node); + + if (!compare(current, threshold)) + continue; + + if(!allowAtBorder && node.atBorder()) + continue; + + tie(arc, nbend) = adjacent_vertices(*node, g); + for (;arc != nbend; ++arc) + if (!compare(current, get(src, *arc))) + break; + + if (arc == nbend) + { + put(dest, *node, marker); + ++count; + } + } + return count; +} + +#endif //#ifdef OPENMP + +} // namespace boost_graph + +namespace lemon_graph { + +#ifdef OPENMP + +//TODO: Get rid of break statement inside the parallel region? +template +unsigned int +localMinMaxGraph(Graph const &g, T1Map const &src, T2Map &dest, - typename property_traits::value_type marker, - typename property_traits::value_type threshold, + typename T2Map::value_type marker, + typename T1Map::value_type threshold, Compare const &compare, bool allowAtBorder = true) { - typedef typename graph_traits::vertex_iterator graph_scanner; - typedef typename graph_traits::adjacency_iterator neighbor_iterator; - - typedef typename property_traits::value_type T1; - - graph_scanner node, srcend; - neighbor_iterator arc, nbend; + typedef typename Graph::NodeIt graph_scanner; + typedef typename Graph::OutArcIt neighbor_iterator; unsigned int count = 0; - tie(node, srcend) = vertices(g); - for (; node != srcend; ++node) - { - const T1 current = get(src, *node); - if (!compare(current, threshold)) - continue; - - if(!allowAtBorder && node.atBorder()) - continue; - - tie(arc, nbend) = adjacent_vertices(*node, g); - for (;arc != nbend; ++arc) - if (!compare(current, get(src, *arc))) - break; - - if (arc == nbend) + #pragma omp parallel shared(count) + #pragma omp single + { + for (graph_scanner node(g); node != INVALID; ++node) { - put(dest, *node, marker); - ++count; + typename T1Map::value_type current = src[*node]; + + if (!compare(current, threshold)) + continue; + + if(!allowAtBorder && node.atBorder()) + continue; + + #pragma omp task firstprivate(node, current) + { + neighbor_iterator arc(g, node); + + for (; arc != INVALID; ++arc) + if (!compare(current, src[g.target(*arc)])) + break; + + if (arc == INVALID) + { + dest[*node] = marker; + #pragma omp atomic + ++count; + } + } } } return count; } -} // namespace boost_graph -namespace lemon_graph { +//TODO: Get rid of break statement inside the parallel region? +template +unsigned int +extendedLocalMinMaxGraph(Graph const &g, + T1Map const &src, + T2Map &dest, + typename T2Map::value_type marker, + typename T1Map::value_type threshold, + Compare const &compare, + Equal const &equal, + bool allowAtBorder = true) +{ + typename Graph::template NodeMap regions(g); + + int max_region_label = labelGraph(g, src, regions, equal); + + // assume that a region is a extremum until the opposite is proved + std::vector isExtremum(max_region_label+1, (unsigned char)1); + + typedef typename Graph::NodeIt graph_scanner; + typedef typename Graph::OutArcIt neighbor_iterator; + + unsigned int count = max_region_label; + + #pragma omp parallel shared(count) + #pragma omp single + { + for (graph_scanner node(g); node != INVALID; ++node) + { + unsigned int label = regions[*node]; + + if(!isExtremum[label]) + continue; + + typename T1Map::value_type current = src[*node]; + + if (!compare(current, threshold) || + (!allowAtBorder && node.atBorder())) + { + isExtremum[label] = 0; + --count; + continue; + } + + #pragma omp task firstprivate(node, label, current) + { + for (neighbor_iterator arc(g, node); arc != INVALID; ++arc) + { + if (label != regions[g.target(*arc)] && compare(src[g.target(*arc)], current)) + { + isExtremum[label] = 0; + #pragma omp atomic + --count; + break; + } + } + } + } + } + + for (graph_scanner node(g); node != INVALID; ++node) + { + if(isExtremum[regions[*node]]) + dest[*node] = marker; + } + return count; +} + +#else //#ifndef OPENMP template unsigned int @@ -137,6 +325,7 @@ localMinMaxGraph(Graph const &g, return count; } + template unsigned int extendedLocalMinMaxGraph(Graph const &g, @@ -194,6 +383,8 @@ extendedLocalMinMaxGraph(Graph const &g, return count; } +#endif //#ifdef OPENMP + } // namespace lemon_graph template void @@ -84,11 +86,42 @@ initMultiArrayImpl(Iterator s, Shape const & shape, Accessor a, VALUETYPE const & v, MetaInt) { Iterator send = s + shape[N]; + + if(N > 2){ + #pragma omp parallel + #pragma omp single + { + for(; s < send; ++s) + { + #pragma omp task shared(shape, a, v), firstprivate(s) + initMultiArrayImpl(s.begin(), shape, a, v, MetaInt()); + } + } + } + else{ + for(; s < send; ++s) + { + initMultiArrayImpl(s.begin(), shape, a, v, MetaInt()); + } + } +} + +#else + +template +void +initMultiArrayImpl(Iterator s, Shape const & shape, Accessor a, + VALUETYPE const & v, MetaInt) +{ + Iterator send = s + shape[N]; for(; s < send; ++s) { initMultiArrayImpl(s.begin(), shape, a, v, MetaInt()); } } + +#endif //#ifdef OPENMP /** \brief Write a value to every element in a multi-dimensional array. @@ -355,6 +388,46 @@ copyMultiArrayImpl(SrcIterator s, SrcShape const & sshape, SrcAccessor src, } } +#ifdef OPENMP + +template +void +copyMultiArrayImpl(SrcIterator s, SrcShape const & sshape, SrcAccessor src, + DestIterator d, DestShape const & dshape, DestAccessor dest, MetaInt) +{ + DestIterator dend = d + dshape[N]; + if(sshape[N] == 1) + { + for(; d < dend; ++d) + { + copyMultiArrayImpl(s.begin(), sshape, src, d.begin(), dshape, dest, MetaInt()); + } + } + else + { + if(N > 2){ + #pragma omp parallel + #pragma omp single + { + for(; d < dend; ++s, ++d) + { + #pragma omp task shared(sshape, dshape, src, dest), firstprivate(s, d) + copyMultiArrayImpl(s.begin(), sshape, src, d.begin(), dshape, dest, MetaInt()); + } + //No need task wait because parent task does not depend on child task's result + } + }else{ + for(; d < dend; ++s, ++d) + { + copyMultiArrayImpl(s.begin(), sshape, src, d.begin(), dshape, dest, MetaInt()); + } + } + } +} + +#else + template void @@ -377,6 +450,7 @@ copyMultiArrayImpl(SrcIterator s, SrcShape const & sshape, SrcAccessor src, } } } +#endif // #ifdef OPENMP /** \brief Copy a multi-dimensional array. @@ -662,7 +736,53 @@ transformMultiArrayExpandImpl(SrcIterator s, SrcShape const & sshape, SrcAccesso transformLine(s, s + sshape[0], src, d, dest, f); } } - + +#ifdef OPENMP + +template +void +transformMultiArrayExpandImpl(SrcIterator s, SrcShape const & sshape, SrcAccessor src, + DestIterator d, DestShape const & dshape, DestAccessor dest, + Functor const & f, MetaInt) +{ + DestIterator dend = d + dshape[N]; + if(sshape[N] == 1) + { + for(; d < dend; ++d) + { + transformMultiArrayExpandImpl(s.begin(), sshape, src, d.begin(), dshape, dest, + f, MetaInt()); + } + } + else + { + if(N > 2){ + #pragma omp parallel + #pragma omp single + { + for(; d < dend; ++s, ++d) + { + #pragma omp task shared(sshape, dshape, src, dest), firstprivate(s, d) + transformMultiArrayExpandImpl(s.begin(), sshape, src, + d.begin(), dshape, dest, + f, MetaInt()); + } + } + }else{ + for(; d < dend; ++s, ++d) + { + transformMultiArrayExpandImpl(s.begin(), sshape, src, + d.begin(), dshape, dest, + f, MetaInt()); + } + } + } +} + +#else + template @@ -689,6 +809,7 @@ transformMultiArrayExpandImpl(SrcIterator s, SrcShape const & sshape, SrcAccesso } } } +#endif //#ifdef OPENMP template +#include +#include +#include +#include +#include +#include + +#ifdef WITH_BOOST_GRAPH +# include +#endif + + +using namespace vigra; + +template +struct GridGraphAlgorithmTests +{ + int max_threads; + + typedef typename MultiArrayShape::type Shape; + + GridGraphAlgorithmTests() +#ifdef OPENMP + : max_threads(std::max(omp_get_num_procs() / 2, 2)) +#else + : max_threads(1) +#endif + { +#ifdef OPENMP + omp_set_num_threads(max_threads); +#endif + } + + template + void testLocalMinMax() + { + typedef GridGraph Graph; + + Graph g(Shape(3), NType); + typename Graph::template NodeMap src(g), dest(g); + + src[Shape(1)] = 1; + + should(1 == boost_graph::localMinMaxGraph(g, src, dest, 1, -9999, std::greater())); + + shouldEqualSequence(src.begin(), src.end(), dest.begin()); + + dest.init(0); + + should(1 == lemon_graph::localMinMaxGraph(g, src, dest, 1, -9999, std::greater())); + + shouldEqualSequence(src.begin(), src.end(), dest.begin()); + } +}; + +template +struct GridgraphTestSuiteN +: public vigra::test_suite +{ + GridgraphTestSuiteN() + : vigra::test_suite((std::string("Gridgraph Test Dimension ") + vigra::asString(N)).c_str()) + { + add(testCase((&GridGraphAlgorithmTests::template testLocalMinMax))); + } +}; + +struct GridgraphTestSuite +: public vigra::test_suite +{ + GridgraphTestSuite() +#ifdef WITH_BOOST_GRAPH + : vigra::test_suite("Gridgraph BGL Test") +#else + : vigra::test_suite("Gridgraph Test") +#endif + { + USETICTOC + + TIC + add(VIGRA_TEST_SUITE(GridgraphTestSuiteN<2>)); + TOC + + TIC + add(VIGRA_TEST_SUITE(GridgraphTestSuiteN<3>)); + TOC + + TIC + add(VIGRA_TEST_SUITE(GridgraphTestSuiteN<4>)); + TOC + } +}; + +int main(int argc, char **argv) +{ + + GridgraphTestSuite gridgraphTest; + + int failed = gridgraphTest.run(vigra::testsToBeExecuted(argc, argv)); + + std::cout << gridgraphTest.report() << std::endl; + + return (failed != 0); +} diff --git a/test/openmp_multi_pointoperators/CMakeLists.txt b/test/openmp_multi_pointoperators/CMakeLists.txt new file mode 100644 index 000000000..de94c41d3 --- /dev/null +++ b/test/openmp_multi_pointoperators/CMakeLists.txt @@ -0,0 +1 @@ +VIGRA_ADD_TEST(test_openmp_multi_pointoperators test.cxx) \ No newline at end of file diff --git a/test/openmp_multi_pointoperators/test.cxx b/test/openmp_multi_pointoperators/test.cxx new file mode 100644 index 000000000..3e63b660d --- /dev/null +++ b/test/openmp_multi_pointoperators/test.cxx @@ -0,0 +1,176 @@ +/************************************************************************/ +/* */ +/* Copyright 2004 by Ullrich Koethe */ +/* */ +/* 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 "vigra/unittest.hxx" +#include "vigra/error.hxx" +#include "vigra/multi_array.hxx" +#include "vigra/multi_pointoperators.hxx" +#include "vigra/multi_math.hxx" +#include "vigra/openmp_vigra.h" +#include "vigra/timing.hxx" + +using namespace std; +using namespace vigra; + +#define DIMENSION_LENGTH 60 + +struct OpenMPMultiArrayPointOperatorsTest { + typedef float PixelType; + typedef MultiArray<5, PixelType> Image5D; + typedef Image5D::difference_type Size5; + Image5D img; + int max_threads; + + OpenMPMultiArrayPointOperatorsTest() : + img(Size5( DIMENSION_LENGTH, DIMENSION_LENGTH, + DIMENSION_LENGTH, DIMENSION_LENGTH, DIMENSION_LENGTH)) +#ifdef OPENMP + , max_threads(std::max(omp_get_num_procs() / 2, 2)) +#else + , max_threads(1) +#endif + { +#ifdef OPENMP + omp_set_num_threads(max_threads); +#endif + int i; + PixelType c = 0.1f; + for (i = 0; i < img.elementCount(); ++i, ++c) + img.data()[i] = c; + } + + +//TODO: test function returns number of threads involved in parallel region +//TODO: Meaningful runtime. + void testInit() + { + Image5D res(img.shape()); + const Image5D::value_type ini = 1.1f; + + should(res.shape() == Size5( DIMENSION_LENGTH, DIMENSION_LENGTH, + DIMENSION_LENGTH, DIMENSION_LENGTH, + DIMENSION_LENGTH )); + + USETICTOC + TIC + initMultiArray(res, ini); + TOC + + using namespace multi_math; + should(all(res == ini)); + + initMultiArray(res, 2.2f); + should(all(res == 2.2f)); + + res = 3.3f; + should(all(res == 3.3f)); + + res.init(4.4f); + should(all(res == 4.4f)); + } + + + void testCopy() { + + Image5D res(img.shape(), 1.0); + + USETICTOC + TIC + vigra::copyMultiArray(img, res); + TOC + + should(img == res); + } + + + void testTransform() + { + using namespace vigra::functor; + + Image5D res(img.shape()); + + USETICTOC + TIC + transformMultiArray(img, res, Arg1() + Arg1()); + TOC + + using namespace multi_math; + should(all(2.0*img == res)); + } + + + void testTransformOuterExpand() + { + using namespace functor; + + Image5D res(img.shape()); + + USETICTOC + TIC + transformMultiArray(img.subarray(Size5(0,0,0,0,0), Size5(DIMENSION_LENGTH,1,1,1,1)), res, + Arg1() + Arg1()); + TOC + + int x,y,z,t,k; + for(k=0; k