diff --git a/Code/FBSD/BSTools/bsdetectionwidget.cpp b/Code/FBSD/BSTools/bsdetectionwidget.cpp
index be62ee7c0dfe34cd96de781ef3e0bf844d2f8867..8ad571e06e81506dc5ec9ba2a82c9bbeb9806cac 100755
--- a/Code/FBSD/BSTools/bsdetectionwidget.cpp
+++ b/Code/FBSD/BSTools/bsdetectionwidget.cpp
@@ -55,7 +55,8 @@ BSDetectionWidget::BSDetectionWidget (QWidget *parent)
   bsHighColor[1] = Qt::blue;
   bsLowColor[1] = Qt::magenta;
   bsHighColor[2] = Qt::yellow;
-  bsLowColor[2] = Qt::red;
+//  bsLowColor[2] = Qt::red;
+  bsLowColor[2] = Qt::black;
   bsHighColor[3] = Qt::white;
   bsLowColor[3] = Qt::gray;
   selectionColor = Qt::red;
@@ -404,7 +405,7 @@ void BSDetectionWidget::mousePressEvent (QMouseEvent *event)
     oldp1.set (p1);
     oldp2.set (p2);
     oldudef = udef;
-    p1 = Pt2i (ex, ey);
+    p1.set (ex, ey);
     if (p1.manhattan (p2) < 10) p1.set (oldp1);
     else if (p1.manhattan (oldp1) < 10) p1.set (p2);
     udef = true;
@@ -418,7 +419,7 @@ void BSDetectionWidget::mouseReleaseEvent (QMouseEvent *event)
   {
     int ex = zoom * (event->pos().x () - xShift);
     int ey = zoom * (event->pos().y () - yShift);
-    p2 = Pt2i (ex, height - 1 - ey);
+    p2.set (ex, height - 1 - ey);
     if (p1.equals (p2))
     {
       p1.set (oldp1);
@@ -443,7 +444,7 @@ void BSDetectionWidget::mouseMoveEvent (QMouseEvent *event)
   {
     int ex = zoom * (event->pos().x () - xShift);
     int ey = zoom * (event->pos().y () - yShift);
-    p2 = Pt2i (ex, height - 1 - ey);
+    p2.set (ex, height - 1 - ey);
     if (verbose) cerr << "(" << p1.x () << ", " << p1.y () << ") ("
                       << p2.x () << ", " << p2.y () << ")" << endl;
     if (p1.manhattan (p2) > 5
@@ -552,11 +553,15 @@ void BSDetectionWidget::keyPressEvent (QKeyEvent *event)
     case Qt::Key_F :
       if (event->modifiers () & Qt::ControlModifier)
       {
-        // Switches initial detection filtering
-        detector.switchFiltering (BSDetector::STEP_INITIAL);
-        cout << "Pre-filtering "
-             << (detector.isFiltering (BSDetector::STEP_INITIAL) ? "on" : "off")
-             << endl;
+        detector.switchNFA ();
+        cout << "NFA-filtering " << (detector.isNFA () ? "on" : "off") << endl;
+        extract ();
+      }
+      else if (detector.isNFA ())
+      {
+        detector.incNfaLengthRatio (
+          event->modifiers () & Qt::ShiftModifier ? -1 : 1);
+        cout << "NFA-lratio: " << detector.nfaLengthRatio () << endl;
         extract ();
       }
       break;
@@ -582,18 +587,6 @@ void BSDetectionWidget::keyPressEvent (QKeyEvent *event)
       }
       break;
 
-    case Qt::Key_H :
-      if (event->modifiers () & Qt::ControlModifier)
-      {
-        // Switches final detection filtering
-        detector.switchFiltering (BSDetector::STEP_FINAL);
-        cout << "Final filtering "
-             << (detector.isFiltering (BSDetector::STEP_FINAL) ? "on" : "off")
-             << endl;
-        extract ();
-      }
-      break;
-
 // DEV IN
     case Qt::Key_J :
       if (event->modifiers () & Qt::ControlModifier)
@@ -1154,6 +1147,7 @@ void BSDetectionWidget::drawBlurredSegment (QPainter &painter,
     if (bsDisplay % 2 == 1) // Points
       drawPoints (painter, bs->getAllPoints (),
                   high ? bsHighColor[bsColor] : bsLowColor[bsColor]);
+
   }
 }
 
@@ -1351,15 +1345,44 @@ void BSDetectionWidget::displayDetectionResult ()
   vector<BlurredSegment *> bss = detector.getBlurredSegments ();
   if (! bss.empty ())
   {
-    if (arlequin != 0) srand (time (NULL));
-    vector<BlurredSegment *>::const_iterator it = bss.begin ();
-    while (it != bss.end ())
+    if (arlequin != 0)
     {
-      if (arlequin != 0) drawArlequinSegment (painter, *it);
-      else
+      srand (time (NULL));
+      vector<BlurredSegment *>::const_iterator it = bss.begin ();
+      while (it != bss.end ())
+      {
+        drawArlequinSegment (painter, *it);
+        it++;
+      }
+    }
+    else if (detector.isNFA ())
+    {
+      vector<BlurredSegment *> vss = detector.getValidSegments ();
+cout << vss.size () << " VALID" << endl;
+      vector<BlurredSegment *>::const_iterator it = vss.begin ();
+      while (it != vss.end ())
+      {
+        drawBlurredSegment (painter, *it, true);
+        it++;
+      }
+      vector<BlurredSegment *> rss = detector.getRejectedSegments ();
+cout << rss.size () << " REJEC" << endl;
+      it = rss.begin ();
+      while (it != rss.end ())
+      {
+        drawBlurredSegment (painter, *it, false);
+        it++;
+      }
+    }
+    else
+    {
+      vector<BlurredSegment *>::const_iterator it = bss.begin ();
+      while (it != bss.end ())
+      {
         drawBlurredSegment (painter, *it,
                  detector.getMaxDetections () == 0 || *it == bss.back ());
-      it++;
+        it++;
+      }
     }
   }
   else drawBlurredSegment (painter, detector.getBlurredSegment ());
@@ -1423,8 +1446,8 @@ int BSDetectionWidget::saveExtractedSegment ()
     {
       ExtractedSegment es;
       es.bs = bs;
-      es.p1 = p1;
-      es.p2 = p2;
+      es.p1.set (p1);
+      es.p2.set (p2);
       extractedSegments.push_back (es);
       detector.preserveFormerBlurredSegment ();
       count = 1;
@@ -1437,8 +1460,8 @@ int BSDetectionWidget::saveExtractedSegment ()
     {
       ExtractedSegment es;
       es.bs = *it++;
-      es.p1 = p1;
-      es.p2 = p2;
+      es.p1.set (p1);
+      es.p2.set (p2);
       extractedSegments.push_back (es);
     }
     detector.preserveFormerBlurredSegments ();
@@ -1761,8 +1784,8 @@ void BSDetectionWidget::localTest ()
     if (i == 4)
     {
       udef = true;
-      p1 = Pt2i (val[0], val[1]);
-      p2 = Pt2i (val[2], val[3]);
+      p1.set (val[0], val[1]);
+      p2.set (val[2], val[3]);
       cout << "Run test on (" << val[0] << ", " << val[1]
            << ") (" << val[2] << ", " << val[3] << ")" << endl;
     }
diff --git a/Code/FBSD/BSTools/bsprofileitem.cpp b/Code/FBSD/BSTools/bsprofileitem.cpp
index 5d32163fcf7f80ac9d0e49139920bd9b48b23d26..b41ddce8c9a3426b25ca58be9a7dacafb7d33df6 100755
--- a/Code/FBSD/BSTools/bsprofileitem.cpp
+++ b/Code/FBSD/BSTools/bsprofileitem.cpp
@@ -115,8 +115,8 @@ void BSProfileItem::setImage (QImage *image, VMap *idata)
 void BSProfileItem::buildScans (Pt2i p1, Pt2i p2)
 {
   // Updates the central scan end points for parallel display
-  this->pt1 = p1;
-  this->pt2 = p2;
+  this->pt1.set (p1);
+  this->pt2.set (p2);
 
   // Resets the profiles
   rightscan.clear ();
diff --git a/Code/FBSD/BSTools/bsrandomtester.cpp b/Code/FBSD/BSTools/bsrandomtester.cpp
index f99e5ff537879e200beca09ffec646b390c9b2df..853cf4d6b4d722689ac37eed57b3028e93868d96 100755
--- a/Code/FBSD/BSTools/bsrandomtester.cpp
+++ b/Code/FBSD/BSTools/bsrandomtester.cpp
@@ -1029,7 +1029,7 @@ void BSRandomTester::generateImage ()
       if (rp1[i].chessboard (rp2[i]) < sminlength) nok = true;
       else
       {
-        rdir[i] = rp1[i].vectorTo (rp2[i]);
+        rdir[i].set (rp1[i].vectorTo (rp2[i]));
         bsc1.set ((rp1[i].x () + rp2[i].x ()) / 2,
                   (rp1[i].y () + rp2[i].y ()) / 2);
         for (int si = 0; (! nok) && si < i; si ++)
@@ -1037,15 +1037,15 @@ void BSRandomTester::generateImage ()
           score1 = rdir[si].squaredScalarProduct (rdir[i])
                           / (double) (rdir[si].norm2 () * rdir[i].norm2 ());
           if (rp1[si].chessboard (bsc1) < sminlength / 2)
-            ali = bsc1.vectorTo (rp2[si]);
-          else ali = rp1[si].vectorTo (bsc1);
+            ali.set (bsc1.vectorTo (rp2[si]));
+          else ali.set (rp1[si].vectorTo (bsc1));
           score2 = rdir[si].squaredScalarProduct (ali)
                    / (double) (rdir[si].norm2 () * ali.norm2 ());
           bsc2.set ((rp1[si].x () + rp2[si].x ()) / 2,
                     (rp1[si].y () + rp2[si].y ()) / 2);
           if (rp1[i].chessboard (bsc2) < sminlength / 2)
-            ali = bsc2.vectorTo (rp2[i]);
-          else ali = rp1[i].vectorTo (bsc2);
+            ali.set (bsc2.vectorTo (rp2[i]));
+          else ali.set (rp1[i].vectorTo (bsc2));
           score3 = rdir[i].squaredScalarProduct (ali)
                    / (double) (rdir[i].norm2 () * ali.norm2 ());
           if (score1 > 0.9 && (score2 > 0.9 || score3 > 0.9)) nok = true;
@@ -1094,11 +1094,11 @@ bool BSRandomTester::reloadImage ()
   {
     input >> x;
     input >> y;
-    rp1[i] = Pt2i (x, y);
+    rp1[i].set (x, y);
     input >> x;
     input >> y;
-    rp2[i] = Pt2i (x, y);
-    rdir[i] = rp1[i].vectorTo (rp2[i]);
+    rp2[i].set (x, y);
+    rdir[i].set (rp1[i].vectorTo (rp2[i]));
     input >> rw[i];
   }
 
diff --git a/Code/FBSD/BSTools/bsstructureitem.cpp b/Code/FBSD/BSTools/bsstructureitem.cpp
index 995c5d7e511b95f3b673ca47f52ba0b63037fd51..5ad949acf1b75bf67f7320a765e2e18ccc00713c 100755
--- a/Code/FBSD/BSTools/bsstructureitem.cpp
+++ b/Code/FBSD/BSTools/bsstructureitem.cpp
@@ -147,18 +147,8 @@ void BSStructureItem::paintBlurredSegment (QPainter *painter, int step)
       addText (painter, QString ("H : pixel lack tolerence = ")
                + QString::number (det->getPixelLackTolerence ()));
       if (step == 0)
-      {
         addText (painter, QString ("D : extension limit = ")
                  + QString::number (det->initialDetectionMaxExtent ()));
-        addText (painter, QString ("P : Prefiltering ")
-                 + (det->isFiltering (step) ?
-                    (QString ("on : ")
-                     + QString::number (det->prefilteringInputSize ())
-                     + QString (" -> ")
-                     + QString::number (det->prefilteringOutputSize ())
-                     + QString (" pixels")) :
-                    QString ("off")));
-      }
     }
   }
 }
@@ -286,21 +276,8 @@ void BSStructureItem::paintScansAndFilter (QPainter *painter, int step)
       paintPixels (painter, stroke, Qt::red);
     }
 
-    // Displays filter output
-    if (det->isFiltering (step))
-    {
-      paintPixels (painter, det->getAccepted (step), Qt::green);
-      paintPixels (painter, det->getRejected (step), Qt::red);
-      BlurredSegment *bs = det->getBlurredSegment (step);
-      if (bs != NULL) paintPixels (painter, bs->getStartPt (), Qt::yellow);
-    }
-
-    // Displays the blurred segment
-    else
-    {
-      BlurredSegment *bs = det->getBlurredSegment (step);
-      if (bs != NULL) paintPixels (painter, bs->getAllPoints (), Qt::green);
-    }
+    BlurredSegment *bs = det->getBlurredSegment (step);
+    if (bs != NULL) paintPixels (painter, bs->getAllPoints (), Qt::green);
   }
 
   if (verbose)
diff --git a/Code/FBSD/BlurredSegment/bsdetector.cpp b/Code/FBSD/BlurredSegment/bsdetector.cpp
index f1cfa36a2fbdeb7e0a5507aeb34b9e623eb050e0..e594d4a3bcee88ede33a72c9f3c8005893956e1b 100755
--- a/Code/FBSD/BlurredSegment/bsdetector.cpp
+++ b/Code/FBSD/BlurredSegment/bsdetector.cpp
@@ -1,5 +1,4 @@
 #include "bsdetector.h"
-//#include "linespacefilter.h"
 
 
 const std::string BSDetector::VERSION = "0.2.2";
@@ -57,12 +56,9 @@ BSDetector::BSDetector ()
       bstStatic->toggleAssignedThicknessControl ();
   }
 
-  prefilteringOn = false;
-  //lsf1 = (prefilteringOn ? new LineSpaceFilter () : NULL);
-  lsf1 = (prefilteringOn ? new BSFilter () : NULL);
-  filteringOn = false;
-  //lsf2 = (filteringOn ? new LineSpaceFilter () : NULL);
-  lsf2 = (filteringOn ? new BSFilter () : NULL);
+  nfaOn = false;
+  nfaf = NULL;
+  nfaf = (nfaOn ? new NFAFilter () : NULL);
 
   oppositeGradientDir = false;   // main edge detection
   initialMinSize = DEFAULT_INITIAL_MIN_SIZE;
@@ -73,7 +69,6 @@ BSDetector::BSDetector ()
   finalSizeTestOn = true;
   finalMinSize = DEFAULT_FINAL_MIN_SIZE;
   finalSparsityTestOn = false;
-  // nbSmallBS = 0;
   multiSelection = false;
   autodet = false;
   autoSweepingStep = DEFAULT_AUTO_SWEEPING_STEP;
@@ -94,8 +89,6 @@ BSDetector::~BSDetector ()
   delete bst1;
   delete bst2;
   
-  if (lsf1 != NULL) delete lsf1;
-  if (lsf2 != NULL) delete lsf2;
   if (bsini != NULL) delete bsini;
   if (bsf != NULL) delete bsf;
   vector <BlurredSegment *>::iterator it = mbsf.begin ();
@@ -105,16 +98,6 @@ BSDetector::~BSDetector ()
 
 void BSDetector::clearAll ()
 {
-  if (lsf1 != NULL)
-  {
-    delete lsf1;
-    lsf1 = NULL;
-  }
-  if (lsf2 != NULL)
-  {
-    delete lsf2;
-    lsf2 = NULL;
-  }
   if (bsini != NULL)
   {
     delete bsini;
@@ -128,6 +111,8 @@ void BSDetector::clearAll ()
   vector <BlurredSegment *>::iterator it = mbsf.begin ();
   while (it != mbsf.end ()) delete (*it++);
   mbsf.clear ();
+  vbsf.clear ();
+  rbsf.clear ();
 }
 
 
@@ -139,6 +124,7 @@ void BSDetector::setGradientMap (VMap *data)
   if (bst1) bst1->setGradientMap (data);
   if (bst2) bst2->setGradientMap (data);
   if (bstStatic) bstStatic->setGradientMap (data);
+  if (nfaf) nfaf->init (data);
 }
 
 
@@ -151,7 +137,6 @@ void BSDetector::detectAll ()
 
   bool isnext = true;
   nbtrials = 0;
-  // nbSmallBS = 0;
   int width = gMap->getWidth ();
   int height = gMap->getHeight ();
   for (int x = width / 2; isnext && x > 0; x -= autoSweepingStep)
@@ -165,7 +150,8 @@ void BSDetector::detectAll ()
        isnext && y < height - 1; y += autoSweepingStep)
     isnext = runMultiDetection (Pt2i (0, y), Pt2i (width - 1, y));
   if (maxtrials > (int) (mbsf.size ())) maxtrials = 0;
-  // cout << nbSmallBS << " petits BS elimines" << endl;
+
+  if (nfaf) nfaf->filter (mbsf, vbsf, rbsf);
 
   gMap->setMasking (false);
 }
@@ -180,7 +166,6 @@ void BSDetector::detectAllWithBalancedXY ()
 
   bool isnext = true;
   nbtrials = 0;
-  // nbSmallBS = 0;
   int width = gMap->getWidth ();
   int height = gMap->getHeight ();
   int xg = width / 2, yb = height / 2;
@@ -214,7 +199,6 @@ void BSDetector::detectAllWithBalancedXY ()
     }
   }
   if (maxtrials > (int) (mbsf.size ())) maxtrials = 0;
-  // cout << nbSmallBS << " petits BS elimines" << endl;
   gMap->setMasking (false);
 }
 
@@ -228,9 +212,7 @@ void BSDetector::detectSelection (const Pt2i &p1, const Pt2i &p2)
     gMap->setMasking (true);
     gMap->clearMask ();
     nbtrials = 0;
-    // nbSmallBS = 0;
     runMultiDetection (p1, p2);
-    // cout << nbSmallBS << " petits BS elimines" << endl;
     if (maxtrials > (int) (mbsf.size ())) maxtrials = 0;
     gMap->setMasking (false);
   }
@@ -309,8 +291,6 @@ int BSDetector::detect (const Pt2i &p1, const Pt2i &p2,
   if (prelimDetectionOn) bst0->clear ();
   bst1->clear ();
   bst2->clear ();
-  if (prefilteringOn) lsf1->clear ();
-  if (filteringOn) lsf2->clear ();
   if (bspre != NULL) delete bspre;
   bspre = NULL;
   if (bsini != NULL) delete bsini;
@@ -376,20 +356,6 @@ int BSDetector::detect (const Pt2i &p1, const Pt2i &p2,
       return RESULT_INITIAL_TOO_SPARSE;
   }
 
-  // Filtering the initial segment
-  //------------------------------
-  if (prefilteringOn)
-  {
-    BlurredSegment *fbs = lsf1->filter (bsini);
-    if (fbs != NULL)
-    {
-      delete bsini;
-      bsini = fbs;
-    }
-    if (bsini->size () < initialMinSize)
-      return RESULT_INITIAL_TOO_MANY_OUTLIERS;
-  }
-
   // Orientation test for automatic extractions
   //-------------------------------------------
   Vr2i bsinidir = bsini->getSupportVector();
@@ -420,10 +386,7 @@ int BSDetector::detect (const Pt2i &p1, const Pt2i &p2,
   {
     // DigitalStraightSegment *dss = bsf->getSegment ();
     if ((int) (bsf->getAllPoints().size ()) < finalMinSize)
-    {
-      // nbSmallBS ++;
       return RESULT_FINAL_TOO_SMALL;
-    }
   }
 
   // Sparsity test
@@ -449,20 +412,6 @@ int BSDetector::detect (const Pt2i &p1, const Pt2i &p2,
     if (bsccp < bssize / 2) return RESULT_FINAL_TOO_FRAGMENTED;
   }
 
-
-  // Final filtering
-  //----------------
-  if (filteringOn)
-  {
-    BlurredSegment *fbsf = lsf2->filter (bsf);
-    if (fbsf != NULL)
-    {
-      delete bsf;
-      bsf = fbsf;
-    }
-    else return RESULT_FINAL_TOO_MANY_OUTLIERS;
-  }
-
   return RESULT_OK;
 }
 
@@ -570,10 +519,7 @@ int BSDetector::staticDetect (const Pt2i &p1, const Pt2i &p2,
   {
     // DigitalStraightSegment *dss = bsf->getSegment ();
     if ((int) (bsf->getAllPoints().size ()) < finalMinSize)
-    {
-      // nbSmallBS ++;
       return RESULT_FINAL_TOO_SMALL;
-    }
   }
 
   // Final sparsity test
@@ -731,31 +677,13 @@ void BSDetector::switchOrthoScans ()
 }
 
 
-vector<Pt2i> BSDetector::getRejected (int step) const
+void BSDetector::switchNFA ()
 {
-  vector<Pt2i> res;
-  if (step == STEP_FINAL)
-  {
-    if (filteringOn) res = lsf2->getRejected ();
-  }
-  else if (prefilteringOn) res = lsf1->getRejected ();
-  return res;
-}
-
-
-void BSDetector::switchFiltering (int step)
-{
-  if (step == STEP_FINAL)
-  {
-    filteringOn = ! filteringOn;
-    //if (filteringOn && lsf2 == NULL) lsf2 = new LineSpaceFilter ();
-    if (filteringOn && lsf2 == NULL) lsf2 = new BSFilter ();
-  }
-  else if (step == STEP_INITIAL)
+  nfaOn = ! nfaOn;
+  if (nfaOn && nfaf == NULL)
   {
-    prefilteringOn = ! prefilteringOn;
-    //if (prefilteringOn && lsf1 == NULL) lsf1 = new LineSpaceFilter ();
-    if (prefilteringOn && lsf1 == NULL) lsf1 = new BSFilter ();
+    nfaf = new NFAFilter ();
+    nfaf->init (gMap);
   }
 }
 
diff --git a/Code/FBSD/BlurredSegment/bsdetector.h b/Code/FBSD/BlurredSegment/bsdetector.h
index 0a53f3e8a060d5426fec4f63c9f9af3490dfb24a..c83967ea424f2c7d3a0f75ebe703b6238a5bc004 100755
--- a/Code/FBSD/BlurredSegment/bsdetector.h
+++ b/Code/FBSD/BlurredSegment/bsdetector.h
@@ -2,7 +2,7 @@
 #define BLURRED_SEGMENT_DETECTOR_H
 
 #include "bstracker.h"
-#include "bsfilter.h"
+#include "nfafilter.h"
 #include <iostream>
 
 using namespace std;
@@ -11,7 +11,6 @@ using namespace std;
 /** 
  * @class BSDetector bsdetector.h
  * \brief Blurred segment detector in grey level images.
- * \author {P. Even}
  */
 class BSDetector
 {
@@ -162,6 +161,18 @@ public:
   inline const vector<BlurredSegment *> getBlurredSegments () const {
     return (mbsf); }
 
+  /**
+   * \brief Returns the list of NFA-validated blurred segments at final step.
+   */
+  inline const vector<BlurredSegment *> getValidSegments () const {
+    return (vbsf); }
+
+  /**
+   * \brief Returns the list of NFA-rejected blurred segments at final step.
+   */
+  inline const vector<BlurredSegment *> getRejectedSegments () const {
+    return (rbsf); }
+
   /**
    * \brief Avoids the deletion of the last extracted blurred segment.
    */
@@ -454,58 +465,27 @@ public:
     bst2->incAssignedThicknessControlDelay (val); }
 
   /**
-   * \brief Returns if one of the filter is activated.
-   */
-  inline bool isOneFilterSet () { return (filteringOn || prefilteringOn); }
-
-  /**
-   * \brief Returns a blurred segment filter.
-   * @param step Initial step addressed if set to 0, final step otherwise.
-   */
-  inline BSFilter *getFilter (int step) const {
-    return (step == STEP_FINAL ? lsf2 : lsf1); }
-
-  /**
-   * \brief Returns the accepted points by the blurred segment filter.
-   * @param step Initial step addressed if set to 0, final step otherwise.
-   */
-  inline vector<Pt2i> getAccepted (int step) const {
-    return (step == STEP_FINAL ?
-      (filteringOn ? lsf2->getAccepted () : bsf->getAllPoints ()) :
-      (prefilteringOn ? lsf1->getAccepted () : bsini->getAllPoints ())); }
-
-  /**
-   * \brief Returns the rejected points by the line space based filter.
-   * @param step Initial step addressed if set to 0, final step otherwise.
+   * \brief Returns if NFA-based filter is activated.
    */
-  vector<Pt2i> getRejected (int step) const;
+  inline bool isNFA () const { return (nfaOn); }
 
   /**
-   * \brief Returns if the given line space based filter is activated.
-   * @param step Initial step addressed if set to 0, final step otherwise.
+   * \brief Toggles the use of the NFA-based filter.
    */
-  inline bool isFiltering (int step) const {
-    if (step == STEP_FINAL) return (filteringOn);
-    else if (step == STEP_INITIAL) return (prefilteringOn);
-    else return false; }
+  void switchNFA ();
 
   /**
-   * \brief Toggles the use of the given line space based filter.
-   * @param step Initial step addressed if set to 0, final step otherwise.
+   * \brief Returns the NFA length ratio parameter value.
    */
-  void switchFiltering (int step);
+  inline double nfaLengthRatio () const {
+    return (nfaOn ? nfaf->lengthRatio () : 0.0); }
 
   /**
-   * \brief Returns the blurred segment size before pre-filtering.
+   * \brief Increments the NFA length ratio parameter value.
+   * @param inc Increment direction.
    */
-  inline int prefilteringInputSize () {
-    return (prefilteringOn ? lsf1->blurredSegmentInitialSize () : 0); }
-
-  /**
-   * \brief Returns the blurred segment size after pre-filtering.
-   */
-  inline int prefilteringOutputSize () {
-    return (prefilteringOn ? lsf1->blurredSegmentFinalSize () : 0); }
+  inline void incNfaLengthRatio (int inc) {
+    if (nfaOn) nfaf->incLengthRatio (inc); }
 
   /**
    * \brief Returns whether the density test at initial step is set.
@@ -742,16 +722,15 @@ private :
   BlurredSegment *bsf;
   /** Detected blurred segments in case of multi-detection (final results). */
   vector<BlurredSegment *> mbsf;
-
-  /** Initial segment filtering modality. */
-  bool prefilteringOn;
-  /** Blurred segment pre-filter. */
-  BSFilter *lsf1;
-
-  /** Detected segment filtering modality. */
-  bool filteringOn;
-  /** Blurred segment post-filter. */
-  BSFilter *lsf2;
+  /** NFA-validated blurred segments. */
+  vector<BlurredSegment *> vbsf;
+  /** NFA-rejected blurred segments. */
+  vector<BlurredSegment *> rbsf;
+
+  /** NFA-based filtering modality. */
+  bool nfaOn;
+  /** NFA-based filter. */
+  NFAFilter *nfaf;
 
 
   /**
diff --git a/Code/FBSD/BlurredSegment/bsfilter.cpp b/Code/FBSD/BlurredSegment/bsfilter.cpp
deleted file mode 100755
index 94695c86e6dbd659018ebaa579e6f7a3a738e850..0000000000000000000000000000000000000000
--- a/Code/FBSD/BlurredSegment/bsfilter.cpp
+++ /dev/null
@@ -1,56 +0,0 @@
-#include "bsfilter.h"
-#include "bsproto.h"
-
-
-
-BSFilter::BSFilter ()
-{
-  bsInitialSize = 0;
-  bsFinalSize = 0;
-}
-
-
-BSFilter::~BSFilter ()
-{
-}
-
-
-void BSFilter::clear ()
-{
-  leftIn.clear ();
-  rightIn.clear ();
-  out.clear ();
-}
-
-
-BlurredSegment *BSFilter::filter (BlurredSegment *bs)
-{
-  leftIn.clear ();
-  rightIn.clear ();
-  vector<Pt2i> *pts = bs->getAllLeft ();
-  vector<Pt2i>::iterator it = pts->end ();
-  while (it-- != pts->begin ()) leftIn.push_back (*it);
-  pts = bs->getAllRight ();
-  it = pts->begin ();
-  while (it != pts->end ()) rightIn.push_back (*it++);
-
-  int mw = 1;
-  AbsRat sw = bs->minimalWidth ();
-  if (sw.denominator () != 0) mw += sw.numerator () / sw.denominator ();
-  BSProto resbs (mw, bs->getCenter (), leftIn, rightIn);
-  bsFinalSize = resbs.size ();
-  return (resbs.endOfBirth ());
-}
-
-
-vector<Pt2i> BSFilter::getAccepted ()
-{
-  vector<Pt2i> res;
-  for (vector<Pt2i>::iterator it = leftIn.begin ();
-       it != leftIn.end (); it ++)
-    res.push_back (*it);
-  for (vector<Pt2i>::iterator it = rightIn.begin ();
-       it != rightIn.end (); it ++)
-    res.push_back (*it);
-  return res;
-}
diff --git a/Code/FBSD/BlurredSegment/bsfilter.h b/Code/FBSD/BlurredSegment/bsfilter.h
deleted file mode 100755
index e7bee9c01a476007674846491ddecc0afb88ed31..0000000000000000000000000000000000000000
--- a/Code/FBSD/BlurredSegment/bsfilter.h
+++ /dev/null
@@ -1,74 +0,0 @@
-#ifndef BLURRED_SEGMENT_FILTER_H
-#define BLURRED_SEGMENT_FILTER_H
-
-#include <vector>
-#include "pt2i.h"
-#include "blurredsegment.h"
-
-using namespace std;
-
-
-/** 
- * @class Blurred segment filter.
- * \author {P. Even}
- */
-class BSFilter
-{
-public:
-
-  /**
-   * \brief Creates a empty blurred segment filter.
-   */
-  BSFilter ();
-
-  /**
-   * \brief Deletes the filter.
-   */
-  virtual ~BSFilter ();
-
-  /**
-   * \brief Clears the stored information.
-   */
-  virtual void clear ();
-
-  /** \brief Filters and returns a blurred segment.
-    * @param bs The blurred segment to be processed.
-    */
-  virtual BlurredSegment *filter (BlurredSegment *bs);
-
-  /**
-   * \brief Returns the accepted points by the line space filter.
-   */
-  vector<Pt2i> getAccepted ();
-
-  /**
-   * \brief Returns the rejected points by the line space filter.
-   */
-  inline vector<Pt2i> getRejected () const { return out; }
-
-  /**
-   * \brief Returns the blurred segment size before filtering.
-   */
-  inline int blurredSegmentInitialSize () const { return bsInitialSize; }
-
-  /**
-   * \brief Returns the blurred segment size after filtering.
-   */
-  inline int blurredSegmentFinalSize () const { return bsFinalSize; }
-
-
-protected :
-
-  /** Blurred segment size before filtering. */
-  int bsInitialSize;
-  /** Blurred segment size after filtering. */
-  int bsFinalSize;
-  /** Accepted left points after line space filtering. */
-  vector<Pt2i> leftIn;
-  /** Accepted right points after line space filtering. */
-  vector<Pt2i> rightIn;
-  /** Rejected points after line space filtering. */
-  vector<Pt2i> out;
-
-};
-#endif
diff --git a/Code/FBSD/BlurredSegment/nfafilter.cpp b/Code/FBSD/BlurredSegment/nfafilter.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..14039e77d677ef855d0fe715eb4de2c222f935a0
--- /dev/null
+++ b/Code/FBSD/BlurredSegment/nfafilter.cpp
@@ -0,0 +1,138 @@
+#include <cmath>
+#include <iostream>
+#include "nfafilter.h"
+
+
+const double NFAFilter::NFA_EPSILON = 1.0;
+const double NFAFilter::DEFAULT_LRATIO = 1.0;
+const int NFAFilter::DEFAULT_MIN_SECTION_LENGTH = 3;
+
+NFAFilter::NFAFilter ()
+{
+  min_section_length = DEFAULT_MIN_SECTION_LENGTH;
+  max_grad2 = 0;
+  gradient_map = NULL;
+  cum_histo = NULL;
+  bs_section_count = 0;
+  lratio = DEFAULT_LRATIO;
+}
+
+
+NFAFilter::~NFAFilter ()
+{
+  delete [] cum_histo;
+}
+
+
+void NFAFilter::init (VMap *gmap)
+{
+  // Stores the reference to the gradient map
+  gradient_map = gmap;
+  int width = gradient_map->getWidth ();
+  int height = gradient_map->getHeight ();
+
+  // Counts the number of non-gradient pixels
+  int m = 0;
+  max_grad2 = 0;
+  for (int j = 0; j < height; j++)
+  {
+    for (int i = 0; i < width; i++)
+    {
+      int gradval = gradient_map->sqNorm (i, j);
+      if (gradval != 0)
+      {
+        m++;
+        if (gradval > max_grad2) max_grad2 = gradval;
+      }
+    }
+  }
+
+  // Rectification 
+  m = (width -2) * (height - 2);
+
+  // Gets gradient histogram
+  int gmax = (int) (sqrt (max_grad2));
+  cum_histo = new double[gmax + 1];
+  for (int j = 0; j < height; j++)
+    for (int i = 0; i < width; i++)
+      cum_histo[(int) (sqrt (gradient_map->sqNorm (i, j)))] ++;
+
+  // Gets cumulated histogram
+  for (int i = gmax; i > 0; i--)
+    cum_histo[i-1] += cum_histo[i];
+
+  // Normalization
+  for (int i = 0; i <= gmax; i++)
+    cum_histo[i] /= m;
+}
+
+
+double NFAFilter::nfaValue (double proba, int length)
+{
+  length = (int) (length / lratio); // Magic number : divForTestSeg
+  double nfa = bs_section_count;
+  for (int i = 0; i < length && nfa > NFA_EPSILON; i++) nfa *= proba;
+  return nfa;
+}
+
+
+bool NFAFilter::filter (const BlurredSegment *bs, int start, int end)
+{
+  int length = end - start;
+  if (length < min_section_length) return false;
+
+  // Gets point with small gradient
+  int gmin = max_grad2;
+  int pmin = -1;
+  vector<Pt2i> pts = bs->getAllPoints ();
+  for (int i = start; i < end; i++)
+  {
+    int gn = (gradient_map->getValue (pts[i])).norm2 ();
+    if (gn < gmin)
+    {
+      gmin = gn;
+      pmin = i;
+    }
+  }
+
+  // Gets NFA and accept or split the segment
+  double nfa = nfaValue (cum_histo[(int) (sqrt (gmin))], length);
+  if (nfa < NFA_EPSILON) return true;
+  return (filter (bs, start, pmin) && filter (bs, pmin + 1, end));
+}
+
+
+void NFAFilter::filter (const vector<BlurredSegment *> &bss,
+                        vector<BlurredSegment *> &vsegs,
+                        vector<BlurredSegment *> &rsegs)
+{
+  vsegs.clear ();
+  rsegs.clear ();
+
+  // Computes Np
+  bs_section_count = 0;
+  vector<BlurredSegment *>::const_iterator it = bss.begin ();
+  while (it != bss.end ())
+  {
+    int length = (*it)->size ();
+    bs_section_count += length * (length - 1) / 2;
+    it ++;
+  }
+
+  // Computes and test each segment NFA
+  it = bss.begin ();
+  while (it != bss.end ())
+  {
+    if (filter (*it, 0, (*it)->size ())) vsegs.push_back (*it);
+    else rsegs.push_back (*it);
+    it ++;
+  }
+}
+
+
+void NFAFilter::incLengthRatio (int inc)
+{
+  lratio += inc * 0.05;
+  if (lratio < 1.0) lratio = 1.0;
+  else if (lratio > 3.0) lratio = 3.0;
+}
diff --git a/Code/FBSD/BlurredSegment/nfafilter.h b/Code/FBSD/BlurredSegment/nfafilter.h
new file mode 100755
index 0000000000000000000000000000000000000000..d112f5f0d96af0f0d35a3805234d852fea666bcc
--- /dev/null
+++ b/Code/FBSD/BlurredSegment/nfafilter.h
@@ -0,0 +1,83 @@
+#ifndef NFA_FILTER_H
+#define NFA_FILTER_H
+
+#include "blurredsegment.h"
+#include "vmap.h"
+
+
+/** 
+ * @class Number of false alarm based filter for selecting blurred segments.
+ */
+class NFAFilter
+{
+public:
+
+  /**
+   * \brief Creates an empty NFA-based filter.
+   */
+  NFAFilter ();
+
+  /**
+   * \brief Deletes the filter.
+   */
+  ~NFAFilter ();
+
+  /**
+   * \brief Initializes the filter before any detection.
+   * @param gmap Reference to used gradient map.
+   */
+  void init (VMap *gmap);
+
+  /**
+    * \brief Filters a set of blurred segments.
+    * @param bss Input set of blurred segments.
+    * @param vbss Output set of valid blurred segments.
+    * @param rbss Output set of rejected blurred segments.
+    */
+  void filter (const vector<BlurredSegment *> &bss,
+               vector<BlurredSegment *> &vbss,
+               vector<BlurredSegment *> &rbss);
+
+  inline double lengthRatio () const { return lratio; }
+  void incLengthRatio (int inc);
+
+private :
+
+  /** Tolered number of false alarm expectation. */
+  static const double NFA_EPSILON;
+  /** Default length ratio for NFA measure. */
+  static const double DEFAULT_LRATIO;
+  /** Default value for the smallest blurred segment section considered. */
+  static const int DEFAULT_MIN_SECTION_LENGTH;
+
+  /** Smallest blurred segment section considered. */
+  int min_section_length;
+  /** Maximal squarred gradient value. */
+  int max_grad2;
+  /** Reference to used gradient map. */
+  VMap *gradient_map;
+  /** Cumulated gradient histogramm (H). */
+  double *cum_histo;
+  /** Count of any blurred segment sections (Np). */
+  int bs_section_count;
+
+  double lratio;
+
+
+  /** 
+    * \brief Computes number of false alarms of a segment section.
+    * @param proba Validation probability associated to min gradient point.
+    * @param length Section length.
+    */
+  double nfaValue (double proba, int length);
+
+  /**
+    * \brief Filters a blurred segment section.
+    * @param bss Reference to input blurred segment.
+    * @param start Index of start point in the section.
+    * @param end Index of first point out of the section.
+    */
+  bool filter (const BlurredSegment *bs, int start, int end);
+
+};
+#endif
diff --git a/Code/FBSD/FBSD.pro b/Code/FBSD/FBSD.pro
index 2a03ddf33c75c4bf9c2048c733aa164c52281aca..55c0a3f953f0bb31dc07d42d85c2a13a5d6e6c37 100644
--- a/Code/FBSD/FBSD.pro
+++ b/Code/FBSD/FBSD.pro
@@ -17,9 +17,9 @@ OBJECTS_DIR = obj
 HEADERS += BlurredSegment/biptlist.h \
            BlurredSegment/blurredsegmentproto.h \
            BlurredSegment/bsdetector.h \
-           BlurredSegment/bsfilter.h \
            BlurredSegment/bsproto.h \
            BlurredSegment/bstracker.h \
+           BlurredSegment/nfafilter.h \
            BSTools/bsdetectionwidget.h \
            BSTools/bsgroundtruthitem.h \
            BSTools/bsgroundtruthview.h \
@@ -60,9 +60,9 @@ SOURCES += main.cpp \
            BlurredSegment/biptlist.cpp \
            BlurredSegment/blurredsegment.cpp \
            BlurredSegment/bsdetector.cpp \
-           BlurredSegment/bsfilter.cpp \
            BlurredSegment/bsproto.cpp \
            BlurredSegment/bstracker.cpp \
+           BlurredSegment/nfafilter.cpp \
            BSTools/bsdetectionwidget.cpp \
            BSTools/bsgroundtruthitem.cpp \
            BSTools/bsgroundtruthview.cpp \
diff --git a/Code/FBSD/ImageTools/pt2i.h b/Code/FBSD/ImageTools/pt2i.h
index c08d4828c309e1bb5f91d9502d68fdcf4cb1a13f..066e88d991a3014122df84a66482f52d0d42b16c 100755
--- a/Code/FBSD/ImageTools/pt2i.h
+++ b/Code/FBSD/ImageTools/pt2i.h
@@ -87,7 +87,7 @@ public:
   /**
    * @fn void set (const Pt2i &p)
    * \brief Sets the point coordinates.
-   * @param point to recopy.
+   * @param p Point to copy.
    */
   inline void set (const Pt2i &p)
   {
diff --git a/Code/FBSD/ImageTools/vmap.cpp b/Code/FBSD/ImageTools/vmap.cpp
index 42cb4a1155dad6c5b73d1101784d5cb2a7fd4bd2..db238265ee72bcc9142c72c73e525604c9d4fa57 100755
--- a/Code/FBSD/ImageTools/vmap.cpp
+++ b/Code/FBSD/ImageTools/vmap.cpp
@@ -3,7 +3,7 @@
 #include "math.h"
 #include <inttypes.h>
 
-using namespace std;
+// using namespace std;
 
 
 const int VMap::TYPE_UNKNOWN = -1;
@@ -394,18 +394,18 @@ void VMap::buildGradientMap (unsigned char *data)
     gm++;
     for (int j = 1; j < width - 1; j++)
     {
-      gm->set (data[(i - 1) * height + j + 1]
-               + 2 * data[i * height + j + 1]
-               + data[(i + 1) * height + j + 1]
-               - data[(i - 1) * height + j - 1]
-               - 2 * data[i * height + j - 1]
-               - data[(i + 1) * height + j - 1],
-               data[(i + 1) * height + j - 1]
-               + 2 * data[(i + 1) * height + j]
-               + data[(i + 1) * height + j + 1]
-               - data[(i - 1) * height + j - 1]
-               - 2 * data[(i - 1) * height + j]
-               - data[(i - 1) * height + j + 1]);
+      gm->set (data[(i - 1) * width + j + 1]
+               + 2 * data[i * width + j + 1]
+               + data[(i + 1) * width + j + 1]
+               - data[(i - 1) * width + j - 1]
+               - 2 * data[i * width + j - 1]
+               - data[(i + 1) * width + j - 1],
+               data[(i + 1) * width + j - 1]
+               + 2 * data[(i + 1) * width + j]
+               + data[(i + 1) * width + j + 1]
+               - data[(i - 1) * width + j - 1]
+               - 2 * data[(i - 1) * width + j]
+               - data[(i - 1) * width + j + 1]);
       gm++;
     }
     gm->set (0, 0);
@@ -435,18 +435,18 @@ void VMap::buildGradientMap (int *data)
     gm++;
     for (int j = 1; j < width - 1; j++)
     {
-      gm->set (data[(i - 1) * height + j + 1]
-               + 2 * data[i * height + j + 1]
-               + data[(i + 1) * height + j + 1]
-               - data[(i - 1) * height + j - 1]
-               - 2 * data[i * height + j - 1]
-               - data[(i + 1) * height + j - 1],
-               data[(i + 1) * height + j - 1]
-               + 2 * data[(i + 1) * height + j]
-               + data[(i + 1) * height + j + 1]
-               - data[(i - 1) * height + j - 1]
-               - 2 * data[(i - 1) * height + j]
-               - data[(i - 1) * height + j + 1]);
+      gm->set (data[(i - 1) * width + j + 1]
+               + 2 * data[i * width + j + 1]
+               + data[(i + 1) * width + j + 1]
+               - data[(i - 1) * width + j - 1]
+               - 2 * data[i * width + j - 1]
+               - data[(i + 1) * width + j - 1],
+               data[(i + 1) * width + j - 1]
+               + 2 * data[(i + 1) * width + j]
+               + data[(i + 1) * width + j + 1]
+               - data[(i - 1) * width + j - 1]
+               - 2 * data[(i - 1) * width + j]
+               - data[(i - 1) * width + j + 1]);
       gm++;
     }
     gm->set (0, 0);
@@ -511,46 +511,46 @@ void VMap::buildSobel5x5Map (unsigned char *data)
     gm++;
     for (int j = 2; j < width - 2; j++)
     {
-      gm->set (5 * data[(i - 2) * height + j + 2]
-                 + 8 * data[(i - 1) * height + j + 2]
-                 + 10 * data[i * height + j + 2]
-                 + 8 * data[(i + 1) * height + j + 2]
-                 + 5 * data[(i + 2) * height + j + 2]
-               + 4 * data[(i - 2) * height + j + 1]
-                 + 10 * data[(i - 1) * height + j + 1]
-                 + 20 * data[i * height + j + 1]
-                 + 10 * data[(i + 1) * height + j + 1]
-                 + 4 * data[(i + 2) * height + j + 1]
-               - 4 * data[(i - 2) * height + j - 1]
-                 - 10 * data[(i - 1) * height + j - 1]
-                 - 20 * data[i * height + j - 1]
-                 - 10 * data[(i + 1) * height + j - 1]
-                 - 4 * data[(i + 2) * height + j - 1] 
-               - 5 * data[(i - 2) * height + j - 2]
-                 - 8 * data[(i - 1) * height + j - 2]
-                 - 10 * data[i * height + j - 2]
-                 - 8 * data[(i + 1) * height + j - 2]
-                 - 5 * data[(i + 2) * height + j - 2],
-               5 * data[(i + 2) * height + j - 2]
-                 + 8 * data[(i + 2) * height + j - 1]
-                 + 10 * data[(i + 2) * height + j]
-                 + 8 * data[(i + 2) * height + j + 1]
-                 + 5 * data[(i + 2) * height + j + 2]
-               + 4 * data[(i + 1) * height + j - 2]
-                 + 10 * data[(i + 1) * height + j - 1]
-                 + 20 * data[(i + 1) * height + j]
-                 + 10 * data[(i + 1) * height + j + 1]
-                 + 4 * data[(i + 1) * height + j + 2]
-               - 4 * data[(i - 1) * height + j - 2]
-                 - 10 * data[(i - 1) * height + j - 1]
-                 - 20 * data[(i - 1) * height + j]
-                 - 10 * data[(i - 1) * height + j + 1]
-                 - 4 * data[(i - 1) * height + j + 2]
-               - 5 * data[(i - 2) * height + j - 2]
-                 - 8 * data[(i - 2) * height + j - 1]
-                 - 10 * data[(i - 2) * height + j]
-                 - 8 * data[(i - 2) * height + j + 1]
-                 - 5 * data[(i - 2) * height + j + 2]);
+      gm->set (5 * data[(i - 2) * width + j + 2]
+                 + 8 * data[(i - 1) * width + j + 2]
+                 + 10 * data[i * width + j + 2]
+                 + 8 * data[(i + 1) * width + j + 2]
+                 + 5 * data[(i + 2) * width + j + 2]
+               + 4 * data[(i - 2) * width + j + 1]
+                 + 10 * data[(i - 1) * width + j + 1]
+                 + 20 * data[i * width + j + 1]
+                 + 10 * data[(i + 1) * width + j + 1]
+                 + 4 * data[(i + 2) * width + j + 1]
+               - 4 * data[(i - 2) * width + j - 1]
+                 - 10 * data[(i - 1) * width + j - 1]
+                 - 20 * data[i * width + j - 1]
+                 - 10 * data[(i + 1) * width + j - 1]
+                 - 4 * data[(i + 2) * width + j - 1] 
+               - 5 * data[(i - 2) * width + j - 2]
+                 - 8 * data[(i - 1) * width + j - 2]
+                 - 10 * data[i * width + j - 2]
+                 - 8 * data[(i + 1) * width + j - 2]
+                 - 5 * data[(i + 2) * width + j - 2],
+               5 * data[(i + 2) * width + j - 2]
+                 + 8 * data[(i + 2) * width + j - 1]
+                 + 10 * data[(i + 2) * width + j]
+                 + 8 * data[(i + 2) * width + j + 1]
+                 + 5 * data[(i + 2) * width + j + 2]
+               + 4 * data[(i + 1) * width + j - 2]
+                 + 10 * data[(i + 1) * width + j - 1]
+                 + 20 * data[(i + 1) * width + j]
+                 + 10 * data[(i + 1) * width + j + 1]
+                 + 4 * data[(i + 1) * width + j + 2]
+               - 4 * data[(i - 1) * width + j - 2]
+                 - 10 * data[(i - 1) * width + j - 1]
+                 - 20 * data[(i - 1) * width + j]
+                 - 10 * data[(i - 1) * width + j + 1]
+                 - 4 * data[(i - 1) * width + j + 2]
+               - 5 * data[(i - 2) * width + j - 2]
+                 - 8 * data[(i - 2) * width + j - 1]
+                 - 10 * data[(i - 2) * width + j]
+                 - 8 * data[(i - 2) * width + j + 1]
+                 - 5 * data[(i - 2) * width + j + 2]);
       gm++;
     }
     gm->set (0, 0);
@@ -584,46 +584,46 @@ void VMap::buildSobel5x5Map (int *data)
     gm++;
     for (int j = 2; j < width - 2; j++)
     {
-      gm->set (5 * data[(i - 2) * height + j + 2]
-                 + 8 * data[(i - 1) * height + j + 2]
-                 + 10 * data[i * height + j + 2]
-                 + 8 * data[(i + 1) * height + j + 2]
-                 + 5 * data[(i + 2) * height + j + 2]
-               + 4 * data[(i - 2) * height + j + 1]
-                 + 10 * data[(i - 1) * height + j + 1]
-                 + 20 * data[i * height + j + 1]
-                 + 10 * data[(i + 1) * height + j + 1]
-                 + 4 * data[(i + 2) * height + j + 1]
-               - 4 * data[(i - 2) * height + j - 1]
-                 - 10 * data[(i - 1) * height + j - 1]
-                 - 20 * data[i * height + j - 1]
-                 - 10 * data[(i + 1) * height + j - 1]
-                 - 4 * data[(i + 2) * height + j - 1] 
-               - 5 * data[(i - 2) * height + j - 2]
-                 - 8 * data[(i - 1) * height + j - 2]
-                 - 10 * data[i * height + j - 2]
-                 - 8 * data[(i + 1) * height + j - 2]
-                 - 5 * data[(i + 2) * height + j - 2],
-               5 * data[(i + 2) * height + j - 2]
-                 + 8 * data[(i + 2) * height + j - 1]
-                 + 10 * data[(i + 2) * height + j]
-                 + 8 * data[(i + 2) * height + j + 1]
-                 + 5 * data[(i + 2) * height + j + 2]
-               + 4 * data[(i + 1) * height + j - 2]
-                 + 10 * data[(i + 1) * height + j - 1]
-                 + 20 * data[(i + 1) * height + j]
-                 + 10 * data[(i + 1) * height + j + 1]
-                 + 4 * data[(i + 1) * height + j + 2]
-               - 4 * data[(i - 1) * height + j - 2]
-                 - 10 * data[(i - 1) * height + j - 1]
-                 - 20 * data[(i - 1) * height + j]
-                 - 10 * data[(i - 1) * height + j + 1]
-                 - 4 * data[(i - 1) * height + j + 2]
-               - 5 * data[(i - 2) * height + j - 2]
-                 - 8 * data[(i - 2) * height + j - 1]
-                 - 10 * data[(i - 2) * height + j]
-                 - 8 * data[(i - 2) * height + j + 1]
-                 - 5 * data[(i - 2) * height + j + 2]);
+      gm->set (5 * data[(i - 2) * width + j + 2]
+                 + 8 * data[(i - 1) * width + j + 2]
+                 + 10 * data[i * width + j + 2]
+                 + 8 * data[(i + 1) * width + j + 2]
+                 + 5 * data[(i + 2) * width + j + 2]
+               + 4 * data[(i - 2) * width + j + 1]
+                 + 10 * data[(i - 1) * width + j + 1]
+                 + 20 * data[i * width + j + 1]
+                 + 10 * data[(i + 1) * width + j + 1]
+                 + 4 * data[(i + 2) * width + j + 1]
+               - 4 * data[(i - 2) * width + j - 1]
+                 - 10 * data[(i - 1) * width + j - 1]
+                 - 20 * data[i * width + j - 1]
+                 - 10 * data[(i + 1) * width + j - 1]
+                 - 4 * data[(i + 2) * width + j - 1] 
+               - 5 * data[(i - 2) * width + j - 2]
+                 - 8 * data[(i - 1) * width + j - 2]
+                 - 10 * data[i * width + j - 2]
+                 - 8 * data[(i + 1) * width + j - 2]
+                 - 5 * data[(i + 2) * width + j - 2],
+               5 * data[(i + 2) * width + j - 2]
+                 + 8 * data[(i + 2) * width + j - 1]
+                 + 10 * data[(i + 2) * width + j]
+                 + 8 * data[(i + 2) * width + j + 1]
+                 + 5 * data[(i + 2) * width + j + 2]
+               + 4 * data[(i + 1) * width + j - 2]
+                 + 10 * data[(i + 1) * width + j - 1]
+                 + 20 * data[(i + 1) * width + j]
+                 + 10 * data[(i + 1) * width + j + 1]
+                 + 4 * data[(i + 1) * width + j + 2]
+               - 4 * data[(i - 1) * width + j - 2]
+                 - 10 * data[(i - 1) * width + j - 1]
+                 - 20 * data[(i - 1) * width + j]
+                 - 10 * data[(i - 1) * width + j + 1]
+                 - 4 * data[(i - 1) * width + j + 2]
+               - 5 * data[(i - 2) * width + j - 2]
+                 - 8 * data[(i - 2) * width + j - 1]
+                 - 10 * data[(i - 2) * width + j]
+                 - 8 * data[(i - 2) * width + j + 1]
+                 - 5 * data[(i - 2) * width + j + 2]);
       gm++;
     }
     gm->set (0, 0);
diff --git a/Code/FBSD/ImageTools/vr2i.h b/Code/FBSD/ImageTools/vr2i.h
index 03168d1158892215089a71944506e8587b295845..ac49a41c6a9c3a2c0d60b5cf3631e8165358d9ef 100755
--- a/Code/FBSD/ImageTools/vr2i.h
+++ b/Code/FBSD/ImageTools/vr2i.h
@@ -58,6 +58,17 @@ public:
     xv = x;
     yv = y; }
 
+  /**
+   * @fn void set (const Vr2i &p)
+   * \brief Copies the vector coordinates.
+   * @param vec Vector to copy.
+   */
+  inline void set (const Vr2i &vec)
+  {
+    xv = vec.xv;
+    yv = vec.yv;
+  }
+
   /**
    * @fn int norm2 ()
    * \brief Returns the squared norm of the vector.