[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/seededregiongrowing3d.hxx
00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn       */
00004 /*                        and Ullrich Koethe                            */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    The VIGRA Website is                                              */
00008 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00034 /*                                                                      */
00035 /************************************************************************/
00036 
00037 #ifndef VIGRA_SEEDEDREGIONGROWING_3D_HXX
00038 #define VIGRA_SEEDEDREGIONGROWING_3D_HXX
00039 
00040 #include <vector>
00041 #include <stack>
00042 #include <queue>
00043 #include "utilities.hxx"
00044 #include "stdimage.hxx"
00045 #include "stdimagefunctions.hxx"
00046 #include "seededregiongrowing.hxx"
00047 #include "multi_pointoperators.hxx"
00048 #include "voxelneighborhood.hxx"
00049 
00050 namespace vigra {
00051 
00052 namespace detail {
00053 
00054 template <class COST, class Diff_type>
00055 class SeedRgVoxel
00056 {
00057 public:
00058     Diff_type location_, nearest_;
00059     COST cost_;
00060     int count_;
00061     int label_;
00062     int dist_;
00063 
00064     SeedRgVoxel()
00065     //: location_(0,0,0), nearest_(0,0,0), cost_(0), count_(0), label_(0)
00066     {
00067         location_ = Diff_type(0,0,0);
00068         nearest_ = Diff_type(0,0,0);
00069         cost_ = 0;
00070         count_ = 0;
00071         label_ = 0;
00072     }
00073 
00074     SeedRgVoxel(Diff_type const & location, Diff_type const & nearest,
00075                 COST const & cost, int const & count, int const & label)
00076     : location_(location), nearest_(nearest),
00077       cost_(cost), count_(count), label_(label)
00078     {
00079         int dx = location_[0] - nearest_[0];
00080         int dy = location_[1] - nearest_[1];
00081         int dz = location_[2] - nearest_[2];
00082         dist_ = dx * dx + dy * dy + dz * dz;
00083     }
00084 
00085     void set(Diff_type const & location, Diff_type const & nearest,
00086              COST const & cost, int const & count, int const & label)
00087     {
00088         location_ = location;
00089         nearest_ = nearest;
00090         cost_ = cost;
00091         count_ = count;
00092         label_ = label;
00093 
00094         int dx = location_[0] - nearest_[0];
00095         int dy = location_[1] - nearest_[1];
00096         int dz = location_[2] - nearest_[2];
00097         dist_ = dx * dx + dy * dy + dz * dz;
00098     }
00099 
00100     struct Compare
00101     {
00102         // must implement > since priority_queue looks for largest element
00103         bool operator()(SeedRgVoxel const & l,
00104                         SeedRgVoxel const & r) const
00105         {
00106             if(r.cost_ == l.cost_)
00107             {
00108                 if(r.dist_ == l.dist_) return r.count_ < l.count_;
00109 
00110                 return r.dist_ < l.dist_;
00111             }
00112 
00113             return r.cost_ < l.cost_;
00114         }
00115         bool operator()(SeedRgVoxel const * l,
00116                         SeedRgVoxel const * r) const
00117         {
00118             if(r->cost_ == l->cost_)
00119             {
00120                 if(r->dist_ == l->dist_) return r->count_ < l->count_;
00121 
00122                 return r->dist_ < l->dist_;
00123             }
00124 
00125             return r->cost_ < l->cost_;
00126         }
00127     };
00128 
00129     struct Allocator
00130     {
00131         ~Allocator()
00132         {
00133             while(!freelist_.empty())
00134             {
00135                 delete freelist_.top();
00136                 freelist_.pop();
00137             }
00138         }
00139 
00140         SeedRgVoxel * create(Diff_type const & location, Diff_type const & nearest,
00141                              COST const & cost, int const & count, int const & label)
00142         {
00143             if(!freelist_.empty())
00144             {
00145                 SeedRgVoxel * res = freelist_.top();
00146                 freelist_.pop();
00147                 res->set(location, nearest, cost, count, label);
00148                 return res;
00149             }
00150 
00151             return new SeedRgVoxel(location, nearest, cost, count, label);
00152         }
00153 
00154         void dismiss(SeedRgVoxel * p)
00155         {
00156             freelist_.push(p);
00157         }
00158 
00159         std::stack<SeedRgVoxel<COST,Diff_type> *> freelist_;
00160     };
00161 };
00162 
00163 } // namespace detail
00164 
00165 /** \addtogroup SeededRegionGrowing
00166     Region segmentation and voronoi tesselation
00167 */
00168 //@{
00169 
00170 /********************************************************/
00171 /*                                                      */
00172 /*                    seededRegionGrowing3D             */
00173 /*                                                      */
00174 /********************************************************/
00175 
00176 /** \brief Three-dimensional Region Segmentation by means of Seeded Region Growing.
00177 
00178     This algorithm implements seeded region growing as described in
00179 
00180     The seed image is a partly segmented multi-dimensional array which contains uniquely
00181     labeled regions (the seeds) and unlabeled voxels (the candidates, label 0).
00182     Seed regions can be as large as you wish and as small as one voxel. If
00183     there are no candidates, the algorithm will simply copy the seed array
00184     into the output array. Otherwise it will aggregate the candidates into
00185     the existing regions so that a cost function is minimized.
00186     Candidates are taken from the neighborhood of the already assigned pixels, 
00187     where the type of neighborhood is determined by parameter <tt>neighborhood</tt>
00188     which can take the values <tt>NeighborCode3DSix()</tt> (the default) 
00189     or <tt>NeighborCode3DTwentySix()</tt>. The algorithm basically works as follows 
00190     (illustrated for 6-neighborhood, but 26-neighborhood works in the same way):
00191 
00192     <ol>
00193 
00194     <li> Find all candidate pixels that are 6-adjacent to a seed region.
00195     Calculate the cost for aggregating each candidate into its adajacent region
00196     and put the candidates into a priority queue.
00197 
00198     <li> While( priority queue is not empty)
00199 
00200         <ol>
00201 
00202         <li> Take the candidate with least cost from the queue. If it has not
00203         already been merged, merge it with it's adjacent region.
00204 
00205         <li> Put all candidates that are 4-adjacent to the pixel just processed
00206         into the priority queue.
00207 
00208         </ol>
00209 
00210     </ol>
00211 
00212     <tt>SRGType</tt> can take the following values:
00213     
00214     <DL>
00215     <DT><tt>CompleteGrow</tt> <DD> produce a complete tesselation of the volume (default).
00216     <DT><tt>KeepContours</tt> <DD> keep a 1-voxel wide unlabeled contour between all regions.
00217     <DT><tt>StopAtThreshold</tt> <DD> stop when the boundary indicator values exceed the 
00218                              threshold given by parameter <tt>max_cost</tt>.
00219     <DT><tt>KeepContours | StopAtThreshold</tt> <DD> keep 1-voxel wide contour and stop at given <tt>max_cost</tt>.
00220     </DL>
00221 
00222     The cost is determined jointly by the source array and the
00223     region statistics functor. The source array contains feature values for each
00224     pixel which will be used by the region statistics functor to calculate and
00225     update statistics for each region and to calculate the cost for each
00226     candidate. The <TT>RegionStatisticsArray</TT> must be compatible to the
00227     \ref ArrayOfRegionStatistics functor and contains an <em> array</em> of
00228     statistics objects for each region. The indices must correspond to the
00229     labels of the seed regions. The statistics for the initial regions must have
00230     been calculated prior to calling <TT>seededRegionGrowing3D()</TT>
00231 
00232     For each candidate
00233     <TT>x</TT> that is adjacent to region <TT>i</TT>, the algorithm will call
00234     <TT>stats[i].cost(as(x))</TT> to get the cost (where <TT>x</TT> is a <TT>SrcImageIterator</TT>
00235     and <TT>as</TT> is
00236     the SrcAccessor). When a candidate has been merged with a region, the
00237     statistics are updated by calling <TT>stats[i].operator()(as(x))</TT>. Since
00238     the <TT>RegionStatisticsArray</TT> is passed by reference, this will overwrite
00239     the original statistics.
00240 
00241     If a candidate could be merged into more than one regions with identical
00242     cost, the algorithm will favour the nearest region. If <tt>StopAtThreshold</tt> is active, 
00243     and the cost of the current candidate at any point in the algorithm exceeds the optional 
00244     <tt>max_cost</tt> value (which defaults to <tt>NumericTraits<double>::max()</tt>), 
00245     region growing is aborted, and all voxels not yet assigned to a region remain unlabeled.
00246 
00247     In some cases, the cost only depends on the feature value of the current
00248     voxel. Then the update operation will simply be a no-op, and the <TT>cost()</TT>
00249     function returns its argument. This behavior is implemented by the
00250     \ref SeedRgDirectValueFunctor.
00251 
00252     <b> Declarations:</b>
00253 
00254     pass arguments explicitly:
00255     \code
00256     namespace vigra {
00257         template <class SrcImageIterator, class Shape, class SrcAccessor,
00258                   class SeedImageIterator, class SeedAccessor,
00259                   class DestImageIterator, class DestAccessor,
00260                   class RegionStatisticsArray, class Neighborhood>
00261         void 
00262         seededRegionGrowing3D(SrcImageIterator srcul, Shape shape, SrcAccessor as,
00263                               SeedImageIterator seedsul, SeedAccessor aseeds,
00264                               DestImageIterator destul, DestAccessor ad,
00265                               RegionStatisticsArray & stats, 
00266                               SRGType srgType = CompleteGrow,
00267                               Neighborhood neighborhood = NeighborCode3DSix(),
00268                               double max_cost = NumericTraits<double>::max());
00269     }
00270     \endcode
00271 
00272     use argument objects in conjunction with \ref ArgumentObjectFactories :
00273     \code
00274     namespace vigra {
00275         template <class SrcImageIterator, class Shape, class SrcAccessor,
00276                   class SeedImageIterator, class SeedAccessor,
00277                   class DestImageIterator, class DestAccessor,
00278                   class RegionStatisticsArray, class Neighborhood>
00279         void
00280         seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> src,
00281                               pair<SeedImageIterator, SeedAccessor> seeds,
00282                               pair<DestImageIterator, DestAccessor> dest,
00283                               RegionStatisticsArray & stats, 
00284                               SRGType srgType = CompleteGrow,
00285                               Neighborhood neighborhood = NeighborCode3DSix(), 
00286                               double max_cost = NumericTraits<double>::max());
00287     }
00288     \endcode
00289 
00290 */
00291 doxygen_overloaded_function(template <...> void seededRegionGrowing3D)
00292 
00293 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
00294           class SeedImageIterator, class SeedAccessor,
00295           class DestImageIterator, class DestAccessor,
00296           class RegionStatisticsArray, class Neighborhood>
00297 void 
00298 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
00299                       SeedImageIterator seedsul, SeedAccessor aseeds,
00300                       DestImageIterator destul, DestAccessor ad,
00301                       RegionStatisticsArray & stats, 
00302                       SRGType srgType,
00303                       Neighborhood,
00304                       double max_cost)
00305 {
00306     SrcImageIterator srclr = srcul + shape;
00307     //int w = srclr.x - srcul.x;
00308     int w = shape[0];
00309     //int h = srclr.y - srcul.y;
00310     int h = shape[1];
00311     //int d = srclr.z - srcul.z;
00312     int d = shape[2];
00313     int count = 0;
00314 
00315     SrcImageIterator isy = srcul, isx = srcul, isz = srcul;  // iterators for the src image
00316 
00317     typedef typename RegionStatisticsArray::value_type RegionStatistics;
00318     typedef typename PromoteTraits<typename RegionStatistics::cost_type, double>::Promote CostType;
00319     typedef detail::SeedRgVoxel<CostType, Diff_type> Voxel;
00320 
00321     typename Voxel::Allocator allocator;
00322 
00323     typedef std::priority_queue< Voxel *,
00324                                  std::vector<Voxel *>,
00325                                  typename Voxel::Compare >  SeedRgVoxelHeap;
00326     typedef MultiArray<3, int> IVolume;
00327 
00328     // copy seed image in an image with border
00329     Diff_type regionshape = shape + Diff_type(2,2,2);
00330     IVolume regions(regionshape);
00331     MultiIterator<3,int> ir = regions.traverser_begin();
00332     ir = ir + Diff_type(1,1,1);
00333     
00334     //IVolume::Iterator iry, irx, irz;
00335     MultiIterator<3,int> iry, irx, irz;
00336 
00337     //initImageBorder(destImageRange(regions), 1, SRGWatershedLabel);
00338     initMultiArrayBorder(destMultiArrayRange(regions), 1, SRGWatershedLabel); 
00339     
00340     copyMultiArray(seedsul, Diff_type(w,h,d), aseeds, ir, AccessorTraits<int>::default_accessor());
00341 
00342     // allocate and init memory for the results
00343 
00344     SeedRgVoxelHeap pheap;
00345     int cneighbor;
00346 
00347     #if 0
00348     static const Diff_type dist[] = { Diff_type(-1, 0, 0), Diff_type( 0,-1, 0),
00349                                       Diff_type( 1, 0, 0), Diff_type( 0, 1, 0),
00350                                       Diff_type( 0, 0,-1), Diff_type( 0, 0, 1) };
00351     #endif
00352 
00353     typedef typename Neighborhood::Direction Direction;
00354     int directionCount = Neighborhood::DirectionCount;
00355 
00356     Diff_type pos(0,0,0);
00357 
00358     for(isz=srcul, irz=ir, pos[2]=0; pos[2]<d;
00359             pos[2]++, isz.dim2()++, irz.dim2()++)
00360     {
00361         //std::cerr << "Z = " << pos[2] << std::endl;
00362 
00363         for(isy=isz, iry=irz, pos[1]=0; pos[1]<h;
00364             pos[1]++, isy.dim1()++, iry.dim1()++)
00365         {
00366             //std::cerr << "Y = " << pos[1] << std::endl;
00367             
00368             for(isx=isy, irx=iry, pos[0]=0; pos[0]<w;
00369                 pos[0]++, isx.dim0()++, irx.dim0()++)
00370             {
00371                 //std::cerr << "X = " << pos[0] << std::endl;
00372                 
00373                 if(*irx == 0)
00374                 {
00375                     // find candidate pixels for growing and fill heap
00376                     for(int i=0; i<directionCount; i++)
00377                     {
00378                         cneighbor = *(irx + Neighborhood::diff((Direction)i));
00379                         if(cneighbor > 0)
00380                         {
00381                             CostType cost = stats[cneighbor].cost(as(isx));
00382 
00383                             Voxel * voxel =
00384                                 allocator.create(pos, pos+Neighborhood::diff((Direction)i), cost, count++, cneighbor);
00385                             pheap.push(voxel);
00386                         }
00387                     }
00388                 }
00389             }
00390         }
00391     }
00392     
00393     // perform region growing
00394     while(pheap.size() != 0)
00395     {
00396         Voxel * voxel = pheap.top();
00397         pheap.pop();
00398 
00399         if((srgType & StopAtThreshold) != 0 && voxel->cost_ > max_cost)
00400             break;
00401 
00402         Diff_type pos = voxel->location_;
00403         Diff_type nearest = voxel->nearest_;
00404         int lab = voxel->label_;
00405 
00406         allocator.dismiss(voxel);
00407 
00408         irx = ir + pos;
00409         isx = srcul + pos;
00410 
00411         if(*irx) // already labelled region / watershed?
00412             continue;
00413 
00414         if((srgType & KeepContours) != 0)
00415         {
00416             for(int i=0; i<directionCount; i++)
00417             {
00418                 cneighbor = * (irx + Neighborhood::diff((Direction)i));
00419                 if((cneighbor>0) && (cneighbor != lab))
00420                 {
00421                     lab = SRGWatershedLabel;
00422                     break;
00423                 }
00424             }
00425         }
00426 
00427         *irx = lab;
00428 
00429         if((srgType & KeepContours) == 0 || lab > 0)
00430         {
00431             // update statistics
00432             stats[*irx](as(isx));
00433 
00434             // search neighborhood
00435             // second pass: find new candidate pixels
00436             for(int i=0; i<directionCount; i++)
00437             {
00438                 if(*(irx + Neighborhood::diff((Direction)i)) == 0)
00439                 {
00440                     CostType cost = stats[lab].cost(as(isx, Neighborhood::diff((Direction)i)));
00441 
00442                     Voxel * new_voxel =
00443                         allocator.create(pos+Neighborhood::diff((Direction)i), nearest, cost, count++, lab);
00444                     pheap.push(new_voxel);
00445                 }
00446             }
00447         }
00448     }
00449 
00450     // write result
00451     transformMultiArray(ir, Diff_type(w,h,d), AccessorTraits<int>::default_accessor(), 
00452                         destul, ad, detail::UnlabelWatersheds());
00453 }
00454 
00455 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
00456           class SeedImageIterator, class SeedAccessor,
00457           class DestImageIterator, class DestAccessor,
00458           class RegionStatisticsArray, class Neighborhood >
00459 inline void
00460 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
00461                       SeedImageIterator seedsul, SeedAccessor aseeds,
00462                       DestImageIterator destul, DestAccessor ad,
00463                       RegionStatisticsArray & stats, SRGType srgType, Neighborhood n)
00464 {
00465     seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 
00466                            destul, ad, stats, srgType, n, NumericTraits<double>::max());
00467 }
00468 
00469 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
00470           class SeedImageIterator, class SeedAccessor,
00471           class DestImageIterator, class DestAccessor,
00472           class RegionStatisticsArray >
00473 inline void
00474 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
00475                       SeedImageIterator seedsul, SeedAccessor aseeds,
00476                       DestImageIterator destul, DestAccessor ad,
00477                       RegionStatisticsArray & stats, SRGType srgType)
00478 {
00479     seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, 
00480                            destul, ad, stats, srgType, NeighborCode3DSix());
00481 }
00482 
00483 template <class SrcImageIterator, class Diff_type, class SrcAccessor,
00484           class SeedImageIterator, class SeedAccessor,
00485           class DestImageIterator, class DestAccessor,
00486           class RegionStatisticsArray >
00487 inline void
00488 seededRegionGrowing3D(SrcImageIterator srcul, Diff_type shape, SrcAccessor as,
00489                       SeedImageIterator seedsul, SeedAccessor aseeds,
00490                       DestImageIterator destul, DestAccessor ad,
00491                       RegionStatisticsArray & stats)
00492 {
00493     seededRegionGrowing3D( srcul, shape, as, seedsul, aseeds, destul, ad, 
00494                            stats, CompleteGrow);
00495 }
00496 
00497 template <class SrcImageIterator, class Shape, class SrcAccessor,
00498           class SeedImageIterator, class SeedAccessor,
00499           class DestImageIterator, class DestAccessor,
00500           class RegionStatisticsArray, class Neighborhood>
00501 inline void
00502 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
00503                       pair<SeedImageIterator, SeedAccessor> img3,
00504                       pair<DestImageIterator, DestAccessor> img4,
00505                       RegionStatisticsArray & stats, 
00506                       SRGType srgType, Neighborhood n, double max_cost)
00507 {
00508     seededRegionGrowing3D(img1.first, img1.second, img1.third,
00509                           img3.first, img3.second,
00510                           img4.first, img4.second,
00511                           stats, srgType, n, max_cost);
00512 }
00513 
00514 template <class SrcImageIterator, class Shape, class SrcAccessor,
00515           class SeedImageIterator, class SeedAccessor,
00516           class DestImageIterator, class DestAccessor,
00517           class RegionStatisticsArray, class Neighborhood>
00518 inline void
00519 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
00520                       pair<SeedImageIterator, SeedAccessor> img3,
00521                       pair<DestImageIterator, DestAccessor> img4,
00522                       RegionStatisticsArray & stats, 
00523                       SRGType srgType, Neighborhood n)
00524 {
00525     seededRegionGrowing3D(img1.first, img1.second, img1.third,
00526                           img3.first, img3.second,
00527                           img4.first, img4.second,
00528                           stats, srgType, n, NumericTraits<double>::max());
00529 }
00530 
00531 template <class SrcImageIterator, class Shape, class SrcAccessor,
00532           class SeedImageIterator, class SeedAccessor,
00533           class DestImageIterator, class DestAccessor,
00534           class RegionStatisticsArray>
00535 inline void
00536 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
00537                       pair<SeedImageIterator, SeedAccessor> img3,
00538                       pair<DestImageIterator, DestAccessor> img4,
00539                       RegionStatisticsArray & stats, SRGType srgType)
00540 {
00541     seededRegionGrowing3D(img1.first, img1.second, img1.third,
00542                           img3.first, img3.second,
00543                           img4.first, img4.second,
00544                           stats, srgType, NeighborCode3DSix());
00545 }
00546 
00547 template <class SrcImageIterator, class Shape, class SrcAccessor,
00548           class SeedImageIterator, class SeedAccessor,
00549           class DestImageIterator, class DestAccessor,
00550           class RegionStatisticsArray>
00551 inline void
00552 seededRegionGrowing3D(triple<SrcImageIterator, Shape, SrcAccessor> img1,
00553                       pair<SeedImageIterator, SeedAccessor> img3,
00554                       pair<DestImageIterator, DestAccessor> img4,
00555                       RegionStatisticsArray & stats)
00556 {
00557     seededRegionGrowing3D(img1.first, img1.second, img1.third,
00558                           img3.first, img3.second,
00559                           img4.first, img4.second,
00560                           stats);
00561 }
00562 
00563 } // namespace vigra
00564 
00565 #endif // VIGRA_SEEDEDREGIONGROWING_HXX
00566 

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)