diff --git a/CGALimageIO/include/CGAL/Image_3.h b/CGALimageIO/include/CGAL/Image_3.h
index 92029144376..bc0f166c030 100644
--- a/CGALimageIO/include/CGAL/Image_3.h
+++ b/CGALimageIO/include/CGAL/Image_3.h
@@ -198,9 +198,10 @@ Image_3::trilinear_interpolation(const Coord_type& x,
   const int dimz = zdim();
   const int dimxy = dimx*dimy;
 
-  const int i1 = (int)(z / image()->vx);
+  // images are indexed by (z,y,x)
+  const int i1 = (int)(z / image()->vz); 
   const int j1 = (int)(y / image()->vy);
-  const int k1 = (int)(x / image()->vz); 
+  const int k1 = (int)(x / image()->vx);
   const int i2 = i1 + 1;
   const int j2 = j1 + 1;
   const int k2 = k1 + 1;
@@ -209,16 +210,17 @@ Image_3::trilinear_interpolation(const Coord_type& x,
    *   a, b, c, d, e, f, g, h are the value of the image at the corresponding
    *   voxels:
    *
-   *     z
-   *     |  e______ h
+   *
+   *     x        z
+   *     |       /
+   *        f___ __ g
    *       /|      /|
-   *     f/_|____g/ |
+   *     e/_|____h/ |
    *     |  |    |  |
-   *     |  |a___|_d|
-   *     | /     | / 
-   *    b|/_____c|/  _y
+   *     |  |b___|_c|
+   *     | /     | /
+   *    a|/_____d|/  _y
    *    
-   *   x/
    *
    * a = val(i1, j1, k1)
    * b = val(i2, j1, k1)
@@ -237,7 +239,7 @@ Image_3::trilinear_interpolation(const Coord_type& x,
      z < 0.f ||
      i1 >= dimx ||
      j1 >= dimy ||
-     k1 >= dimy)
+     k1 >= dimz)
   {
     return outside;
   }
@@ -314,12 +316,12 @@ Image_3::trilinear_interpolation(const Coord_type& x,
     }
   }
 
-  const Target_word_type di2 = i2 - x;
-  const Target_word_type di1 = x - i1;
+  const Target_word_type di2 = i2 - z;
+  const Target_word_type di1 = z - i1;
   const Target_word_type dj2 = j2 - y;
   const Target_word_type dj1 = y - j1;
-  const Target_word_type dk2 = k2 - z;
-  const Target_word_type dk1 = z - k1;
+  const Target_word_type dk2 = k2 - x;
+  const Target_word_type dk1 = x - k1;
 //   std::cerr << di2 << " " << di1 << "\n";
 //   std::cerr << dj2 << " " << dj1 << "\n";
 //   std::cerr << dk2 << " " << dk1 << "\n";
diff --git a/CGALimageIO/test/CGALimageIO/test_trilinear_interpolation.cpp b/CGALimageIO/test/CGALimageIO/test_trilinear_interpolation.cpp
index 98db29b22a1..0977788eb93 100644
--- a/CGALimageIO/test/CGALimageIO/test_trilinear_interpolation.cpp
+++ b/CGALimageIO/test/CGALimageIO/test_trilinear_interpolation.cpp
@@ -15,40 +15,14 @@ int main() {
 
   ImageIO_free(image.data());
   image.set_data(&data[0]);
-  _image* im = image.image();
 
-//   im->xdim = 2;
-//   im->ydim = 2;
-//   im->zdim = 2;
-//   im->vdim = 1;
-//   im->vx = 1;
-//   im->vy = 1;
-//   im->vz = 1;
-
-//   im->cx = im->cy = im->cz = 0;
-//   im->tx = im->ty = im->tz = 0.0;
-//   im->rx = im->ry = im->rz = 0.0;
-
-
-//   im->fd = NULL;
-//   im->openMode = OM_CLOSE;
-//   im->endianness = END_UNKNOWN;
-
-//   im->dataMode = DM_BINARY;
-
-//   // word type (unsigned char)
-//   im->wdim = 1;
-//   im->wordKind = WK_FIXED;
-//   im->vectMode = VM_SCALAR;
-//   im->sign = SGN_UNSIGNED;
-//   im->imageFormat = NULL;
   std::cerr << std::setprecision(2) << std::fixed;
 
   for(int x = 0; x <= 1; ++x)
     for(int y = 0; y <= 1; ++y)
       for(int z = 0; z <= 1; ++z)
       {
-	data[x * 4 + y * 2 + z] = 1;
+	data[z * 4 + y * 2 + x] = 1;
 
 	std::cerr << "#### data"
 		  << "[" << x << "]"
@@ -63,10 +37,10 @@ int main() {
 	      assert((int)d_y == 0);
 	      assert((int)d_z == 0);
 	      const double value = 
-		image.trilinear_interpolation<Word, double, float>(d_x,
-								   d_y,
-								   d_z,
-								   255);
+		image.trilinear_interpolation<Word, double, double>(d_x,
+								    d_y,
+								    d_z,
+								    255);
 	      std::cerr << "val(" << d_x << ", " << d_y << " , " << d_z << ") = "
 			<< value << std::endl;
 	      const double sq_dist = 
@@ -83,9 +57,27 @@ int main() {
 		assert(value >= 0.12);
 	      else
 		assert(value <= 0.001);
+
+	      const float value2 = 
+		image.trilinear_interpolation<Word, float, float>(d_x,
+								  d_y,
+								  d_z,
+								  0);
+
+	      const float value3 = triLinInterp(image.image(),
+						d_x, 
+						d_y,
+						d_z,
+						0.f);
+	      std::cerr << "tri(" << d_x << ", " << d_y << " , " << d_z << ") = "
+			<< value3 << std::endl;
+	      if(value2 != value3)
+		std::cerr << std::setprecision(30)
+			  << "   " << value2 << "\n!= " << value3 << std::endl;
+	      assert(value2 == value3);
 	    }
 
-	data[x * 4 + y * 2 + z] = 0;
+	data[z * 4 + y * 2 + x] = 0;
       }
   image.set_data(0);
 }