diff --git a/app/paint/gimpheal.c b/app/paint/gimpheal.c
index 15a0108..f25b4e4 100644
--- a/app/paint/gimpheal.c
+++ b/app/paint/gimpheal.c
@@ -51,6 +51,20 @@
  * terrible and I can't see a use for it right now.
  */
 
+/* NOTES v2.6.10+
+ * I have changed the algorithm quite a bit as the proposed relighting invariant
+ * method was not giving proper results. The method used here is similar 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 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.
+ *
+ * Jean-Yves Couleaud cjyves@free.fr
+ */
 
 static gboolean     gimp_heal_start              (GimpPaintCore    *paint_core,
                                                   GimpDrawable     *drawable,
@@ -166,6 +180,9 @@ gimp_heal_start (GimpPaintCore     *paint_core,
  * 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)
+ *
+ * JYC: we don't use this function anymore
+ *
  */
 static void
 gimp_heal_substitute_0_for_1 (PixelRegion *pr)
@@ -201,6 +218,9 @@ gimp_heal_substitute_0_for_1 (PixelRegion *pr)
 
 /*
  * Divide topPR by bottomPR and store the result as a double
+ *
+ * JYC: Actually now, we substract.
+ *
  */
 static void
 gimp_heal_divide (PixelRegion *topPR,
@@ -231,7 +251,7 @@ 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;
@@ -246,6 +266,9 @@ gimp_heal_divide (PixelRegion *topPR,
 
 /*
  * multiply first by secondPR and store the result as a PixelRegion
+ *
+ * JYC: Actually now, we add.
+ *
  */
 static void
 gimp_heal_multiply (gdouble     *first,
@@ -276,7 +299,7 @@ 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 +325,96 @@ 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 + k) = 0.25 *
-                                    ( *(matrix - depth + k)       /* west */
-                                    + *(matrix + depth + k)       /* east */
-                                    + *(matrix - rowstride + k)   /* north */
-                                    + *(matrix + rowstride + k)); /* south */
+                  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 + k) - tmp;
+                  diff = solution[off + k] - tmp;
                   err += diff*diff;
                 }
             }
+        }
+    }
 
-          /* advance pointers to data */
-          matrix += depth;
-          solution += depth;
-          mask++;
+  /* do blacks */
+  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;
+
+          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;
+                }
+            }
         }
     }
 
@@ -359,7 +432,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 +448,18 @@ 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/
+ *
+ * The new algorithm is also from T. Georgiev.
  */
 static PixelRegion *
 gimp_heal_region (PixelRegion *tempPR,
@@ -398,7 +473,9 @@ gimp_heal_region (PixelRegion *tempPR,
   /* substitute 0's for 1's for the division and multiplication operations that
    * come later
    */
+#if 0
   gimp_heal_substitute_0_for_1 (srcPR);
+#endif
 
   /* divide tempPR by srcPR and store the result as a double in i_1 */
   gimp_heal_divide (tempPR, srcPR, i_1);
