diff --git a/app/paint/gimpheal.c b/app/paint/gimpheal.c
index 15a0108..d24f67e 100644
--- a/app/paint/gimpheal.c
+++ b/app/paint/gimpheal.c
@@ -45,26 +45,33 @@
 #include "gimp-intl.h"
 
 
-/* NOTES:
+
+/* NOTES
+ * The method used here is similar to the lighting invariant correctin method
+ * but slightly different: we do not divide the RGB components, but substract
+ * them I2 = I0 - I1, where I0 is the sample image to be corrected, I1 is the
+ * reference pattern. Then we solve DeltaI=0 (Laplace) with I2 Dirichlet
+ * conditions at the borders of the mask. The solver is a unoptimized red/black
+ * checker Gauss-Siedel with an over-relaxation factor of 1.8. It can benefit
+ * from a multi-grid evaluation of an initial solution before the main
+ * iteration loop.
+ * I reduced the convergence criteria to .1% (0.001) as we are dealing here
+ * with RGB integer components, more is overkill.
  *
- * I had the code working for healing from a pattern, but the results look
- * terrible and I can't see a use for it right now.
+ * Jean-Yves Couleaud cjyves@free.fr
  */
 
-
 static gboolean     gimp_heal_start              (GimpPaintCore    *paint_core,
                                                   GimpDrawable     *drawable,
                                                   GimpPaintOptions *paint_options,
                                                   const GimpCoords *coords,
                                                   GError          **error);
 
-static void         gimp_heal_substitute_0_for_1 (PixelRegion      *pr);
-
-static void         gimp_heal_divide             (PixelRegion      *topPR,
+static void         gimp_heal_sub                (PixelRegion      *topPR,
                                                   PixelRegion      *bottomPR,
                                                   gdouble          *result);
 
-static void         gimp_heal_multiply           (gdouble          *first,
+static void         gimp_heal_add                (gdouble          *first,
                                                   PixelRegion      *secondPR,
                                                   PixelRegion      *resultPR);
 
@@ -154,7 +161,7 @@ gimp_heal_start (GimpPaintCore     *paint_core,
   if (! source_core->set_source && gimp_drawable_is_indexed (drawable))
     {
       g_set_error_literal (error, GIMP_ERROR, GIMP_FAILED,
-			   _("Healing does not operate on indexed layers."));
+                           _("Healing does not operate on indexed layers."));
       return FALSE;
     }
 
@@ -162,50 +169,13 @@ gimp_heal_start (GimpPaintCore     *paint_core,
 }
 
 /*
- * Substitute any zeros in the input PixelRegion for ones.  This is needed by
- * the algorithm to avoid division by zero, and to get a more realistic image
- * since multiplying by zero is often incorrect (i.e., healing from a dark to
- * a light region will have incorrect spots of zero)
- */
-static void
-gimp_heal_substitute_0_for_1 (PixelRegion *pr)
-{
-  gint     i, j, k;
-
-  gint     height = pr->h;
-  gint     width  = pr->w;
-  gint     depth  = pr->bytes;
-
-  guchar  *pr_data = pr->data;
-
-  guchar  *p;
-
-  for (i = 0; i < height; i++)
-    {
-      p = pr_data;
-
-      for (j = 0; j < width; j++)
-        {
-          for (k = 0; k < depth; k++)
-            {
-              if (p[k] == 0)
-                p[k] = 1;
-            }
-
-          p += depth;
-        }
-
-      pr_data += pr->rowstride;
-    }
-}
-
-/*
- * Divide topPR by bottomPR and store the result as a double
+ * Substract bottomPR to topPR and store the result as a double
+ *
  */
 static void
-gimp_heal_divide (PixelRegion *topPR,
-                  PixelRegion *bottomPR,
-                  gdouble     *result)
+gimp_heal_sub (PixelRegion *topPR,
+               PixelRegion *bottomPR,
+               gdouble     *result)
 {
   gint     i, j, k;
 
@@ -231,9 +201,8 @@ gimp_heal_divide (PixelRegion *topPR,
         {
           for (k = 0; k < depth; k++)
             {
-              r[k] = (gdouble) (t[k]) / (gdouble) (b[k]);
+              r[k] = (gdouble) (t[k]) - (gdouble) (b[k]);
             }
-
           t += depth;
           b += depth;
           r += depth;
@@ -245,12 +214,13 @@ gimp_heal_divide (PixelRegion *topPR,
 }
 
 /*
- * multiply first by secondPR and store the result as a PixelRegion
+ * add first to secondPR and store the result as a PixelRegion
+ *
  */
 static void
-gimp_heal_multiply (gdouble     *first,
-                    PixelRegion *secondPR,
-                    PixelRegion *resultPR)
+gimp_heal_add (gdouble     *first,
+               PixelRegion *secondPR,
+               PixelRegion *resultPR)
 {
   gint     i, j, k;
 
@@ -276,7 +246,8 @@ gimp_heal_multiply (gdouble     *first,
         {
           for (k = 0; k < depth; k++)
             {
-              r[k] = (guchar) CLAMP0255 (ROUND (((gdouble) (f[k])) * ((gdouble) (s[k]))));
+              r[k] = (guchar) CLAMP0255 (ROUND (((gdouble) (f[k])) +
+                     ((gdouble) (s[k]))));
             }
 
           f += depth;
@@ -302,46 +273,95 @@ gimp_heal_laplace_iteration (gdouble *matrix,
                              gdouble *solution,
                              guchar  *mask)
 {
-  const gint  rowstride = width * depth;
-  gint        i, j, k;
-  gdouble     tmp, diff;
-  gdouble     err       = 0.0;
-
+  const gint    rowstride = width * depth;
+  gint          i, j, k, off, offm, offm0, off0;
+  gdouble       tmp, diff;
+  gdouble       err       = 0.0;
+  const gdouble w         = 1.80 * 0.25; /* Over-relaxation = 1.8 */
+
+/* we use a red/black checker model of the discretization grid */
+/* do reds */
   for (i = 0; i < height; i++)
     {
-      for (j = 0; j < width; j++)
+      off0  = i * rowstride;
+      offm0 = i * width;
+
+      for (j = i % 2; j < width; j += 2)
         {
-          if ((0 == *mask) ||
+          off  = off0 + j * depth;
+          offm = offm0 + j;
+
+          if ((0 == mask[offm]) ||
               (i == 0) || (i == (height - 1)) ||
               (j == 0) || (j == (width - 1)))
             {
               /* do nothing at the boundary or outside mask */
               for (k = 0; k < depth; k++)
-                *(solution + k) = *(matrix + k);
+                solution[off + k] = matrix[off + k];
             }
           else
             {
-              /* calculate the mean of four neighbours for all color channels */
-              /* v[i][j] = 0.45 * (v[i][j-1]+v[i][j+1]+v[i-1][j]+v[i+1][j]) */
+              /* Use Gauss Siedel to get the correction factor then over-relax
+               * it */
               for (k = 0; k < depth; k++)
                 {
-                  tmp = *(solution + k);
+                  tmp = solution[off + k];
+                  solution[off + k] = matrix[off + k] +
+                                      w *
+                                      ( matrix[off - depth + k]     /* west */
+                                      + matrix[off + depth + k]     /* east */
+                                      + matrix[off - rowstride + k] /* north */
+                                      + matrix[off + rowstride + k] - 4.0 *
+                                        matrix[off+k]);             /* south */
+
+                  diff = solution[off + k] - tmp;
+                  err += diff * diff;
+                }
+            }
+        }
+    }
 
-                  *(solution + k) = 0.25 *
-                                    ( *(matrix - depth + k)       /* west */
-                                    + *(matrix + depth + k)       /* east */
-                                    + *(matrix - rowstride + k)   /* north */
-                                    + *(matrix + rowstride + k)); /* south */
+/* Do blacks
+ * As we've done the reds earlier, we can use them right now to accelerate the
+ * convergence. So we have "solution" in the solver instead of "matrix" above
+ */
+  for (i = 0; i < height; i++)
+    {
+      off0 =  i * rowstride;
+      offm0 = i * width;
+
+      for (j = (i % 2) + 1; j < width; j += 2)
+        {
+          off = off0 + j * depth;
+          offm = offm0 + j;
 
-                  diff = *(solution + k) - tmp;
+          if ((0 == mask[offm]) ||
+              (i == 0) || (i == (height - 1)) ||
+              (j == 0) || (j == (width - 1)))
+            {
+              /* do nothing at the boundary or outside mask */
+              for (k = 0; k < depth; k++)
+                solution[off + k] = matrix[off + k];
+            }
+          else
+            {
+              /* Use Gauss Siedel to get the correction factor then over-relax
+               * it */
+              for (k = 0; k < depth; k++)
+                {
+                  tmp = solution[off + k];
+                  solution[off + k] = matrix[off + k] +
+                                      w *
+                                      ( solution[off - depth + k]     /* west */
+                                      + solution[off + depth + k]     /* east */
+                                      + solution[off - rowstride + k] /* north */
+                                      + solution[off + rowstride + k] - 4.0 *
+                                        matrix[off+k]);               /* south */
+
+                  diff = solution[off + k] - tmp;
                   err += diff*diff;
                 }
             }
-
-          /* advance pointers to data */
-          matrix += depth;
-          solution += depth;
-          mask++;
         }
     }
 
@@ -359,7 +379,7 @@ gimp_heal_laplace_loop (gdouble *matrix,
                         gdouble *solution,
                         guchar  *mask)
 {
-#define EPSILON   0.0001
+#define EPSILON   0.001
 #define MAX_ITER  500
   gint i;
 
@@ -375,16 +395,16 @@ gimp_heal_laplace_loop (gdouble *matrix,
       /* copy solution to matrix */
       memcpy (matrix, solution, width * height * depth * sizeof (double));
 
-      if (sqr_err < SQR (EPSILON))
+      if (sqr_err < EPSILON)
         break;
     }
 }
 
 /*
- * Algorithm Design:
+ * Original Algorithm Design:
  *
- * T. Georgiev, "Image Reconstruction Invariant to Relighting", EUROGRAPHICS
- * 2005, http://www.tgeorgiev.net/
+ * T. Georgiev, "Photoshop Healing Brush: a Tool for Seamless Cloning
+ * http://www.tgeorgiev.net/Photoshop_Healing.pdf
  */
 static PixelRegion *
 gimp_heal_region (PixelRegion *tempPR,
@@ -395,20 +415,16 @@ gimp_heal_region (PixelRegion *tempPR,
   gdouble *i_2  = g_new (gdouble, tempPR->h * tempPR->bytes * tempPR->w);
   guchar  *mask = temp_buf_get_data (mask_buf);
 
-  /* substitute 0's for 1's for the division and multiplication operations that
-   * come later
-   */
-  gimp_heal_substitute_0_for_1 (srcPR);
-
-  /* divide tempPR by srcPR and store the result as a double in i_1 */
-  gimp_heal_divide (tempPR, srcPR, i_1);
+  /* substract pattern to image and store the result as a double in i_1 */
+  gimp_heal_sub (tempPR, srcPR, i_1);
 
   /* FIXME: is a faster implementation needed? */
   gimp_heal_laplace_loop (i_1, tempPR->h, tempPR->bytes, tempPR->w, i_2, mask);
 
-  /* multiply a double by srcPR and store in tempPR */
-  gimp_heal_multiply (i_2, srcPR, tempPR);
+  /* add solution to original image and store in tempPR */
+  gimp_heal_add (i_2, srcPR, tempPR);
 
+  /* clean up */
   g_free (i_1);
   g_free (i_2);
 
