|| || || || || || || ||

---

QbWatershed.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_Trait_QbWatershed_h
#define Impala_Core_Array_Trait_QbWatershed_h

#include "Core/Geometry/PointZ.h"
#include "Core/Array/Pattern/FcvArray.h"
#include "Core/Array/Pattern/QCollection.h"
#include "Core/Array/Pattern/Array2D.h"
#include "Core/Array/Arrays.h"

namespace Impala
{
namespace Core
{
namespace Array
{
namespace Trait
{


const int CONTOURVAL    = -3;
const int MASKVAL       = -2;
const int INITVAL       = -1;
const int WSHEDVAL      =  0;


template<class DstArrayT, class Src1ArrayT, class Src2ArrayT>
class QbWatershed {

    typedef typename DstArrayT::ArithType DstArithT;
    typedef typename Src1ArrayT::ArithType Src1ArithT;
    typedef typename Src2ArrayT::ArithType Src2ArithT;

public:

    typedef struct structPointValue {
        DstArithT arith; 
        Src1ArithT img; 
        Src1ArithT mask; 
        Geometry::PointZ point;
        int cost;
        int lastDist;

        bool
        less(const struct structPointValue &rhs) const 
        { 
            if (cost<rhs.cost)
                return true;
            if (cost>rhs.cost)
                return false;
            // for equal values be a FIFO=QUEUE
            return localOrderCounter < rhs.localOrderCounter; 
        }

        bool
        operator < (const struct structPointValue &rhs) const
        { 
//          return !less(rhs); // for priority_queue returning highest value
            return less(rhs); // for set/multiset returning lowest value
        }
 
        unsigned localOrderCounter;
    } PointValueT;

    typedef Pattern::Qset<PointValueT> QT;

    typedef struct {
        Geometry::PointZ point;
        DstArithT arith;
        Src1ArithT img;
        Src1ArithT mask;
    } NeighborsT;
 
    typedef Pattern::FcvArray<NeighborsT> VecNeighborsT;

    QbWatershed(int width, int height, int minVal, int maxVal,
                Array2dScalarInt32* exportDistImage = 0)
        : costImage(width, height, 0), distImage(width, height, 0)
    { 
        labelcode=0;
        globalOrderCounter=0;
        curdist=1;
        h = hmin = minVal;
        hmax = maxVal;
        for(int i=0;i<256;i++)
            hist[i]=false;
        exportDistIm = exportDistImage;
    }

    bool
    GlobalPixelInit(PointValueT&vp, DstArithT& arith, Src1ArithT img, Src1ArithT)
    {
        arith = INITVAL;
        return false;
    }

    bool
    WantAnotherLoop()
    { 
        h++; 
        if (h<=hmax) 
            return true; 
        else {
            //exports distImage if that was required by the user
            if (exportDistIm) {
                for (int y=0 ; y<distImage.getDim(1) ; y++)
                    for (int x=0 ; x<distImage.getDim(0) ; x++)
                        *Pattern::ArrayCPB(exportDistIm, x, y) = distImage(x, y);
            }
            return false;
        }
    }

    bool
    WantFreshStartLocalPixelInit()
    {
        return true; //this is the real action of this function
    }

    bool
    LocalPixelInit(PointValueT& vp, DstArithT& arith, Src1ArithT img,
                   Src1ArithT mask, bool& continua)
    { 
        vp.arith=arith;
        vp.img=img;
        vp.mask=mask;
        continua = true;
        bool q=false;

//when I'm first time in this function (for level hmin) I fill the histogram
        if (h==hmin) {
            hist[img]=true;
        }
//then, chech if the current level h exist in the image
        if (h!=hmin && !hist[h] ) {
            continua = false;
            return q;
        }

//init for searching new local minima (with this one starts the algorithm) (point 4)
        if (img == h && arith == INITVAL && arith != CONTOURVAL ) {
            q=true; // do queue
            distImage(vp.point.mX, vp.point.mY)=0; // this I don't have to do, it
                                                 // was done in globalinit
            costImage(vp.point.mX, vp.point.mY)=1;
            vp.arith = MASKVAL;
            vp.lastDist = 0;
        }
//init for the part of propagating "contour" pixels (point 1)
        if (img==h && (arith == CONTOURVAL || arith == WSHEDVAL) ) {
            distImage(vp.point.mX, vp.point.mY)= 1;
            costImage(vp.point.mX, vp.point.mY)= -2;
            curdist=1;
            vp.lastDist=1;
            q=true; // do queue
        }

        if (q) { 
            vp.localOrderCounter=globalOrderCounter++; // for equal values
                                                       // be a FIFO=QUEUE
            vp.cost=costImage(vp.point.mX, vp.point.mY);
        }
        arith=vp.arith;
        return q;
    }

    // first from queue
    void
    First(const PointValueT& val, DstArithT arith, Src1ArithT img, Src1ArithT)
    { 
        center=val; 
        center.arith=arith; //take  the current value from the image, not from queue
    }

    // sequentially read point/value and action to be done; 
    // return false if no more data
    Pattern::QAction
    Result()
    {
        if (index>=neighboursout.size())
            return Pattern::STOP; // quit asking 
        resultVP=neighboursout[index];
        Pattern::QAction action=actionsout[index];
        index++;
        return action;
    }

    void
    GetItemToQueue(PointValueT& vp)
    {
        vp=resultVP;
    }

    // speed up searching in set<less..>, provide low/up bound where
    // resultVP is never out
    void
    GetItemToRemove(PointValueT& lower, PointValueT& upper)
    {
        lower.cost=upper.cost=resultVP.cost;
        lower.localOrderCounter=0; // unsigned int
        upper.localOrderCounter=4000000000;
    }

    bool
    KillThisOne(const PointValueT& vp)
    {
        if (resultVP.cost==vp.cost && resultVP.point==vp.point) 
            return true;
        return false; 
    }

    void
    GetItemToWrite(Geometry::PointZ& point, DstArithT& arith)
    {
        point=resultVP.point;
        arith=resultVP.arith;
    }

    void
    Calculate(VecNeighborsT& neighboursin)
    {
        index=0; // for result()
        neighboursout.clear();
        actionsout.clear();

        int i;
        bool step2=false;
//here starts the propagation part (in algorithm point 2 to 4)
        if (!neighboursin.empty() && (center.arith == CONTOURVAL)) {
            // here we increment curDist, every time a new layer of pixels is
            // gonna be processed.
            if (center.lastDist > curdist) {
                curdist++;
            }

            step2=true;
            for (i=0; i<neighboursin.size(); i++) {
                NeighborsT val=neighboursin[i];
                PointValueT neighbor;
                neighbor.point=val.point; 
                neighbor.arith=val.arith; 
                neighbor.img=val.img; 
                neighbor.mask=val.mask;
                neighbor.cost=costImage(neighbor.point.mX, neighbor.point.mY);
                int dn= distImage(neighbor.point.mX, neighbor.point.mY);
                if ((neighbor.arith > 0 || neighbor.arith == WSHEDVAL
                       || neighbor.arith == CONTOURVAL)) { //not labeled yet
                    if (neighbor.arith > 0) {
                        if (center.arith == WSHEDVAL
                                 || center.arith == CONTOURVAL  ) {
                            center.arith = neighbor.arith;
                            neighboursout.push_back(center);
                            actionsout.push_back(Pattern::WRITE);
                        } else if (center.arith != neighbor.arith) {
                            center.arith = WSHEDVAL;
                            neighboursout.push_back(center);
                            actionsout.push_back(Pattern::WRITE);
                        }
                    } else if (center.arith < 0 && center.img == h) {
                        center.arith = WSHEDVAL;
                        neighboursout.push_back(center);
                        actionsout.push_back(Pattern::WRITE);
                    }
                }
                else if((neighbor.arith <0 && neighbor.img == h) 
                            && distImage(neighbor.point.mX, neighbor.point.mY) == 0) {
                    // nu aici trebuie incrementat current distance, sau in orice
                    // caz nu de fiecare data cumva trebuie sa semnalizez cind
                    // iau de pe stiva, primul pixel pus aici
                    // pot folosi de exemplu cimputl cost...
                    distImage(neighbor.point.mX, neighbor.point.mY)=curdist+1; 
                    neighbor.lastDist=distImage(neighbor.point.mX, neighbor.point.mY);
                    costImage(neighbor.point.mX, neighbor.point.mY) = -1;// curdist+2;
                    neighbor.localOrderCounter=globalOrderCounter++; // for equal
                    // values be a FIFO=QUEUE
                    neighbor.cost=costImage(neighbor.point.mX, neighbor.point.mY);
                    neighboursout.push_back(neighbor);
                    actionsout.push_back(Pattern::QUEUE);

                    neighbor.arith = CONTOURVAL;
                    neighboursout.push_back(neighbor);
                    actionsout.push_back(Pattern::WRITE);

                }
                //leon a adaugat asta, ca sa markez pixeli de contur
                //else
                if (neighbor.img!=center.img && center.img==h 
                        && neighbor.arith<0) {
                    // we mark this pixel, to be treated later, but we do not put
                    // it on Queue
                    neighbor.arith=CONTOURVAL; // line 9
                    distImage(neighbor.point.mX, neighbor.point.mY)= 1; 
                    costImage(neighbor.point.mX, neighbor.point.mY)= -2;
                    neighbor.localOrderCounter=globalOrderCounter++; // for equal
                    // values be a FIFO=QUEUE
                    neighbor.cost=costImage(neighbor.point.mX, neighbor.point.mY);
                    neighboursout.push_back(neighbor);
                    actionsout.push_back(Pattern::WRITE);
                }

            }
        }


//the next LOOP(2) should start only AFTER ALL pixels from the queue were tried
// to be processed with LOOP1
//here starts the propagation part searching for local minima (in algorithm point 4)

        if (!step2) {
            //if not labeled yet
            if(!neighboursin.empty() && center.arith == MASKVAL) {
                labelcode++;
                center.arith = labelcode;
                neighboursout.push_back(center);
                actionsout.push_back(Pattern::WRITE);
            } //end if(!neighboursin.empty() && center.arith == MASKVAL)

            if(!neighboursin.empty() && center.arith > 0) {
                for (i=0; i<neighboursin.size(); i++) {
                    NeighborsT val=neighboursin[i];
                    PointValueT neighbor;
                    neighbor.point=val.point; 
                    neighbor.arith=val.arith; 
                    neighbor.img=val.img; 
                    neighbor.mask=val.mask;
                    neighbor.cost=costImage(neighbor.point.mX, neighbor.point.mY);

                    if(neighbor.arith == MASKVAL 
                        || neighbor.arith == INITVAL )
//if we uncomment this, the watershed lines will be removed
//                      || neighbor.arith == WSHEDVAL
//                      || neighbor.arith == CONTOURVAL) //not labeled yet
                    {
                        if(neighbor.img==center.img && center.img==h ) {
                            neighbor.arith=center.arith; // line 9
                            costImage(neighbor.point.mX, neighbor.point.mY)= 0;
                            neighbor.localOrderCounter=globalOrderCounter++; // for
                            // equal values be a FIFO=QUEUE
                            neighbor.cost=costImage(neighbor.point.mX,
                                                    neighbor.point.mY);
                            neighboursout.push_back(neighbor);
                            actionsout.push_back(Pattern::WRITE);
                            neighboursout.push_back(neighbor);
                            actionsout.push_back(Pattern::QUEUE);
                        }
                        else if(neighbor.img!=center.img && neighbor.arith<0) {
                            // this is contour pixel
                            // we mark this pixel, to be treated by the other loop,
                            // but we do not put it on Queue
                            neighbor.arith=CONTOURVAL; // line 9
                            distImage(neighbor.point.mX, neighbor.point.mY)= 1; 
                            costImage(neighbor.point.mX, neighbor.point.mY)= -2;
                            neighbor.localOrderCounter=globalOrderCounter++; // for
                            // equal values be a FIFO=QUEUE
                            neighbor.cost=costImage(neighbor.point.mX,
                                                    neighbor.point.mY);
                            neighboursout.push_back(neighbor);
                            actionsout.push_back(Pattern::WRITE);
                        }
                    }
                }
            }
        }
//here ends the propagation part searching for local minima (in algorithm point 4)
        
    } //end calculate

private:

    PointValueT center;
    Pattern::FcvArray<PointValueT> neighboursout;
    Pattern::FcvArray<Pattern::QAction> actionsout;
    int index;
    int labelcode;
    int curdist;
    PointValueT resultVP;

    int h;  //current gray value
    int hmin, hmax; //the min and max grey value in the image

    //we asume the input image has only 256 levels of gray
    bool hist[256];

    unsigned globalOrderCounter;
    Pattern::Array2D<int> costImage;
    Pattern::Array2D<int> distImage; 
    Array2dScalarInt32* exportDistIm;
};

} // namespace Trait
} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---