From 9557d6f5054bf6a54ef22ad172201f157ab123f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Markus=20D=C3=B6ring?= Date: Wed, 26 Mar 2014 13:39:04 +0100 Subject: [PATCH] Squashed commit of the following: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit commit 57cada32d220e3743a61b80e89770bf66c36981b Author: Markus Döring Date: Wed Mar 26 13:37:33 2014 +0100 correct path to python unittests commit f13569068647a5322876b67e1de850546328b55e Author: Markus Döring Date: Wed Mar 26 13:29:19 2014 +0100 updated travis instructions to include python tests commit 66a161e4b22acd09d0f10e40646fd6e72015c22a Author: Markus Döring Date: Wed Mar 26 13:21:38 2014 +0100 disabled old style lazyflow tests commit 382bc2c68fd20b63cc2195e25d947ae424386f60 Author: Markus Döring Date: Wed Mar 26 13:03:53 2014 +0100 reintroduced const where applicable commit 039444693f23942833797478e427d29b86e2f828 Author: Markus Döring Date: Wed Mar 26 12:53:06 2014 +0100 adhere to coding style commit b57bee19db786e4569cee4f02b0654531c660fef Author: Markus Döring Date: Tue Mar 25 19:15:15 2014 +0100 travis: update to newer h5py version commit 6fbcf3126ffd08b831a31608064f8e1fdfb8e06b Author: Markus Döring Date: Tue Mar 25 18:54:05 2014 +0100 revert unnecessary file changes commit fcdf0010b1bc1313d2a9d4d3ce2f438775645cf1 Author: Markus Döring Date: Tue Mar 25 18:05:56 2014 +0100 introduced ROI interface to adapters, better roi interface commit 83207e554f1aba3fdca051a579088b498fa9e302 Author: Markus Döring Date: Sat Mar 22 15:51:51 2014 +0100 uint32 support commit 903995e098d8f00c5532907ae22ccffffd48bba0 Author: Markus Döring Date: Sat Mar 22 15:21:13 2014 +0100 wrote labeling operator commit fcbb833cd72471ed7ff65afcd2fe8f25a7689ea6 Author: Markus Döring Date: Fri Mar 21 17:53:01 2014 +0100 all adapters exported commit cfc07813f994b946f426432eb2620908de07e797 Author: Markus Döring Date: Thu Mar 20 20:45:40 2014 +0100 removed unnecessary const attribute in sink commit 89da0e38b9a53723dd407582e1a52a58ca29f92f Author: Markus Döring Date: Thu Mar 20 19:03:10 2014 +0100 provide more types for standard source/sink commit 51d7093594ad6ce833ab0e741aef7fb60663f1c2 Author: Markus Döring Date: Wed Jan 22 13:48:07 2014 +0100 more doc on pyABC, tuple instead of ndarray commit c1c1a24c2fb99e19992d0dc92a8d914f9213e231 Author: Markus Döring Date: Mon Jan 20 19:22:42 2014 +0100 ubuntu does not know how to convert ndarray to tinyvector commit f808dc4fc2e6bef8ac4865baa30898d0ee39e829 Author: Markus Döring Date: Mon Jan 20 12:20:33 2014 +0100 removed output to stdout commit f4ca707be73716efde9db472cf2e08176ae4a77c Author: Markus Döring Date: Mon Jan 6 13:51:37 2014 +0100 python interface for cc complete commit 51ed35628460a5dbeb4c2ce9c483c3268bbb370f Author: Markus Döring Date: Sun Jan 5 17:00:04 2014 +0100 source and sink working in python commit ab3263bd2fb9c956cf99ec3156e694582f9e678a Author: Markus Döring Date: Sat Jan 4 18:48:04 2014 +0100 getting there commit 3ead3ee78f3992461414688ba15aa5ab8be4ecb7 Author: Markus Döring Date: Fri Jan 3 19:03:27 2014 +0100 try again, fail better commit b12f5141e639e8076d7be0257ee4a75f7f7e1997 Author: Markus Döring Date: Thu Jan 2 16:34:14 2014 +0100 started developmnet --- .travis.yml | 2 + .../requirements/development-stage2.txt | 3 +- blockedarray/.gitignore | 1 + blockedarray/CMakeLists.txt | 3 +- blockedarray/__init__.py | 2 +- blockedarray/adapters.py | 148 ++++++++++++++++ blockedarray/adapters_py.cxx | 140 +++++++++++++++ blockedarray/adapters_py.h | 7 + blockedarray/blockwisecc_py.cxx | 160 ++++++++++-------- .../lazyflow_cpp/testOpArrayCacheCpp.py | 38 +---- blockedarray/module_py.cxx | 4 +- blockedarray/opBlockedConnectedComponents.py | 97 +++++++++++ .../testOpBlockedConnectedComponents.py | 91 ++++++++++ blockedarray/test_adapters.py | 59 +++++++ blockedarray/test_blockedarray.py | 1 + blockedarray/test_connectedcomponents.py | 92 ++++++++++ include/bw/connectedcomponents.h | 70 +++++--- include/bw/roi.h | 16 ++ include/bw/sink.h | 10 +- include/bw/source.h | 8 +- 20 files changed, 810 insertions(+), 142 deletions(-) create mode 100644 blockedarray/adapters.py create mode 100644 blockedarray/adapters_py.cxx create mode 100644 blockedarray/adapters_py.h create mode 100644 blockedarray/opBlockedConnectedComponents.py create mode 100644 blockedarray/testOpBlockedConnectedComponents.py create mode 100644 blockedarray/test_adapters.py create mode 100644 blockedarray/test_connectedcomponents.py diff --git a/.travis.yml b/.travis.yml index 9c00c87..66cd33f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -14,3 +14,5 @@ before_script: - "export PYTHONPATH=$VIRTUAL_ENV/lib/python2.7/site-packages:$PYTHONPATH" script: - make test + - cd .. + - nosetests blockedarray diff --git a/.travis_scripts/requirements/development-stage2.txt b/.travis_scripts/requirements/development-stage2.txt index 80035be..ab103d7 100644 --- a/.travis_scripts/requirements/development-stage2.txt +++ b/.travis_scripts/requirements/development-stage2.txt @@ -1 +1,2 @@ -h5py==2.1.3 +h5py==2.2.1 +nose diff --git a/blockedarray/.gitignore b/blockedarray/.gitignore index 76e8f19..4d83ce1 100644 --- a/blockedarray/.gitignore +++ b/blockedarray/.gitignore @@ -1,2 +1,3 @@ _blockedarray.so *.pyc +*.h5 diff --git a/blockedarray/CMakeLists.txt b/blockedarray/CMakeLists.txt index 2f746c5..4b9ecb5 100644 --- a/blockedarray/CMakeLists.txt +++ b/blockedarray/CMakeLists.txt @@ -24,6 +24,7 @@ add_library(_blockedarray SHARED module_py.cxx blockedarray_py.cxx blockwisecc_py.cxx + adapters_py.cxx ${EXTRA_SRCS} ) @@ -34,9 +35,9 @@ else() endif() target_link_libraries(_blockedarray + ${Boost_PYTHON_LIBRARIES} snappy ${PYTHON_LIBRARY} - ${Boost_PYTHON_LIBRARIES} ${HDF5_LIBRARY} ${HDF5_HL_LIBRARY} ${VIGRA_IMPEX_LIBRARY} diff --git a/blockedarray/__init__.py b/blockedarray/__init__.py index a9c83d9..a706a1c 100644 --- a/blockedarray/__init__.py +++ b/blockedarray/__init__.py @@ -1,3 +1,3 @@ from _blockedarray import * - +from opBlockedConnectedComponents import OpBlockedConnectedComponents diff --git a/blockedarray/adapters.py b/blockedarray/adapters.py new file mode 100644 index 0000000..e97a0c7 --- /dev/null +++ b/blockedarray/adapters.py @@ -0,0 +1,148 @@ + +import vigra + +import numpy as np +from _blockedarray import * + + +# This file contains examples for using the Source and Sink classes from +# blockedarray. These should be viewed as guidelines on how to use the exported +# interfaces +# dim[2|3].Source[U|S][8|16|32|64] +# dim[2|3].Sink[U|S][8|16|32|64] +# The ABC classes are for reference on what methods are exposed, the Example +# classes show how the base classes can be used. There is a complete workflow +# for blockwise connected components in the test file +# test_connectedcomponents.py + + +# examples deal with 8bit image input and 32bit image output +_Source = dim3.PySourceU8 +_Sink = dim3.PySinkU32 + + +## Source Interface +# +# This class provides the python interface to the C++ class `BW::Source`. If +# you inherit from this class, be sure to implement *all* abstract methods. +class SourceABC(_Source): + + ## Constructor + # + # Be sure to call super().__init__ if you override the constructor! + def __init__(self): + super(SourceABC, self).__init__() + + ## Set a custom ROI + # selects only the region of interest given from the + # underlying data source. When readBlock() is used, the coordinates + # are relative to roi[0] + # + # @param roi a tuple containing 2 lists of length 3 (incl. start, excl. stop) + def pySetRoi(self, roi): + raise NotImplementedError + #pass + + ## Volume shape getter + # + # @return must be a tuple of (python) integers + def pyShape(self): + raise NotImplementedError + #return (10, 10, 10) + + ## Block getter + # + # read a block of data into a 3d numpy array + # + # @param roi a tuple containing 2 lists of length 3 (incl. start, excl. stop) + # @param output a 3d numpy.ndarray with shape roi[1]-roi[0] and dtype uint8 + def pyReadBlock(self, roi, output): + raise NotImplementedError + #return True + + +## Sink Interface +# +# This class provides the python interface to the C++ class `BW::Sink`. If +# you inherit from this class, be sure to implement *all* abstract methods. +class SinkABC(_Sink): + + ## Constructor + # + # Be sure to call super().__init__ if you override the constructor! + def __init__(self): + super(SinkABC, self).__init__() + + ## Write a block of data + # + # @param roi a tuple containing 2 lists of length 3 (incl. start, excl. stop) + # @param block a 3d numpy.ndarray with shape roi[1]-roi[0] and dtype int32 + def pyWriteBlock(self, roi, block): + raise NotImplementedError + #return True + + +class DummySource(SourceABC): + def __init__(self): + super(DummySource, self).__init__() + + def pySetRoi(self, roi): + pass + + def pyShape(self): + return (100, 100, 10) + + def pyReadBlock(self, roi, output): + output[...] = 0 + return True + + +class ExampleSource(SourceABC): + def __init__(self, vol): + # the call to super() is super-important! + super(ExampleSource, self).__init__() + self._vol = vol + self._p = np.zeros((len(vol.shape),), dtype=np.long) + self._q = np.asarray(vol.shape, dtype=np.long) + + def pySetRoi(self, roi): + self._p = np.asarray(roi.p, dtype=np.long) + self._q = np.asarray(roi.q, dtype=np.long) + + def pyShape(self): + return self._vol.shape + + def pyReadBlock(self, roi, output): + roiP = np.asarray(roi.p) + roiQ = np.asarray(roi.q) + p = self._p + roiP + q = p + roiQ - roiP + if np.any(q > self._q): + raise IndexError("Requested roi is too large for selected " + "roi (previous call to setRoi)") + s = _roi2slice(p, q) + output[...] = self._vol[s] + return True + + +class ExampleSink(SinkABC): + def __init__(self): + super(ExampleSink, self).__init__() + self.vol = None + + def pyWriteBlock(self, roi, block): + if self.vol is None: + shape = self.shape + shape = _v2tup(shape) + self.vol = np.zeros(shape, dtype=np.uint8) + s = _roi2slice(roi.p, roi.q) + self.vol[s] = block + + +def _v2tup(v, d=3): + return tuple([v[i] for i in range(d)]) + + +def _roi2slice(p, q): + s = [slice(p[i], q[i]) for i in range(len(p))] + return tuple(s) diff --git a/blockedarray/adapters_py.cxx b/blockedarray/adapters_py.cxx new file mode 100644 index 0000000..d2761f6 --- /dev/null +++ b/blockedarray/adapters_py.cxx @@ -0,0 +1,140 @@ +#define PY_ARRAY_UNIQUE_SYMBOL blockedarray_PyArray_API +#define NO_IMPORT_ARRAY + +#include "Python.h" + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include "adapters_py.h" + +#include + +using namespace BW; + + +template +struct PySourceABC : Source, boost::python::wrapper > +{ +public: + + typedef typename Source::V TinyVec; + + PySourceABC() {} + virtual ~PySourceABC() {} + + void setRoi(Roi roi) { + this->pySetRoi(roi); + } + + TinyVec shape() const { + return this->pyShape(); + } + + bool readBlock(Roi roi, vigra::MultiArrayView& block) const { + bool ret; + + //temporary NumpyArray, because MultiArrayView is not convertible to python + vigra::NumpyArray tempArray(roi.q-roi.p); + + if (!this->pyReadBlock(roi, tempArray)) + return false; + + block.copy(tempArray); + + return true; + } + + void pySetRoi(Roi roi) { + this->get_override("pySetRoi")(roi); + } + + TinyVec pyShape() const { + return this->get_override("pyShape")(); + }; + + bool pyReadBlock(Roi roi, vigra::NumpyArray& block) const { + return this->get_override("pyReadBlock")(roi, block); + }; +}; + +template +struct PySinkABC : Sink, boost::python::wrapper > { + public: + typedef typename Roi::V V; + + PySinkABC() {}; + virtual ~PySinkABC() {}; + + bool writeBlock(Roi roi, const vigra::MultiArrayView& block) { + vigra::NumpyArray tempArray(block); + return this->pyWriteBlock(roi, tempArray); + } + + bool pyWriteBlock(Roi roi, const vigra::NumpyArray& block) { + return this->get_override("pyWriteBlock")(roi, block); + }; + + + V getShape() const { + return this->shape_; + } + + void setShape(V shape) { + this->shape_ = shape; + } + + V getBlockShape() const { + return this->blockShape_; + } + + void setBlockShape(V shape) { + this->blockShape_ = shape; + } +}; + + + +template +void exportSpecificSourceAdapter(std::string suffix) { + using namespace boost::python; + + class_, boost::noncopyable>(("PySource" + suffix).c_str()) + .def("pySetRoi", pure_virtual(&PySourceABC::pySetRoi)) + .def("pyShape", pure_virtual(&PySourceABC::pyShape)) + .def("pyReadBlock", pure_virtual(&PySourceABC::pyReadBlock)) + ; + + class_, boost::noncopyable>(("PySink" + suffix).c_str()) + .def("pyWriteBlock", pure_virtual(&PySinkABC::pyWriteBlock)) + .add_property("shape", &PySinkABC::getShape, &PySinkABC::setShape) + .add_property("blockShape", &PySinkABC::getBlockShape, &PySinkABC::setBlockShape) + ; +} + +template +void exportSourceAdaptersForDim() { + + exportSpecificSourceAdapter("U8"); + //exportSpecificSourceAdapter("U16"); + exportSpecificSourceAdapter("U32"); + + exportSpecificSourceAdapter("S8"); + //exportSpecificSourceAdapter("S16"); + exportSpecificSourceAdapter("S32"); + + exportSpecificSourceAdapter("F"); + exportSpecificSourceAdapter("D"); +} + +template void exportSourceAdaptersForDim<2>(); +template void exportSourceAdaptersForDim<3>(); + diff --git a/blockedarray/adapters_py.h b/blockedarray/adapters_py.h new file mode 100644 index 0000000..38c30b4 --- /dev/null +++ b/blockedarray/adapters_py.h @@ -0,0 +1,7 @@ +#ifndef ADAPTERS_PY_H +#define ADAPTERS_PY_H + +template +void exportSourceAdaptersForDim(); + +#endif /* ADAPTERS_PY_H */ diff --git a/blockedarray/blockwisecc_py.cxx b/blockedarray/blockwisecc_py.cxx index fbd7a98..8b9e518 100644 --- a/blockedarray/blockwisecc_py.cxx +++ b/blockedarray/blockwisecc_py.cxx @@ -43,98 +43,108 @@ #include #include "blockwisecc_py.h" +#include "adapters_py.h" #include using namespace BW; -template -struct PyConnectedComponents { - typedef ConnectedComponents BCC; -}; - -template -struct ExportV { - static void export_(); -}; - -template -struct ExportV<2, V> { - static void export_() { - using namespace boost::python; - class_("V", init()); - } -}; - - -template -struct ExportV<3, V> { - static void export_() { - using namespace boost::python; - class_("V", init()); - } -}; - -template -struct ExportV<4, V> { - static void export_() { - using namespace boost::python; - class_("V", init()); - } -}; - -template -void blockwiseCC() { +template +void exportSpecificCC(std::string suffix) { using namespace boost::python; - typedef ConnectedComponents BCC; - typedef Thresholding BWT; - typedef ChannelSelector BWCS; - typedef PyConnectedComponents PyBCC; - typedef SourceHDF5 HDF5BP_T; + using namespace BW; + typedef ConnectedComponents BCC; + + class_(("ConnectedComponents"+suffix).c_str(), + init*, typename BCC::V>()) + .def("writeResult", &BCC::writeResult, + (arg("hdf5file"), arg("hdf5group"), arg("compression")=1)) + .def("writeToSink", &BCC::writeToSink, + (arg("sink"))) + ; +} - std::stringstream n; n << N; +/* CC conversion */ +template +void exportCCForDim() { + exportSpecificCC("U8"); + exportSpecificCC("U32"); +} - std::stringstream fullModname; fullModname << "_blockedarray.dim" << N; - std::stringstream modName; modName << "dim" << N; +/* ROI conversion */ +template +void exportRoiForDim() { + using namespace boost::python; + using namespace BW; + typedef typename Roi::V V; + + class_< Roi >("Roi", init()) + //.def_readwrite("p", &Roi::p) + //.def_readwrite("q", &Roi::p) + //.def("getP", &Roi::getP) + //.def("getQ", &Roi::getQ) + .add_property("p", &Roi::getP, &Roi::setP) + .add_property("q", &Roi::getQ, &Roi::setQ) + ; +} - //see: http://isolation-nation.blogspot.de/2008/09/packages-in-python-extension-modules.html - object module(handle<>(borrowed(PyImport_AddModule(fullModname.str().c_str())))); - scope().attr(modName.str().c_str()) = module; - scope s = module; - ExportV::export_(); +template +void exportSpecificSource(std::string suffix) { + + using namespace boost::python; + const char *source = ("Source"+suffix).c_str(); + const char *sink = ("Sink"+suffix).c_str(); + class_ >(source, no_init); + class_ >(sink, no_init); +} - class_ >("Source", no_init); - class_ >("Sink", no_init); +template +void exportSourceForDim() { + exportSpecificSource("U8"); + //exportSpecificSource("U16"); + exportSpecificSource("U32"); - class_, bases > >("SourceHDF5", - init()) - ; - class_, bases > >("SinkHDF5", - init()) - ; + exportSpecificSource("S8"); + //exportSpecificSource("S16"); + exportSpecificSource("S32"); - class_("Thresholding", - init*, typename BWT::V>()) - .def("run", vigra::registerConverters(&BWT::run), - (arg("threshold"), arg("ifLower"), arg("ifHigher"), arg("sink"))) - ; + exportSpecificSource("F"); + exportSpecificSource("D"); +} - class_("ChannelSelector", - init*, typename BWCS::V>()) - .def("run", vigra::registerConverters(&BWCS::run), - (arg("dimension"), arg("channel"), arg("sink"))) - ; - class_("ConnectedComponents", - init*, typename BCC::V>()) - .def("writeResult", &BCC::writeResult, - (arg("hdf5file"), arg("hdf5group"), arg("compression")=1)) - ; +template +void exportAllForDim() { + + using namespace boost::python; + + // set the correct module + std::stringstream n; n << N; + + std::stringstream fullModname; fullModname << "_blockedarray.dim" << N; + std::stringstream modName; modName << "dim" << N; + + //see: http://isolation-nation.blogspot.de/2008/09/packages-in-python-extension-modules.html + object module(handle<>(borrowed(PyImport_AddModule(fullModname.str().c_str())))); + scope().attr(modName.str().c_str()) = module; + scope s = module; + + + // export source and sink + exportSourceForDim(); + + exportCCForDim(); + + exportRoiForDim(); + + exportSourceAdaptersForDim(); + } + void export_blockwiseCC() { - blockwiseCC<2, float>(); - blockwiseCC<3, float>(); + exportAllForDim<2>(); + exportAllForDim<3>(); } diff --git a/blockedarray/lazyflow_cpp/testOpArrayCacheCpp.py b/blockedarray/lazyflow_cpp/testOpArrayCacheCpp.py index 2ab1766..8bcce50 100644 --- a/blockedarray/lazyflow_cpp/testOpArrayCacheCpp.py +++ b/blockedarray/lazyflow_cpp/testOpArrayCacheCpp.py @@ -1,10 +1,13 @@ import threading +import unittest import numpy import vigra from lazyflow.graph import Graph from lazyflow.roi import sliceToRoi, roiToSlice from lazyflow.operators import OpArrayPiper -from lazyflow.operators.opArrayCacheCpp import OpArrayCacheCpp + +# does not exist anymore!! +#from lazyflow.operators.opArrayCacheCpp import OpArrayCacheCpp class KeyMaker(): def __getitem__(self, *args): @@ -26,7 +29,8 @@ def execute(self, slot, subindex, roi, result): super(OpArrayPiperWithAccessCount, self).execute(slot, subindex, roi, result) -class TestOpArrayCache(object): +@unittest.skip("OpArrayCacheCpp is outdated") +class TestOpArrayCache(unittest.TestCase): def setUp(self): print "RRRRRRRRRRRRRRRRRRRRRRRRRR" @@ -285,31 +289,5 @@ def handleDirty(slot, roi): sys.argv.append("--nologcapture") # Don't set the logging level to DEBUG. Leave it alone. ret = nose.run(defaultTest=__file__) - - - - - - - - - - - - - - - - - - - - - - - - - - - - if not ret: sys.exit(1) + if not ret: + sys.exit(1) diff --git a/blockedarray/module_py.cxx b/blockedarray/module_py.cxx index 42ce85d..022570e 100644 --- a/blockedarray/module_py.cxx +++ b/blockedarray/module_py.cxx @@ -34,12 +34,14 @@ #include #include +#include #include "blockedarray_py.h" #include "blockwisecc_py.h" + BOOST_PYTHON_MODULE_INIT(_blockedarray) { - _import_array(); + //_import_array(); vigra::import_vigranumpy(); boost::python::numeric::array::set_module_and_type("numpy", "ndarray"); diff --git a/blockedarray/opBlockedConnectedComponents.py b/blockedarray/opBlockedConnectedComponents.py new file mode 100644 index 0000000..5f1469f --- /dev/null +++ b/blockedarray/opBlockedConnectedComponents.py @@ -0,0 +1,97 @@ + +import numpy as np +import vigra + +from lazyflow.operator import Operator, InputSlot, OutputSlot +from lazyflow.operators import OpCompressedCache +from lazyflow.operators.opLabelVolume import OpNonLazyCC +from lazyflow.rtype import SubRegion + +from _blockedarray import dim2, dim3 + + +## compute connected components blockwise +# +# Input must be 'xyzct' so this operator can be wrapped with +# lazyflow.OpLabelVolume +class OpBlockedConnectedComponents(OpNonLazyCC): + name = "OpBlockedConnectedComponents" + supportedDtypes = [np.uint8, np.uint32] + + def __init__(self, *args, **kwargs): + super(OpBlockedConnectedComponents, self).__init__(*args, **kwargs) + self._sourceMap3 = dict(zip(self.supportedDtypes, + [dim3.PySourceU8, dim3.PySourceU32])) + self._sinkMap3 = {np.uint32: dim3.PySinkU32} + self._ccMap3 = dict(zip(self.supportedDtypes, + [dim3.ConnectedComponentsU8, + dim3.ConnectedComponentsU32])) + + def _updateSlice(self, c, t, bg): + bg = self.Background[0, 0, 0, c, t].wait() + assert bg == 0, "Non-zero background not supported" + source = self._getSource(c, t) + sink = self._getSink(c, t) + #TODO enable 2d + blockShape = _v2tup(self._cache.BlockShape.value, dim=3) + CC = self._ccMap3[self.Input.meta.dtype] + cc = CC(source, blockShape) + + cc.writeToSink(sink) + + ## factory for dynamic source creation + def _getSource(self, c, t): + SourceClass = self._sourceMap3[self.Input.meta.dtype] + shape = tuple([int(s) for s in self.Input.meta.shape[:3]]) + + class TempSource(SourceClass): + def __init__(self, *args, **kwargs): + self._slot = kwargs['slot'] + del kwargs['slot'] + super(TempSource, self).__init__(*args, **kwargs) + def pySetRoi(self, roi): + assert False, "Not supported" + def pyShape(self): + return shape + def pyReadBlock(self, roi, block): + start = roi.p + (c, t) + stop = roi.q + (c+1, t+1) + subr = SubRegion(self._slot, start=start, stop=stop) + block[:] = self._slot.get(subr).wait()[..., 0, 0] + return True + + return TempSource(slot=self.Input) + + ## factory for dynamic sink creation + def _getSink(self, c, t): + SinkClass = self._sinkMap3[self.Output.meta.dtype] + + class TempSink(SinkClass): + def __init__(self, *args, **kwargs): + self._cache = kwargs['cache'] + del kwargs['cache'] + super(TempSink, self).__init__(*args, **kwargs) + self.shape = _v2tup(self._cache.Input.meta.shape) + self.blockShape = _v2tup(self._cache.BlockShape.value) + def pyWriteBlock(self, roi, block): + block = vigra.taggedView(block, axistags='xyz') + block = block.withAxes(*'xyzct') + start = roi.p + (c, t) + stop = roi.q + (c+1, t+1) + subr = SubRegion(self._cache.Input, start=start, stop=stop) + self._cache.setInSlot(self._cache.Input, (), subr, block) + return True + + return TempSink(cache=self._cache) + + def setupOutputs(self): + super(OpBlockedConnectedComponents, self).setupOutputs() + assert len(self.Input.meta.shape) == 5, "Input must be 5d" + if self.Input.meta.axistags: + # if axistags are given, they must be xyzct + s = "".join(self.Input.meta.getAxisKeys()) + assert s == "xyzct", "Input must be in xyzct order, if any" + + +def _v2tup(v, dim=3): + return tuple([int(v[i]) for i in range(dim)]) diff --git a/blockedarray/testOpBlockedConnectedComponents.py b/blockedarray/testOpBlockedConnectedComponents.py new file mode 100644 index 0000000..fc9550f --- /dev/null +++ b/blockedarray/testOpBlockedConnectedComponents.py @@ -0,0 +1,91 @@ +import numpy as np +import vigra + +import unittest + +from lazyflow.graph import Graph + +from numpy.testing import assert_array_equal, assert_array_almost_equal + +#import blockedarray.OpBlockedConnectedComponents +from blockedarray import OpBlockedConnectedComponents + + + +class TestOpBlockedConnectedComponents(unittest.TestCase): + + def setUp(self): + bg = np.zeros((1,1,1,1,1), dtype=np.uint8) + bg = vigra.taggedView(bg, axistags='xyzct') + self.bg = bg + + def testSimpleUsage(self): + vol = np.random.randint(255, size=(1000, 100, 10)) + vol = vol.astype(np.uint8) + vol = vigra.taggedView(vol, axistags='xyz') + vol = vol.withAxes(*'xyzct') + + op = OpBlockedConnectedComponents(graph=Graph()) + op.Input.setValue(vol) + op.Background.setValue(self.bg) + + out = op.Output[...].wait() + + assert_array_equal(vol.shape, out.shape) + + def testCorrectLabeling(self): + vol = np.zeros((1000, 100, 10)) + vol = vol.astype(np.uint8) + vol = vigra.taggedView(vol, axistags='xyz') + + vol[20:40, 10:30, 2:4] = 1 + vol = vol.withAxes(*'xyzct') + + op = OpBlockedConnectedComponents(graph=Graph()) + op.Input.setValue(vol) + op.Background.setValue(self.bg) + + out = op.Output[...].wait() + tags = op.Output.meta.getTaggedShape() + out = vigra.taggedView(out, axistags=op.Output.meta.axistags) + + assert_array_equal(vol, out) + + def testCorrectLabelingInt(self): + vol = np.zeros((1000, 100, 10)) + vol = vol.astype(np.uint32) + vol = vigra.taggedView(vol, axistags='xyz') + + vol[20:40, 10:30, 2:4] = 1 + vol = vol.withAxes(*'xyzct') + + op = OpBlockedConnectedComponents(graph=Graph()) + op.Input.setValue(vol) + op.Background.setValue(self.bg) + + out = op.Output[...].wait() + tags = op.Output.meta.getTaggedShape() + out = vigra.taggedView(out, axistags=op.Output.meta.axistags) + + assert_array_equal(vol, out) + + +class TestSimpleThings(unittest.TestCase): + def testRoi(self): + from blockedarray import dim3 + roi = dim3.Roi((0,0,0), (2,3,4)) + p = roi.p + assert isinstance(p, tuple) + assert len(p) ==3 + q = roi.q + assert isinstance(q, tuple) + assert len(q) ==3 + + roi.p = q + roi.q = p + + tempQ = roi.p + assert isinstance(tempQ, tuple) + assert tempQ == q + + diff --git a/blockedarray/test_adapters.py b/blockedarray/test_adapters.py new file mode 100644 index 0000000..cc494e2 --- /dev/null +++ b/blockedarray/test_adapters.py @@ -0,0 +1,59 @@ +import numpy as np +import unittest + +from adapters import ExampleSource, ExampleSink +from _blockedarray import dim3 + + +class TestSource(unittest.TestCase): + + def setUp(self): + self.vol = np.random.randint(0, 2**16-1, (100, 100, 10)).astype(np.uint32) + + def testExampleSource(self): + vol = self.vol + s = ExampleSource(vol) + + roi = dim3.Roi((0,)*len(vol.shape), + tuple(vol.shape)) + newVol = np.zeros(vol.shape, dtype=np.uint32) + assert s.pyReadBlock(roi, newVol) + assert np.all(vol == newVol) + + def testRoi(self): + vol = self.vol + s = ExampleSource(vol) + + reqRoi = [(50, 50, 2), (70, 70, 4)] + reqRoi = dim3.Roi(reqRoi[0], reqRoi[1]) + s.pySetRoi(reqRoi) + roi = [(0, 0, 0), (20, 20, 2)] + roi = dim3.Roi(roi[0], roi[1]) + newVol = np.zeros((20, 20, 2), dtype=np.uint32) + assert s.pyReadBlock(roi, newVol) + assert np.all(vol[50:70, 50:70, 2:4] == newVol) + + def testRoi2(self): + vol = self.vol + s = ExampleSource(vol) + roi = [(0, 0, 2), (20, 20, 4)] + roi = dim3.Roi(roi[0], roi[1]) + newVol = np.zeros((20, 20, 2), dtype=np.uint32) + assert s.pyReadBlock(roi, newVol) + + +class TestSink(unittest.TestCase): + + def setUp(self): + pass + + def testExampleSink(self): + s = ExampleSink() + s.shape = (100, 100, 10) + s.blockShape = (10, 10, 10) + roi = [(15, 20, 2), (30, 30, 6)] + roi = dim3.Roi(roi[0], roi[1]) + s.pyWriteBlock(roi, + np.ones((15, 10, 4), dtype=np.uint8)) + + assert np.all(s.vol[15:30, 20:30, 2:6] == 1) diff --git a/blockedarray/test_blockedarray.py b/blockedarray/test_blockedarray.py index aaf4a8f..42ee75e 100644 --- a/blockedarray/test_blockedarray.py +++ b/blockedarray/test_blockedarray.py @@ -74,6 +74,7 @@ def test1(): ba.setCompressionEnabled(True) ba.setCompressionEnabled(False) + if __name__ == "__main__": test1() print "success" diff --git a/blockedarray/test_connectedcomponents.py b/blockedarray/test_connectedcomponents.py new file mode 100644 index 0000000..bc46da1 --- /dev/null +++ b/blockedarray/test_connectedcomponents.py @@ -0,0 +1,92 @@ +import unittest +import tempfile +import os + +import numpy as np +from numpy.testing import assert_almost_equal +import vigra + +from adapters import DummySource, ExampleSource, ExampleSink +from _blockedarray import dim3 + + +CC = dim3.ConnectedComponentsU8 + +class TestConnectedComponents(unittest.TestCase): + + def setUp(self): + shape = (100, 110, 120) + self.blockShape = (10, 10, 10) + vol = np.zeros(shape, dtype=np.uint8) + vol[7:12, 17:22, 27:32] = 42 + vol[:, :, 72:75] = 111 + vol[:, 72:75, :] = 111 + self.vol = vol + + def testBlockShape(self): + s = DummySource() + + v = (10, 10, 10) + cc = CC(s, v) + + ''' + # these are not expected to work + v = np.asarray((10, 10, 10)) + cc = CC(s, v) + + v = np.asarray((10, 10, 10), dtype=np.long) + cc = CC(s, v) + + v = np.asarray((10, 10, 10), dtype=np.int) + cc = CC(s, v) + ''' + + def testCC(self): + tmp = tempfile.NamedTemporaryFile(suffix='.h5', delete=False) + try: + h5f = (tmp.name, "/data") + s = ExampleSource(self.vol) + cc = CC(s, self.blockShape) + + cc.writeResult(h5f[0], h5f[1]) + labels = vigra.readHDF5(h5f[0], h5f[1]) + except Exception: + raise + finally: + pass + os.remove(h5f[0]) + + # vigra and hdf5 have incompatible axis order, see + # http://ukoethe.github.io/vigra/doc/vigra/classvigra_1_1HDF5File.html#a638be2d51070009b02fed0f234f4b0bf + labels = np.transpose(labels, (2, 1, 0)) + + self._check(labels) + + def testCCsink(self): + source = ExampleSource(self.vol) + sink = ExampleSink() + cc = CC(source, self.blockShape) + cc.writeToSink(sink) + + self._check(sink.vol) + + def _check(self, labels): + comp1 = np.concatenate((labels[:, :, 72:75].flatten(), + labels[:, 72:75, :].flatten())) + comp2 = labels[7:12, 17:22, 27:32].flatten() + + # assert that each component has exactly one positive label + assert_almost_equal(comp1.var(), 0) + assert np.all(comp1 > 0) + assert_almost_equal(comp2.var(), 0) + assert np.all(comp2 > 0) + + # assert that each component is labeled with a unique label + assert comp1[0] != comp2[0] + + # assert that no other pixel is labeled + labels[:, :, 72:75] = 0 + labels[:, 72:75, :] = 0 + labels[7:12, 17:22, 27:32] = 0 + assert_almost_equal(labels.sum(), 0) + diff --git a/include/bw/connectedcomponents.h b/include/bw/connectedcomponents.h index bc8c704..53129a8 100644 --- a/include/bw/connectedcomponents.h +++ b/include/bw/connectedcomponents.h @@ -51,12 +51,14 @@ #include namespace BW { + +typedef unsigned int LabelType; using vigra::detail::UnionFindArray; template void blockMerge( - UnionFindArray& ufd, + UnionFindArray& ufd, Roi roi1, T offset1, const vigra::MultiArrayView& block1, @@ -76,8 +78,8 @@ void blockMerge( typename ArrayView::iterator it1 = ov1.begin(); typename ArrayView::iterator it2 = ov2.begin(); for(; it1 != ov1.end(); ++it1, ++it2) { - int x = *it1+offset1; - int y = *it2+offset2; + LabelType x = *it1+offset1; + LabelType y = *it2+offset2; //if(*it1 == 0) ufd.makeUnion(0,y); //if(*it2 == 0) ufd.makeUnion(x,0); ufd.makeUnion(x, y); @@ -110,16 +112,14 @@ struct ConnectedComponentsComputer<3,T,LabelType> { /** * Compute connected components block-wise (less limited to RAM) */ -template +template class ConnectedComponents { public: typedef typename Roi::V V; - typedef int LabelType; - typedef CompressedArray Compressed; - ConnectedComponents(Source* blockProvider, V blockShape) + ConnectedComponents(Source* blockProvider, V blockShape) : blockProvider_(blockProvider) , blockShape_(blockShape) { @@ -139,13 +139,13 @@ class ConnectedComponents { vector offsets(blockRois_.size()+1); ccBlocks_.resize(blockRois_.size()); - std::cout << "* run connected components on each of " << blockRois_.size() << " blocks separately" << std::endl; +// std::cout << "* run connected components on each of " << blockRois_.size() << " blocks separately" << std::endl; size_t sizeBytes = 0; size_t sizeBytesUncompressed = 0; for(size_t i=0; i& block = blockRois_[i].second; //std::cout << " - CC on block " << block << std::endl; @@ -153,12 +153,12 @@ class ConnectedComponents { // vigra::srcMultiArrayRange( labels_.subarray(block.p, block.q) ), // vigra::destMultiArray( ccBlocks[i] ), vigra::NeighborCode3DSix(), 0); - MultiArray inBlock(block.q-block.p); + MultiArray inBlock(block.q-block.p); blockProvider_->readBlock(block, inBlock); MultiArray cc(block.q-block.p); - LabelType maxLabel = ConnectedComponentsComputer::compute(inBlock, cc); + LabelType maxLabel = ConnectedComponentsComputer::compute(inBlock, cc); //LabelType maxLabel = vigra::labelImageWithBackground( // inBlock, // cc, @@ -175,18 +175,18 @@ class ConnectedComponents { offsets[i+1] = offsets[i] + maxLabel + 1; } - std::cout << std::endl; - std::cout << " " << sizeBytes/(1024.0*1024.0) << " MB (vs. " << sizeBytesUncompressed/(1024.0*1024.0) << "MB uncompressed)" << std::endl; +// std::cout << std::endl; +// std::cout << " " << sizeBytes/(1024.0*1024.0) << " MB (vs. " << sizeBytesUncompressed/(1024.0*1024.0) << "MB uncompressed)" << std::endl; LabelType maxOffset = offsets[offsets.size()-1]; - std::cout << " initialize union find datastructure with maxOffset = " << maxOffset << std::endl; +// std::cout << " initialize union find datastructure with maxOffset = " << maxOffset << std::endl; UnionFindArray ufd(maxOffset); // // merge adjacent blocks // - std::cout << "* merge adjacent blocks" << std::endl; +// std::cout << "* merge adjacent blocks" << std::endl; //FIXME: use adjacency relation of blocks std::vector< std::pair > adjBlocks; @@ -202,7 +202,7 @@ class ConnectedComponents { } for(vector >::iterator it=adjBlocks.begin(); it!=adjBlocks.end(); ++it) { - std::cout << " pair " << std::distance(adjBlocks.begin(), it)+1 << "/" << adjBlocks.size() << " \r" << std::flush; +// std::cout << " pair " << std::distance(adjBlocks.begin(), it)+1 << "/" << adjBlocks.size() << " \r" << std::flush; size_t i = it->first; size_t j = it->second; @@ -220,15 +220,15 @@ class ConnectedComponents { block1, offsets[i], cc1, block2, offsets[j], cc2); } - std::cout << std::endl; +// std::cout << std::endl; LabelType maxLabel = ufd.makeContiguous(); - std::cout << "* relabel" << std::endl; +// std::cout << "* relabel" << std::endl; sizeBytes = 0; sizeBytesUncompressed = 0; for(size_t i=0; i& roi = blockRois_[i].second; //MultiArray cc = ccBlocks_[i]; @@ -244,8 +244,8 @@ class ConnectedComponents { sizeBytes += ccBlocks_[i].currentSizeBytes(); sizeBytesUncompressed += ccBlocks_[i].uncompressedSizeBytes(); } - std::cout << std::endl; - std::cout << " " << sizeBytes/(1024.0*1024.0) << " MB (vs. " << sizeBytesUncompressed/(1024.0*1024.0) << "MB uncompressed)" << std::endl; +// std::cout << std::endl; +// std::cout << " " << sizeBytes/(1024.0*1024.0) << " MB (vs. " << sizeBytesUncompressed/(1024.0*1024.0) << "MB uncompressed)" << std::endl; } @@ -254,11 +254,12 @@ class ConnectedComponents { vigra_precondition(compression >= 1 && compression <= 9, "compression must be >= 1 and <= 9"); - std::cout << "* write " << hdf5file << "/" << hdf5group << std::endl; +// std::cout << "* write " << hdf5file << "/" << hdf5group << std::endl; + HDF5File out(hdf5file, HDF5File::Open); out.createDataset(hdf5group, blockProvider_->shape(), 0, blockShape_, compression); for(size_t i=0; i& roi = blockRois_[i].second; //cc = ccBlocks_[i]; @@ -268,12 +269,31 @@ class ConnectedComponents { out.writeBlock(hdf5group, roi.p, cc); } - std::cout << std::endl; +// std::cout << std::endl; out.close(); } + + void writeToSink(Sink* sink) + { +// std::cout << "* writing to sink" << std::endl; + + sink->setShape(blockProvider_->shape()); + sink->setBlockShape(blockShape_); + + for(size_t i=0; i& roi = blockRois_[i].second; + + vigra::MultiArray cc(roi.q-roi.p); + ccBlocks_[i].readArray(cc); + + sink->writeBlock(roi, cc); + } +// std::cout << std::endl; + } private: - Source* blockProvider_; + Source* blockProvider_; V blockShape_; std::vector< std::pair > > blockRois_; diff --git a/include/bw/roi.h b/include/bw/roi.h index 7eb2dca..294076a 100644 --- a/include/bw/roi.h +++ b/include/bw/roi.h @@ -182,6 +182,22 @@ class Roi { } return ret; } + + V getP() const { + return p; + } + + V getQ() const { + return q; + } + + void setP(V p) { + this->p = p; + } + + void setQ(V q) { + this->q = q; + } V p; V q; diff --git a/include/bw/sink.h b/include/bw/sink.h index 1308aa6..40c1e85 100644 --- a/include/bw/sink.h +++ b/include/bw/sink.h @@ -45,20 +45,22 @@ class Sink { typedef typename Roi::V V; Sink() {} - virtual ~Sink() {}; + virtual ~Sink() {} + /* has to be called before any calls to writeBlock */ void setShape(V shape) { shape_ = shape; } + /* has to be called before any calls to writeBlock */ void setBlockShape(V shape) { blockShape_ = shape; } - V shape() const { return shape_; }; - V blockShape() const { return blockShape_; }; + V shape() const { return shape_; } + V blockShape() const { return blockShape_; } - virtual bool writeBlock(Roi roi, const vigra::MultiArrayView& block) { return true; }; + virtual bool writeBlock(Roi roi, const vigra::MultiArrayView& block) { return true; } protected: V shape_; diff --git a/include/bw/source.h b/include/bw/source.h index 0f8a92c..6e13837 100644 --- a/include/bw/source.h +++ b/include/bw/source.h @@ -45,17 +45,17 @@ class Source { typedef typename Roi::V V; Source() {} - virtual ~Source() {}; + virtual ~Source() {} /** * selects only the region of interest given from the * underlying data source. When readBlock() is used, the coordinates * are relative to roi.q */ - virtual void setRoi(Roi roi) {}; + virtual void setRoi(Roi roi) {} - virtual V shape() const { return V(); }; - virtual bool readBlock(Roi roi, vigra::MultiArrayView& block) const { return true; }; + virtual V shape() const { return V(); } + virtual bool readBlock(Roi roi, vigra::MultiArrayView& block) const { return true; } }; } /* namespace BW */