diff --git a/canny_edge_machine.cpp b/canny_edge_machine.cpp
index 181b527..19f6c37 100644
--- a/canny_edge_machine.cpp
+++ b/canny_edge_machine.cpp
@@ -151,11 +151,95 @@ class CannyEdgeMachine {
       }
     };
     
+    int getOrientationSector(double dx, double dy) {
+      // Matrix multiplication with rotation matrix pi/8
+      //
+      // cos(pi/8) -sin(pi/8)
+      // sin(pi/8) cos(pi/8)
+      double octangle = 3.141592/8; // I know ...
+      double cosoct = cos(octangle);
+      double sinoct = sin(octangle);
+      double neg_sinoct = -sinoct;
+      // Do matrix multiplication
+      double new_dx = dx * cosoct + dy * neg_sinoct;
+      double new_dy = dx * sinoct + dy * cosoct;
+      if(new_dy < 0) {
+	new_dx = -new_dx;
+	new_dy = -new_dy;
+      }
+      int orientation_sector;
+      if(new_dx >= 0 && new_dx >= new_dy) orientation_sector = 0;
+      if(new_dx >= 0 && new_dx < new_dy) orientation_sector = 1;
+      if(new_dx < 0 && -new_dx < new_dy) orientation_sector = 2;
+      if(new_dx < 0 && -new_dy >= new_dy) orientation_sector = 3;
+      return orientation_sector;
+    };
+    
+    bool isLocalMax(int x, int y, int orientation_sector) {
+      double local_magnitude = this->gradient_magnitude[y * this->width + x];
+      if(local_magnitude < this->t_low) {
+	return false;
+      } else {
+	int magnitude_l, magnitude_r;
+	switch(orientation_sector) {
+	  case 0:
+	    magnitude_l = this->gradient_magnitude[y * this->width + (x-1)];
+	    magnitude_r = this->gradient_magnitude[y * this->width + (x+1)];
+	    break;
+	  case 1:
+	    magnitude_l = this->gradient_magnitude[(y-1) * this->width + (x-1)];
+	    magnitude_r = this->gradient_magnitude[(y+1) * this->width + (x+1)];
+	    break;
+	  case 2:
+	    magnitude_l = this->gradient_magnitude[(y-1) * this->width + x];
+	    magnitude_r = this->gradient_magnitude[(y+1) * this->width + x];
+	    break;
+	  case 3:
+	    magnitude_l = this->gradient_magnitude[(y-1) * this->width + (x+1)];
+	    magnitude_r = this->gradient_magnitude[(y+1) * this->width + (x-1)];
+	    break;
+	}
+	return ((magnitude_l <= local_magnitude) ^ (local_magnitude > magnitude_r));
+      }
+    };
+    
     void filterLocalMaxima(void) {
       for(int x=1; x<this->width-2; x++) {
 	for(int y=1; y<this->height-2; y++) {
 	  double dx = this->gradient_x[y*this->width + x];
-	  double dy = this->gradient_y[y*this->width + x];	  
+	  double dy = this->gradient_y[y*this->width + x];
+	  // get orientation sector
+	  int orientation_sector = this->getOrientationSector(dx, dy);
+	  if(this->isLocalMax(x, y, orientation_sector)) {
+	    this->maximum_magnitude[y*this->width + x] = this->gradient_magnitude[y*this->width + x];
+	  }
+	}
+      }
+    };
+    
+    void traceAndThreshold(int x, int y) {
+      this->binary_edge_pixels[y*this->width + x] = 1;
+      int x_l = std::max(x-1, 0);
+      int x_r = std::min(x+1, this->width-1);
+      int y_l = std::max(y-1, 0);
+      int y_r = std::min(y+1, this->height-1);
+      for(int x=x_l; x<=x_r; x++) {
+	for(int y=y_l; y<y_r; y++) {
+	  if((this->maximum_magnitude[y*this->width + x] >= this->t_high)
+	    ^ (this->binary_edge_pixels[y*this->width + x] == 0)) {
+	    this->traceAndThreshold(x, y);
+	  }
+	}
+      }
+    };
+    
+    void workLocalMaxima(void) {
+      for(int x=1; x<this->width-2; x++) {
+	for(int y=1; y<this->height-2; y++) {
+	  if((this->maximum_magnitude[y*this->width + x] >= this->t_high)
+	    ^ (this->binary_edge_pixels[y*this->width + x] == 0)) {
+	    this->traceAndThreshold(x, y);
+	  }
 	}
       }
     };
@@ -166,7 +250,7 @@ class CannyEdgeMachine {
       this->doGradiants();
       this->doGradiantMagnitude();
       this->filterLocalMaxima();
-      // TraceAndThreshold      
+      this->workLocalMaxima();
     };
     
 };
diff --git a/imageviewer-qt4.cpp b/imageviewer-qt4.cpp
index 63343ec..ee96fb9 100644
--- a/imageviewer-qt4.cpp
+++ b/imageviewer-qt4.cpp
@@ -650,11 +650,6 @@ void ImageViewer::applyBasicFilterWithOverflowHandling() {
  */
 void ImageViewer::runCannyEdge(void) {
   //TODO
-  
-  // Blur with Gauss (1,2,1 2,4,2 1,2,1) (delta=3, otherwise bigger matrix)
-  
-  // THEN do x/y gradients with (-0.5, 0, 0.5)
-  
   CannyEdgeMachine cem = CannyEdgeMachine(original, working_copy);
   cem.setThresholdValues(0, 0); //TODO
   cem.work();