Index: /trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- /trunk/psLib/src/imageops/psImageInterpolate.c	(revision 12775)
+++ /trunk/psLib/src/imageops/psImageInterpolate.c	(revision 12776)
@@ -7,6 +7,6 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-05 23:57:42 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-10 03:59:13 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -132,11 +132,18 @@
 
     int xNum, yNum;                     // Number of interpolation kernel pixels
+    int xCentral, yCentral;             // Central pixel of the convolution
     switch (options->mode) {
       case PS_INTERPOLATE_BILINEAR:
         xNum = yNum = 2;
+        // Central pixel is the pixel below the point of interest
+        xCentral = floor(x - 0.5 + FLT_EPSILON);
+        yCentral = floor(y - 0.5 + FLT_EPSILON);
         break;
       case PS_INTERPOLATE_BICUBE:
       case PS_INTERPOLATE_GAUSS:
         xNum = yNum = 3;
+        // Central pixel is the closest pixel to the point of interest
+        xCentral = x;
+        yCentral = y;
         break;
       case PS_INTERPOLATE_FLAT:
@@ -149,11 +156,9 @@
 
     const psImage *image = options->image; // Image of interest
-    int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x
-    int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y
     int xLast = image->numCols - 1;     // Last pixel in x
     int yLast = image->numRows - 1;     // Last pixel in y
 
-    if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast ||
-        yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) {
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
         // At least one pixel of the interpolation kernel is off the image
         if (imageValue) {
@@ -171,6 +176,6 @@
     switch (options->mode) {
       case PS_INTERPOLATE_BILINEAR: {
-          double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
-          double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
+          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
+          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
           kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac);
           kernel[0][1] = xFrac * (1.0 - yFrac);
@@ -180,6 +185,6 @@
       }
       case PS_INTERPOLATE_BICUBE: {
-          double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
-          double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
+          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
+          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
           // Calculation variables
           double xxFrac = xFrac * xFrac / 6.0;
@@ -200,13 +205,18 @@
       }
       case PS_INTERPOLATE_GAUSS: {
-          double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
-          double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
-          double xGaussFrac = 2.0 * erf((double)xNum / 4.0) - 1.0; // Fraction of Gaussian in x
-          double yGaussFrac = 2.0 * erf((double)yNum / 4.0) - 1.0; // Fraction of Gaussian in x
-          double norm = 1.0 / (double)xNum / (double)yNum *
-              (1.0 - xGaussFrac) / xGaussFrac * (1.0 - yGaussFrac) / yGaussFrac; // Normalisation
+          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
+          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
+          double sigma = 0.5;           // Gaussian sigma
+          double norm = 0.0;            // Normalisation
           for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {
               for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {
-                  kernel[j][i] = norm * exp(-0.5 * (PS_SQR(xPos - xFrac) + PS_SQR(yPos - yFrac)));
+                  norm += kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
+                                                                     PS_SQR(yPos - yFrac)));
+              }
+          }
+          norm = 1.0 / norm;
+          for (int j = 0; j < yNum; j++) {
+              for (int i = 0; i < xNum; i++) {
+                  kernel[j][i] *= norm;
               }
           }
@@ -224,6 +234,6 @@
     #define KERNEL_IMAGE_CASE(TYPE) \
         case PS_TYPE_##TYPE: { \
-            for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
-                for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
+            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
+                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                     value += values[j][i] = kernel[j][i] * image->data.TYPE[yPix][xPix]; \
                 } \
@@ -262,6 +272,6 @@
         int badPix = 0;                 // Number of bad pixels
         *maskValue = 0;
-        for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) {
-            for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) {
+        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
+            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
                     badValue += values[j][i];
@@ -289,6 +299,6 @@
         #define KERNEL_VARIANCE_CASE(TYPE) \
             case PS_TYPE_##TYPE: { \
-                for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
-                    for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
+                for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
+                    for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                         *varianceValue += PS_SQR(kernel[j][i]) * variance->data.TYPE[yPix][xPix]; \
                     } \
@@ -337,9 +347,14 @@
         }
     } else {
+        double norm3 = 0.0;             // Normalisation
         double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
-        double norm2 = 2.0 / (num / 2);  // Normalisation for sinc functions
-        double pos = frac - (num - 1)/2;  // Position of interest
+        double norm2 = 4.0 / (float)num; // 2.0 / (num / 2);  // Normalisation for sinc functions
+        double pos = - (num - 1)/2 - frac;  // Position of interest
         for (int i = 0; i < num; i++, pos += 1.0) {
-            values[i] = norm1 * sin(M_PI * pos * norm2) * sin(M_PI_2 * pos * norm2) / PS_SQR(pos);
+            norm3 += values[i] = norm1 * sin(M_PI * pos * norm2) * sin(M_PI_2 * pos * norm2) / PS_SQR(pos);
+        }
+        norm3 = 1.0 / norm3;
+        for (int i = 0; i < num; i++, pos += 1.0) {
+            values[i] *= norm3;
         }
     }
@@ -373,12 +388,12 @@
     }
 
+    // Central pixel is the pixel below the point of interest
+    int xCentral = floor(x - 0.5), yCentral = floor(y - 0.5); // Central pixel of the convolution
     const psImage *image = options->image; // Image of interest
-    int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x
-    int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y
     int xLast = image->numCols - 1;     // Last pixel in x
     int yLast = image->numRows - 1;     // Last pixel in y
 
-    if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast ||
-        yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) {
+    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
+        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
         // At least one pixel of the interpolation kernel is off the image
         if (imageValue) {
@@ -400,5 +415,5 @@
       case PS_INTERPOLATE_LANCZOS3:
       case PS_INTERPOLATE_LANCZOS4: {
-          double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x
+          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
 #if 0
           if (fabs(xFrac) < DBL_EPSILON) {
@@ -411,5 +426,5 @@
           }
 #endif
-          double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y
+          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
 #if 0
           if (fabs(yFrac) < DBL_EPSILON) {
@@ -435,7 +450,7 @@
     #define SEPARATE_IMAGE_CASE(TYPE) \
       case PS_TYPE_##TYPE: { \
-        for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
+        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
             double xInterpValue = 0.0; /* Interpolation in x */ \
-            for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
+            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                 xInterpValue += values[j][i] = xKernel[i] * image->data.TYPE[yPix][xPix]; \
             } \
@@ -475,8 +490,8 @@
         int badPix = 0;                 // Number of bad pixels
         *maskValue = 0;
-        for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) {
+        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
             // Interpolation in x
             double xInterpValue = 0.0;
-            for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) {
+            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
                     xInterpValue += values[j][i];
@@ -506,7 +521,7 @@
         #define SEPARATE_VARIANCE_CASE(TYPE) \
           case PS_TYPE_##TYPE: { \
-            for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
+            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                 double xInterpValue = 0.0; /* Interpolation in x */ \
-                for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
+                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                     xInterpValue += PS_SQR(xKernel[i]) * variance->data.TYPE[yPix][xPix]; \
                 } \