|| ||           || Doxygen's  

---

HxWatershedSlow.h File Reference

More...

#include "HxSF.h"

Go to the source code of this file.

Functions
HxImageRep L_HXIMAGEREP	HxWatershedSlow (HxImageRep im, HxSF sf, HxString linereg)
	Watershed function creates the result image by detecting the domain of the catchment basins of "im" indicated by the "immersion definition", according to the connectivity defined by "sf". More...

---

Detailed Description

---

Function Documentation

HxImageRep L_HXIMAGEREP HxWatershedSlow (  HxImageRep    im, HxSF    sf, HxString    linereg )

Watershed function creates the result image by detecting the domain of the catchment basins of "im" indicated by the "immersion definition", according to the connectivity defined by "sf". According to the flag "linereg" the result image will be a labeled image of the catchment basins domain or just a binary image that presents the watershed lines. Implementation based on Vincent algorithm (PAMI) Evaluation using SDC Matlab toolbox (www.morph.com) im is the input image (HxImageRep) mark is the marker image (HxImageRep) sf is the structuring function; default type is crossSF (HxSF) linereg HxString. Possible values 'LINES' or ' REGIONS'. Default: "LINES". Returns: the watershed Image, ImageRep * references: L. Vincent and P. Soille, Watersheds in digital spaces: an efficient algorithm based on immersion simulations, IEEE Transactions on Pattern Analysis and Machine Intelligence. 13:583-598, 1991. Remarks: this queue based implementation is based on sorted lists of pixels for color images there is no partial ordering therefore we cannot use this for color { HxSizes sz = im.sizes(); HxImageSignature sig = im.signature(); //because there is no ordering in 3D space of color images, here we consider gray level images only // if( sig.pixelDimensionality() == 3 ) { im = HxRGB2Intensity(im); im = HxImageAsInt(im); HX_COUT << "this function is not implemented for color images!" << STD_ENDL; } int conn = sf.getConnectivity(); HxImageSignature sigOut; //default constructor, make a singature 2DINT //initialize the resltImage with -1 HxImageRep im0 = HxMakeFromValue(sigOut, sz, -1); HxImageRep imd = HxMakeFromValue(sigOut, sz, 0); int currentLabel = 0; int currentDist; HxPoint2 p; HxPoint2 p1(-1,-1); HxImgAsListIterator p1lIter, p2lIter; //fifo std::list< HxPoint2> fifo; int mask = -2; //the initial value of a threshold level int wshed = -3; //0; //value of pixels belonging to the watersheds int h, hmin,hmax; hmin = HxPixMin(im).HxScalarIntValue().x(); //what about color image? hmax = HxPixMax(im).HxScalarIntValue().x(); //sort the pixels in im in the increasing order of their gray values //use HxGetPoints, HxGetValues, and sort these lists HxValueList vl; HxPointList pl; HxPointListIter plIter; HxValueListIter vlIter; HxTimer timer; timer.start(); HxGetPoints(im, std::back_inserter(pl)); HxGetValues(im, pl.begin(), pl.end(), std::back_inserter(vl)); HxImgAsList iml; // printf("\nSTART imlsize=%d\n", iml.size()); for(plIter = pl.begin(), vlIter = vl.begin(); plIter!= pl.end() && vlIter!= vl.end(); plIter++, vlIter++) { HxPointAndValue pv((*plIter), (*vlIter).HxScalarIntValue()); iml.push_back(pv); } //this is not possible because HxValue does not implement the < operator //vl.sort(); // printf("\nimlsize=%d\n", iml.size()); // iml.sort(); //sort the vector, not the list std::sort(iml.begin(), iml.end()); //printf("\n"); timer.stop(); show(timer); timer.start(); // printf("\nimlsize=%d\n", iml.size()); // printf("hmin=%d hmax=%d\n", hmin, hmax); HxImgAsListIterator imlIter,imlIter2, startIter, startIterMin ; startIter = iml.begin(); startIterMin = iml.begin(); //<<1>> for(h = hmin; h<=hmax; h++) { bool hIsNotPresent = false; bool processed = false; //geodesic SKIZ of level h-1 inside level h for(imlIter=startIter; imlIter!=iml.end(); imlIter++) { if( (*imlIter).v == h) //LT here you can speedup using a set instead of list of values { processed = true; p = (*imlIter).p; //im0(p) <- mask im0.setAt(mask, p.x(), p.y() ); //now check the vecinity of pixel p in im0 if( HxCompareNeighborPixels(im0, p, conn, GT, 0) || //if one of the nighbor is already labeled. HxCompareNeighborPixels(im0, p, conn, EQ, wshed) ) { imd.setAt(1, p.x(), p.y() ); //and put 'p' in fifo list fifo.push_back(p); } } else //because the image is sorted by pixelvalues, you can speed up by breaking this loop { if((*imlIter).v > h && !processed) //this level doesn't exist in the image, so break the main loop { hIsNotPresent = true; break; } if(processed) { startIter = imlIter; break; } } } if(hIsNotPresent) continue; // printf("h=%3d: fifo size=%7d\n",h, fifo.size()); //<<2>> currentDist=1; HxPoint2 fictionPoint(-1,-1); fifo.push_back(fictionPoint); while(true) { p = *fifo.begin(); fifo.erase(fifo.begin()); if ( p==fictionPoint) { if(fifo.empty()) break; else { fifo.push_back(fictionPoint); currentDist++; p = *fifo.begin(); fifo.erase(fifo.begin()); } } // printf("currentDist=%d\n",currentDist); //<<3>> //hey! what if point p is not correct? HxImgAsList p1l=HxGetNeighborPixels(im0, p, conn); // if(!p1l.empty()) for(p1lIter=p1l.begin(); p1lIter!=p1l.end(); p1lIter++) { p1 = (*p1lIter).p; if( (imd.getAt(p1.x(),p1.y()).HxScalarIntValue() < currentDist ) && ((im0.getAt(p1.x(),p1.y()).HxScalarIntValue() > 0 ) || (im0.getAt(p1.x(),p1.y()).HxScalarIntValue() == wshed )) ) { //i.e. p' belongs to an already labeled basin or to the watersheds //then p get the value of p1 or 'wshed' value if( im0.getAt(p1.x(),p1.y()).HxScalarIntValue() > 0 ) //a basin point { if( (im0.getAt(p.x(),p.y()).HxScalarIntValue() == mask) || (im0.getAt(p.x(),p.y()).HxScalarIntValue() == wshed) ) im0.setAt(im0.getAt(p1.x(),p1.y()).HxScalarIntValue(),p.x(),p.y()); else if( im0.getAt(p.x(),p.y()).HxScalarIntValue() != im0.getAt(p1.x(),p1.y()).HxScalarIntValue() ) //im0(p)!=im0(p') im0.setAt(wshed,p.x(),p.y()); } else //i.e. p' is a wshed point //then check the value of p if(im0.getAt(p.x(),p.y()).HxScalarIntValue() == mask ) im0.setAt(wshed,p.x(),p.y()); } else if( (im0.getAt(p1.x(),p1.y()).HxScalarIntValue() == mask) && (imd.getAt(p1.x(),p1.y()).HxScalarIntValue() == 0) ) { imd.setAt(currentDist+1,p1.x(),p1.y()); fifo.push_back(p1); } } }//end while(true) -> endless loop; break when fifo is empty //lt for debug //export the distance image // HxImageRep imdE = HxImageAsByte(imd); // HxWriteFile(imdE, "dist.tif"); //<<4>> // check if new minima have been discovered bool minFound=false; for(imlIter2 = startIterMin; imlIter2 != iml.end(); imlIter2++) { if( (*imlIter2).v.x() == h) { minFound = true; p = (*imlIter2).p; //the distance associated with p is reset to 0 imd.setAt(0, p.x(),p.y()); if(im0.getAt(p.x(),p.y()).HxScalarIntValue().x() == mask ) { currentLabel++; fifo.push_back(p); im0.setAt(currentLabel, p.x(),p.y()); while(!fifo.empty()) { p1 = *fifo.begin(); fifo.erase(fifo.begin()); HxPoint2 p2; HxImgAsList p2l = HxGetNeighborPixels(im0, p1, conn); for(p2lIter=p2l.begin(); p2lIter!=p2l.end(); p2lIter++) { p2 = (*p2lIter).p; if(im0.getAt(p2.x(),p2.y()).HxScalarIntValue().x() == mask) { fifo.push_back(p2); im0.setAt(currentLabel, p2.x(),p2.y()); } } } } } else { if(minFound) { startIterMin=imlIter2; break; } } } // HxImageRep im0E = HxImageAsByte(im0); // HxWriteFile(im0E, "im0before4.tif"); }//end for every h level in [hmin,hmax] timer.stop(); show(timer); /*************************/ //end lucVincent algorithm, now make the tesselation image /************************/ //implement this with neighbourhood operators timer.start(); HxImageRep imres = im0; pl.clear(); vl.clear(); //take care HxGetPoints doesn't include the zero pixels! //but here I changed whsed value from 0 to -3 HxGetPoints(im0, std::back_inserter(pl)); HxGetValues(im0, pl.begin(), pl.end(), std::back_inserter(vl)); // printf("plsize=%d, vlsize=%d\n",pl.size(), vl.size()); STD_COUT << "imres : " << imres.signature() << ", sizes: " << imres.sizes() << STD_ENDL; for(plIter = pl.begin(), vlIter = vl.begin(); plIter!= pl.end() && vlIter!= vl.end(); plIter++, vlIter++) { HxPoint2 pc = (*plIter); int v = (*vlIter).HxScalarIntValue().x(); if( v == wshed) { //leon it is better to give to wshed pixel the label of //the smallest label among neighbors - this is just to be consistent HxImgAsList p0l=HxGetNeighborPixels(im0, pc, conn); int maxValNg = 0; for(p1lIter=p0l.begin(); p1lIter!=p0l.end(); p1lIter++) { if(( *p1lIter).v.x() > maxValNg) maxValNg = (*p1lIter).v.x(); } STD_COUT << " in loop before setAt" << STD_ENDL; if(maxValNg > 0 ) { STD_COUT << "values : " << maxValNg << "," << pc.x() << "," << pc.y() << STD_ENDL; HxValue maxValNgV(maxValNg); int x = pc.x(); int y = pc.y(); STD_COUT << "values local : " << maxValNgV << "," << x << "," << y << STD_ENDL; HxValue orgV = imres.getAt(x, y); STD_COUT << "org value : " << orgV << STD_ENDL; //imres.setAt(maxValNg, pc.x(),pc.y()); //imres.setAt(maxValNgV, x, y); imres.setAt(orgV, x, y); } STD_COUT << " in loop after setAt" << STD_ENDL; } } im0 = imres; //to separate catchment basines if(linereg=="LINES"){ //im0 = HxLWshed(im0,conn,wshed); HxTagList tags; HxAddTag(tags, "conn", conn); HxAddTag(tags, "wshedval", wshed); im0 = im0.neighbourhoodOp("lwshed", tags); //this can be implemented iwth LUT, some unary pixel operation anyway... pl.clear(); vl.clear(); HxGetPoints(im0, std::back_inserter(pl)); HxGetValues(im0, pl.begin(), pl.end(), std::back_inserter(vl)); for(plIter = pl.begin(), vlIter = vl.begin(); plIter!= pl.end() && vlIter!= vl.end(); plIter++, vlIter++) { HxPoint2 pc = (*plIter); int v = (*vlIter).HxScalarIntValue().x(); if(v == wshed ) imres.setAt(255, pc.x(),pc.y()); else imres.setAt(0, pc.x(),pc.y()); } im0=imres; } if(linereg=="REGIONS"){ //implemented with neighbourhood operator //im0 = HxLWshed(im0,conn,wshed); HxTagList tags; HxAddTag(tags, "conn", conn); HxAddTag(tags, "wshedval", wshed); im0 = im0.neighbourhoodOp("lwshed", tags); pl.clear(); vl.clear(); HxGetPoints(im0, std::back_inserter(pl)); HxGetValues(im0, pl.begin(), pl.end(), std::back_inserter(vl)); for(plIter = pl.begin(), vlIter = vl.begin(); plIter!= pl.end() && vlIter!= vl.end(); plIter++, vlIter++) { HxPoint2 pc = (*plIter); int v = (*vlIter).HxScalarIntValue().x(); if(v == wshed ) imres.setAt(0, pc.x(),pc.y()); } im0=imres; } timer.stop(); show(timer); return im0; }

---