Update to latest code from website

CMakeLists.txt

@@ -97,6 +97,7 @@ set(GC_SOURCES
     Gc/System/Log.cpp    
     Gc/System/Collection/Tuple.h
     Gc/System/Collection/Array.h
+    Gc/System/Collection/Pair.h
     Gc/System/Collection/ArrayBuilder.h
     Gc/System/Collection/IArrayMask.h
     Gc/System/Collection/BoolArrayMask.h    

Doc/ChangeLog.txt

@@ -1,3 +1,19 @@
+2012-03-05 Ondrej Danek <ondrej.danek@gmail.com>
+    * Matlab examples: CreateNeighbourhood function removed as it is no longer necessary
+        when the Gc::Energy::Neighbourhood object supports also N98 neighbourhood in the
+        Common method.
+    * Gc/Flow/Grid/DanekLabels.h: Small documentation improvements.
+
+2012-03-03 Ondrej Danek <ondrej.danek@gmail.com>
+    * Gc/System/Collection/Pair.h: Pair container added.
+    * Gc/System/Algo/Sort/Heap.h: Simultaneous sorting removed (now pairs of values
+        are sorted to achieve the same effect) and predicate option added.
+    * Gc/Energy/Potential/Metric/RiemannianDanek.*: Slightly faster specialization of the
+        approximation in 2D added.
+
+2012-02-10 Ondrej Danek <ondrej.danek@gmail.com>
+    * Gc/Flow/Grid/DanekLabels.cpp: Bug fixed.
+
 2011-11-15 Ondrej Danek <ondrej.danek@gmail.com>
     * Gc/Energy/Neighbourhood.h: Added 98-neighbourhood to the Common() method. Exception
         is thrown when unexpected neighbourhood size is supplied.

Examples/Matlab/ChanVese/GcChanVese.cpp

@@ -48,8 +48,7 @@ static void Segment(const mxArray *mx_in, const mxArray *mx_mask, T l1, T l2, T
     }
 
     // Create neighbourhood object
-    Gc::Energy::Neighbourhood<N,Gc::Int32> nb;
-    Gc::Examples::Matlab::CreateNeighbourhood(str_nb, nb);
+    Gc::Energy::Neighbourhood<N,Gc::Int32> nb((Gc::Size)atoi(str_nb + 1), false);
 
     // Segment
     Gc::System::Collection::Array<N,bool> seg;

Examples/Matlab/ChanVese/GcChanVeseLab.cpp

@@ -59,8 +59,7 @@ static void Segment(const mxArray *mx_in, const mxArray *mx_interface, const mxA
     }
 
     // Create neighbourhood object
-    Gc::Energy::Neighbourhood<N,Gc::Int32> nb;
-    Gc::Examples::Matlab::CreateNeighbourhood(str_nb, nb);
+    Gc::Energy::Neighbourhood<N,Gc::Int32> nb((Gc::Size)atoi(str_nb + 1), false);
 
     // Segment
     Gc::System::Collection::Array<N,bool> seg;

Examples/Matlab/ChanVese/GcChanVeseTp.cpp

@@ -58,8 +58,7 @@ static void Segment(const mxArray *mx_in, const mxArray *mx_interface, const mxA
     }
 
     // Create neighbourhood object
-    Gc::Energy::Neighbourhood<N,Gc::Int32> nb;
-    Gc::Examples::Matlab::CreateNeighbourhood(str_nb, nb);
+    Gc::Energy::Neighbourhood<N,Gc::Int32> nb((Gc::Size)atoi(str_nb + 1), false);
 
     // Segment
     Gc::System::Collection::Array<N,bool> seg;

Examples/Matlab/ChanVese/gc_chan_vese_lab.m

@@ -1,5 +1,6 @@
-% gc_chan_vese_lab   Calculate label preserving Chan-Vese segmentation via 
-%                    graph cuts.
+% gc_chan_vese_lab   Calculate graph cut based Chan-Vese segmentation
+%                    integrated with object indication function.
+%                    
 %
 % SYNOPSIS:
 %  [seg energy iter fc1 fc2] = gc_chan_vese_lab(img, ilab, mu, lambda1,
@@ -14,10 +15,10 @@
 %  [seg2 energy iter c1 c2] = gc_chan_vese_lab(s2n, seg1, 1, 10, 10);
 %
 % PARAMETERS:
-%  ilab: Initial labeling of the image.
+%  ilab: Initial object indicators.
 %  See gc_chan_vese for the documentation of the remaining input parameters.
 %
-%  seg: Segmentation mask.
+%  seg: Final segmentation (object indicators).
 %  energy: Energy of the solution.
 %  iter: Number of iterations performed.
 %  fc1: Final background mean intensity.
@@ -57,8 +58,8 @@
 %  c1, c2     Initialized using gc_weighted_kmeans function
 %
 % DESCRIPTION:
-%  Computes a label preserving Chan-Vese segmentation of the input image 
-%  with given weights via graph-cut based optimization. TODO.
+%  Computes the Chan-Vese segmentation of the input image integrated with
+%  object indication function via graph-cut based optimization. TODO.
 %
 % NOTES:
 %  See gc_chan_vese.

Examples/Matlab/GcMatlabTools.h

@@ -21,43 +21,6 @@ namespace Gc
         {
             /*************************************************************************************/
 
-            // Create neighbourhood corresponding to a string identifier
-            // Supported neighbourhoods are: 
-            //  for 2D: N4, N8, N16, N32
-            //  for 3D: N6, N18, N26, N98
-            template <Gc::Size N>
-            void CreateNeighbourhood(const char *str_nb, Gc::Energy::Neighbourhood<N,Gc::Int32> &nb)
-            {
-                if (N == 2)
-                {
-                    if (!strcmp(str_nb, "N4") || !strcmp(str_nb, "N8") || 
-                        !strcmp(str_nb, "N16") || !strcmp(str_nb, "N32"))
-                    {
-                        nb.Common((Gc::Size)atoi(str_nb + 1), false);
-                        return;
-                    }
-                }
-
-                if (N == 3)
-                {
-                    if (!strcmp(str_nb, "N6") || !strcmp(str_nb, "N18") || 
-                        !strcmp(str_nb, "N26"))
-                    {
-                        nb.Common((Gc::Size)atoi(str_nb + 1), false);
-                        return;
-                    }
-                    if (!strcmp(str_nb, "N98"))
-                    {
-                        nb.Box(Gc::Math::Algebra::Vector<N,Gc::Size>(2), true, false);
-                        return;
-                    }
-                }
-
-                mexErrMsgTxt("Unsupported neighbourhood type for given image dimensionality.");
-            }
-
-            /*************************************************************************************/
-
             // Create general max-flow algorithm corresponging to given string identifier
             template <class T>
             static Gc::Flow::IMaxFlow<T,T,T> *CreateGeneralMaxFlow(const char *name)

Examples/Matlab/MumfordShah/GcMumfordShah.cpp

@@ -37,8 +37,7 @@ static void Segment(const mxArray *mx_in, Gc::Size k, T lambda, T conv,
     img.SetSpacing(Gc::Math::Algebra::Vector<N,T>(1));
 
     // Create neighbourhood object
-    Gc::Energy::Neighbourhood<N,Gc::Int32> nb;
-    Gc::Examples::Matlab::CreateNeighbourhood(str_nb, nb);
+    Gc::Energy::Neighbourhood<N,Gc::Int32> nb((Gc::Size)atoi(str_nb + 1), false);
 
     // Weights - same weight lambda for all partitions
     Gc::System::Collection::Array<1,T> l(k, lambda);

Examples/Matlab/RoussonDeriche/GcRoussonDeriche.cpp

@@ -35,8 +35,7 @@ static void Segment(const mxArray *mx_in, T lambda, T conv, Gc::Size max_iter,
     img.SetSpacing(Gc::Math::Algebra::Vector<N,T>(1));
 
     // Create neighbourhood object
-    Gc::Energy::Neighbourhood<N,Gc::Int32> nb;
-    Gc::Examples::Matlab::CreateNeighbourhood(str_nb, nb);
+    Gc::Energy::Neighbourhood<N,Gc::Int32> nb((Gc::Size)atoi(str_nb + 1), false);
 
     // Initialization
     Gc::Algo::Segmentation::RoussonDeriche::Params<T> rdpar;

Gc/Energy/Potential/Metric/RiemannianDanek.cpp

@@ -28,6 +28,7 @@
 
 #include "../../../Math/Constant.h"
 #include "../../../Math/Geometry/Voronoi.h"
+#include "../../../System/Algo/Sort/Heap.h"
 #include "RiemannianDanek.h"
 
 namespace Gc
@@ -84,11 +85,11 @@ namespace Gc
 
                     // Calculate delta rho and capacities
                     T coef = CauchyCroftonCoef<N,T>();
-                    T cell_area = mt.Determinant();
+                    T gridCellArea = mt.Determinant();
 
                     for (Size i = 0; i < m_nb.Elements(); i++)
                     {
-                        T drho = cell_area / mt.Mul(m_nb[i]).Length();
+                        T drho = gridCellArea / mt.Mul(m_nb[i]).Length();
                         m_rw[i] = (dphi[i] * drho) / coef;
                     }
 
@@ -106,6 +107,85 @@ namespace Gc
                 /** @endcond */
 
                 /******************************************************************************/
+
+                template <class T, class SZ>
+                struct OrientationPred
+                {
+                    bool operator()(const System::Collection::Pair<Math::Algebra::Vector<2,T>,SZ>& v1, 
+                        const System::Collection::Pair<Math::Algebra::Vector<2,T>,SZ>& v2) const
+                    {
+                        return (Math::Algebra::AngularOrientation(v1.m_first) < 
+                            Math::Algebra::AngularOrientation(v2.m_first));
+                    }
+                };
+
+                /******************************************************************************/
+
+                template <class T>
+                RiemannianDanek2D<T>::RiemannianDanek2D(const Neighbourhood<2,T> &n)
+                {
+                    // Init members
+                    m_rw.Resize(n.Elements());
+
+                    // Sort neighbourhood according to the angular orientation
+                    // but remember indexes to the original array
+                    m_nb.Resize(n.Elements());
+
+                    for (Size i = 0; i < n.Elements(); i++)
+                    {
+                        m_nb[i].m_first = n[i];
+                        m_nb[i].m_second = (Uint8)i;
+                    }
+
+                    // Sort angular orientations and the indexes
+                    System::Algo::Sort::Heap(m_nb.Begin(), m_nb.End(), OrientationPred<T,Uint8>());
+                }
+
+                /******************************************************************************/
+
+                template <class T>
+                RiemannianDanek2D<T>& RiemannianDanek2D<T>::SetTransformationMatrix
+                    (const Math::Algebra::SquareMatrix<2,T> &mt)
+                {                    
+                    Math::Algebra::Vector<2,T> vecPrev, vecCur, vecNext;
+                    T gridCellArea = mt.Determinant();
+
+                    // Working set of transformed vectors
+                    vecPrev = mt.Mul(m_nb[m_nb.Elements() - 1].m_first);
+                    vecCur = mt.Mul(m_nb[0].m_first);
+
+                    for (Size i = 0; i < m_nb.Elements(); i++)
+                    {
+                        // Following edge
+                        Size nextIdx = (i == m_nb.Elements() - 1) ? 0 : (i + 1);
+                        vecNext = mt.Mul(m_nb[nextIdx].m_first);
+
+                        // Delta phi
+                        T dphi = Math::Algebra::ClockwiseAngle(vecNext, vecPrev) / T(2);
+
+                        // Delta rho
+                        T drho = gridCellArea / vecCur.Length();
+
+                        // Edge weight
+                        m_rw[m_nb[i].m_second] = (dphi * drho) / T(2.0);
+
+                        // Shift working set
+                        vecPrev = vecCur;
+                        vecCur = vecNext;                        
+                    }
+
+                    return *this;
+                }
+
+                /******************************************************************************/
+
+                // Explicit instantiations
+                /** @cond */
+                template class GC_DLL_EXPORT RiemannianDanek2D<Float32>;
+                template class GC_DLL_EXPORT RiemannianDanek2D<Float64>;
+                /** @endcond */
+
+                /******************************************************************************/
             }
         }
     }

Gc/Energy/Potential/Metric/RiemannianDanek.h

@@ -32,6 +32,7 @@
 #include "../../../Core.h"
 #include "../../../Type.h"
 #include "../../../Math/Algebra/SquareMatrix.h"
+#include "../../../System/Collection/Pair.h"
 #include "../../Neighbourhood.h"
 
 namespace Gc
@@ -150,18 +151,60 @@ namespace Gc
                     RiemannianDanek<N,T>& SetTransformationMatrix
                         (const Math::Algebra::SquareMatrix<N,T> &mt);
 
-                    /** Get edge weights approximating the Euclidean metric. */
+                    /** Get edge weights approximating the Riemannian metric. */
                     const System::Collection::Array<1,T>& EdgeWeights() const
                     {
                         return m_rw;
                     }
 
-                    /** Get edge weight for i-th neighbour. */
+                    /** Get edge weight for i-th neighbourhood vector. */
                     T operator[] (Size i) const
                     {
                         return m_rw[i];
                     }
                 };
+
+                /** Faster specialization of RiemannianDanek for a 2D space. 
+
+                    This class does the same as RiemannianDanek. It offers faster recomputation
+                    when the transformation matrix changes, but supports only 2D space and requires 
+                    slightly more memory.
+
+                    @see RiemannianDanek.
+
+                    @tparam T Precision.
+                */
+                template <class T>
+                class GC_DLL_EXPORT RiemannianDanek2D
+                {
+                protected:
+                    /** Neighbourhood sorted according to the angular orientation and cross indexes. */
+                    System::Collection::Array<1,System::Collection::Pair<Math::Algebra::Vector<2,T>,Uint8> > m_nb;
+                    /** Calculated edge weights. */
+                    System::Collection::Array<1,T> m_rw;
+
+                public:
+                    /** Constructor. */
+                    RiemannianDanek2D(const Neighbourhood<2,T> &n);
+
+                    /** Specify the transformation of the Riemannian space. */
+                    RiemannianDanek2D<T>& SetTransformationMatrix
+                        (const Math::Algebra::SquareMatrix<2,T> &mt);
+
+                    /** Get edge weights approximating the Riemannian metric. */
+                    const System::Collection::Array<1,T>& EdgeWeights() const
+                    {
+                        return m_rw;
+                    }
+
+                    /** Get edge weight for i-th neighbourhood vector. */
+                    T operator[](Size i) const
+                    {
+                        return m_rw[i];
+                    }
+
+                protected:
+                };
             }
         }
     }

Gc/Flow/Grid/DanekLabels.cpp

@@ -338,15 +338,15 @@ namespace Gc
                         {
                             if (head->m_tag >= 8) // Free node
                             {
-                                head->m_parent_arc = Sister(i);
-                                head->m_tag = n->m_tag;
-                                head->m_label = n->m_label;
-
-                                if ((n->m_tag & 1) != (head->m_tag & 1))
+                                if (n->m_label != head->m_label)
                                 {
                                     head->m_dist = n->m_dist - 1;
                                 }
 
+                                head->m_parent_arc = Sister(i);
+                                head->m_tag = n->m_tag;
+                                head->m_label = n->m_label;
+
                                 EnqueueNode(head);
                             }
                             else if ((n->m_tag & 2) == (head->m_tag & 2)) // Same tree