//
// 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 <Misc.hpp>
#include <SystemProcessLoad.hpp>
#include <QtXml>
#include <QRegExp>
#include <MGXViewer/qglviewer.h>
#include <Progress.hpp>
#include <GraphUtils.hpp>
#include <Information.hpp>

using namespace qglviewer;

namespace mgx 
{
  typedef Vector<1, float> Point1f;
  typedef Vector<1, int> Point1i;
  typedef Vector<1, uint> Point1u;
  typedef Vector<3, vertex> Point3v;
  
  struct FieldFormat {
    enum Type { BINARY, ASCII, APPENDED } type;
    size_t offset;
  
    FieldFormat() : type(ASCII), offset(0) {}
  
    FieldFormat(Type t, int o) : type(t), offset(o) {}
  
    FieldFormat& operator=(const FieldFormat& other)
    {
      type = other.type;
      offset = other.offset;
      return *this;
    }
  
    static FieldFormat fromElement(QDomElement& elm)
    {
      Type type;
      int offset = 0;
      QString format = elm.attribute("format").toLower();
      if(format == "ascii")
        type = ASCII;
      else if(format == "binary")
        type = BINARY;
      else if(format == "appended") {
        type = APPENDED;
        bool ok;
        QString offset_s = elm.attribute("offset", "no offset attribute");
        offset = offset_s.toUInt(&ok);
        if(!ok)
          throw QString("loadMesh::Error:invalid offset for appended format: %1").arg(offset_s);
      } else
        throw QString("loadMesh::Error:invalid format '%1': should be 'ascii', 'binary' or 'appended'.")
              .arg(format);
      return FieldFormat(type, offset);
    }
  };
  
  struct Types {
    enum Type { Int, UInt, Float };
  
    Types(Type t = Int, int s = 1) : size(s) , type(t) {}
  
    static Types fromString(const QString& s)
    {
      if(s == "Int8")
        return Types(Int, 1);
      else if(s == "Int16")
        return Types(Int, 2);
      else if(s == "Int32")
        return Types(Int, 4);
      else if(s == "Int64")
        return Types(Int, 8);
      else if(s == "UInt8")
        return Types(UInt, 1);
      else if(s == "UInt16")
        return Types(UInt, 2);
      else if(s == "UInt32")
        return Types(UInt, 4);
      else if(s == "UInt64")
        return Types(UInt, 8);
      else if(s == "Float32")
        return Types(Float, 4);
      else if(s == "Float64")
        return Types(Float, 8);
      else
        throw QString(
                "loadMesh::Error:type must follow the pattern {Int|UInt}{8|16|32|64} or be Float32 or Float64.");
    }
  
    int size;
    Type type;
  };
  
  struct AppendedData {
    enum Encoding { INVALID, BASE64, RAW };
  
    AppendedData() : encoding(INVALID) {}
  
    AppendedData(const QByteArray& d, Encoding enc) : data(d) , encoding(enc) {}
  
    QByteArray data;
    Encoding encoding;
  
    static int checkBeforeIgnoreSpaces(const QByteArray& data, char value, int idx)
    {
      for(int j = idx; j >= 0; ++j) {
        if(data[j] == value)
          return j;
        if(data[j] != ' ' and data[j] != '\t' and data[j] != '\n')
          break;
      }
      return -1;
    }
  
    static int checkAfterIgnoreSpaces(const QByteArray& data, char value, int idx)
    {
      for(int j = idx; j < data.size(); ++j) {
        if(data[j] == value)
          return j;
        if(data[j] != ' ' and data[j] != '\t' and data[j] != '\n')
          break;
      }
      return -1;
    }
  
    // Retrieve and remove the AppendedData from the ba
    static AppendedData fromByteArray(QByteArray& ba)
    {
      // First, find the first occurence of <\s*AppendedData\s*encoding=".*"\s*>.*_
      AppendedData result;
      int pos_appended_data = 0;
      Encoding encoding;
      while(true) {
        pos_appended_data = ba.indexOf("AppendedData", pos_appended_data);
        if(pos_appended_data == -1)
          break;
        if(checkBeforeIgnoreSpaces(ba, '<', pos_appended_data - 1) == -1)
          continue;
        int pos_e = checkAfterIgnoreSpaces(ba, 'e', pos_appended_data + 12);
        if(pos_e == -1)
          continue;
        if(ba.mid(pos_e, 8) != "encoding")
          continue;
        int pos_equal = checkAfterIgnoreSpaces(ba, '=', pos_e + 8);
        int pos_quote = checkAfterIgnoreSpaces(ba, '"', pos_equal + 1);
        int pos_end_quote = ba.indexOf('"', pos_quote + 1);
        QByteArray enc = ba.mid(pos_quote + 1, pos_end_quote - pos_quote - 1).trimmed().toLower();
        if(enc == "raw")
          encoding = RAW;
        else if(enc == "base64")
          encoding = BASE64;
        else
          throw QString("loadMesh::Error:Invalid AppendedData encoding: '%1'").arg(QString::fromLocal8Bit(enc));
        // Now, find the last </AppendedData>
        int pos_end_tag = ba.indexOf('>', pos_end_quote);
        if(pos_end_tag == -1)
          continue;
  
        int pos_underscore = ba.indexOf('_', pos_end_tag);
        if(pos_underscore == -1)
          continue;
  
        int last_pos = -1;
  
        while(true) {
          last_pos = ba.lastIndexOf("AppendedData", last_pos);
          if(last_pos == -1)
            break;
          int pos_slash = checkBeforeIgnoreSpaces(ba, '/', last_pos - 1);
          if(pos_slash == -1)
            continue;
          if(checkAfterIgnoreSpaces(ba, '>', last_pos + 12) == -1)
            continue;
          int pos_bracket = checkBeforeIgnoreSpaces(ba, '<', pos_slash - 1);
          if(pos_bracket == -1)
            continue;
          last_pos = pos_bracket;
          break;
        }
  
        if(last_pos != -1) {
          result.encoding = encoding;
          result.data = ba.mid(pos_underscore + 1, last_pos - pos_underscore - 1);
          ba.remove(pos_end_tag + 1, last_pos - pos_end_tag - 1);
          break;
        }
      }
      return result;
    }
  };
  
  template <typename T, typename EncodedType>
  static void vtkRawDecode(std::vector<T>& data, const QByteArray& text, size_t offset, size_t nb_comp)
  {
    // It has to be ascii
    QByteArray size = text.mid(offset, 4);
    int s = *reinterpret_cast<qint32*>(size.data());
    if(data.empty()) {
      if(s % sizeof(EncodedType) != 0)
        throw QString("loadMesh::Error:%1 is an invalid size for an array of %2").arg(s).arg(
                qdemangle(typeid(EncodedType).name()));
      int nb_et = s / sizeof(EncodedType);
      if(nb_et % nb_comp != 0)
        throw QString("loadMesh::Error:%1 is an invalid number of elements for an array of %2 components")
              .arg(nb_et)
              .arg(nb_comp);
      data.resize(nb_et / nb_comp);
    } else if((size_t)s != sizeof(EncodedType) * nb_comp * data.size())
      throw QString("loadMesh::Error:not enough data for the field: expected '%1' and got '%2'").arg(data.size()).arg(
              s / sizeof(EncodedType));
    const EncodedType* array = reinterpret_cast<const EncodedType*>(text.data() + offset + 4);
    for(size_t i = 0, k = 0; i < data.size(); ++i)
      for(size_t j = 0; j < nb_comp; ++j, ++k)
        data[i][j] = (typename T::value_type)array[k];
  }
  
  template <typename T, typename EncodedType>
  static void vtkBase64Decode(std::vector<T>& data, const QByteArray& text, size_t offset, size_t nb_comp)
  {
    // It has to be ascii
    QByteArray size = text.mid(offset, 6);
    QByteArray ba, decoded;
    size = QByteArray::fromBase64(size);
    if(size.size() != 4)
      throw QString("loadMesh::Error:invalid size specification in base64 text");
    int s = *reinterpret_cast<qint32*>(size.data());
    int encoded_size;
    if(text[int(offset + 6)] == '=') {
      while(text[int(offset + 6)] == '=')
        ++offset;
      encoded_size = s * 4;
      if(encoded_size % 3 != 0)
        encoded_size = encoded_size / 3 + 1;
      else
        encoded_size /= 3;
      if(int(encoded_size + offset) > text.size())
        throw QString("loadMesh::Error:Not enough data to read: %1 bytes expected, %2 bytes available.")
              .arg(encoded_size)
              .arg(text.size() - offset);
      ba = text.mid(offset + 6, encoded_size);
      decoded = QByteArray::fromBase64(ba);
    } else {
      encoded_size = (s + 4) * 4;
      if(encoded_size % 3 != 0)
        encoded_size = encoded_size / 3 + 1;
      else
        encoded_size /= 3;
      if(int(encoded_size + offset) > text.size())
        throw QString("loadMesh::Error:Not enough data to read: %1 bytes expected, %2 bytes available.")
              .arg(encoded_size)
              .arg(text.size() - offset);
      ba = text.mid(offset, encoded_size);
      decoded = QByteArray::fromBase64(ba).mid(4);
    }
    if(s != decoded.size()) {
      if(encoded_size < 20)
        Information::out << QString("Read %1 bytes from offset %2: '%3' when %4 bytes where expected")
          .arg(encoded_size).arg(offset).arg(QString::fromLocal8Bit(ba)).arg(s) << endl;
      else
        Information::out << QString("Read %1 bytes from offset %2 when %3 bytes where expected")
          .arg(encoded_size).arg(offset).arg(s) << endl;
      throw QString("loadMesh::Error:size specification (%1) doesn't correspond to actual data size (%2).")
            .arg(s)
            .arg(decoded.size());
    }
    if(data.empty()) {
      if(s % sizeof(EncodedType) != 0)
        throw QString("loadMesh::Error:%1 is an invalid size for an array of %2").arg(s).arg(
                qdemangle(typeid(EncodedType).name()));
      int nb_et = s / sizeof(EncodedType);
      if(nb_et % nb_comp != 0)
        throw QString("loadMesh::Error:%1 is an invalid number of elements for an array of %2 components")
              .arg(nb_et)
              .arg(nb_comp);
      data.resize(nb_et / nb_comp);
    } else if((size_t)s != sizeof(EncodedType) * nb_comp * data.size())
      throw QString("loadMesh::Error:not enough data for the field: expected '%1' and got '%2'").arg(data.size()).arg(
              s / sizeof(EncodedType));
    const EncodedType* array = reinterpret_cast<const EncodedType*>(decoded.data());
    for(size_t i = 0, k = 0; i < data.size(); ++i)
      for(size_t j = 0; j < nb_comp; ++j, ++k)
        data[i][j] = (typename T::value_type)array[k];
  }
  
  template <typename T, typename EncodedType>
  static void vtkAppendedDecode(std::vector<T>& data, const AppendedData& appended, size_t offset, size_t nb_comp)
  {
    switch(appended.encoding) {
    case AppendedData::BASE64:
      vtkBase64Decode<T, EncodedType>(data, appended.data, offset, nb_comp);
      break;
    case AppendedData::RAW:
      vtkRawDecode<T, EncodedType>(data, appended.data, offset, nb_comp);
      break;
    case AppendedData::INVALID:
      throw QString("loadMesh::Error:Tag content refers to non-existing AppendedData");
    }
  }
  
  template <typename T, typename EncodedType>
  static void vtkBinaryDecode(std::vector<T>& data, const QString& text, size_t nb_comp)
  {
    // It has to be ascii
    QByteArray ba = text.toLatin1();
    vtkBase64Decode<T, EncodedType>(data, ba, 0, nb_comp);
  }
  
  template <typename EncodedType> struct TypeAsciiReader;
  
  struct TypeIntAsciiReader {
    void setValue(int& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = int(value);
    }
  
    void setValue(uint& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = uint(value);
    }
  
    void setValue(long& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = value;
    }
  
    void setValue(ulong& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = ulong(value);
    }
  
    void setValue(long long& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = (long long)(value);
    }
  
    void setValue(unsigned long long& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = (unsigned long long)(value);
    }
  
    void setValue(float& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = float(value);
    }
  
    void setValue(double& v, const QString& s)
    {
      bool ok;
      long long value = s.toLongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = double(value);
    }
  };
  
  template <> struct TypeAsciiReader<qint8> : public TypeIntAsciiReader {
  };
  template <> struct TypeAsciiReader<qint16> : public TypeIntAsciiReader {
  };
  template <> struct TypeAsciiReader<qint32> : public TypeIntAsciiReader {
  };
  template <> struct TypeAsciiReader<qint64> : public TypeIntAsciiReader {
  };
  
  struct TypeUIntAsciiReader {
    void setValue(int& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = int(value);
    }
  
    void setValue(uint& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = uint(value);
    }
  
    void setValue(long& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = long(value);
    }
  
    void setValue(ulong& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = ulong(value);
    }
  
    void setValue(long long& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = (long long)(value);
    }
  
    void setValue(unsigned long long& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = (unsigned long long)(value);
    }
  
    void setValue(float& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = float(value);
    }
  
    void setValue(double& v, const QString& s)
    {
      bool ok;
      unsigned long long value = s.toULongLong(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = double(value);
    }
  };
  
  template <> struct TypeAsciiReader<quint8> : public TypeUIntAsciiReader {
  };
  template <> struct TypeAsciiReader<quint16> : public TypeUIntAsciiReader {
  };
  template <> struct TypeAsciiReader<quint32> : public TypeUIntAsciiReader {
  };
  template <> struct TypeAsciiReader<quint64> : public TypeUIntAsciiReader {
  };
  
  struct TypeFloatAsciiReader {
    void setValue(int& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = int(value);
    }
  
    void setValue(uint& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = uint(value);
    }
  
    void setValue(long& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = long(value);
    }
  
    void setValue(ulong& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = ulong(value);
    }
  
    void setValue(long long& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid integer");
      v = (long long)(value);
    }
  
    void setValue(unsigned long long& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      if(!ok)
        throw QString("loadMesh::Error:Invalid unsigned integer");
      v = (unsigned long long)(value);
    }
  
    void setValue(float& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      v = float(value);
    }
  
    void setValue(double& v, const QString& s)
    {
      bool ok;
      double value = s.toDouble(&ok);
      v = double(value);
    }
  };
  
  template <> struct TypeAsciiReader<float> : public TypeFloatAsciiReader {
  };
  
  template <> struct TypeAsciiReader<double> : public TypeFloatAsciiReader {
  };
  
  template <typename T, typename EncodedType>
  static void vtkAsciiDecode(std::vector<T>& data, const QString& text, size_t nb_comp)
  {
    TypeAsciiReader<EncodedType> decoder;
    QStringList values = text.split(QRegExp("[\\s\\n]+"));
    if(data.empty()) {
      if(data.size() % nb_comp != 0)
        throw QString("loadMesh::Error:invalid number of elements '%1' as there must be a multiple of '%2'")
              .arg(data.size())
              .arg(nb_comp);
      data.resize(values.size() / nb_comp);
    } else if(data.size() * nb_comp != (size_t)values.size())
      throw QString("loadMesh::Error:Error, expected %1 values, got %2 values").arg(data.size() * nb_comp).arg(
              values.size());
    for(size_t i = 0, k = 0; i < data.size(); ++i)
      for(size_t j = 0; j < nb_comp; ++j, ++k)
        decoder.setValue(data[i][j], values[k]);
  }
  
  QString textFromTag(QDomElement& tag)
  {
    if(tag.childNodes().length() != 1)
      throw QString("loadMesh::Error:%1 must contain just a text.").arg(tag.tagName());
    QDomText tag_text = tag.firstChild().toText();
    if(tag_text.isNull())
      throw QString("loadMesh::Error:%1 must contain just a text.").arg(tag.tagName());
    QString text = tag_text.data().trimmed();
    return text;
  }
  
  template <typename T, typename EncodedType>
  static void vtkDecode(std::vector<T>& data, QDomElement& tag, size_t nb_comp, const FieldFormat& format,
                        const AppendedData& appended)
  {
    switch(format.type) {
    case FieldFormat::BINARY:
      vtkBinaryDecode<T, EncodedType>(data, textFromTag(tag), nb_comp);
      break;
    case FieldFormat::ASCII: {
      vtkAsciiDecode<T, EncodedType>(data, textFromTag(tag), nb_comp);
    } break;
    case FieldFormat::APPENDED:
      vtkAppendedDecode<T, EncodedType>(data, appended, format.offset, nb_comp);
      break;
    }
  }
  
  template <typename T>
  static void vtkDecode(std::vector<T>& data, QDomElement& tag, size_t nb_comp, const FieldFormat& format,
                        const AppendedData& appended, Types type)
  {
    switch(type.type) {
    case Types::Int:
      switch(type.size) {
      case 1:
        vtkDecode<T, qint8>(data, tag, nb_comp, format, appended);
        return;
      case 2:
        vtkDecode<T, qint16>(data, tag, nb_comp, format, appended);
        return;
      case 4:
        vtkDecode<T, qint32>(data, tag, nb_comp, format, appended);
        return;
      case 8:
        vtkDecode<T, qint64>(data, tag, nb_comp, format, appended);
        return;
      default:
        throw QString("loadMesh::Error:Error, invalid type Int%1").arg(8 * type.size);
      }
      break;
    case Types::UInt:
      switch(type.size) {
      case 1:
        vtkDecode<T, quint8>(data, tag, nb_comp, format, appended);
        return;
      case 2:
        vtkDecode<T, quint16>(data, tag, nb_comp, format, appended);
        return;
      case 4:
        vtkDecode<T, quint32>(data, tag, nb_comp, format, appended);
        return;
      case 8:
        vtkDecode<T, quint64>(data, tag, nb_comp, format, appended);
        return;
      default:
        throw QString("loadMesh::Error:Error, invalid type UInt%1").arg(8 * type.size);
      }
      break;
    case Types::Float:
      switch(type.size) {
      case 4:
        vtkDecode<T, float>(data, tag, nb_comp, format, appended);
        return;
      case 8:
        vtkDecode<T, double>(data, tag, nb_comp, format, appended);
        return;
      default:
        throw QString("loadMesh::Error:Error, invalid type Float%1").arg(8 * type.size);
      }
    }
  }
  
  static Point3f realPos(Point3f pos, bool scale, bool transform, const Mesh* mesh)
  {
    if(!scale)
      pos = mesh->stack()->abstractToWorld(pos);
    if(transform)
      pos = Point3f(mesh->stack()->frame().coordinatesOf(Vec(pos)));
    return pos;
  }
  
  bool insertEdges(vvGraph& S, vertex v1, vertex v2, vertex v3)
  {
    if(S.empty(v1)) {
      S.insertEdge(v1, v2);
      S.insertEdge(v1, v3);
      return true;
    } else if(S.edge(v1, v2)) {
      if(!S.edge(v1, v3))
        S.spliceAfter(v1, v2, v3);
      else if(!S.nextTo(v1, v2) == v3)
        throw QString("loadMesh::Error:inconsistent triangulation");
      return true;
    } else if(S.edge(v1, v3)) {
      if(!S.edge(v1, v2))
        S.spliceBefore(v1, v3, v2);
      else if(!S.prevTo(v1, v3) == v2)
        throw QString("loadMesh::Error:inconsistent triangulation");
      return true;
    }
    return false;
  }
  
  static bool addTriangle(vvGraph& S, vertex v1, vertex v2, vertex v3, std::list<Point3v>& to_process)
  {
    if(!insertEdges(S, v1, v2, v3) or !insertEdges(S, v2, v3, v1) or !insertEdges(S, v3, v1, v2)) {
      to_process.push_back(Point3v(v1, v2, v3));
      return false;
    }
    return true;
  }
  
  bool MeshImport::loadMeshVTK(Mesh* mesh, const QString& filename, bool& scale, bool transform, bool add)
  {
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
      throw QString("loadMesh::Error:Cannot open input file: %1").arg(filename);
  
    QByteArray content = file.readAll();
  
    file.close();
  
    // First, find and remove appended data if any
  
    // Then, find if there is some Appended section
    AppendedData appended = AppendedData::fromByteArray(content);
  
    QDomDocument doc("VTKFile");
    if(!doc.setContent(content))
      throw QString("loadMesh::Error:File '%1' is not a valid XML file").arg(filename);
  
    QDomElement root = doc.documentElement();
    if(root.tagName() != "VTKFile")
      throw QString("loadMesh::Error:File '%1' is not a VTKFile: the root tag is '%2'").arg(filename).arg(
              root.tagName());
  
    if(root.attribute("type", "") != "UnstructuredGrid")
      throw QString("loadMesh::Error:File '%1' doesn't contain an UnstructuredGrid. (%2)").arg(filename).arg(
              root.attribute("type", ""));
  
    QDomNodeList grids = root.elementsByTagName("UnstructuredGrid");
    if(grids.isEmpty())
      throw QString("loadMesh::Error:File '%1' doesn't contain any UnstructuredGrid").arg(filename);
  
    vvGraph& S = mesh->graph();
  
    if(!add)
      S.clear();
  
    FieldFormat format;
    Types type;
  
    // First, find the total number of elements
    size_t nb_phases = 0, cur_phase = 0;
    for(int grid_id = 0; grid_id < grids.length(); ++grid_id) {
      QDomElement grid = grids.at(grid_id).toElement();
      QDomNodeList pieces = grid.elementsByTagName("Piece");
      if(pieces.isEmpty())
        throw QString("loadMesh::Error::File '%1' contains a grid without any piece.").arg(filename);
      nb_phases += 3 * pieces.size();
    }
  
    progressStart(QString("Loading mesh from file '%1'").arg(filename), nb_phases);
  
    for(int grid_id = 0; grid_id < grids.length(); ++grid_id) {
      QDomElement grid = grids.at(grid_id).toElement();
  
      // First, read global fields, if any, as they may contain things like scaling
      float scale_factor = 1.0f;
      QDomElement field_data = grid.firstChildElement("FieldData");
      if(!field_data.isNull()) {
        if(!field_data.nextSiblingElement("FieldData").isNull())
          throw QString("loadMesh::Error:more than one FieldData tag in a grid");
        QDomElement field_elm = field_data.firstChildElement();
        while(!field_elm.isNull()) {
          if(field_elm.tagName() != "DataArray")
            throw QString("loadMesh::Error:FieldData can only contain DataArray tags.");
          QString name = field_elm.attribute("Name");
          if(name == "Scale") {
            FieldFormat format = FieldFormat::fromElement(field_elm);
            Types type = Types::fromString(field_elm.attribute("type"));
            bool ok;
            QString nb_tuples_s = field_elm.attribute("NumberOfTuples", "1");
            uint nb_tuples = nb_tuples_s.toUInt(&ok);
            if(!ok)
              throw QString("loadMesh::Error:Invalid number of tuples: '%1'").arg(nb_tuples_s);
            if(nb_tuples == 1) {
              std::vector<Point1f> scales(1);
              vtkDecode(scales, field_elm, 1, format, appended, type);
              scale_factor = scales[0][0] * 1e6;             // in micro-meter for us!
              scale = true;
            }
          }
  
          field_elm = field_elm.nextSiblingElement();
        }
      }
  
      // Now, find the pieces
      QDomNodeList pieces = grid.elementsByTagName("Piece");
      if(pieces.isEmpty())
        throw QString("loadMesh::Error::File '%1' contains a grid without any piece.").arg(filename);
  
      for(int piece_id = 0; piece_id < pieces.length(); ++piece_id) {
        QDomElement piece = pieces.at(piece_id).toElement();
  
        QString nb_cells_str = piece.attribute("NumberOfCells");
        QString nb_points_str = piece.attribute("NumberOfPoints");
  
        if(nb_cells_str.isEmpty())
          throw QString("loadMesh::Error:Piece tag doesn't have NumberOfCells attribute");
  
        if(nb_points_str.isEmpty())
          throw QString("loadMesh::Error:Piece tag doesn't have NumberOfPoints attribute");
  
        bool ok;
        int nb_cells = nb_cells_str.toInt(&ok);
        if(!ok)
          throw QString("loadMesh::Error:NumberOfCells attribute of tag 'Piece' is not an integer");
        int nb_points = nb_points_str.toInt(&ok);
        if(!ok)
          throw QString("loadMesh::Error:NumberOfPoints attribute of tag 'Piece' is not an integer");
  
        SETSTATUS("Found piece with " << nb_points << " points and " << nb_cells << " cells.");
  
        SETSTATUS("Read points and points data");
  
        std::vector<Point3f> points_data(nb_points);
        std::vector<Point3f> points_normal(nb_points);
        std::vector<Point1f> points_signal(nb_points);
        std::vector<Point1f> points_color(nb_points, Point1f(1.0f));
        std::vector<Point1i> points_label(nb_points);
        std::vector<Point1u> cells_connectivity;
        std::vector<Point1u> cells_offsets(nb_cells);
        std::vector<Point1u> cells_types(nb_cells);
  
        bool has_normals = false;
  
        QDomNodeList points_elmts = piece.elementsByTagName("Points");
        if(points_elmts.length() != 1)
          throw QString("loadMesh::Error:Piece needs exactly one list of points");
        QDomElement points_elm = points_elmts.at(0).toElement();
  
        // Now, read the points
        QDomNodeList points_data_elms = points_elm.elementsByTagName("DataArray");
        if(points_data_elms.length() != 1)
          throw QString("loadMesh::Error:Points tag must contain a single DataArray");
        QDomElement points_data_elm = points_data_elms.at(0).toElement();
  
        QString nb_comp_str = points_data_elm.attribute("NumberOfComponents");
        if(nb_comp_str.isEmpty())
          throw QString("loadMesh::Error:NumberOfComponents missing in points data array");
        size_t nb_comp = nb_comp_str.toUInt(&ok);
        if(!ok)
          throw QString("loadMesh::Error:NumberOfComponents is not an integer");
        if(nb_comp > 3)
          throw QString("loadMesh::Error:MorphoGraphX can only handle up to 3 dimensions. PointsData has %1.")
                .arg(nb_comp);
        format = FieldFormat::fromElement(points_data_elm);
        QString s_type = points_data_elm.attribute("type");
        if(s_type.isEmpty())
          throw QString("loadMesh::Error:missing type in points data array");
        type = Types::fromString(s_type);
        vtkDecode(points_data, points_data_elm, nb_comp, format, appended, type);
  
        // And the data attached to the points
  
        QDomNodeList pointdata_elmts = piece.elementsByTagName("PointData");
        if(pointdata_elmts.length() > 1)
          throw QString("loadMesh::Error:The piece cannot have more than ont PointData tag.");
  
        QDomElement pointdata_elm = pointdata_elmts.at(0).toElement();
        QString normal_array = pointdata_elm.attribute("Normals", "Normals");
        QDomElement childElement = pointdata_elm.firstChildElement();
        while(!childElement.isNull()) {
          SETSTATUS("Found point element " << childElement.tagName() << " named "
                                           << childElement.attribute("Name"));
          if(childElement.tagName() == "DataArray") {
            QString array_name = childElement.attribute("Name", "##Not an Array##");
            if(array_name == "Signal" or array_name == "Color" or array_name == "Label"
               or array_name == normal_array) {
              bool ok;
              QString nb_comp_s = childElement.attribute("NumberOfComponents", "1");
              size_t nb_comp = nb_comp_s.toUInt(&ok);
              if(!ok)
                throw QString("loadMesh::Error::NumberOfComponents attribute '%1' is not an integer")
                      .arg(nb_comp_s);
              FieldFormat format = FieldFormat::fromElement(childElement);
              Types type = Types::fromString(childElement.attribute("type"));
              if(array_name == "Signal" and nb_comp == 1)
                vtkDecode(points_signal, childElement, 1, format, appended, type);
              else if(array_name == "Color" and nb_comp == 1)
                vtkDecode(points_color, childElement, 1, format, appended, type);
              else if(array_name == "Label" and nb_comp == 1)
                vtkDecode(points_label, childElement, 1, format, appended, type);
              else if(array_name == normal_array and nb_comp == 3) {
                has_normals = true;
                vtkDecode(points_normal, childElement, 3, format, appended, type);
              }
            }
          }
          childElement = childElement.nextSiblingElement();
        }
  
        if(!progressAdvance(++cur_phase))
          userCancel();
  
        SETSTATUS("Read cells and cells data");
  
        // Now, read the cells
        QDomNodeList cells_elmts = piece.elementsByTagName("Cells");
        if(points_elmts.length() != 1)
          throw QString("loadMesh::Error:Piece needs exactly one list of cells");
        QDomElement cells_elm = cells_elmts.at(0).toElement();
  
        QDomNodeList cells_desc_elmts = cells_elm.elementsByTagName("DataArray");
        bool has_connectivity = false, has_offsets = false, has_types = false;
        QDomElement cells_desc = cells_elm.firstChildElement();
        while(!cells_desc.isNull()) {
          if(cells_desc.tagName() != "DataArray")
            throw QString("loadMesh::Error:Cells tag can only contain DataArray tags, not '%1'")
                  .arg(cells_desc.tagName());
          QString attribute_name = cells_desc.attribute("Name");
          if(attribute_name != "connectivity" and attribute_name != "offsets" and attribute_name != "types")
            throw QString("loadMesh::Error:Data arrays of Cells must be one of 'connectivity', 'offsets' or "
                          "'types', not '%1'").arg(attribute_name);
          FieldFormat format = FieldFormat::fromElement(cells_desc);
          Types type = Types::fromString(cells_desc.attribute("type"));
          if(attribute_name == "connectivity") {
            if(has_connectivity)
              throw QString("loadMesh::Error:connectivity of Cells defined more than once");
            vtkDecode(cells_connectivity, cells_desc, 1, format, appended, type);
          } else if(attribute_name == "offsets") {
            if(has_offsets)
              throw QString("loadMesh::Error:offset of Cells defined more than once");
            vtkDecode(cells_offsets, cells_desc, 1, format, appended, type);
          } else if(attribute_name == "types") {
            if(has_types)
              throw QString("loadMesh::Error:types of Cells defined more than once");
            vtkDecode(cells_types, cells_desc, 1, format, appended, type);
          }
          cells_desc = cells_desc.nextSiblingElement();
        }
        if(!progressAdvance(++cur_phase))
          userCancel();
  
        SETSTATUS("Create the graph");
  
        // Now, create the vertices
        std::vector<vertex> vertices(nb_points, vertex(0));
        for(int i = 0; i < nb_points; ++i) {
          vertex v;
          v->signal = points_signal[i][0];
          v->pos = Point3d(scale_factor * realPos(points_data[i], scale, transform, mesh));
          v->label = points_label[i][0];
          if(has_normals)
            v->nrml = Point3d(points_normal[i]);
          if(v->label > mesh->viewLabel())
            mesh->setLabel(v->label);
          v->saveId = i;
          vertices[i] = v;
          S.insert(v);
        }
  
        std::list<Point3v> to_process;
        // Now, create the cells
        size_t idx = 0;
        for(int i = 0; i < nb_cells; ++i) {
          std::vector<vertex> connectivity(cells_offsets[i][0] - idx);
          for(int j = 0; idx < cells_offsets[i][0]; ++idx, ++j) {
            size_t vid = cells_connectivity[idx][0];
            if(vid >= vertices.size())
              throw QString("loadMesh::Error:invalid vertex id '%1'").arg(vid);
            connectivity[j] = vertices[vid];
          }
          switch(cells_types[i][0]) {
          case 5:         // Triangle
          {
            if(connectivity.size() != 3)
              throw QString("loadMesh::Error:triangle with %1 sides").arg(connectivity.size());
            addTriangle(S, connectivity[0], connectivity[1], connectivity[2], to_process);
          } break;
          case 6:         // Triangle strip
          {
            if(connectivity.size() < 3)
              throw QString("loadMesh::Error:triangle strip with %1 sides").arg(connectivity.size());
            for(size_t i = 2; i < connectivity.size(); ++i) {
              if(i % 2 == 0)
                addTriangle(S, connectivity[i - 2], connectivity[i - 1], connectivity[i], to_process);
              else
                addTriangle(S, connectivity[i - 2], connectivity[i], connectivity[i - 1], to_process);
            }
          } break;
          case 8:         // 'pixel', i.e. quad but as '0,1,3,2'
          {
            if(connectivity.size() != 4)
              throw QString("loadMesh::Error:pixel with %1 sides").arg(connectivity.size());
            double d1 = norm(connectivity[3]->pos - connectivity[0]->pos);
            double d2 = norm(connectivity[2]->pos - connectivity[1]->pos);
            if(d2 < d1) {
              addTriangle(S, connectivity[0], connectivity[1], connectivity[2], to_process);
              addTriangle(S, connectivity[1], connectivity[3], connectivity[2], to_process);
            } else {
              addTriangle(S, connectivity[0], connectivity[1], connectivity[3], to_process);
              addTriangle(S, connectivity[0], connectivity[3], connectivity[2], to_process);
            }
          } break;
          case 9:         // quad
          {
            if(connectivity.size() != 4)
              throw QString("loadMesh::Error:quad with %1 sides").arg(connectivity.size());
            double d1 = norm(connectivity[2]->pos - connectivity[0]->pos);
            double d2 = norm(connectivity[3]->pos - connectivity[1]->pos);
            if(d2 < d1) {
              addTriangle(S, connectivity[0], connectivity[1], connectivity[3], to_process);
              addTriangle(S, connectivity[1], connectivity[2], connectivity[3], to_process);
            } else {
              addTriangle(S, connectivity[0], connectivity[1], connectivity[2], to_process);
              addTriangle(S, connectivity[0], connectivity[2], connectivity[3], to_process);
            }
          } break;
          default:
            break;
          }
        }
  
        while(!to_process.empty()) {
          bool advanced = false;
          std::list<Point3v> processing;
          std::swap(processing, to_process);
          forall(const Point3v& tr, processing) {
            if(addTriangle(S, tr[0], tr[1], tr[2], to_process))
              advanced = true;
          }
          if(!advanced)
            throw QString("loadMesh::Error:cannot recontruct mesh from triangulation. Make sure triangulation "
                          "is of a 2-Manifold.");
        }
  
        if(!has_normals)
          setNormals(mesh->graph());
  
        if(!progressAdvance(++cur_phase))
          userCancel();
      }
    }
  
    // Now, remove empty vertices
    std::list<vertex> to_remove;
    forall(const vertex& v, S) {
      if(S.empty(v))
        to_remove.push_back(v);
    }
    forall(const vertex& v, to_remove)
      S.erase(v);
  
    if(S.empty())
      throw QString("loadMesh::Error:This file doesn't contain any valid surface");
  
    mesh->setMeshType("MGXM");
    mesh->clearImgTex();
  
    SETSTATUS("Created mesh with " << mesh->graph().size() << " vertices");
  
    return true;
  }
}
