//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2016 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include "Image.hpp"

#include "Information.hpp"
#include <QImageReader>
#include <QImageWriter>
#include <math.h>
#include "Progress.hpp"
#include "Process.hpp"

#include <tiffio.h>
#include <tiff.h>
#include <QDomDocument>
#include <QDomElement>
#include <CImg.h>

namespace mgx 
{
  using namespace cimg_library;
  typedef CImg<ushort> CImgUS;
  
  QStringList supportedImageReadFormats()
  {
    static QStringList formats;
    if(formats.empty()) {
      formats << "dlm" << "tif" << "cimg" << "CImg Auto";
      // foreach(QByteArray ba, QImageReader::supportedImageFormats())
      //{
      //  formats << QString::fromLocal8Bit(ba).toLower();
      //}
      // formats.removeDuplicates();
      formats.sort();
    }
    return formats;
  }
  
  QStringList supportedImageWriteFormats()
  {
    static QStringList formats;
    if(formats.empty()) {
      formats << "dlm" << "tif" << "cimg" << "CImg Auto";
      // foreach(QByteArray ba, QImageWriter::supportedImageFormats())
      //{
      //  formats << QString::fromLocal8Bit(ba).toLower();
      //}
      // formats.removeDuplicates();
      formats.sort();
    }
    return formats;
  }
  
  static void save_dlm(const QString& fn, const CImgUS& img) { img.save_dlm(fn.toLocal8Bit()); }
  static void save_cimg(const QString& fn, const CImgUS& img) { img.save_cimg(fn.toLocal8Bit()); }
  static void save_tiff(const QString& fn, const CImgUS& img) { img.save_tiff(fn.toLocal8Bit()); }
  static void save_auto(const QString& fn, const CImgUS& img) { img.save(fn.toLocal8Bit()); }
  
  bool saveImage(QString filename, const Image3D& data, QString type, unsigned int nb_digits)
  {
    if(type.isEmpty()) {
      if(filename.contains("."))
        type = filename.mid(filename.indexOf('.') + 1);
      else
        throw QString("Either the type is set, or the filename must specify an extension.");
    }
  
    const Point3u& size = data.size;
  
    unsigned int min_digits = (unsigned int)ceil(log10((float)size.z()));
  
    if(nb_digits > 0 and nb_digits < min_digits)
      throw QString("Error, you specified %1 digits for %2 images. This is not enough.").arg(nb_digits).arg(size.z());
  
    if(type != "CImg Auto") {
      if(filename.contains(".")) {
        QString ext = filename.mid(filename.lastIndexOf(".") + 1);
        if(ext.toLower() != type)
          filename += "." + type;
      } else {
        filename += "." + type;
      }
    } else if(nb_digits == 0)
      nb_digits = min_digits;
  
    void (*save)(const QString& fn, const CImgUS& img) = 0;
    if(type == "dlm")
      save = &save_dlm;
    else if(type == "cimg")
      save = &save_cimg;
    else if(type == "tif")
      save = &save_tiff;
    else
      save = &save_auto;
    if(nb_digits == 0) {
      CImgUS image(&(*data.data)[0], size.x(), size.y(), size.z(), 1, true);
      save(filename, image);
    } else {
      int index_dot = filename.lastIndexOf(".");
      QString prefix = filename.left(index_dot);
      QString suffix = filename.mid(index_dot);
      unsigned int shift = size.x() * size.y();
      for(unsigned int z = 0; z < size.z(); ++z) {
        CImgUS image(&(*data.data)[z * shift], size.x(), size.y(), 1, 1, true);
        QString fn = QString("%1%2%3").arg(prefix).arg(z, nb_digits, 10, QLatin1Char('0')).arg(suffix);
        save(fn, image);
      }
    }
    return true;
  }
  
  HVecUS resize(const HVecUS& data, const Point3i& before, const Point3i& after, bool center)
  {
    HVecUS result;
    size_t tot_before = size_t(before.x()) * size_t(before.y()) * size_t(before.z());
    size_t tot_after = size_t(after.x()) * size_t(after.y()) * size_t(after.z());
    if(tot_before != data.size())
      return result;
    result.resize(tot_after, 0);
    uint length_x = std::min(before.x(), after.x());
    uint length_y = std::min(before.y(), after.y());
    uint length_z = std::min(before.z(), after.z());
    uint before_shift_x = 0;
    uint before_shift_y = 0;
    uint before_shift_z = 0;
    uint after_shift_x = 0;
    uint after_shift_y = 0;
    uint after_shift_z = 0;
    uint before_inter_y = before.x() * (before.y() - length_y);
    uint after_inter_y = after.x() * (after.y() - length_y);
    if(center) {
      if(length_x == (uint)before.x())
        after_shift_x = (after.x() - length_x) / 2;
      else
        before_shift_x = (before.x() - length_x) / 2;
      if(length_y == (uint)before.y())
        after_shift_y = (after.y() - length_y) / 2;
      else
        before_shift_y = (before.y() - length_y) / 2;
      if(length_z == (uint)before.z())
        after_shift_z = (after.z() - length_z) / 2;
      else
        before_shift_z = (before.z() - length_z) / 2;
    }
    const ushort *b = &data[0] + (before_shift_z * before.y() + before_shift_y) * before.x() + before_shift_x;
    ushort* a = &result[(after_shift_z * after.y() + after_shift_y) * after.x() + after_shift_x];
    for(uint z = 0; z < length_z; ++z) {
      for(uint y = 0; y < length_y; ++y) {
        memcpy(a, b, sizeof(ushort) * length_x);
        b += before.x();
        a += after.x();
      }
      b += before_inter_y;
      a += after_inter_y;
    }
    return result;
  }
  
  Image3D::Image3D() : data(0), size(0u), step(1.0f), minc(0), maxc(0), brightness(1.0f),
                       plane(-1), labels(false), allocated(false) {}
  
  Image3D::Image3D(const Image3D& copy) : data(copy.data), size(copy.size), step(copy.step),
     minc(0), maxc(0), brightness(1.0f), plane(-1), labels(copy.labels), allocated(false) {}
  
  Image3D::Image3D(HVecUS& data, const Point3u& size, const Point3f& step, bool l)
    : data(&data), size(size), step(step), minc(0), maxc(0), brightness(1.0f), 
      plane(-1), labels(l), allocated(false) {}
  
  void Image3D::allocate(const Point3u& s)
  {
    size = s;
    allocate();
  }
  
  void Image3D::allocate()
  {
    minc = 0;
    maxc = 0;
    size_t s = size_t(size.x()) * size_t(size.y()) * size_t(size.z());
    if(s > 0) {
      if(data) {
        data->resize(s);
      } else {
        allocated = true;
        data = new HVecUS(s);
      }
    } else {
      if(allocated) {
        delete data;
        allocated = false;
        data = 0;
      } else
        data->clear();
    }
  }
  
  Image3D::~Image3D()
  {
    if(data and allocated)
      delete data;
    allocated = false;
    data = 0;
  }
  
  static void tiffErrorHandler(const char* module, const char* fmt, va_list ap)
  {
    Information::out << "TIFF Error in module " << module << endl;
    vfprintf(stderr, fmt, ap);
  }
  
  static void tiffWarningHandler(const char* module, const char* fmt, va_list ap)
  {
    Information::out << "TIFF warning in module " << module << endl;
    vfprintf(stderr, fmt, ap);
  }

  struct StripDataIncrement {
    int start_row;
    int start_column;
    int shift_row;
    int shift_column;
    int line_end;
  };
  
  template <typename T> struct StripBuffer {
    StripBuffer(T* buf)
      : buffer(buf)
    {
    }
    T* buffer;
    StripBuffer& operator++()
    {
      ++buffer;
      return *this;
    }
    void newLine() {
    }
    uint operator*() const;
  };
  
  template <> uint StripBuffer<signed char>::operator*() const
  {
    int val = *buffer;
    if(val < 0)
      return 0;
    return uint(val << 8);
  }
  
  template <> uint StripBuffer<unsigned char>::operator*() const 
  {
    return uint(*buffer) << 8;
  }
  
  template <> uint StripBuffer<short>::operator*() const
  {
    int val = *buffer;
    if(val < 0)
      return 0;
    return uint(val);
  }
  
  template <> uint StripBuffer<ushort>::operator*() const 
  {
    return uint(*buffer);
  }
  
  template <> uint StripBuffer<int>::operator*() const
  {
    int val = *buffer;
    if(val < 0)
      return 0;
    return uint(val >> 16);
  }
  
  template <> uint StripBuffer<uint>::operator*() const
  {
    int val = *buffer;
    return (val >> 16);
  }
  
  template <> uint StripBuffer<float>::operator*() const
  {
    float val = *buffer;
    if(val < 0)
      return 0;
    return uint(floor(val));
  }
  
  template <typename T> StripBuffer<T> stripBuffer(T* buf) 
  {
    return StripBuffer<T>(buf);
  }
  
  /*
     struct BitStripBuffer
     {
     BitStripBuffer(unsigned char* buf, int nb)
      : buffer(buf)
      , cur_part(0)
      , nb_bits(nb)
      , nb_parts(8/nb_bits)
     {
      if(nb_bits == 1)
        bit_mask = 0x8000;
      else if(nb_bits == 2)
        bit_mask = 0xB000;
      else if(nb_bits == 4)
        bit_mask = 0xF000;
      else
        bit_mask = 0xFFFF;
     }
  
     unsigned char* buffer;
     int cur_part;
     int nb_bits, nb_parts, bit_mask;
     BitStripBuffer& operator++()
     {
     ++cur_part;
      if(cur_part >= nb_parts)
      {
     ++buffer;
        cur_part = 0;
      }
      return *this;
     }
     void newLine() { if(cur_part != 0) { ++buffer; cur_part = 0; } }
     uint operator*() const
     {
      return ((uint(*buffer) << (16 - nb_bits*(nb_parts-cur_part))) & bit_mask);
     }
     };
   */
  
  template <typename Buffer>
  static void readStrip(Image3D& data, ushort*& p, const StripDataIncrement& incr, Buffer buf, size_t nb_pixels,
                        uint16 spp)
  {
    int col_count = 0;
    for(size_t i = 0; i < nb_pixels; ++i) {
      uint val = 0;
      for(size_t j = 0; j < spp; ++j, ++buf)
        val += *buf;
      val *= data.brightness;
      if(val > 0xFFFF)
        *p = 0xFFFF;
      else
        *p = ushort(val);
      p += incr.shift_column;
      if(++col_count >= incr.line_end) {
        col_count = 0;
        p += incr.shift_row;
        buf.newLine();
      }
      if(data.minc > val)
        data.minc = val;
      if(data.maxc < val)
        data.maxc = val;
    }
  }
  
  struct TIFFHandler 
  {
    typedef TIFF* pointer;
    TIFFHandler(TIFF* t)
      : tif(t)
    {
    }
  
    ~TIFFHandler()
    {
      if(tif)
        TIFFClose(tif);
      tif = 0;
    }
  
    operator bool() const { return bool(tif); }
  
    operator pointer() {
      return tif;
    }
  
    TIFF* tif;
  };
  
  struct TIFFAllocated 
  {
    TIFFAllocated(tdata_t t) : tif(t) {}
  
    ~TIFFAllocated()
    {
      if(tif)
        _TIFFfree(tif);
      tif = 0;
    }
  
    operator bool() const { return bool(tif); }
  
    operator tdata_t() {
      return tif;
    }
  
    tdata_t& pointer() {
      return tif;
    }
  
    tdata_t tif;
  };
  
  #define NB_UNITS 22
  
  static const QString units[NB_UNITS] = {
    "ym",    // yocto - 10^-24
    "zm",    // zepto - 10^-21
    "am",    // atto  - 10^-18
    "fm",    // femto - 10^-15
    "pm",    // pico  - 10^-12
    "nm",    // nano  - 10^-9
    "um",    // micro - 10^-6
    UM,      // micro - 10^-6
    "mm",    // milli - 10^-3
    "cm",    // centi - 10^-2
    "dm",    // deci  - 10^-1
    "m",     // unit  - 10^0
    "dam",   // deca  - 10^1
    "hm",    // hecto - 10^2
    "km",    // kilo  - 10^3
    "Mm",    // mega  - 10^6
    "Gm",    // giga  - 10^9
    "Tm",    // tera  - 10^12
    "Pm",    // peta  - 10^15
    "Em",    // exa   - 10^18
    "Zm",    // zetta - 10^21
    "Ym"     // yota  - 10^24
  };
  
  static const float unit_value[NB_UNITS] = { 1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-9, 1e-6, 1e-6, 1e-3, 1e-2, 1e-1,
                                              1,     1e1,   1e2,   1e3,   1e6,   1e9,  1e12, 1e15, 1e18, 1e21, 1e24 };
  
  static void readResolution(TIFF* tif, Image3D& data)
  {
    // First, check if this is ImageJ-like
    float xres = data.step.x(), yres = data.step.y(), zres = data.step.z();
    float unit = 1e-6;   // um by default
    uint16 resunit;
    char* desc;
    bool has_xres = TIFFGetField(tif, TIFFTAG_XRESOLUTION, &xres);
    bool has_yres = TIFFGetField(tif, TIFFTAG_YRESOLUTION, &yres);
    bool has_zres = false;
    bool has_desc = TIFFGetField(tif, TIFFTAG_IMAGEDESCRIPTION, &desc);
    bool has_resunit = TIFFGetFieldDefaulted(tif, TIFFTAG_RESOLUTIONUNIT, &resunit);
    bool labels = false;
    if(has_xres and xres > 0 and has_yres and yres > 0) {
      xres = 1.0 / xres;
      yres = 1.0 / yres;
    }
    if(has_resunit) {
      if(resunit == RESUNIT_INCH) {
        unit = 0.0254;
      } else if(resunit == RESUNIT_CENTIMETER) {
        unit = 1e-2;
      }
    }
    QString description = (has_desc ? QString::fromLocal8Bit(desc) : "");
    if(description.startsWith("ImageJ")) {
      Information::out << "**** Detected ImageJ TIFF:" << description << endl;
      QStringList fields = description.split("\n");
      foreach(QString f, fields) {
        QStringList entry = f.split("=");
        if(entry.size() == 2) {
          QString key = entry[0].trimmed();
          QString value = entry[1].trimmed();
          if(key == "spacing") {
            Information::out << "Found spacing = " << value << endl;
            has_zres = true;
            bool ok;
            zres = value.toFloat(&ok);
            if(!ok or isNan(zres) or zres <= 0)
              zres = 1.0;
          } else if(key == "unit") {
            for(int i = 0; i < NB_UNITS; ++i) {
              if(value == units[i]) {
                unit = unit_value[i];
                Information::out << "Found unit = " << value << " => " << unit << "m" << endl;
                break;
              }
            }
          } else if(key == "labels") {
            Information::out << "Found labels = " << value << endl;
            static QStringList t = QStringList() << "t" << "true" << "on" << "yes" << "y" << "1";
            labels = t.contains(value.toLower());
          }
        }
      }
    } else if(description.startsWith("<?xml"))   // This may not be OME TIFF, but it is at least XML
    {
      QDomDocument doc;
      Information::out << "**** Detected XML-TIFF" << endl;
      if(doc.setContent(description)) {
        QDomElement root = doc.documentElement();
        if(root.tagName() == "OME") {
          Information::out << "**** Detected OME-TIFF" << endl;
          // We have an OME-TIFF: find the resolution
          QDomElement image = root.firstChildElement("Image");
          if(!image.isNull()) {
            QDomElement pixels = image.firstChildElement("Pixels");
            if(!pixels.isNull()) {
              QString val = pixels.attribute("PhysicalSizeZ");
              bool ok;
              if(!val.isEmpty()) {
                zres = val.toDouble(&ok) * 1e-6 / unit;
                Information::out << "zres = " << zres << endl;
                has_zres = ok;
              }
              /*
               * val = pixels.attribute("SizeZ");
               * if(!val.isEmpty())
               *{
               *  int nb_slices = val.toInt(&ok);
               *  if(ok)
               *    zres /= float(nb_slices);
               *}
               */
              /*
                 val = pixels.attribute("PhysicalSizeY");
                 double r;
                 if(!val.isEmpty())
                 {
                 r = val.toDouble(&ok) * 1e-6;
                 if(ok) yres = r;
                 }
                 val = pixels.attribute("PhysicalSizeX");
                 if(!val.isEmpty())
                 {
                 r = val.toDouble(&ok) * 1e-6;
                 if(ok) xres = r;
                 }
               */
            }
          }
        }
      }
    } else
      Information::out << "*** TIFF of unknown source" << endl;

    unit *= 1e6;
    if(not has_zres) {
      zres = 1; // Make the z-step one micron
      Information::out << QString::fromWCharArray(L"Warning, no z resolution found, using default: %1 %2")
         .arg(unit*zres).arg(UM) << endl;
    }
    data.step = unit * Point3f(xres, yres, zres);
    data.labels = labels;
  }
  
  bool loadTIFFImage(QString filename, Image3D& data, bool allocate_data)
  {
    QByteArray ba = filename.toLocal8Bit();
    TIFFSetErrorHandler(tiffErrorHandler);
#ifdef WIN32
    // the default warning handler on Windows pops up dialogs
    // redirect them to the console window
    TIFFSetWarningHandler(tiffWarningHandler);
#endif
    TIFFHandler tif = TIFFOpen(ba.data(), "r");
    uint32 depth = 0;
    data.step = Point3f(1, 1, 1);
    if(tif) {
      tdir_t first_dir = TIFFCurrentDirectory(tif);
      // First, scan the number of directories
      do {
        depth++;
      } while(TIFFReadDirectory(tif));
      TIFFSetDirectory(tif, first_dir);
      readResolution(tif, data);
      Information::out << depth << " images in this TIFF file" << endl;
      do {
        uint32 w, h, rps;
        uint32* bc;
        uint16 planar_config, bps, spp, format, orientation;
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_IMAGEWIDTH, &w) != 1) {
          Information::out << "Error, no width" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_IMAGELENGTH, &h) != 1) {
          Information::out << "Error, no height" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_BITSPERSAMPLE, &bps) != 1) {
          Information::out << "Error, no bits per sample" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLESPERPIXEL, &spp) != 1) {
          Information::out << "Error, no samples per pixel" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_PLANARCONFIG, &planar_config) != 1) {
          Information::out << "Error, no planar configuration" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_STRIPBYTECOUNTS, &bc) != 1) {
          Information::out << "Error, no strip bute counts" << endl;
          return false;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &format) != 1) {
          format = 1;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_ORIENTATION, &orientation) != 1) {
          orientation = 1;
        }
        if(TIFFGetFieldDefaulted(tif, TIFFTAG_ROWSPERSTRIP, &rps) != 1) {
          Information::out << "Error, no rows per strip" << endl;
          return false;
        }
        if(bps != 8 and bps != 16 and bps != 32) {
          Information::out << "Error, this reader handles only 8, 16 and 32 bits per sample" << endl;
          return false;
        }
        StripDataIncrement incr;
        switch(orientation) {
        case ORIENTATION_TOPRIGHT:       // line = top, column = right
          incr.shift_column = -1;
          incr.shift_row = 0;
          incr.start_column = w - 1;
          incr.start_row = h - 1;
          incr.line_end = w;
          break;
        case ORIENTATION_BOTRIGHT:       // line = bottom column = right
          incr.shift_column = -1;
          incr.shift_row = 2 * w;
          incr.start_column = w - 1;
          incr.start_row = 0;
          incr.line_end = w;
          break;
        case ORIENTATION_BOTLEFT:       // line = bottom, column = left
          incr.shift_column = 1;
          incr.shift_row = 0;
          incr.start_column = 0;
          incr.start_row = 0;
          incr.line_end = w;
          break;
        case ORIENTATION_LEFTTOP:       // line = left, column = top
          std::swap(w, h);
          incr.shift_column = -w;
          incr.shift_row = h * w + 1;
          incr.start_column = 0;
          incr.start_row = h - 1;
          incr.line_end = h;
          break;
        case ORIENTATION_RIGHTTOP:       // line = right, column = top
          std::swap(w, h);
          incr.shift_column = -w;
          incr.shift_row = h * w - 1;
          incr.start_column = w - 1;
          incr.start_row = h - 1;
          incr.line_end = h;
          break;
        case ORIENTATION_RIGHTBOT:       // line = right, column = bottom
          std::swap(w, h);
          incr.shift_column = w;
          incr.shift_row = -h * w - 1;
          incr.start_column = w - 1;
          incr.start_row = 0;
          incr.line_end = h;
          break;
        case ORIENTATION_LEFTBOT:       // line = left, column = bottom
          std::swap(w, h);
          incr.shift_column = w;
          incr.shift_row = -h * w + 1;
          incr.start_column = 0;
          incr.start_row = 0;
          incr.line_end = h;
          break;
        case ORIENTATION_TOPLEFT:       // line = top, column = left
        default:
          incr.shift_column = 1;
          incr.shift_row = -2 * w;
          incr.start_column = 0;
          incr.start_row = h - 1;
          incr.line_end = w;
          break;
        }
        if(allocate_data) {
          allocate_data = false;
          if(data.size.x() != w or data.size.y() != h or (data.plane == -1 and data.size.z() != depth)) {
            Point3u size = Point3u(w, h, depth);
            data.allocate(size);
          }
        } else if(data.data) {
          if(data.plane == -1) {
            if(data.size.x() != w or data.size.y() != h or data.size.z() != depth) {
              Information::out << QString("Error, writing an image of size %1x%2x%3, but an image of size %4x%5x%6 "
                                "has been allocated")
                        .arg(w).arg(h).arg(depth).arg(data.size.x()).arg(data.size.y()).arg(data.size.z()) << endl;
              return false;
            }
          } else {
            if(data.size.x() != w or data.size.y() != h or data.plane < 0 or (size_t) data.plane >= data.size.z()) {
              Information::out << QString("Error, writing on plane %3 an image of size %1x%2, but an image of size "
                                "%4x%5x%6 has been allocated")
                        .arg(w).arg(h).arg(data.plane).arg(data.size.x()).arg(data.size.y()).arg(data.size.z()) << endl;
              return false;
            }
          }
        }
        if(data.plane == -1)
          data.plane = 0;
        if(data.data) {
          size_t tot_pixels = 0, max_pixels = size_t(data.size.x()) * size_t(data.size.y());
          ushort *p = &data[0] + size_t(data.size.x()) * data.size.y() * data.plane + incr.start_column + incr.start_row * w;
          switch(planar_config) {
          case PLANARCONFIG_SEPARATE: {
            Information::out << "Error, PLANARCONFIG_SEPARATE is not yet handled." << endl;
          } break;
          case PLANARCONFIG_CONTIG: {
            uint32 strip_size = TIFFStripSize(tif);
            size_t nb_pixels = rps * w;
            TIFFAllocated alloc = _TIFFmalloc(strip_size + 1);
            for(size_t strip = 0; strip < TIFFNumberOfStrips(tif); ++strip) {
              size_t size_buf = TIFFReadEncodedStrip(tif, strip, alloc.pointer(), (size_t)-1);
              if(size_buf == 0)
                Information::err << "Warning: empty buffer!" << endl;
              // Now, convert the data
              if((tot_pixels + nb_pixels) > max_pixels)
                nb_pixels = max_pixels - tot_pixels;
              tot_pixels += nb_pixels;
              switch(bps) {
              /*
                 case 1:
                 {
                  readStrip(data, p, incr, BitStripBuffer((unsigned char*)alloc.pointer(), 1), nb_pixels,
                 spp);
                 }
                 break;
                 case 2:
                 {
                  readStrip(data, p, incr, BitStripBuffer((unsigned char*)alloc.pointer(), 2), nb_pixels,
                 spp);
                 }
                 break;
                 case 4:
                 {
                  readStrip(data, p, incr, BitStripBuffer((unsigned char*)alloc.pointer(), 4), nb_pixels,
                 spp);
                 }
                 break;
               */
              case 8: {
                switch(format) {
                case SAMPLEFORMAT_UINT:
                  readStrip(data, p, incr, stripBuffer((unsigned char*)alloc.pointer()), nb_pixels, spp);
                  break;
                case SAMPLEFORMAT_INT:
                  readStrip(data, p, incr, stripBuffer((signed char*)alloc.pointer()), nb_pixels, spp);
                  break;
                default:
                  Information::out << QString("Error, cannot handle 8 bits data format: %1").arg(format) << endl;
                  return false;
                }
              } break;
              case 16: {
                switch(format) {
                case SAMPLEFORMAT_UINT:
                  readStrip(data, p, incr, stripBuffer((unsigned short*)alloc.pointer()), nb_pixels, spp);
                  break;
                case SAMPLEFORMAT_INT:
                  readStrip(data, p, incr, stripBuffer((signed short*)alloc.pointer()), nb_pixels, spp);
                  break;
                default:
                  Information::out << QString("Error, cannot handle 16 bits data format: %1").arg(format) << endl;
                  return false;
                }
              } break;
              case 32: {
                switch(format) {
                case SAMPLEFORMAT_UINT:
                  readStrip(data, p, incr, stripBuffer((uint*)alloc.pointer()), nb_pixels, spp);
                  break;
                case SAMPLEFORMAT_INT:
                  readStrip(data, p, incr, stripBuffer((int*)alloc.pointer()), nb_pixels, spp);
                  break;
                case SAMPLEFORMAT_IEEEFP:
                  readStrip(data, p, incr, stripBuffer((float*)alloc.pointer()), nb_pixels, spp);
                  break;
                default:
                  Information::out << QString("Error, cannot handle 32 bits data format: %1").arg(format) << endl;
                  return false;
                }
              } break;
              default:
                Information::out << QString("Error, cannot handle %1 bits per sample").arg(bps) << endl;
                return false;
              }
            }
          } break;
          default:
            Information::out << "Error, unknown planar configuration: " << planar_config << endl;
            return false;
          }
          if(tot_pixels != max_pixels)
            Information::out << "Warning, read " << tot_pixels << " when " << max_pixels << " where expected" << endl;
        } else {
          data.size.x() = w;
          data.size.y() = h;
          data.size.z() = depth;
          break;
        }
        if(!progressAdvance())
          throw UserCancelException();
        data.plane++;
      } while(TIFFReadDirectory(tif));
    } else {
      Information::out << "Cannot open image file " << filename << endl;
      return false;
    }
    return true;
  }
  
  bool saveTIFFImage(QString filename, const Image3D& data)
  {
    QByteArray ba = filename.toLocal8Bit();
    TIFFSetErrorHandler(tiffErrorHandler);
    TIFFHandler tif = TIFFOpen(ba.data(), "w");
    if(tif) {
      tsize_t image_size = data.size.x() * data.size.y();
      for(uint z = 0; z < data.size.z(); ++z) {
        uint32 w = data.size.x(), h = data.size.y(), rps = data.size.y();
        if(TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, 16) != 1) {
          Information::out << "Error writing bits per sample" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, 1) != 1) {
          Information::out << "Error writing samples per pixel" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG) != 1) {
          Information::out << "Error writing planar configuration" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_UINT) != 1) {
          Information::out << "Error writing sample format" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT) != 1) {
          Information::out << "Error writing orientation" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, w) != 1) {
          Information::out << "Error writing width" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_IMAGELENGTH, h) != 1) {
          Information::out << "Error writing height" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, rps) != 1) {
          Information::out << "Error writing rows per strip" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_COMPRESSION, COMPRESSION_LZW) != 1) {
          Information::out << "Error setting LZW compression scheme" << endl;
          return false;
        }
  
        float xres = 1 / data.step.x(), yres = 1 / data.step.y(), zres = data.step.z();
        uint16 resunit = RESUNIT_NONE;
  
        if(TIFFSetField(tif, TIFFTAG_XRESOLUTION, xres) != 1) {
          Information::out << "Error writing x resolution" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_YRESOLUTION, yres) != 1) {
          Information::out << "Error writing y resolution" << endl;
          return false;
        }
        QString description
          = QString("ImageJ=1.45l\nimages=%1\nslices=%1\nunit=um\nspacing=%2\nloop=false\nlabels=%3\n")
            .arg(data.size.z())
            .arg(zres)
            .arg(data.labels ? "true" : "false");
        QByteArray desc = description.toLocal8Bit();
        if(TIFFSetField(tif, TIFFTAG_IMAGEDESCRIPTION, desc.data()) != 1) {
          Information::out << "Error writing description" << endl;
          return false;
        }
        if(TIFFSetField(tif, TIFFTAG_RESOLUTIONUNIT, resunit) != 1) {
          Information::out << "Error writing resolution unit" << endl;
          return false;
        }
  
        if(TIFFWriteEncodedStrip(tif, 0, (tdata_t) &(*data.data)[image_size * z], 2 * image_size) == -1) {
          Information::out << "Error writing image data" << endl;
          return false;
        }
  
        TIFFWriteDirectory(tif);
      }
    }
    return true;
  }
  
  bool loadImage(QString filename, Image3D& data, bool allocate_data)
  {
    if(!(filename.endsWith(".tif", Qt::CaseInsensitive) or filename.endsWith(".tiff", Qt::CaseInsensitive))) {
      CImgUS image(filename.toLocal8Bit().data());
      if(allocate_data) {
        if((int)data.size.x() != image.width() or (int) data.size.y() != image.height()
           or (data.plane == -1 and (int) data.size.z() != image.depth())) {
          Point3u size = Point3u(image.width(), image.height(), 1);
          data.allocate(size);
        }
      }
      if(data.data) {
        if(int(data.size.x()) != image.width() or int(data.size.y()) != image.height()) {
          Information::out << QString("Error: Bad image file: %1 size: %2x%3").arg(filename).arg(image.width()).arg(
                      image.height()) << endl;
          return false;
        }
        ushort* p = 0;
        if(data.plane == -1) {
          p = &data[0];
        } else {
          p = &data[data.size.x() * data.size.y() * data.plane];
        }
        for(uint y = 0; y < data.size.y(); y++)
          for(uint x = 0; x < data.size.x(); x++) {
            float voxel = 0;
            for(uint s = 0; s < (uint)image.spectrum(); s++)
              voxel += image(x, y, 0, s);
            voxel *= data.brightness;
            *p++ = ushort(voxel);
            if(voxel > data.maxc)
              data.maxc = voxel;
            if(voxel < data.minc)
              data.minc = voxel;
          }
      } else {
        data.size.x() = image.width();
        data.size.y() = image.height();
        data.size.z() = image.depth();
      }
      return true;
    } else   // Use libtiff!
      return loadTIFFImage(filename, data, allocate_data);
  }
}
