// Copyright (c) 2005-2008  INRIA Sophia-Antipolis (France).
//               2008 GeometryFactory
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// Author(s)     : Laurent Rineau, Pierre Alliez

#ifndef CGAL_IMAGE_3_H
#define CGAL_IMAGE_3_H

#include <CGAL/basic.h>

#include <boost/shared_ptr.hpp>

#include <boost/format.hpp>
#include <CGAL/ImageIO.h>
#include <CGAL/function_objects.h>

#include <limits>

class vtkImageData;

namespace CGAL {

class Image_3
{
  struct Image_deleter {
    void operator()(_image* image)
    {
      ::_freeImage(image);
    }
  };
public:
  typedef boost::shared_ptr<_image> Image_shared_ptr;
  typedef Image_shared_ptr Pointer;

protected:
  Image_shared_ptr image_ptr;

   // implementation in src/CGALimageIO/Image_3.cpp
  bool private_read(_image* im);

public:
  Image_3()
    : image_ptr()
  {
  }

  Image_3(const Image_3& bi)
    : image_ptr(bi.image_ptr)
  {
    std::cerr << "Image_3::copy_constructor\n";
  }

  Image_3(_image* im) 
  {
    private_read(im);
  }

  ~Image_3()
  {
  }

  const _image* image() const
  {
    return image_ptr.get();
  }

  _image* image()
  {
    return image_ptr.get();
  }

  const void* data() const
  {
    return image()->data;
  }

  void* data()
  {
    return image()->data;
  }

  void set_data(void* d)
  {
    image()->data = d;
  }

  unsigned int xdim() const { return image_ptr->xdim; }
  unsigned int ydim() const { return image_ptr->ydim; }
  unsigned int zdim() const { return image_ptr->zdim; }

  unsigned int size() const { return xdim() * ydim() * zdim(); }

  double vx() const { return image_ptr->vx; }
  double vy() const { return image_ptr->vy; }
  double vz() const { return image_ptr->vz; }

  float value(const unsigned int i,
              const unsigned int j,
              const unsigned int k) const
  {
    return ::evaluate(image(),i,j,k);
  }

public:

  bool read(const char* file)
  {
    return private_read(::_readImage(file));
  }

  bool read_raw(const char* file,
                const unsigned int rx,
                const unsigned int ry,
                const unsigned int rz,
                const double vx = 1,
                const double vy = 1,
                const double vz = 1,
		const unsigned int offset = 0)
  {
    return private_read(::_readImage_raw(file,
                                         rx,ry,rz,
                                         vx,vy,vz,offset));
  }

#ifdef CGAL_USE_VTK
  bool read_vtk_image_data(vtkImageData*);
#endif // CGAL_USE_VTK

  // implementation in src/CGALimageIO/Image_3.cpp
  void gl_draw(const float point_size,
               const unsigned char r,
               const unsigned char g,
               const unsigned char b);

  // implementation in src/CGALimageIO/Image_3.cpp
  void gl_draw_bbox(const float line_width,
                    const unsigned char red,
                    const unsigned char green,
                    const unsigned char blue);

public:
  template <typename Image_word_type,
	    typename Target_word_type,
	    typename Coord_type,
	    class Image_transform>
  Target_word_type 
  trilinear_interpolation(const Coord_type&x, 
			  const Coord_type&y, 
			  const Coord_type&z,
			  const Image_word_type& value_outside = 
			    Image_word_type(),
			  Image_transform transform = 
			    Image_transform() ) const;

  // default Image_transform = CGAL::Identity
  template <typename Image_word_type,
	    typename Target_word_type,
	    typename Coord_type>
  Target_word_type 
  trilinear_interpolation(const Coord_type&x, 
			  const Coord_type&y, 
			  const Coord_type&z,
			  const Image_word_type& value_outside = 
			  Image_word_type()) const 
  {
    return trilinear_interpolation<
      Image_word_type,
      Target_word_type>(x, y, z, value_outside,
			  CGAL::Identity<Image_word_type>());
  }
}; // end Image_3

template <typename Image_word_type,
	  typename Target_word_type,
	  typename Coord_type,
	  class Image_transform>
Target_word_type 
Image_3::trilinear_interpolation(const Coord_type& x, 
				 const Coord_type& y, 
				 const Coord_type& z,
				 const Image_word_type& value_outside,
				 Image_transform transform) const 
{
  const int dimx = xdim();
  const int dimy = ydim();
  const int dimz = zdim();
  const int dimxy = dimx*dimy;

  // 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()->vx);
  const int i2 = i1 + 1;
  const int j2 = j1 + 1;
  const int k2 = k1 + 1;

  /*   We assume (x,y,z) lies in the following cube.
   *   a, b, c, d, e, f, g, h are the value of the image at the corresponding
   *   voxels:
   *
   *
   *     x        z
   *     |       /
   *        f___ __ g
   *       /|      /|
   *     e/_|____h/ |
   *     |  |    |  |
   *     |  |b___|_c|
   *     | /     | /
   *    a|/_____d|/  _y
   *    
   *
   * a = val(i1, j1, k1)
   * b = val(i2, j1, k1)
   * c = val(i2, j2, k1)
   * d = val(i1, j2, k1)
   * e = val(i1, j1, k2)
   * f = val(i2, j1, k2)
   * g = val(i2, j2, k2)
   * h = val(i1, j2, k2)
   */

  const Target_word_type outside = transform(value_outside);

  if(x < 0.f ||
     y < 0.f ||
     z < 0.f ||
     i1 >= dimx ||
     j1 >= dimy ||
     k1 >= dimz)
  {
    return outside;
  }

  Target_word_type a, b, c, d, e, f, g, h; 

  if(k1 < 0) {
    a = b = c = d = outside;
  }
  else {
    if(j1 < 0) {
      a = b = outside;
    }
    else {
      if(i1 < 0)
	a = outside;
      else
	a = ((Image_word_type*)image()->data)[i1 * dimxy + j1 * dimx + k1];

      if(i2 >= dimx)
	b = outside;
      else 
	b = ((Image_word_type*)image()->data)[i2 * dimxy + j1 * dimx + k1];
    }

    if(j2 >= dimy) {
      c = d = outside;
    }
    else {
      if(i1 < 0)
	d = outside;
      else
	d = ((Image_word_type*)image()->data)[i1 * dimxy + j2 * dimx + k1];

      if(i2 >= dimx)
	c = outside;
      else
	c = ((Image_word_type*)image()->data)[i2 * dimxy + j2 * dimx + k1];
    }
  }

  if(k2 >= dimz) {
    e = f = g = h = outside;
  }
  else {
    if(j1 < 0) {
      e = f = outside;
    }
    else {
      if(i1 < 0)
	e = outside;
      else
	e = ((Image_word_type*)image()->data)[i1 * dimxy + j1 * dimx + k2];

      if(i2 >= dimx)
	f = outside;
      else 
	f = ((Image_word_type*)image()->data)[i2 * dimxy + j1 * dimx + k2];
    }

    if(j2 >= dimy) {
      g = h = outside;
    }
    else {
      if(i1 < 0)
	h = outside;
      else
	h = ((Image_word_type*)image()->data)[i1 * dimxy + j2 * dimx + k2];

      if(i2 >= dimx)
	g = outside;
      else
	g = ((Image_word_type*)image()->data)[i2 * dimxy + j2 * dimx + k2];
    }
  }

  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 - x;
  const Target_word_type dk1 = x - k1;
//   std::cerr << di2 << " " << di1 << "\n";
//   std::cerr << dj2 << " " << dj1 << "\n";
//   std::cerr << dk2 << " " << dk1 << "\n";

  return ( (  ( a * di2 + b * di1 ) * dj2 + 
	      ( d * di2 + c * di1 ) * dj1   ) * dk2 +
	   (  ( e * di2 + f * di1 ) * dj2 + 
	      ( h * di2 + g * di1 ) * dj1   ) * dk1 );
}


} // end namespace CGAL

 
#endif // CGAL_IMAGE_3_H