Utility class to convert between anatomical coordinate systems

vtkAnatomicalOrientation.cpp

#include <sstream>
#include <algorithm>

// Utility methods to be able to convert between various anatomical coordinate systems
namespace vtkAnatomicalOrientation {
  enum class Axis { L, R, P, A, S, I, None };
  static const Axis ValidAxes[6] = { Axis::L, Axis::R, Axis::P, Axis::A, Axis::S, Axis::I };

  static std::string AxisToString(Axis axis) {
    switch (axis) {
      case Axis::L: return "Left";
      case Axis::R: return "Right";
      case Axis::P: return "Posterior";
      case Axis::A: return "Anterior";
      case Axis::S: return "Superior";
      case Axis::I: return "Inferior";
      case Axis::None: return "";
    }
  }

  static std::string AxisToLowercaseString(Axis axis) {
    std::string name = AxisToString(axis);
    std::transform(name.begin(), name.end(), name.begin(), ::tolower);
    return name;
  }

  static Axis AxisFromString(std::string name) {
    for (const Axis& axis : ValidAxes) {
      if (name == AxisToString(axis) ||
          name == AxisToLowercaseString(axis))
        return axis;
    }
    return Axis::None;
  }

  static const char AxisToChar(Axis dir) {
    switch (dir) {
      case Axis::L: return 'L';
      case Axis::R: return 'R';
      case Axis::P: return 'P';
      case Axis::A: return 'A';
      case Axis::S: return 'S';
      case Axis::I: return 'I';
      case Axis::None: return 0;
    }
  }

  static Axis AxisFromChar(const char letter) {
    for (const Axis& axis : ValidAxes) {
      if (letter == AxisToChar(axis))
        return axis;
    }
    return Axis::None;
  }

  static Axis AxisInverse(Axis dir) {
    switch (dir) {
      case Axis::L: return Axis::R;
      case Axis::R: return Axis::L;
      case Axis::P: return Axis::A;
      case Axis::A: return Axis::P;
      case Axis::S: return Axis::I;
      case Axis::I: return Axis::S;
      case Axis::None: return Axis::None;
    }
  }

  struct CoordinateSystem {
    Axis X, Y, Z;

    bool operator==(const CoordinateSystem& rhs) {
      return X == rhs.X && Y == rhs.X && Z == rhs.Z;
    }

    bool operator!=(const CoordinateSystem& rhs) {
      return !(*this == rhs);
    }

    Axis& operator[](int i) {
      if (i == 0) return X;
      if (i == 1) return Y;
      if (i == 2) return Z;
      else {
        cerr << "Out of bounds of the coordinate system" << endl;
        exit(0);
      }
    }

    std::string GetAsAcronym() {
      std::string acronym = "";
      acronym += AxisToChar(X);
      acronym += AxisToChar(Y);
      acronym += AxisToChar(Z);
      return acronym;
    }

    void SetForAcronym(const std::string acronym) {
      try {
        X = AxisFromChar(acronym.at(0));
      } catch (const std::exception&) {
        X = Axis::None;
      }
      try {
        Y = AxisFromChar(acronym.at(1));
      } catch (const std::exception&) {
        Y = Axis::None;
      }
      try {
        Z = AxisFromChar(acronym.at(2));
      } catch (const std::exception&) {
        Z = Axis::None;
      }
    }

    std::string GetAsString(std::string separator) {
      std::string XStr = AxisToLowercaseString(X);
      std::string YStr = AxisToLowercaseString(Y);
      std::string ZStr = AxisToLowercaseString(Z);
      std::string str = XStr;
      if (YStr.length() > 0) {
        if (str.length() > 0) {
          str += separator;
        }
        str += YStr;
      }
      if (ZStr.length() > 0) {
        if (str.length() > 0) {
          str += separator;
        }
        str += ZStr;
      }
      return str;
    }

    void SetForString(const std::string str, const char delim) {
      std::string XYZStr[3];
      int i = 0;
      std::stringstream ss(str);
      std::string token;
      while (std::getline(ss, token, delim) && i < 3) {
        XYZStr[i] = token;
        i++;
      }
      X = AxisFromString(XYZStr[0]);
      Y = AxisFromString(XYZStr[1]);
      Z = AxisFromString(XYZStr[2]);
    }

    void GetTransformTo(CoordinateSystem newSystem, double* transform) {
      for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
          if ((*this)[i] == newSystem[j])
            transform[i*3+j] = 1;
          else if ((*this)[i] == AxisInverse(newSystem[j]))
            transform[i*3+j] = -1;
          else
            transform[i*3+j] = 0;
        }
      }
    }

    CoordinateSystem() : X(Axis::None), Y(Axis::None), Z(Axis::None) {}
    CoordinateSystem(Axis X, Axis Y, Axis Z) : X(X), Y(Y), Z(Z) {}
    CoordinateSystem(std::string str) {
      SetForAcronym(str);
      if (X == Axis::None || Y == Axis::None) {
        SetForString(str, '-');
      }
    }
  };

  static const CoordinateSystem LPS("LPS");
};

// Example: transformation from ILA to LPS
{
  vtkAnatomicalOrientation::CoordinateSystem inputCoordinateSystem("ILA");
  double ILAtoLPS[9]={0.0};
  inputCoordinateSystem.GetTransformTo(vtkAnatomicalOrientation::LPS, ILAtoLPS);
  // Result: ILAtoLPS = [0, 0, -1, 1, 0, 0, 0, -1, 0]
}