diff --git a/CMakeLists.txt b/CMakeLists.txt
index 221218bc73..7401d833fd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -118,6 +118,7 @@ set_property(GLOBAL PROPERTY USE_FOLDERS ON)
 set(FETCH_FILE_PATH "${simplnx_BINARY_DIR}/simplnx_fetch_remote_files.cmake")
 set_property(GLOBAL PROPERTY FETCH_FILE_PATH "${FETCH_FILE_PATH}")
 file(WRITE ${FETCH_FILE_PATH} "# -----------------------------------------------------------------------
+# This file is programmatically generated from the CMake file ${CMAKE_CURRENT_LIST_FILE}
 # This file has the commands to download each of the test files.
 # The WORKING_DIRECTORY is set to the following CMake code:
 # ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}/$<$<CONFIG:Debug>:Debug>$<$<CONFIG:Release>:Release>
diff --git a/src/Plugins/OrientationAnalysis/CMakeLists.txt b/src/Plugins/OrientationAnalysis/CMakeLists.txt
index 8a61f045ec..f2d6825b79 100644
--- a/src/Plugins/OrientationAnalysis/CMakeLists.txt
+++ b/src/Plugins/OrientationAnalysis/CMakeLists.txt
@@ -243,6 +243,9 @@ set(PLUGIN_EXTRA_SOURCES
   "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/SIMPLConversion.hpp"
   "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/TiffWriter.hpp"
   "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/TiffWriter.cpp"
+  "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/delaunator.cpp"
+  "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/delaunator.h"
+  "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities/IntersectionUtilities.hpp"
 )
 target_sources(${PLUGIN_NAME} PRIVATE ${PLUGIN_EXTRA_SOURCES})
 source_group(TREE "${${PLUGIN_NAME}_SOURCE_DIR}/src/${PLUGIN_NAME}/utilities" PREFIX ${PLUGIN_NAME} FILES ${PLUGIN_EXTRA_SOURCES})
diff --git a/src/Plugins/OrientationAnalysis/docs/WritePoleFigureFilter.md b/src/Plugins/OrientationAnalysis/docs/WritePoleFigureFilter.md
index 366d855b27..46d15cb023 100644
--- a/src/Plugins/OrientationAnalysis/docs/WritePoleFigureFilter.md
+++ b/src/Plugins/OrientationAnalysis/docs/WritePoleFigureFilter.md
@@ -6,9 +6,9 @@ IO (Output)
 
 ## Description
 
-This **Filter** creates a standard crystallographic pole figure image for each **Ensemble** (phase)in a selected **Data Container**. The **Filter** uses Euler angles in radians and requires the crystal structures and material names for each **Ensemble** array and the corresponding **Ensemble** Ids on the **Cells**. The **Filter** also optionally can use a *mask* array to determine which **Cells** are valid for the pole figure computation.
+This **Filter** creates a standard crystallographic pole figure image for each **Ensemble** (phase) in a selected **Data Container**. The **Filter** uses Euler angles in radians and requires the crystal structures and material names for each **Ensemble** array and the corresponding **Ensemble** Ids on the **Cells**. The **Filter** also optionally can use a *mask* array to determine which angles are valid for the pole figure computation.
 
-In a practicale sense, this means that the following information is available to the filter:
+In a practical sense, this means that the following information is available to the filter:
 
 - Cell Level
 
@@ -23,17 +23,36 @@ In a practicale sense, this means that the following information is available to
 
 ### Algorithm Choice
 
-1: The pole figure algorithm uses a *modified Lambert square* to perform the interpolations onto the unit circle. This is an alternate type of interpolation that the EBSD OEMs do not perform which may make the output from DREAM.3D look slightly different than output obtained from the OEM programs.
+1: The pole figure algorithm uses a *modified Lambert square* to perform the interpolations onto the unit circle. This is an alternate type
+of interpolation that the EBSD OEMs do not perform which may make the output from DREAM.3D look slightly different than output
+obtained from the OEM programs.
 
-**Only an advanced user with intimate knowledge of the modified Lambert projection should attempt to change the value for the "Lambert Image Size (Pixels)" input parameter.**
+**Only an advanced user with intimate knowledge of the modified Lambert projection should attempt to change the value for
+the "Lambert Image Size (Pixels)" input parameter.**
 
 2: Discrete Pole figure. The algorithm will simply mark each pixel that had at least 1 count as a black pixel.
 
-### Layout
+## Output Options
+
+### Write Image to Disk
+
+The user can select to have the combined set of pole figures written to disk as a tiff image file
+
+### Save Pole Figure as Image Geometry
+
+The combined pole figure image will be saved into the DataStructure as an Image Geometry
+
+### Save Raw Intensity Data
+
+The normalized count data is saved for each pole figure into a Data Array that is stored inside an Image Geometry. This allows
+the user to select their own color plots. The Image Geometry will also have a string DataArray that lists the pertinent
+data that went into the creation: Number of points, which hemisphere, Phase Name, etc.
+
+### Image Layout
 
 The 3 pole figures can be laid out in a Square, Horizontal row or vertical column. Supporting information (including the color bar legend for color pole figures) will also be printed on the image.
 
-| Lambert Projection | Discrete |
+| Colorized Intensity | Discrete |
 |--------------------|----------|
 | ![Example Pole Figure Using Square Layout](Images/PoleFigure_Example.png) | ![Example Pole Figure Using Square Layout](Images/Pole_Figure_Discrete_Example.png) |
 
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.cpp
index 72c9b5625b..72e67555e2 100644
--- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.cpp
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.cpp
@@ -2,16 +2,22 @@
 
 #include "OrientationAnalysis/utilities/FiraSansRegular.hpp"
 #include "OrientationAnalysis/utilities/Fonts.hpp"
+#include "OrientationAnalysis/utilities/IntersectionUtilities.hpp"
 #include "OrientationAnalysis/utilities/LatoBold.hpp"
 #include "OrientationAnalysis/utilities/LatoRegular.hpp"
 #include "OrientationAnalysis/utilities/TiffWriter.hpp"
+#include "OrientationAnalysis/utilities/delaunator.h"
 
 #include "simplnx/Common/Constants.hpp"
 #include "simplnx/Common/RgbColor.hpp"
 #include "simplnx/DataStructure/DataArray.hpp"
 #include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
+#include "simplnx/DataStructure/Geometry/TriangleGeom.hpp"
+#include "simplnx/Pipeline/Pipeline.hpp"
 #include "simplnx/Utilities/DataArrayUtilities.hpp"
 #include "simplnx/Utilities/ParallelTaskAlgorithm.hpp"
+#include "simplnx/Utilities/Parsing/DREAM3D/Dream3dIO.hpp"
+#include "simplnx/Utilities/RTree.hpp"
 #include "simplnx/Utilities/StringUtilities.hpp"
 
 #include "EbsdLib/Core/EbsdLibConstants.h"
@@ -26,7 +32,12 @@
 #include "EbsdLib/LaueOps/TriclinicOps.h"
 #include "EbsdLib/LaueOps/TrigonalLowOps.h"
 #include "EbsdLib/LaueOps/TrigonalOps.h"
-#include "EbsdLib/Utilities/ComputeStereographicProjection.h"
+#include "EbsdLib/Utilities/LambertUtilities.h"
+#include "EbsdLib/Utilities/ModifiedLambertProjection.h"
+
+#include "H5Support/H5Lite.h"
+#include "H5Support/H5ScopedSentinel.h"
+#include "H5Support/H5Utilities.h"
 
 #define CANVAS_ITY_IMPLEMENTATION
 #include <canvas_ity.hpp>
@@ -37,6 +48,363 @@ namespace
 {
 const bool k_UseDiscreteHeatMap = false;
 
+class ComputeIntensityStereographicProjection
+{
+public:
+  ComputeIntensityStereographicProjection(EbsdLib::FloatArrayType* xyzCoords, PoleFigureConfiguration_t* config, EbsdLib::DoubleArrayType* intensity, bool normalizeToMRD)
+  : m_XYZCoords(xyzCoords)
+  , m_Config(config)
+  , m_Intensity(intensity)
+  , m_NormalizeToMRD(normalizeToMRD)
+  {
+  }
+
+  void operator()() const
+  {
+    m_Intensity->resizeTuples(m_Config->imageDim * m_Config->imageDim);
+    m_Intensity->initializeWithZeros();
+
+    if(m_Config->discrete)
+    {
+      int halfDim = m_Config->imageDim / 2;
+      double* intensity = m_Intensity->getPointer(0);
+      size_t numCoords = m_XYZCoords->getNumberOfTuples();
+      float* xyzPtr = m_XYZCoords->getPointer(0);
+      for(size_t i = 0; i < numCoords; i++)
+      {
+        if(xyzPtr[i * 3 + 2] < 0.0f)
+        {
+          xyzPtr[i * 3 + 0] *= -1.0f;
+          xyzPtr[i * 3 + 1] *= -1.0f;
+          xyzPtr[i * 3 + 2] *= -1.0f;
+        }
+        float x = xyzPtr[i * 3] / (1 + xyzPtr[i * 3 + 2]);
+        float y = xyzPtr[i * 3 + 1] / (1 + xyzPtr[i * 3 + 2]);
+
+        int xCoord = static_cast<int>(x * static_cast<float32>(halfDim - 1)) + halfDim;
+        int yCoord = static_cast<int>(y * static_cast<float32>(halfDim - 1)) + halfDim;
+
+        size_t index = (yCoord * m_Config->imageDim) + xCoord;
+
+        intensity[index]++;
+      }
+    }
+    else
+    {
+      ModifiedLambertProjection::Pointer lambert = ModifiedLambertProjection::LambertBallToSquare(m_XYZCoords, m_Config->lambertDim, m_Config->sphereRadius);
+      if(m_NormalizeToMRD)
+      {
+        lambert->normalizeSquaresToMRD();
+      }
+      lambert->createStereographicProjection(m_Config->imageDim, *m_Intensity);
+
+      // This next function (writeLambertData) is experimental but should be left in to
+      // make sure it will still compile. At some point I will get back to the development
+      // of this alternate pole figure generation.
+      // writeLambertData(lambert);
+    }
+  }
+
+  template <typename V, typename T, typename W>
+  void unstructuredGridInterpolator(nx::core::IFilter* filter, nx::core::TriangleGeom* delaunayGeom, std::vector<V>& xPositionsPtr, std::vector<V>& yPositionsPtr, T* xyValues,
+                                    typename std::vector<W>& outputValues) const
+  {
+    using Vec3f = nx::core::Vec3<float>;
+    using RTreeType = RTree<size_t, float, 2, float>;
+
+    // filter->notifyStatusMessage(QString("Starting Interpolation...."));
+    nx::core::IGeometry::SharedFaceList& delTriangles = delaunayGeom->getFacesRef();
+    size_t numTriangles = delaunayGeom->getNumberOfFaces();
+    int percent = 0;
+    int counter = xPositionsPtr.size() / 100;
+    RTreeType m_RTree;
+    // Populate the RTree
+
+    size_t numTris = delaunayGeom->getNumberOfFaces();
+    for(size_t tIndex = 0; tIndex < numTris; tIndex++)
+    {
+      std::array<float, 6> boundBox = nx::IntersectionUtilities::GetBoundingBoxAtTri(*delaunayGeom, tIndex);
+      m_RTree.Insert(boundBox.data(), boundBox.data() + 3, tIndex); // Note, all values including zero are fine in this version
+    }
+
+    for(size_t vertIndex = 0; vertIndex < xPositionsPtr.size(); vertIndex++)
+    {
+      Vec3f rayOrigin(xPositionsPtr[vertIndex], yPositionsPtr[vertIndex], 1.0F);
+      Vec3f rayDirection(0.0F, 0.0F, -1.0F);
+      Vec3f barycentricCoord(0.0F, 0.0F, 0.0F);
+      //      int xPos = xPositionsPtr[vertIndex];
+      //      int yPos = yPositionsPtr[vertIndex];
+      if(counter != 0)
+      {
+
+        if(vertIndex % counter == 0)
+        {
+          //          QString ss = QObject::tr("Interpolating || %1% Complete").arg(percent);
+          //          filter->notifyStatusMessage(ss);
+          percent += 1;
+        }
+      }
+
+      // Create these reusable variables to save the reallocation each time through the loop
+
+      std::vector<size_t> hitTriangleIds;
+      std::function<bool(size_t)> func = [&](size_t id) {
+        hitTriangleIds.push_back(id);
+        return true; // keep going
+      };
+
+      int nhits = m_RTree.Search(rayOrigin.data(), rayOrigin.data(), func);
+      for(auto triIndex : hitTriangleIds)
+      {
+        barycentricCoord = {0.0F, 0.0F, 0.0F};
+        std::array<size_t, 3> triVertIndices;
+        // Get the Vertex Coordinates for each of the 3 vertices
+        std::array<nx::core::Point3Df, 3> verts;
+        delaunayGeom->getFaceCoordinates(triIndex, verts);
+        Vec3f v0 = verts[0];
+        Vec3f v1 = verts[1];
+        Vec3f v2 = verts[2];
+
+        // Get the vertex Indices from the triangle
+        delaunayGeom->getFacePointIds(triIndex, triVertIndices);
+        bool inTriangle = nx::IntersectionUtilities::RayTriangleIntersect2(rayOrigin, rayDirection, v0, v1, v2, barycentricCoord);
+        if(inTriangle)
+        {
+          // Linear Interpolate dx and dy values using the barycentric coordinates
+          delaunayGeom->getFaceCoordinates(triIndex, verts);
+          float f0 = xyValues[triVertIndices[0]];
+          float f1 = xyValues[triVertIndices[1]];
+          float f2 = xyValues[triVertIndices[2]];
+
+          float interpolatedVal = (barycentricCoord[0] * f0) + (barycentricCoord[1] * f1) + (barycentricCoord[2] * f2);
+
+          outputValues[vertIndex] = interpolatedVal;
+
+          break;
+        }
+      }
+    }
+  }
+
+  int writeLambertData(ModifiedLambertProjection::Pointer lambert) const
+  {
+    int err = -1;
+
+    int m_Dimension = lambert->getDimension();
+    EbsdLib::DoubleArrayType::Pointer m_NorthSquare = lambert->getNorthSquare();
+    EbsdLib::DoubleArrayType::Pointer m_SouthSquare = lambert->getSouthSquare();
+
+    // We want half the sphere area for each square because each square represents a hemisphere.
+    const float32 sphereRadius = 1.0f;
+    float halfSphereArea = 4.0f * EbsdLib::Constants::k_PiF * sphereRadius * sphereRadius / 2.0f;
+    // The length of a side of the square is the square root of the area
+    float squareEdge = std::sqrt(halfSphereArea);
+    float32 m_StepSize = squareEdge / static_cast<float>(m_Dimension);
+
+    float32 m_MaxCoord = squareEdge / 2.0f;
+    float32 m_MinCoord = -squareEdge / 2.0f;
+    std::array<float, 3> vert = {0.0f, 0.0f, 0.0f};
+
+    std::vector<float> squareCoords(m_Dimension * m_Dimension * 3);
+
+    // Northern Hemisphere Coordinates
+    std::vector<float> northSphereCoords(m_Dimension * m_Dimension * 3);
+    std::vector<float> northStereoCoords(m_Dimension * m_Dimension * 3);
+
+    // Southern Hemisphere Coordinates
+    std::vector<float> southSphereCoords(m_Dimension * m_Dimension * 3);
+    std::vector<float> southStereoCoords(m_Dimension * m_Dimension * 3);
+
+    size_t index = 0;
+
+    const float origin = m_MinCoord + (m_StepSize / 2.0f);
+    for(int32_t y = 0; y < m_Dimension; ++y)
+    {
+      for(int x = 0; x < m_Dimension; ++x)
+      {
+        vert[0] = origin + (static_cast<float>(x) * m_StepSize);
+        vert[1] = origin + (static_cast<float>(y) * m_StepSize);
+
+        squareCoords[index * 3] = vert[0];
+        squareCoords[index * 3 + 1] = vert[1];
+        squareCoords[index * 3 + 2] = 0.0f;
+
+        LambertUtilities::LambertSquareVertToSphereVert(vert.data(), LambertUtilities::Hemisphere::North);
+
+        northSphereCoords[index * 3] = vert[0];
+        northSphereCoords[index * 3 + 1] = vert[1];
+        northSphereCoords[index * 3 + 2] = vert[2];
+
+        northStereoCoords[index * 3] = vert[0] / (1.0f + vert[2]);
+        northStereoCoords[index * 3 + 1] = vert[1] / (1.0f + vert[2]);
+        northStereoCoords[index * 3 + 2] = 0.0f;
+
+        // Reset the Lambert Square Coord
+        vert[0] = origin + (static_cast<float>(x) * m_StepSize);
+        vert[1] = origin + (static_cast<float>(y) * m_StepSize);
+        LambertUtilities::LambertSquareVertToSphereVert(vert.data(), LambertUtilities::Hemisphere::South);
+
+        southSphereCoords[index * 3] = vert[0];
+        southSphereCoords[index * 3 + 1] = vert[1];
+        southSphereCoords[index * 3 + 2] = vert[2];
+
+        southStereoCoords[index * 3] = vert[0] / (1.0f + vert[2]);
+        southStereoCoords[index * 3 + 1] = vert[1] / (1.0f + vert[2]);
+        southStereoCoords[index * 3 + 2] = 0.0f;
+
+        index++;
+      }
+    }
+
+    //**********************************************************************************************************************
+    // Triangulate the stereo coordinates
+
+    DataStructure dataStructure;
+
+    auto* triangleGeomPtr = TriangleGeom::Create(dataStructure, "Delaunay");
+    auto& triangleGeom = *triangleGeomPtr;
+
+    DataPath sharedFaceListPath({"Delaunay", "SharedFaceList"});
+
+    usize numPts = northStereoCoords.size() / 3;
+    // Create the default DataArray that will hold the FaceList and Vertices. We
+    // size these to 1 because the Csv parser will resize them to the appropriate number of tuples
+    using DimensionType = std::vector<size_t>;
+
+    DimensionType faceTupleShape = {0};
+    Result result = CreateArray<IGeometry::MeshIndexType>(dataStructure, faceTupleShape, {3ULL}, sharedFaceListPath, IDataAction::Mode::Execute);
+    if(result.invalid())
+    {
+      return -1;
+      // return MergeResults(result, MakeErrorResult(-5509, fmt::format("{}CreateGeometry2DAction: Could not allocate SharedTriList '{}'", prefix, trianglesPath.toString())));
+    }
+    auto& sharedFaceListRef = dataStructure.getDataRefAs<IGeometry::MeshIndexArrayType>(sharedFaceListPath);
+    triangleGeom.setFaceList(sharedFaceListRef);
+
+    // Create the Vertex Array with a component size of 3
+    DataPath vertexPath({"Delaunay", "SharedVertexList"});
+
+    DimensionType vertexTupleShape = {0};
+    result = CreateArray<float>(dataStructure, vertexTupleShape, {3}, vertexPath, IDataAction::Mode::Execute);
+    if(result.invalid())
+    {
+      return -2;
+      // return MergeResults(result, MakeErrorResult(-5510, fmt::format("{}CreateGeometry2DAction: Could not allocate SharedVertList '{}'", prefix, vertexPath.toString())));
+    }
+    Float32Array* vertexArray = ArrayFromPath<float>(dataStructure, vertexPath);
+    triangleGeom.setVertices(*vertexArray);
+    triangleGeom.resizeVertexList(numPts);
+
+    auto* vertexAttributeMatrix = AttributeMatrix::Create(dataStructure, "VertexData", {numPts}, triangleGeom.getId());
+    if(vertexAttributeMatrix == nullptr)
+    {
+      return -3;
+      // return MakeErrorResult(-5512, fmt::format("CreateGeometry2DAction: Failed to create attribute matrix: '{}'", prefix, vertexDataPath.toString()));
+    }
+    triangleGeom.setVertexAttributeMatrix(*vertexAttributeMatrix);
+
+    if(nullptr != triangleGeom.getVertexAttributeMatrix())
+    {
+      triangleGeom.getVertexAttributeMatrix()->resizeTuples({numPts});
+    }
+    auto vertexCoordsPtr = triangleGeom.getVerticesRef();
+
+    // Create Coords for the Delaunator Algorithm
+    std::vector<float64> coords(2 * numPts, 0.0);
+    for(usize i = 0; i < numPts; i++)
+    {
+      coords[i * 2] = northStereoCoords[i * 3];
+      coords[i * 2 + 1] = northStereoCoords[i * 3 + 1];
+      vertexCoordsPtr[i * 3] = northStereoCoords[i * 3];
+      vertexCoordsPtr[i * 3 + 1] = northStereoCoords[i * 3 + 1];
+      vertexCoordsPtr[i * 3 + 2] = northStereoCoords[i * 3 + 2];
+    }
+
+    // Perform the triangulation
+    nx::delaunator::Delaunator d(coords);
+
+    usize numTriangles = d.triangles.size();
+    triangleGeom.resizeFaceList(numTriangles / 3);
+    auto sharedTriListPtr = triangleGeom.getFacesRef();
+    // usize triangleIndex = 0;
+    for(usize i = 0; i < numTriangles; i += 3)
+    {
+      std::array<usize, 3> triIDs = {d.triangles[i], d.triangles[i + 1], d.triangles[i + 2]};
+      sharedTriListPtr[i] = d.triangles[i];
+      sharedTriListPtr[i + 1] = d.triangles[i + 1];
+      sharedTriListPtr[i + 2] = d.triangles[i + 2];
+    }
+
+    // Create the vertex and face AttributeMatrix
+    auto* faceAttributeMatrix = AttributeMatrix::Create(dataStructure, "Face Data", {numTriangles / 3}, triangleGeom.getId());
+    if(faceAttributeMatrix == nullptr)
+    {
+      return -4;
+      // return MakeErrorResult(-5511, fmt::format("{}CreateGeometry2DAction: Failed to create attribute matrix: '{}'", prefix, faceDataPath.toString()));
+    }
+    triangleGeom.setFaceAttributeMatrix(*faceAttributeMatrix);
+
+    Pipeline pipeline;
+    const Result<> result2 = DREAM3D::WriteFile(fmt::format("/tmp/delaunay_triangulation.dream3d"), dataStructure, pipeline, true);
+    //******************************************************************************************************************************
+
+    //******************************************************************************************************************************
+    // Perform a Bi-linear Interpolation
+    // Generate a regular grid of XY points
+    size_t numSteps = 1024;
+    float32 stepInc = 2.0f / static_cast<float>(numSteps);
+    std::vector<float32> xcoords(numSteps * numSteps);
+    std::vector<float32> ycoords(numSteps * numSteps);
+    for(size_t y = 0; y < numSteps; ++y)
+    {
+      for(size_t x = 0; x < numSteps; x++)
+      {
+        size_t idx = y * numSteps + x;
+        xcoords[idx] = -1.0f + static_cast<float32>(x) * stepInc;
+        ycoords[idx] = -1.0f + static_cast<float32>(y) * stepInc;
+      }
+    }
+
+    std::vector<double> outputValues(numSteps * numSteps);
+    unstructuredGridInterpolator<float32, double, double>(nullptr, &triangleGeom, xcoords, ycoords, m_NorthSquare->data(), outputValues);
+
+    //******************************************************************************************************************************
+
+    //******************************************************************************************************************************
+    // Write out all the data
+    {
+      hid_t groupId = H5Support::H5Utilities::createFile("/tmp/lambert_data.h5");
+      H5Support::H5ScopedFileSentinel fileSentinel(groupId, false);
+
+      std::vector<hsize_t> dims = {static_cast<hsize_t>(m_Dimension), static_cast<hsize_t>(m_Dimension)};
+      err = H5Support::H5Lite::writePointerDataset(groupId, m_NorthSquare->getName(), 2, dims.data(), m_NorthSquare->data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, m_SouthSquare->getName(), 2, dims.data(), m_SouthSquare->data());
+      dims[0] = m_Dimension * m_Dimension;
+      dims[1] = 3ULL;
+
+      err = H5Support::H5Lite::writePointerDataset(groupId, "Lambert Square Coords", 2, dims.data(), squareCoords.data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, "North Sphere Coords", 2, dims.data(), northSphereCoords.data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, "North Stereo Coords", 2, dims.data(), northStereoCoords.data());
+
+      err = H5Support::H5Lite::writePointerDataset(groupId, "South Sphere Coords", 2, dims.data(), southSphereCoords.data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, "South Stereo Coords", 2, dims.data(), southStereoCoords.data());
+
+      dims[0] = numSteps * numSteps;
+      err = H5Support::H5Lite::writePointerDataset(groupId, "X Coords", 1, dims.data(), xcoords.data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, "Y Coords", 1, dims.data(), ycoords.data());
+      err = H5Support::H5Lite::writePointerDataset(groupId, "Interpolated Values", 1, dims.data(), outputValues.data());
+    }
+
+    return err;
+  }
+
+private:
+  EbsdLib::FloatArrayType* m_XYZCoords = nullptr;
+  PoleFigureConfiguration_t* m_Config = nullptr;
+  EbsdLib::DoubleArrayType* m_Intensity = nullptr;
+  bool m_NormalizeToMRD = false;
+};
+
 // -----------------------------------------------------------------------------
 template <typename Ops>
 std::vector<EbsdLib::UInt8ArrayType::Pointer> makePoleFigures(PoleFigureConfiguration_t& config)
@@ -46,7 +414,7 @@ std::vector<EbsdLib::UInt8ArrayType::Pointer> makePoleFigures(PoleFigureConfigur
 }
 
 template <typename OpsType>
-std::vector<EbsdLib::DoubleArrayType ::Pointer> createIntensityPoleFigures(PoleFigureConfiguration_t& config)
+std::vector<EbsdLib::DoubleArrayType::Pointer> createIntensityPoleFigures(PoleFigureConfiguration_t& config, bool normalizeToMRD)
 {
   OpsType ops;
   std::string label0 = std::string("<001>");
@@ -92,69 +460,21 @@ std::vector<EbsdLib::DoubleArrayType ::Pointer> createIntensityPoleFigures(PoleF
   // Compute the Stereographic Data in parallel
   ParallelTaskAlgorithm taskRunner;
   taskRunner.setParallelizationEnabled(true);
-  taskRunner.execute(ops.ComputeStereographicProjection(xyz001.get(), &config, intensity001.get()));
-  taskRunner.execute(ComputeStereographicProjection(xyz011.get(), &config, intensity011.get()));
-  taskRunner.execute(ComputeStereographicProjection(xyz111.get(), &config, intensity111.get()));
+  taskRunner.execute(ComputeIntensityStereographicProjection(xyz001.get(), &config, intensity001.get(), normalizeToMRD));
+  taskRunner.execute(ComputeIntensityStereographicProjection(xyz011.get(), &config, intensity011.get(), normalizeToMRD));
+  taskRunner.execute(ComputeIntensityStereographicProjection(xyz111.get(), &config, intensity111.get(), normalizeToMRD));
   taskRunner.wait(); // This will spill over if the number of DataArrays to process does not divide evenly by the number of threads.
 
-  // Find the Max and Min values based on ALL 3 arrays. We can color scale them all the same
-  double max = std::numeric_limits<double>::min();
-  double min = std::numeric_limits<double>::max();
-
-  double* dPtr = intensity001->getPointer(0);
-  size_t count = intensity001->getNumberOfTuples();
-  for(size_t i = 0; i < count; ++i)
-  {
-    if(dPtr[i] > max)
-    {
-      max = dPtr[i];
-    }
-    if(dPtr[i] < min)
-    {
-      min = dPtr[i];
-    }
-  }
-
-  dPtr = intensity011->getPointer(0);
-  count = intensity011->getNumberOfTuples();
-  for(size_t i = 0; i < count; ++i)
-  {
-    if(dPtr[i] > max)
-    {
-      max = dPtr[i];
-    }
-    if(dPtr[i] < min)
-    {
-      min = dPtr[i];
-    }
-  }
-
-  dPtr = intensity111->getPointer(0);
-  count = intensity111->getNumberOfTuples();
-  for(size_t i = 0; i < count; ++i)
-  {
-    if(dPtr[i] > max)
-    {
-      max = dPtr[i];
-    }
-    if(dPtr[i] < min)
-    {
-      min = dPtr[i];
-    }
-  }
-
-  config.minScale = min;
-  config.maxScale = max;
-
   return {intensity001, intensity011, intensity111};
 }
 
 // -----------------------------------------------------------------------------
-EbsdLib::UInt8ArrayType::Pointer flipAndMirrorPoleFigure(EbsdLib::UInt8ArrayType* src, const PoleFigureConfiguration_t& config)
+template <typename T>
+typename EbsdDataArray<T>::Pointer flipAndMirrorPoleFigure(EbsdDataArray<T>* src, const PoleFigureConfiguration_t& config)
 {
-  EbsdLib::UInt8ArrayType::Pointer converted = EbsdLib::UInt8ArrayType::CreateArray(static_cast<size_t>(config.imageDim * config.imageDim), {4}, src->getName(), true);
+  typename EbsdDataArray<T>::Pointer converted = EbsdDataArray<T>::CreateArray(config.imageDim * config.imageDim, src->getComponentDimensions(), src->getName(), true);
   // We need to flip the image "vertically", which means the bottom row becomes
-  // the top row and convert from BGRA to RGB ordering (This is a Little Endian code)
+  // the top row and convert from BGRA to RGBA ordering (This is a Little Endian code)
   // If this is ever compiled on a BIG ENDIAN machine the colors will be off.
   for(int y = 0; y < config.imageDim; y++)
   {
@@ -164,18 +484,32 @@ EbsdLib::UInt8ArrayType::Pointer flipAndMirrorPoleFigure(EbsdLib::UInt8ArrayType
       const size_t indexSrc = y * config.imageDim + x;
       const size_t indexDest = destY * config.imageDim + x;
 
-      uint8_t* argbPtr = src->getTuplePointer(indexSrc);
-      uint8_t* destPtr = converted->getTuplePointer(indexDest);
-
-      destPtr[0] = argbPtr[2];
-      destPtr[1] = argbPtr[1];
-      destPtr[2] = argbPtr[0];
-      destPtr[3] = argbPtr[3];
+      T* argbPtr = src->getTuplePointer(indexSrc);
+      converted->setTuple(indexDest, argbPtr);
     }
   }
   return converted;
 }
 
+template <typename T>
+typename EbsdDataArray<T>::Pointer convertColorOrder(EbsdDataArray<T>* src, const PoleFigureConfiguration_t& config)
+{
+  typename EbsdDataArray<T>::Pointer converted = EbsdDataArray<T>::CreateArray(config.imageDim * config.imageDim, src->getComponentDimensions(), src->getName(), true);
+  // BGRA to RGBA ordering (This is a Little Endian code)
+  // If this is ever compiled on a BIG ENDIAN machine the colors will be off.
+  size_t numTuples = src->getNumberOfTuples();
+  for(size_t tIdx = 0; tIdx < numTuples; tIdx++)
+  {
+    T* argbPtr = src->getTuplePointer(tIdx);
+    T* destPtr = converted->getTuplePointer(tIdx);
+    destPtr[0] = argbPtr[2];
+    destPtr[1] = argbPtr[1];
+    destPtr[2] = argbPtr[0];
+    destPtr[3] = argbPtr[3];
+  }
+  return converted;
+}
+
 // -----------------------------------------------------------------------------
 void drawInformationBlock(canvas_ity::canvas& context, const PoleFigureConfiguration_t& config, const std::pair<float32, float32>& position, float margins, float fontPtSize, int32_t phaseNum,
                           std::vector<unsigned char>& fontData, const std::string& laueGroupName, const std::string& materialName)
@@ -187,7 +521,7 @@ void drawInformationBlock(canvas_ity::canvas& context, const PoleFigureConfigura
   const float colorHeight = (static_cast<float>(imageHeight)) / static_cast<float>(config.numColors);
   //
   using RectFType = std::pair<float, float>;
-  const RectFType rect = std::make_pair(imageWidth * scaleBarRelativeWidth, colorHeight * 1.00000f);
+  const RectFType rect = std::make_pair(static_cast<float>(imageWidth) * scaleBarRelativeWidth, colorHeight * 1.00000f);
   //
   const std::array<canvas_ity::baseline_style, 6> baselines = {canvas_ity::alphabetic, canvas_ity::top, canvas_ity::middle, canvas_ity::bottom, canvas_ity::hanging, canvas_ity::ideographic};
 
@@ -211,7 +545,7 @@ void drawInformationBlock(canvas_ity::canvas& context, const PoleFigureConfigura
     context.set_font(fontData.data(), static_cast<int>(fontData.size()), fontPtSize);
     context.set_color(canvas_ity::fill_style, 0.0f, 0.0f, 0.0f, 1.0f);
     context.text_baseline = baselines[0];
-    context.fill_text(label.c_str(), position.first + margins + rect.first + margins, position.second + margins + (imageHeight / 3.0) + (heightInc * fontPtSize));
+    context.fill_text(label.c_str(), position.first + margins + rect.first + margins, position.second + margins + (static_cast<float>(imageHeight) / 3.0f) + (heightInc * fontPtSize));
     context.close_path();
     heightInc++;
   }
@@ -236,7 +570,7 @@ void drawScalarBar(canvas_ity::canvas& context, const PoleFigureConfiguration_t&
     r = colors[3 * i];
     g = colors[3 * i + 1];
     b = colors[3 * i + 2];
-    colorTable[i] = RgbColor::dRgb(r * 255, g * 255, b * 255, 255);
+    colorTable[i] = RgbColor::dRgb(static_cast<uint8>(r * 255.0f), static_cast<uint8>(g * 255.0f), static_cast<uint8>(b * 255.0f), 255);
   }
 
   // Now start from the bottom and draw colored lines up the scale bar
@@ -249,7 +583,7 @@ void drawScalarBar(canvas_ity::canvas& context, const PoleFigureConfiguration_t&
 
   using RectFType = std::pair<float, float>;
 
-  const RectFType rect = std::make_pair(imageWidth * scaleBarRelativeWidth, colorHeight * 1.00000f);
+  const RectFType rect = std::make_pair(static_cast<float32>(imageWidth) * scaleBarRelativeWidth, colorHeight * 1.00000f);
 
   const std::array<canvas_ity::baseline_style, 6> baselines = {canvas_ity::alphabetic, canvas_ity::top, canvas_ity::middle, canvas_ity::bottom, canvas_ity::hanging, canvas_ity::ideographic};
 
@@ -268,7 +602,7 @@ void drawScalarBar(canvas_ity::canvas& context, const PoleFigureConfiguration_t&
   context.set_font(fontData.data(), static_cast<int>(fontData.size()), fontPtSize);
   context.set_color(canvas_ity::fill_style, 0.0f, 0.0f, 0.0f, 1.0f);
   context.text_baseline = baselines[0];
-  context.fill_text(minStr.c_str(), position.first + 2.0F * margins + rect.first, position.second + (2 * margins) + (2 * fontPtSize) + ((numColors)*colorHeight));
+  context.fill_text(minStr.c_str(), position.first + 2.0F * margins + rect.first, position.second + (2 * margins) + (2 * fontPtSize) + (static_cast<float32>(numColors) * colorHeight));
   context.close_path();
 
   // Draw the color bar
@@ -278,7 +612,7 @@ void drawScalarBar(canvas_ity::canvas& context, const PoleFigureConfiguration_t&
     std::tie(r, g, b) = RgbColor::fRgb(c);
 
     const float32 x = position.first + margins;
-    const float32 y = position.second + (2 * margins) + (2 * fontPtSize) + (i * colorHeight);
+    const float32 y = position.second + (2 * margins) + (2 * fontPtSize) + (static_cast<float32>(i) * colorHeight);
 
     context.begin_path();
     context.set_color(canvas_ity::fill_style, r, g, b, 1.0f);
@@ -336,11 +670,11 @@ Result<> WritePoleFigure::operator()()
 
   const std::vector<LaueOps::Pointer> orientationOps = LaueOps::GetAllOrientationOps();
 
-  const nx::core::Float32Array& eulers = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->CellEulerAnglesArrayPath);
-  nx::core::Int32Array& phases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->CellPhasesArrayPath);
+  const nx::core::Float32Array& eulerAngles = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->CellEulerAnglesArrayPath);
+  auto& phases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->CellPhasesArrayPath);
 
-  nx::core::UInt32Array& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
-  nx::core::StringArray& materialNames = m_DataStructure.getDataRefAs<StringArray>(m_InputValues->MaterialNameArrayPath);
+  auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
+  auto& materialNames = m_DataStructure.getDataRefAs<StringArray>(m_InputValues->MaterialNameArrayPath);
 
   std::unique_ptr<MaskCompare> maskCompare = nullptr;
   if(m_InputValues->UseMask)
@@ -358,7 +692,7 @@ Result<> WritePoleFigure::operator()()
 
   // Find the total number of angles we have based on the number of Tuples of the
   // Euler Angles array
-  const size_t numPoints = eulers.getNumberOfTuples();
+  const size_t numPoints = eulerAngles.getNumberOfTuples();
   // Find how many phases we have by getting the number of Crystal Structures
   const size_t numPhases = crystalStructures.getNumberOfTuples();
 
@@ -371,7 +705,7 @@ Result<> WritePoleFigure::operator()()
   imageGeom.setDimensions({static_cast<usize>(m_InputValues->ImageSize), static_cast<usize>(m_InputValues->ImageSize), 1});
   imageGeom.getCellData()->resizeTuples(tupleShape);
 
-  // Loop over all the voxels gathering the Eulers for a specific phase into an array
+  // Loop over all the voxels gathering the Euler angles for a specific phase into an array
   for(size_t phase = 1; phase < numPhases; ++phase)
   {
     size_t count = 0;
@@ -388,11 +722,11 @@ Result<> WritePoleFigure::operator()()
       }
     }
     const std::vector<size_t> eulerCompDim = {3};
-    const EbsdLib::FloatArrayType::Pointer subEulersPtr = EbsdLib::FloatArrayType::CreateArray(count, eulerCompDim, "Eulers_Per_Phase", true);
-    subEulersPtr->initializeWithValue(std::numeric_limits<float>::signaling_NaN());
-    EbsdLib::FloatArrayType& subEulers = *subEulersPtr;
+    const EbsdLib::FloatArrayType::Pointer subEulerAnglesPtr = EbsdLib::FloatArrayType::CreateArray(count, eulerCompDim, "Euler_Angles_Per_Phase", true);
+    subEulerAnglesPtr->initializeWithValue(std::numeric_limits<float>::signaling_NaN());
+    EbsdLib::FloatArrayType& subEulerAngles = *subEulerAnglesPtr;
 
-    // Now loop through the eulers again and this time add them to the sub-Eulers Array
+    // Now loop through the Euler angles again and this time add them to the sub-Euler angle Array
     count = 0;
     for(size_t i = 0; i < numPoints; ++i)
     {
@@ -400,22 +734,23 @@ Result<> WritePoleFigure::operator()()
       {
         if(!m_InputValues->UseMask || maskCompare->isTrue(i))
         {
-          subEulers[count * 3] = eulers[i * 3];
-          subEulers[count * 3 + 1] = eulers[i * 3 + 1];
-          subEulers[count * 3 + 2] = eulers[i * 3 + 2];
+          subEulerAngles[count * 3] = eulerAngles[i * 3];
+          subEulerAngles[count * 3 + 1] = eulerAngles[i * 3 + 1];
+          subEulerAngles[count * 3 + 2] = eulerAngles[i * 3 + 2];
           count++;
         }
       }
     }
-    if(subEulersPtr->getNumberOfTuples() == 0)
+    if(subEulerAnglesPtr->getNumberOfTuples() == 0)
     {
       continue;
     } // Skip because we have no Pole Figure data
 
     std::vector<EbsdLib::UInt8ArrayType::Pointer> figures;
+    std::vector<EbsdLib::DoubleArrayType::Pointer> intensityImages;
 
     PoleFigureConfiguration_t config;
-    config.eulers = subEulersPtr.get();
+    config.eulers = subEulerAnglesPtr.get();
     config.imageDim = m_InputValues->ImageSize;
     config.lambertDim = m_InputValues->LambertSize;
     config.numColors = m_InputValues->NumColors;
@@ -428,54 +763,111 @@ Result<> WritePoleFigure::operator()()
     {
     case EbsdLib::CrystalStructure::Cubic_High:
       figures = makePoleFigures<CubicOps>(config);
+      intensityImages = createIntensityPoleFigures<CubicOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Cubic_Low:
       figures = makePoleFigures<CubicLowOps>(config);
+      intensityImages = createIntensityPoleFigures<CubicLowOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Hexagonal_High:
       figures = makePoleFigures<HexagonalOps>(config);
+      intensityImages = createIntensityPoleFigures<HexagonalOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Hexagonal_Low:
       figures = makePoleFigures<HexagonalLowOps>(config);
+      intensityImages = createIntensityPoleFigures<HexagonalLowOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Trigonal_High:
       figures = makePoleFigures<TrigonalOps>(config);
+      intensityImages = createIntensityPoleFigures<TrigonalOps>(config, m_InputValues->NormalizeToMRD);
       //   setWarningCondition(-1010, "Trigonal High Symmetry is not supported for Pole figures. This phase will be omitted from results");
       break;
     case EbsdLib::CrystalStructure::Trigonal_Low:
       figures = makePoleFigures<TrigonalLowOps>(config);
+      intensityImages = createIntensityPoleFigures<TrigonalLowOps>(config, m_InputValues->NormalizeToMRD);
       //  setWarningCondition(-1010, "Trigonal Low Symmetry is not supported for Pole figures. This phase will be omitted from results");
       break;
     case EbsdLib::CrystalStructure::Tetragonal_High:
       figures = makePoleFigures<TetragonalOps>(config);
+      intensityImages = createIntensityPoleFigures<TetragonalOps>(config, m_InputValues->NormalizeToMRD);
       //  setWarningCondition(-1010, "Tetragonal High Symmetry is not supported for Pole figures. This phase will be omitted from results");
       break;
     case EbsdLib::CrystalStructure::Tetragonal_Low:
       figures = makePoleFigures<TetragonalLowOps>(config);
+      intensityImages = createIntensityPoleFigures<TetragonalLowOps>(config, m_InputValues->NormalizeToMRD);
       // setWarningCondition(-1010, "Tetragonal Low Symmetry is not supported for Pole figures. This phase will be omitted from results");
       break;
     case EbsdLib::CrystalStructure::OrthoRhombic:
       figures = makePoleFigures<OrthoRhombicOps>(config);
+      intensityImages = createIntensityPoleFigures<OrthoRhombicOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Monoclinic:
       figures = makePoleFigures<MonoclinicOps>(config);
+      intensityImages = createIntensityPoleFigures<MonoclinicOps>(config, m_InputValues->NormalizeToMRD);
       break;
     case EbsdLib::CrystalStructure::Triclinic:
       figures = makePoleFigures<TriclinicOps>(config);
+      intensityImages = createIntensityPoleFigures<TriclinicOps>(config, m_InputValues->NormalizeToMRD);
       break;
     default:
       break;
     }
 
+    if(m_InputValues->SaveIntensityData && intensityImages.size() == 3)
+    {
+      DataPath amPath = m_InputValues->IntensityGeometryDataPath.createChildPath(write_pole_figure::k_ImageAttrMatName);
+      // If there is more than a single phase we will need to add more arrays to the DataStructure
+      if(phase > 1)
+      {
+        const std::vector<size_t> intensityImageDims = {static_cast<usize>(config.imageDim), static_cast<usize>(config.imageDim), 1ULL};
+        DataPath arrayDataPath = amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot1Name));
+        Result<> result = CreateArray<float64>(m_DataStructure, intensityImageDims, {1ULL}, arrayDataPath, IDataAction::Mode::Execute);
+
+        arrayDataPath = amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot2Name));
+        result = CreateArray<float64>(m_DataStructure, intensityImageDims, {1ULL}, arrayDataPath, IDataAction::Mode::Execute);
+
+        arrayDataPath = amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot3Name));
+        result = CreateArray<float64>(m_DataStructure, intensityImageDims, {1ULL}, arrayDataPath, IDataAction::Mode::Execute);
+      }
+
+      auto intensityPlot1Array = m_DataStructure.getDataRefAs<Float64Array>(amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot1Name)));
+      auto intensityPlot2Array = m_DataStructure.getDataRefAs<Float64Array>(amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot2Name)));
+      auto intensityPlot3Array = m_DataStructure.getDataRefAs<Float64Array>(amPath.createChildPath(fmt::format("Phase_{}_{}", phase, m_InputValues->IntensityPlot3Name)));
+
+      std::vector<size_t> compDims = {1ULL};
+      for(int imageIndex = 0; imageIndex < figures.size(); imageIndex++)
+      {
+        intensityImages[imageIndex] = flipAndMirrorPoleFigure<double>(intensityImages[imageIndex].get(), config);
+      }
+
+      std::copy(intensityImages[0]->begin(), intensityImages[0]->end(), intensityPlot1Array.begin());
+      std::copy(intensityImages[1]->begin(), intensityImages[1]->end(), intensityPlot2Array.begin());
+      std::copy(intensityImages[2]->begin(), intensityImages[2]->end(), intensityPlot3Array.begin());
+
+      DataPath metaDataPath = m_InputValues->IntensityGeometryDataPath.createChildPath(write_pole_figure::k_MetaDataName);
+      auto metaDataArrayRef = m_DataStructure.getDataRefAs<StringArray>(metaDataPath);
+      if(metaDataArrayRef.getNumberOfTuples() != numPhases)
+      {
+        metaDataArrayRef.resizeTuples(std::vector<usize>{numPhases});
+      }
+
+      std::vector<std::string> laueNames = LaueOps::GetLaueNames();
+      const uint32_t laueIndex = crystalStructures[phase];
+      const std::string materialName = materialNames[phase];
+
+      metaDataArrayRef[phase] = fmt::format("Phase Num: {}\nMaterial Name: {}\nLaue Group: {}\nHemisphere: Northern\nSamples: {}\nLambert Square Dim: {}", phase, materialName, laueNames[laueIndex],
+                                            config.eulers->getNumberOfTuples(), config.lambertDim);
+    }
+
     if(figures.size() == 3)
     {
-      const uint32 imageWidth = static_cast<uint32>(config.imageDim);
-      const uint32 imageHeight = static_cast<uint32>(config.imageDim);
-      const float32 fontPtSize = imageHeight / 16.0f;
-      const float32 margins = imageHeight / 32.0f;
+      const auto imageWidth = static_cast<int32>(config.imageDim);
+      const auto imageHeight = static_cast<int32>(config.imageDim);
+      const float32 fontPtSize = static_cast<float>(imageHeight) / 16.0f;
+      const float32 margins = static_cast<float>(imageHeight) / 32.0f;
 
       int32 pageWidth = 0;
-      int32 pageHeight = margins + fontPtSize;
+      auto pageHeight = static_cast<int32>(margins + fontPtSize);
 
       float32 xCharWidth = 0.0f;
       {
@@ -485,14 +877,14 @@ Result<> WritePoleFigure::operator()()
         xCharWidth = tempContext.measure_text(buf.data());
       }
       // Each Pole Figure gets its own Square mini canvas to draw into.
-      const float32 subCanvasWidth = margins + imageWidth + xCharWidth + margins;
-      const float32 subCanvasHeight = margins + fontPtSize + imageHeight + fontPtSize * 2 + margins * 2;
+      const float32 subCanvasWidth = margins + static_cast<float32>(imageWidth) + xCharWidth + margins;
+      const float32 subCanvasHeight = margins + fontPtSize + static_cast<float32>(imageHeight) + fontPtSize * 2 + margins * 2;
 
       std::vector<std::pair<float32, float32>> globalImageOrigins(4);
       if(static_cast<WritePoleFigure::LayoutType>(m_InputValues->ImageLayout) == WritePoleFigure::LayoutType::Horizontal)
       {
-        pageWidth = subCanvasWidth * 4;
-        pageHeight = pageHeight + subCanvasHeight;
+        pageWidth = static_cast<int32>(subCanvasWidth) * 4;
+        pageHeight = pageHeight + static_cast<int32>(subCanvasHeight);
         globalImageOrigins[0] = std::make_pair(0.0f, static_cast<float>(pageHeight) - subCanvasHeight);
         globalImageOrigins[1] = std::make_pair(subCanvasWidth, static_cast<float>(pageHeight) - subCanvasHeight);
         globalImageOrigins[2] = std::make_pair(subCanvasWidth * 2.0f, static_cast<float>(pageHeight) - subCanvasHeight);
@@ -500,8 +892,8 @@ Result<> WritePoleFigure::operator()()
       }
       else if(static_cast<WritePoleFigure::LayoutType>(m_InputValues->ImageLayout) == WritePoleFigure::LayoutType::Vertical)
       {
-        pageWidth = subCanvasWidth;
-        pageHeight = pageHeight + subCanvasHeight * 4.0f;
+        pageWidth = static_cast<int32>(subCanvasWidth);
+        pageHeight = pageHeight + static_cast<int32>(subCanvasHeight) * 4;
         globalImageOrigins[0] = std::make_pair(0.0f, margins + fontPtSize);
         globalImageOrigins[1] = std::make_pair(0.0f, margins + fontPtSize + subCanvasHeight * 1.0f);
         globalImageOrigins[2] = std::make_pair(0.0f, margins + fontPtSize + subCanvasHeight * 2.0f);
@@ -509,8 +901,8 @@ Result<> WritePoleFigure::operator()()
       }
       else if(static_cast<WritePoleFigure::LayoutType>(m_InputValues->ImageLayout) == nx::core::WritePoleFigure::LayoutType::Square)
       {
-        pageWidth = subCanvasWidth * 2.0f;
-        pageHeight = pageHeight + subCanvasHeight * 2.0f;
+        pageWidth = static_cast<int32>(subCanvasWidth) * 2;
+        pageHeight = pageHeight + static_cast<int32>(subCanvasHeight) * 2;
         globalImageOrigins[0] = std::make_pair(0.0f, (static_cast<float>(pageHeight) - 2.0f * subCanvasHeight));           // Upper Left
         globalImageOrigins[1] = std::make_pair(subCanvasWidth, (static_cast<float>(pageHeight) - 2.0f * subCanvasHeight)); // Upper Right
         globalImageOrigins[2] = std::make_pair(0.0f, (static_cast<float>(pageHeight) - subCanvasHeight));                  // Lower Left
@@ -535,23 +927,27 @@ Result<> WritePoleFigure::operator()()
       context.set_color(canvas_ity::fill_style, 1.0f, 1.0f, 1.0f, 1.0f);
       context.fill();
 
+      std::vector<size_t> compDims = {4ULL};
       for(int imageIndex = 0; imageIndex < figures.size(); imageIndex++)
       {
         figures[imageIndex] = flipAndMirrorPoleFigure(figures[imageIndex].get(), config);
+        figures[imageIndex] = convertColorOrder(figures[imageIndex].get(), config);
       }
 
       for(int i = 0; i < 3; i++)
       {
         std::array<float, 2> figureOrigin = {0.0f, 0.0f};
         std::tie(figureOrigin[0], figureOrigin[1]) = globalImageOrigins[i];
-        context.draw_image(figures[i]->getPointer(0), imageWidth, imageHeight, imageWidth * 4, figureOrigin[0] + margins, figureOrigin[1] + fontPtSize * 2 + margins * 2, imageWidth, imageHeight);
+        context.draw_image(figures[i]->getPointer(0), imageWidth, imageHeight, imageWidth * figures[i]->getNumberOfComponents(), figureOrigin[0] + margins,
+                           figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f, static_cast<float32>(imageWidth), static_cast<float32>(imageHeight));
 
         // Draw an outline on the figure
         context.begin_path();
         context.line_cap = canvas_ity::circle;
         context.set_line_width(3.0f);
         context.set_color(canvas_ity::stroke_style, 0.0f, 0.0f, 0.0f, 1.0f);
-        context.arc(figureOrigin[0] + margins + m_InputValues->ImageSize / 2.0f, figureOrigin[1] + fontPtSize * 2 + margins * 2 + m_InputValues->ImageSize / 2.0f, m_InputValues->ImageSize / 2.0f, 0,
+        context.arc(figureOrigin[0] + margins + static_cast<float32>(m_InputValues->ImageSize) / 2.0f,
+                    figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize) / 2.0f, static_cast<float32>(m_InputValues->ImageSize) / 2.0f, 0,
                     nx::core::Constants::k_2Pi<float>);
         context.stroke();
         context.close_path();
@@ -561,8 +957,9 @@ Result<> WritePoleFigure::operator()()
         context.line_cap = canvas_ity::square;
         context.set_line_width(2.0f);
         context.set_color(canvas_ity::stroke_style, 0.0f, 0.0f, 0.0f, 1.0f);
-        context.move_to(figureOrigin[0] + margins, figureOrigin[1] + fontPtSize * 2 + margins * 2 + m_InputValues->ImageSize / 2.0f);
-        context.line_to(figureOrigin[0] + margins + m_InputValues->ImageSize, figureOrigin[1] + fontPtSize * 2 + margins * 2 + m_InputValues->ImageSize / 2.0f);
+        context.move_to(figureOrigin[0] + margins, figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize) / 2.0f);
+        context.line_to(figureOrigin[0] + margins + static_cast<float32>(m_InputValues->ImageSize),
+                        figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize) / 2.0f);
         context.stroke();
         context.close_path();
 
@@ -571,8 +968,9 @@ Result<> WritePoleFigure::operator()()
         context.line_cap = canvas_ity::square;
         context.set_line_width(2.0f);
         context.set_color(canvas_ity::stroke_style, 0.0f, 0.0f, 0.0f, 1.0f);
-        context.move_to(figureOrigin[0] + margins + m_InputValues->ImageSize / 2.0f, figureOrigin[1] + fontPtSize * 2 + margins * 2);
-        context.line_to(figureOrigin[0] + margins + m_InputValues->ImageSize / 2.0f, figureOrigin[1] + fontPtSize * 2 + margins * 2 + m_InputValues->ImageSize);
+        context.move_to(figureOrigin[0] + margins + static_cast<float32>(m_InputValues->ImageSize) / 2.0f, figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f);
+        context.line_to(figureOrigin[0] + margins + static_cast<float32>(m_InputValues->ImageSize) / 2.0f,
+                        figureOrigin[1] + fontPtSize * 2.0f + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize));
         context.stroke();
         context.close_path();
 
@@ -581,7 +979,8 @@ Result<> WritePoleFigure::operator()()
         context.set_font(m_LatoBold.data(), static_cast<int>(m_LatoBold.size()), fontPtSize);
         context.set_color(canvas_ity::fill_style, 0.0f, 0.0f, 0.0f, 1.0f);
         context.text_baseline = baselines[0];
-        context.fill_text("X", figureOrigin[0] + margins * 2.0f + m_InputValues->ImageSize, figureOrigin[1] + fontPtSize * 2.25 + margins * 2 + m_InputValues->ImageSize / 2.0f);
+        context.fill_text("X", figureOrigin[0] + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize),
+                          figureOrigin[1] + fontPtSize * 2.25f + margins * 2.0f + static_cast<float32>(m_InputValues->ImageSize) / 2.0f);
         context.close_path();
 
         // Draw Y Axis Label
@@ -590,7 +989,7 @@ Result<> WritePoleFigure::operator()()
         context.set_color(canvas_ity::fill_style, 0.0f, 0.0f, 0.0f, 1.0f);
         context.text_baseline = baselines[0];
         const float yFontWidth = context.measure_text("Y");
-        context.fill_text("Y", figureOrigin[0] + margins - (0.5f * yFontWidth) + m_InputValues->ImageSize / 2.0f, figureOrigin[1] + fontPtSize * 2 + margins);
+        context.fill_text("Y", figureOrigin[0] + margins - (0.5f * yFontWidth) + static_cast<float32>(m_InputValues->ImageSize) / 2.0f, figureOrigin[1] + fontPtSize * 2.0f + margins);
         context.close_path();
 
         // Draw the figure subtitle. This is usually the direction or plane family
@@ -651,11 +1050,12 @@ Result<> WritePoleFigure::operator()()
       // Now draw the Color Scalar Bar if needed.
       if(config.discrete)
       {
-        drawInformationBlock(context, config, globalImageOrigins[3], margins, imageHeight / 20.0f, static_cast<int32_t>(phase), m_LatoRegular, laueNames[laueIndex], materialName);
+        drawInformationBlock(context, config, globalImageOrigins[3], margins, static_cast<float32>(imageHeight) / 20.0f, static_cast<int32_t>(phase), m_LatoRegular, laueNames[laueIndex],
+                             materialName);
       }
       else
       {
-        drawScalarBar(context, config, globalImageOrigins[3], margins, imageHeight / 20.0f, static_cast<int32_t>(phase), m_LatoRegular, laueNames[laueIndex], materialName);
+        drawScalarBar(context, config, globalImageOrigins[3], margins, static_cast<float32>(imageHeight) / 20.0f, static_cast<int32_t>(phase), m_LatoRegular, laueNames[laueIndex], materialName);
       }
 
       // Fetch the rendered RGBA pixels from the entire canvas.
@@ -670,7 +1070,7 @@ Result<> WritePoleFigure::operator()()
         tupleShape[1] = pageHeight;
         tupleShape[2] = pageWidth;
         // Create an output array to hold the RGB formatted color image
-        auto imageArrayPath = cellAttrMatPath.createChildPath(fmt::format("{}{}", m_InputValues->ImagePrefix, phase));
+        auto imageArrayPath = cellAttrMatPath.createChildPath(fmt::format("Phase_{}", phase));
         auto arrayCreationResult = nx::core::CreateArray<uint8>(m_DataStructure, tupleShape, {3ULL}, imageArrayPath, IDataAction::Mode::Execute);
         if(arrayCreationResult.invalid())
         {
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.hpp
index e42f854cee..57dde873d3 100644
--- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.hpp
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/WritePoleFigure.hpp
@@ -17,7 +17,10 @@ namespace write_pole_figure
 {
 const std::string k_ImageAttrMatName("Cell Data");
 const std::string k_ImageDataName("Image");
+const std::string k_MetaDataName("MetaData");
+
 } // namespace write_pole_figure
+
 struct ORIENTATIONANALYSIS_EXPORT WritePoleFigureInputValues
 {
   StringParameter::ValueType Title;
@@ -39,6 +42,13 @@ struct ORIENTATIONANALYSIS_EXPORT WritePoleFigureInputValues
   bool SaveAsImageGeometry;
   bool WriteImageToDisk;
   DataPath OutputImageGeometryPath;
+
+  bool SaveIntensityData;
+  bool NormalizeToMRD;
+  DataPath IntensityGeometryDataPath;
+  std::string IntensityPlot1Name;
+  std::string IntensityPlot2Name;
+  std::string IntensityPlot3Name;
 };
 
 /**
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.cpp
index c10473cec7..aeae9c61b2 100644
--- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.cpp
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.cpp
@@ -9,7 +9,9 @@
 #include "simplnx/DataStructure/StringArray.hpp"
 #include "simplnx/Filter/Actions/CreateArrayAction.hpp"
 #include "simplnx/Filter/Actions/CreateImageGeometryAction.hpp"
+#include "simplnx/Filter/Actions/CreateStringArrayAction.hpp"
 #include "simplnx/Filter/Actions/DeleteDataAction.hpp"
+#include "simplnx/Parameters/ArrayCreationParameter.hpp"
 #include "simplnx/Parameters/ArraySelectionParameter.hpp"
 #include "simplnx/Parameters/BoolParameter.hpp"
 #include "simplnx/Parameters/ChoicesParameter.hpp"
@@ -122,6 +124,14 @@ Parameters WritePoleFigureFilter::parameters() const
   params.insertLinkableParameter(std::make_unique<BoolParameter>(k_SaveAsImageGeometry_Key, "Save Output as Image Geometry", "Save the generated pole figure as an ImageGeometry", true));
   params.insert(std::make_unique<DataGroupCreationParameter>(k_ImageGeometryPath_Key, "Output Image Geometry", "The path to the created Image Geometry", DataPath({"PoleFigure"})));
 
+  params.insertSeparator(Parameters::Separator{"Output Count Data Arrays"});
+  params.insertLinkableParameter(std::make_unique<BoolParameter>(k_SaveIntensityDataArrays, "Save Count Images", "Save the Count Plots (x3)", true));
+  params.insert(std::make_unique<DataGroupCreationParameter>(k_IntensityGeometryPath, "Output Count Image Geometry", "The path to the created Count Image Geometries", DataPath({"Intensity Data"})));
+  params.insert(std::make_unique<BoolParameter>(k_NormalizeToMRD, "Normalize Count Data to MRD", "The Pole Figure data should be normalized to MRD values", true));
+  params.insert(std::make_unique<DataObjectNameParameter>(k_IntensityPlot1Name, "Count Plot 1", "The counts data for the plot", "<001>"));
+  params.insert(std::make_unique<DataObjectNameParameter>(k_IntensityPlot2Name, "Count Plot 2", "The counts data for the plot", "<011>"));
+  params.insert(std::make_unique<DataObjectNameParameter>(k_IntensityPlot3Name, "Count Plot 3", "The counts data for the plot", "<111>"));
+
   // Associate the Linkable Parameter(s) to the children parameters that they control
   params.linkParameters(k_UseMask_Key, k_MaskArrayPath_Key, true);
 
@@ -132,6 +142,12 @@ Parameters WritePoleFigureFilter::parameters() const
   params.linkParameters(k_WriteImageToDisk, k_OutputPath_Key, true);
   params.linkParameters(k_WriteImageToDisk, k_ImagePrefix_Key, true);
 
+  params.linkParameters(k_SaveIntensityDataArrays, k_ImageGeometryPath_Key, true);
+  params.linkParameters(k_SaveIntensityDataArrays, k_IntensityPlot1Name, true);
+  params.linkParameters(k_SaveIntensityDataArrays, k_IntensityPlot2Name, true);
+  params.linkParameters(k_SaveIntensityDataArrays, k_IntensityPlot3Name, true);
+  params.linkParameters(k_SaveIntensityDataArrays, k_NormalizeToMRD, true);
+
   return params;
 }
 
@@ -161,6 +177,12 @@ IFilter::PreflightResult WritePoleFigureFilter::preflightImpl(const DataStructur
   auto pWriteImageToDisk = filterArgs.value<bool>(k_WriteImageToDisk);
   auto pOutputImageGeometryPath = filterArgs.value<DataPath>(k_ImageGeometryPath_Key);
 
+  auto pSaveIntensityPlot = filterArgs.value<bool>(k_SaveIntensityDataArrays);
+  auto pOutputIntensityGeometryPath = filterArgs.value<DataPath>(k_IntensityGeometryPath);
+  auto pIntensityPlot1Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot1Name);
+  auto pIntensityPlot2Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot2Name);
+  auto pIntensityPlot3Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot3Name);
+
   PreflightResult preflightResult;
 
   const auto* materialNamePtr = dataStructure.getDataAs<StringArray>(pMaterialNameArrayPathValue);
@@ -198,11 +220,12 @@ IFilter::PreflightResult WritePoleFigureFilter::preflightImpl(const DataStructur
     pageWidth = subCanvasWidth * 2.0f;
     pageHeight = pageHeight + subCanvasHeight * 2.0f;
   }
-  const std::vector<size_t> dims = {static_cast<usize>(pageWidth), static_cast<usize>(pageHeight), 1ULL};
-  auto createImageGeometryAction =
-      std::make_unique<CreateImageGeometryAction>(pOutputImageGeometryPath, dims, std::vector<float>{0.0f, 0.0f, 0.0f}, std::vector<float>{1.0f, 1.0f, 1.0f}, write_pole_figure::k_ImageAttrMatName);
-  resultOutputActions.value().appendAction(std::move(createImageGeometryAction));
-
+  {
+    const std::vector<size_t> dims = {static_cast<usize>(pageWidth), static_cast<usize>(pageHeight), 1ULL};
+    auto createImageGeometryAction =
+        std::make_unique<CreateImageGeometryAction>(pOutputImageGeometryPath, dims, std::vector<float>{0.0f, 0.0f, 0.0f}, std::vector<float>{1.0f, 1.0f, 1.0f}, write_pole_figure::k_ImageAttrMatName);
+    resultOutputActions.value().appendAction(std::move(createImageGeometryAction));
+  }
   if(!pSaveAsImageGeometry)
   {
     // After the execute function has been done, delete the original image geometry
@@ -215,6 +238,31 @@ IFilter::PreflightResult WritePoleFigureFilter::preflightImpl(const DataStructur
     preflightUpdatedValues.push_back({"Example Output File.", fmt::format("{}/{}Phase_1.tiff", pOutputPathValue.string(), pImagePrefixValue)});
   }
 
+  if(pSaveIntensityPlot)
+  {
+    const std::vector<size_t> intensityImageDims = {static_cast<usize>(pImageSizeValue), static_cast<usize>(pImageSizeValue), 1ULL};
+    auto createImageGeometryAction = std::make_unique<CreateImageGeometryAction>(pOutputIntensityGeometryPath, intensityImageDims, std::vector<float>{0.0f, 0.0f, 0.0f},
+                                                                                 std::vector<float>{1.0f, 1.0f, 1.0f}, write_pole_figure::k_ImageAttrMatName);
+    resultOutputActions.value().appendAction(std::move(createImageGeometryAction));
+
+    std::vector<size_t> cDims = {1ULL};
+    DataPath path = pOutputIntensityGeometryPath.createChildPath(write_pole_figure::k_ImageAttrMatName).createChildPath(fmt::format("Phase_{}_{}", 1, pIntensityPlot1Name));
+    auto createArray1 = std::make_unique<CreateArrayAction>(DataType::float64, intensityImageDims, cDims, path);
+    resultOutputActions.value().appendAction(std::move(createArray1));
+
+    path = pOutputIntensityGeometryPath.createChildPath(write_pole_figure::k_ImageAttrMatName).createChildPath(fmt::format("Phase_{}_{}", 1, pIntensityPlot2Name));
+    auto createArray2 = std::make_unique<CreateArrayAction>(DataType::float64, intensityImageDims, cDims, path);
+    resultOutputActions.value().appendAction(std::move(createArray2));
+
+    path = pOutputIntensityGeometryPath.createChildPath(write_pole_figure::k_ImageAttrMatName).createChildPath(fmt::format("Phase_{}_{}", 1, pIntensityPlot3Name));
+    auto createArray3 = std::make_unique<CreateArrayAction>(DataType::float64, intensityImageDims, cDims, path);
+    resultOutputActions.value().appendAction(std::move(createArray3));
+
+    path = pOutputIntensityGeometryPath.createChildPath(write_pole_figure::k_MetaDataName);
+    auto createMetaData = std::make_unique<CreateStringArrayAction>(std::vector<usize>{0ULL}, path);
+    resultOutputActions.value().appendAction(std::move(createMetaData));
+  }
+
   // Return both the resultOutputActions and the preflightUpdatedValues via std::move()
   return {std::move(resultOutputActions), std::move(preflightUpdatedValues)};
 }
@@ -244,6 +292,13 @@ Result<> WritePoleFigureFilter::executeImpl(DataStructure& dataStructure, const
   inputValues.WriteImageToDisk = filterArgs.value<bool>(k_WriteImageToDisk);
   inputValues.OutputImageGeometryPath = filterArgs.value<DataPath>(k_ImageGeometryPath_Key);
 
+  inputValues.SaveIntensityData = filterArgs.value<bool>(k_SaveIntensityDataArrays);
+  inputValues.NormalizeToMRD = filterArgs.value<bool>(k_NormalizeToMRD);
+  inputValues.IntensityGeometryDataPath = filterArgs.value<DataPath>(k_IntensityGeometryPath);
+  inputValues.IntensityPlot1Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot1Name);
+  inputValues.IntensityPlot2Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot2Name);
+  inputValues.IntensityPlot3Name = filterArgs.value<DataObjectNameParameter::ValueType>(k_IntensityPlot3Name);
+
   return WritePoleFigure(dataStructure, messageHandler, shouldCancel, &inputValues)();
 }
 
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.hpp
index 6ddbf31cf0..806309745a 100644
--- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.hpp
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/WritePoleFigureFilter.hpp
@@ -42,6 +42,12 @@ class ORIENTATIONANALYSIS_EXPORT WritePoleFigureFilter : public IFilter
   static inline constexpr StringLiteral k_SaveAsImageGeometry_Key = "save_as_image_geometry";
   static inline constexpr StringLiteral k_WriteImageToDisk = "write_image_to_disk";
   static inline constexpr StringLiteral k_ImageGeometryPath_Key = "output_image_geometry_path";
+  static inline constexpr StringLiteral k_SaveIntensityDataArrays = "save_intensity_plots";
+  static inline constexpr StringLiteral k_NormalizeToMRD = "normalize_to_mrd";
+  static inline constexpr StringLiteral k_IntensityGeometryPath = "intensity_geometry_path";
+  static inline constexpr StringLiteral k_IntensityPlot1Name = "intensity_plot_1_name";
+  static inline constexpr StringLiteral k_IntensityPlot2Name = "intensity_plot_2_name";
+  static inline constexpr StringLiteral k_IntensityPlot3Name = "intensity_plot_3_name";
 
   /**
    * @brief Reads SIMPL json and converts it simplnx Arguments.
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/IntersectionUtilities.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/IntersectionUtilities.hpp
new file mode 100644
index 0000000000..2ae27de9db
--- /dev/null
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/IntersectionUtilities.hpp
@@ -0,0 +1,222 @@
+#pragma once
+
+#include "simplnx/Common/Constants.hpp"
+#include "simplnx/DataStructure/Geometry/IGeometry.hpp"
+#include "simplnx/DataStructure/Geometry/TriangleGeom.hpp"
+
+#include <array>
+namespace nx
+{
+namespace IntersectionUtilities
+{
+inline constexpr float k_Epsilon = 1e-8;
+inline constexpr uint32_t width = 310;
+inline constexpr uint32_t height = 150;
+
+/**
+ *@brief Bilinear Interpolator implements function to execute a 2D interpolation
+ */
+using Vec3f = nx::core::Vec3<float>;
+
+inline float deg2rad(float deg)
+{
+  return deg * nx::core::Constants::k_PiOver180<float>;
+}
+
+inline float clamp(float lo, float hi, float v)
+{
+  return std::max(lo, std::min(hi, v));
+}
+
+inline std::array<float, 6> GetBoundingBoxAtTri(nx::core::TriangleGeom& tri, size_t triId)
+{
+  nx::core::IGeometry::SharedTriList& triList = tri.getFacesRef();
+  // nx::core::IGeometry::SharedVertexList& vertList = tri.getVerticesRef();
+
+  nx::core::uint64* Tri = &triList[triId];
+  std::array<nx::core::Point3Df, 3> triCoords;
+  tri.getFaceCoordinates(triId, triCoords);
+  nx::core::Point3Df vert1 = triCoords[0];
+  nx::core::Point3Df vert2 = triCoords[1];
+  nx::core::Point3Df vert3 = triCoords[2];
+  // nx::core::float32* vert1 = &vertList[Tri[0] * 3];
+  // nx::core::float32* vert2 = &vertList[Tri[1]*3+1];
+  // nx::core::float32* vert3 = &vertList[Tri[2]*3+2];
+  auto xMinMax = std::minmax({vert1[0], vert2[0], vert3[0]});
+  auto yMinMax = std::minmax({vert1[1], vert2[1], vert3[1]});
+  auto zMinMax = std::minmax({vert1[2], vert2[2], vert3[2]});
+  return {xMinMax.first, yMinMax.first, zMinMax.first, xMinMax.second, yMinMax.second, zMinMax.second};
+}
+
+/**
+ * @brief This version of the rayTriangleIntersect algorithm uses a more traditional algorithm to
+ * determine if a point is inside a triangle. This version should be more computationally
+ * efficient than the other version
+ * @param rayOrigin
+ * @param rayDirection
+ * @param v0 First Vertex of Triangle
+ * @param v1  Second Vertex of Triangle
+ * @param v2  Third Vertex of Triangle
+ * @param t t part of Barycentric Coord
+ * @param u u part of Barycentric Coord
+ * @param v v part of Barycentric Coord
+ * @return
+ */
+inline bool RayTriangleIntersect(const Vec3f& rayOrigin, const Vec3f& rayDirection, const Vec3f& v0, const Vec3f& v1, const Vec3f& v2, Vec3f& barycentricCoord)
+{
+  // compute plane's normal
+  Vec3f v0v1 = v1 - v0;
+  Vec3f v0v2 = v2 - v0;
+  // no need to normalize
+  Vec3f N = v0v1.cross(v0v2); // N
+  float denom = N.dot(N);
+
+  // Step 1: finding P
+
+  // check if ray and plane are parallel ?
+  float NdotRayDirection = N.dot(rayDirection);
+
+  if(fabs(NdotRayDirection) < k_Epsilon) // almost 0
+  {
+    return false;
+  } // they are parallel so they don't intersect !
+
+  // compute d parameter using equation 2
+  float d = -N.dot(v0);
+
+  // compute t (equation 3)
+  barycentricCoord[2] = -(N.dot(rayOrigin) + d) / NdotRayDirection;
+
+  // check if the triangle is in behind the ray
+  if(barycentricCoord[2] < 0)
+  {
+    return false;
+  } // the triangle is behind
+
+  // compute the intersection point using equation 1
+  Vec3f P = rayOrigin + (rayDirection * barycentricCoord[2]);
+
+  // Step 2: inside-outside test
+  Vec3f C; // vector perpendicular to triangle's plane
+
+  // edge 0
+  Vec3f edge0 = v1 - v0;
+  Vec3f vp0 = P - v0;
+  C = edge0.cross(vp0);
+  if(N.dot(C) < 0)
+  {
+    return false;
+  } // P is on the right side
+
+  // edge 1
+  Vec3f edge1 = v2 - v1;
+  Vec3f vp1 = P - v1;
+  C = edge1.cross(vp1);
+  if((barycentricCoord[0] = N.dot(C)) < 0)
+  {
+    return false;
+  } // P is on the right side
+
+  // edge 2
+  Vec3f edge2 = v0 - v2;
+  Vec3f vp2 = P - v2;
+  C = edge2.cross(vp2);
+  if((barycentricCoord[1] = N.dot(C)) < 0)
+  {
+    return false; // P is on the right side;
+  }
+
+  barycentricCoord[0] /= denom;
+  barycentricCoord[1] /= denom;
+  barycentricCoord[2] = 1 - barycentricCoord[0] - barycentricCoord[1];
+
+  return true; // this ray hits the triangle
+}
+
+/* Adapted from https://github.com/erich666/jgt-code/blob/master/Volume_02/Number_1/Moller1997a/raytri.c */
+/* code rewritten to do tests on the sign of the determinant */
+/* the division is before the test of the sign of the det    */
+/* and one CROSS has been moved out from the if-else if-else
+ * @param rayOrigin
+ * @param rayDirection
+ * @param v0 First Vertex of Triangle
+ * @param v1  Second Vertex of Triangle
+ * @param v2  Third Vertex of Triangle
+ * @param t part of Barycentric Coord
+ * @param u part of Barycentric Coord
+ * @param v part of Barycentric Coord
+ * @return If the point is within the triangle
+ */
+inline bool RayTriangleIntersect2(const Vec3f& orig, const Vec3f& dir, const Vec3f& vert0, const Vec3f& vert1, const Vec3f& vert2, Vec3f& bcoord)
+{
+  Vec3f edge1, edge2, tvec, pvec, qvec;
+  double det, inv_det;
+
+  /* find vectors for two edges sharing vert0 */
+  edge1 = vert1 - vert0;
+  edge2 = vert2 - vert0;
+
+  /* begin calculating determinant - also used to calculate U parameter */
+  pvec = dir.cross(edge2);
+
+  /* if determinant is near zero, ray lies in plane of triangle */
+  det = edge1.dot(pvec);
+
+  /* calculate distance from vert0 to ray origin */
+  tvec = orig - vert0;
+
+  inv_det = 1.0 / det;
+
+  qvec = tvec.cross(edge1);
+
+  if(det > k_Epsilon)
+  {
+    bcoord[0] = tvec.dot(pvec);
+    if(bcoord[0] < 0.0 || bcoord[0] > det)
+    {
+      return false;
+    }
+
+    /* calculate V parameter and test bounds */
+    // bcoord[1] = DOT(dir, qvec);
+    bcoord[1] = dir.dot(qvec);
+    if(bcoord[1] < 0.0 || bcoord[0] + bcoord[1] > det)
+    {
+      return false;
+    }
+  }
+  else if(det < -k_Epsilon)
+  {
+    /* calculate U parameter and test bounds */
+    bcoord[0] = tvec.dot(pvec);
+    if(bcoord[0] > 0.0 || bcoord[0] < det)
+    {
+      return false;
+    }
+
+    /* calculate V parameter and test bounds */
+    bcoord[1] = dir.dot(qvec);
+    if(bcoord[1] > 0.0 || bcoord[0] + bcoord[1] < det)
+    {
+      return false;
+    }
+  }
+  else
+  {
+    return false;
+  } /* ray is parallel to the plane of the triangle */
+
+  // float t = edge2.dotProduct(qvec) * inv_det;
+  float u = bcoord[0] * inv_det;
+  float v = bcoord[1] * inv_det;
+
+  bcoord[0] = 1 - u - v;
+  bcoord[1] = u;
+  bcoord[2] = v;
+
+  return true;
+}
+
+} // namespace IntersectionUtilities
+
+} // namespace nx
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.cpp
new file mode 100644
index 0000000000..3be9c96596
--- /dev/null
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.cpp
@@ -0,0 +1,662 @@
+/*
+* MIT License
+
+Copyright (c) 2018 Volodymyr Bilonenko
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+                                                             copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in all
+   copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+       SOFTWARE.
+           */
+#include "delaunator.h"
+
+#include <iostream>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <numeric>
+#include <stdexcept>
+#include <tuple>
+#include <vector>
+
+namespace nx::delaunator
+{
+
+//@see https://stackoverflow.com/questions/33333363/built-in-mod-vs-custom-mod-function-improve-the-performance-of-modulus-op/33333636#33333636
+inline size_t fast_mod(const size_t i, const size_t c)
+{
+  return i >= c ? i % c : i;
+}
+
+// Kahan and Babuska summation, Neumaier variant; accumulates less FP error
+inline double sum(const std::vector<double>& x)
+{
+  double sum = x[0];
+  double err = 0.0;
+
+  for(size_t i = 1; i < x.size(); i++)
+  {
+    const double k = x[i];
+    const double m = sum + k;
+    err += std::fabs(sum) >= std::fabs(k) ? sum - m + k : k - m + sum;
+    sum = m;
+  }
+  return sum + err;
+}
+
+inline double dist(const double ax, const double ay, const double bx, const double by)
+{
+  const double dx = ax - bx;
+  const double dy = ay - by;
+  return dx * dx + dy * dy;
+}
+
+inline double circumradius(const Point& p1, const Point& p2, const Point& p3)
+{
+  Point d = Point::vector(p1, p2);
+  Point e = Point::vector(p1, p3);
+
+  const double bl = d.magnitude2();
+  const double cl = e.magnitude2();
+  const double det = Point::determinant(d, e);
+
+  Point radius((e.y() * bl - d.y() * cl) * 0.5 / det, (d.x() * cl - e.x() * bl) * 0.5 / det);
+
+  if((bl > 0.0 || bl < 0.0) && (cl > 0.0 || cl < 0.0) && (det > 0.0 || det < 0.0))
+    return radius.magnitude2();
+  return (std::numeric_limits<double>::max)();
+}
+
+inline double circumradius(const double ax, const double ay, const double bx, const double by, const double cx, const double cy)
+{
+  const double dx = bx - ax;
+  const double dy = by - ay;
+  const double ex = cx - ax;
+  const double ey = cy - ay;
+
+  const double bl = dx * dx + dy * dy;
+  const double cl = ex * ex + ey * ey;
+  const double d = dx * ey - dy * ex;
+
+  const double x = (ey * bl - dy * cl) * 0.5 / d;
+  const double y = (dx * cl - ex * bl) * 0.5 / d;
+
+  if((bl > 0.0 || bl < 0.0) && (cl > 0.0 || cl < 0.0) && (d > 0.0 || d < 0.0))
+  {
+    return x * x + y * y;
+  }
+  else
+  {
+    return (std::numeric_limits<double>::max)();
+  }
+}
+
+inline bool clockwise(const Point& p0, const Point& p1, const Point& p2)
+{
+  Point v0 = Point::vector(p0, p1);
+  Point v1 = Point::vector(p0, p2);
+  double det = Point::determinant(v0, v1);
+  double dist = v0.magnitude2() + v1.magnitude2();
+  double dist2 = Point::dist2(v0, v1);
+  if(det == 0)
+  {
+    return false;
+  }
+  double reldet = std::abs(dist / det);
+  if(reldet > 1e14)
+    return false;
+  return det < 0;
+}
+
+inline bool clockwise(double px, double py, double qx, double qy, double rx, double ry)
+{
+  Point p0(px, py);
+  Point p1(qx, qy);
+  Point p2(rx, ry);
+  return clockwise(p0, p1, p2);
+}
+
+inline bool counterclockwise(const Point& p0, const Point& p1, const Point& p2)
+{
+  Point v0 = Point::vector(p0, p1);
+  Point v1 = Point::vector(p0, p2);
+  double det = Point::determinant(v0, v1);
+  double dist = v0.magnitude2() + v1.magnitude2();
+  double dist2 = Point::dist2(v0, v1);
+  if(det == 0)
+    return false;
+  double reldet = std::abs(dist / det);
+  if(reldet > 1e14)
+    return false;
+  return det > 0;
+}
+
+inline bool counterclockwise(double px, double py, double qx, double qy, double rx, double ry)
+{
+  Point p0(px, py);
+  Point p1(qx, qy);
+  Point p2(rx, ry);
+  return counterclockwise(p0, p1, p2);
+}
+
+inline Point circumcenter(const double ax, const double ay, const double bx, const double by, const double cx, const double cy)
+{
+  const double dx = bx - ax;
+  const double dy = by - ay;
+  const double ex = cx - ax;
+  const double ey = cy - ay;
+
+  const double bl = dx * dx + dy * dy;
+  const double cl = ex * ex + ey * ey;
+  // ABELL - This is suspect for div-by-0.
+  const double d = dx * ey - dy * ex;
+
+  const double x = ax + (ey * bl - dy * cl) * 0.5 / d;
+  const double y = ay + (dx * cl - ex * bl) * 0.5 / d;
+
+  return Point(x, y);
+}
+
+inline bool in_circle(const double ax, const double ay, const double bx, const double by, const double cx, const double cy, const double px, const double py)
+{
+  const double dx = ax - px;
+  const double dy = ay - py;
+  const double ex = bx - px;
+  const double ey = by - py;
+  const double fx = cx - px;
+  const double fy = cy - py;
+
+  const double ap = dx * dx + dy * dy;
+  const double bp = ex * ex + ey * ey;
+  const double cp = fx * fx + fy * fy;
+
+  return (dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx)) < 0.0;
+}
+
+constexpr double EPSILON = std::numeric_limits<double>::epsilon();
+
+inline bool check_pts_equal(double x1, double y1, double x2, double y2)
+{
+  return std::fabs(x1 - x2) <= EPSILON && std::fabs(y1 - y2) <= EPSILON;
+}
+
+// monotonically increases with real angle, but doesn't need expensive trigonometry
+inline double pseudo_angle(const double dx, const double dy)
+{
+  const double p = dx / (std::abs(dx) + std::abs(dy));
+  return (dy > 0.0 ? 3.0 - p : 1.0 + p) / 4.0; // [0..1)
+}
+
+Delaunator::Delaunator(std::vector<double> const& in_coords)
+: coords(in_coords)
+, m_points(in_coords)
+{
+  std::size_t n = coords.size() >> 1;
+
+  std::vector<std::size_t> ids(n);
+  std::iota(ids.begin(), ids.end(), 0);
+
+  double max_x = std::numeric_limits<double>::lowest();
+  double max_y = std::numeric_limits<double>::lowest();
+  double min_x = (std::numeric_limits<double>::max)();
+  double min_y = (std::numeric_limits<double>::max)();
+  for(const Point& p : m_points)
+  {
+    min_x = std::min(p.x(), min_x);
+    min_y = std::min(p.y(), min_y);
+    max_x = std::max(p.x(), max_x);
+    max_y = std::max(p.y(), max_y);
+  }
+  double width = max_x - min_x;
+  double height = max_y - min_y;
+  double span = width * width + height * height; // Everything is square dist.
+
+  Point center((min_x + max_x) / 2, (min_y + max_y) / 2);
+
+  std::size_t i0 = INVALID_INDEX;
+  std::size_t i1 = INVALID_INDEX;
+  std::size_t i2 = INVALID_INDEX;
+
+  // pick a seed point close to the centroid
+  double min_dist = (std::numeric_limits<double>::max)();
+  for(size_t i = 0; i < m_points.size(); ++i)
+  {
+    const Point& p = m_points[i];
+    const double d = Point::dist2(center, p);
+    if(d < min_dist)
+    {
+      i0 = i;
+      min_dist = d;
+    }
+  }
+
+  const Point& p0 = m_points[i0];
+
+  min_dist = (std::numeric_limits<double>::max)();
+
+  // find the point closest to the seed
+  for(std::size_t i = 0; i < n; i++)
+  {
+    if(i == i0)
+      continue;
+    const double d = Point::dist2(p0, m_points[i]);
+    if(d < min_dist && d > 0.0)
+    {
+      i1 = i;
+      min_dist = d;
+    }
+  }
+
+  const Point& p1 = m_points[i1];
+
+  double min_radius = (std::numeric_limits<double>::max)();
+
+  // find the third point which forms the smallest circumcircle
+  // with the first two
+  for(std::size_t i = 0; i < n; i++)
+  {
+    if(i == i0 || i == i1)
+      continue;
+
+    const double r = circumradius(p0, p1, m_points[i]);
+    if(r < min_radius)
+    {
+      i2 = i;
+      min_radius = r;
+    }
+  }
+
+  if(!(min_radius < (std::numeric_limits<double>::max)()))
+  {
+    throw std::runtime_error("not triangulation");
+  }
+
+  const Point& p2 = m_points[i2];
+
+  if(counterclockwise(p0, p1, p2))
+    std::swap(i1, i2);
+
+  double i0x = p0.x();
+  double i0y = p0.y();
+  double i1x = m_points[i1].x();
+  double i1y = m_points[i1].y();
+  double i2x = m_points[i2].x();
+  double i2y = m_points[i2].y();
+
+  m_center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);
+
+  // Calculate the distances from the center once to avoid having to
+  // calculate for each compare.  This used to be done in the comparator,
+  // but GCC 7.5+ would copy the comparator to iterators used in the
+  // sort, and this was excruciatingly slow when there were many points
+  // because you had to copy the vector of distances.
+  std::vector<double> dists;
+  dists.reserve(m_points.size());
+  for(const Point& p : m_points)
+    dists.push_back(dist(p.x(), p.y(), m_center.x(), m_center.y()));
+
+  // sort the points by distance from the seed triangle circumcenter
+  std::sort(ids.begin(), ids.end(), [&dists](std::size_t i, std::size_t j) { return dists[i] < dists[j]; });
+
+  // initialize a hash table for storing edges of the advancing convex hull
+  m_hash_size = static_cast<std::size_t>(std::ceil(std::sqrt(n)));
+  m_hash.resize(m_hash_size);
+  std::fill(m_hash.begin(), m_hash.end(), INVALID_INDEX);
+
+  // initialize arrays for tracking the edges of the advancing convex hull
+  hull_prev.resize(n);
+  hull_next.resize(n);
+  hull_tri.resize(n);
+
+  hull_start = i0;
+
+  size_t hull_size = 3;
+
+  hull_next[i0] = hull_prev[i2] = i1;
+  hull_next[i1] = hull_prev[i0] = i2;
+  hull_next[i2] = hull_prev[i1] = i0;
+
+  hull_tri[i0] = 0;
+  hull_tri[i1] = 1;
+  hull_tri[i2] = 2;
+
+  m_hash[hash_key(i0x, i0y)] = i0;
+  m_hash[hash_key(i1x, i1y)] = i1;
+  m_hash[hash_key(i2x, i2y)] = i2;
+
+  // ABELL - Why are we doing this is n < 3?  There is no triangulation if
+  //  there is no triangle.
+
+  std::size_t max_triangles = n < 3 ? 1 : 2 * n - 5;
+  triangles.reserve(max_triangles * 3);
+  halfedges.reserve(max_triangles * 3);
+  add_triangle(i0, i1, i2, INVALID_INDEX, INVALID_INDEX, INVALID_INDEX);
+  double xp = std::numeric_limits<double>::quiet_NaN();
+  double yp = std::numeric_limits<double>::quiet_NaN();
+
+  // Go through points based on distance from the center.
+  for(std::size_t k = 0; k < n; k++)
+  {
+    const std::size_t i = ids[k];
+    const double x = coords[2 * i];
+    const double y = coords[2 * i + 1];
+
+    // skip near-duplicate points
+    if(k > 0 && check_pts_equal(x, y, xp, yp))
+      continue;
+    xp = x;
+    yp = y;
+
+    // ABELL - This is dumb.  We have the indices.  Use them.
+    // skip seed triangle points
+    if(check_pts_equal(x, y, i0x, i0y) || check_pts_equal(x, y, i1x, i1y) || check_pts_equal(x, y, i2x, i2y))
+      continue;
+
+    // find a visible edge on the convex hull using edge hash
+    std::size_t start = 0;
+
+    size_t key = hash_key(x, y);
+    for(size_t j = 0; j < m_hash_size; j++)
+    {
+      start = m_hash[fast_mod(key + j, m_hash_size)];
+
+      // ABELL - Not sure how hull_next[start] could ever equal start
+      // I *think* hull_next is just a representation of the hull in one
+      // direction.
+      if(start != INVALID_INDEX && start != hull_next[start])
+        break;
+    }
+
+    // ABELL
+    // Make sure what we found is on the hull.
+    assert(hull_prev[start] != start);
+    assert(hull_prev[start] != INVALID_INDEX);
+
+    start = hull_prev[start];
+    size_t e = start;
+    size_t q;
+
+    // Advance until we find a place in the hull where our current point
+    // can be added.
+    while(true)
+    {
+      q = hull_next[e];
+      if(Point::equal(m_points[i], m_points[e], span) || Point::equal(m_points[i], m_points[q], span))
+      {
+        e = INVALID_INDEX;
+        break;
+      }
+      if(counterclockwise(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1]))
+        break;
+      e = q;
+      if(e == start)
+      {
+        e = INVALID_INDEX;
+        break;
+      }
+    }
+
+    // ABELL
+    // This seems wrong.  Perhaps we should check what's going on?
+    if(e == INVALID_INDEX) // likely a near-duplicate point; skip it
+      continue;
+
+    // add the first triangle from the point
+    std::size_t t = add_triangle(e, i, hull_next[e], INVALID_INDEX, INVALID_INDEX, hull_tri[e]);
+
+    hull_tri[i] = legalize(t + 2); // Legalize the triangle we just added.
+    hull_tri[e] = t;
+    hull_size++;
+
+    // walk forward through the hull, adding more triangles and
+    // flipping recursively
+    std::size_t next = hull_next[e];
+    while(true)
+    {
+      q = hull_next[next];
+      if(!counterclockwise(x, y, coords[2 * next], coords[2 * next + 1], coords[2 * q], coords[2 * q + 1]))
+        break;
+      t = add_triangle(next, i, q, hull_tri[i], INVALID_INDEX, hull_tri[next]);
+      hull_tri[i] = legalize(t + 2);
+      hull_next[next] = next; // mark as removed
+      hull_size--;
+      next = q;
+    }
+
+    // walk backward from the other side, adding more triangles and flipping
+    if(e == start)
+    {
+      while(true)
+      {
+        q = hull_prev[e];
+        if(!counterclockwise(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1]))
+          break;
+        t = add_triangle(q, i, e, INVALID_INDEX, hull_tri[e], hull_tri[q]);
+        legalize(t + 2);
+        hull_tri[q] = t;
+        hull_next[e] = e; // mark as removed
+        hull_size--;
+        e = q;
+      }
+    }
+
+    // update the hull indices
+    hull_prev[i] = e;
+    hull_start = e;
+    hull_prev[next] = i;
+    hull_next[e] = i;
+    hull_next[i] = next;
+
+    m_hash[hash_key(x, y)] = i;
+    m_hash[hash_key(coords[2 * e], coords[2 * e + 1])] = e;
+  }
+}
+
+double Delaunator::get_hull_area()
+{
+  std::vector<double> hull_area;
+  size_t e = hull_start;
+  size_t cnt = 1;
+  do
+  {
+    hull_area.push_back((coords[2 * e] - coords[2 * hull_prev[e]]) * (coords[2 * e + 1] + coords[2 * hull_prev[e] + 1]));
+    cnt++;
+    e = hull_next[e];
+  } while(e != hull_start);
+  return sum(hull_area);
+}
+
+double Delaunator::get_triangle_area()
+{
+  std::vector<double> vals;
+  for(size_t i = 0; i < triangles.size(); i += 3)
+  {
+    const double ax = coords[2 * triangles[i]];
+    const double ay = coords[2 * triangles[i] + 1];
+    const double bx = coords[2 * triangles[i + 1]];
+    const double by = coords[2 * triangles[i + 1] + 1];
+    const double cx = coords[2 * triangles[i + 2]];
+    const double cy = coords[2 * triangles[i + 2] + 1];
+    double val = std::fabs((by - ay) * (cx - bx) - (bx - ax) * (cy - by));
+    vals.push_back(val);
+  }
+  return sum(vals);
+}
+
+std::size_t Delaunator::legalize(std::size_t a)
+{
+  std::size_t i = 0;
+  std::size_t ar = 0;
+  m_edge_stack.clear();
+
+  // recursion eliminated with a fixed-size stack
+  while(true)
+  {
+    const size_t b = halfedges[a];
+
+    /* if the pair of triangles doesn't satisfy the Delaunay condition
+     * (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
+     * then do the same check/flip recursively for the new pair of triangles
+     *
+     *           pl                    pl
+     *          /||\                  /  \
+     *       al/ || \bl            al/    \a
+     *        /  ||  \              /      \
+     *       /  a||b  \    flip    /___ar___\
+     *     p0\   ||   /p1   =>   p0\---bl---/p1
+     *        \  ||  /              \      /
+     *       ar\ || /br             b\    /br
+     *          \||/                  \  /
+     *           pr                    pr
+     */
+    const size_t a0 = 3 * (a / 3);
+    ar = a0 + (a + 2) % 3;
+
+    if(b == INVALID_INDEX)
+    {
+      if(i > 0)
+      {
+        i--;
+        a = m_edge_stack[i];
+        continue;
+      }
+      else
+      {
+        // i = INVALID_INDEX;
+        break;
+      }
+    }
+
+    const size_t b0 = 3 * (b / 3);
+    const size_t al = a0 + (a + 1) % 3;
+    const size_t bl = b0 + (b + 2) % 3;
+
+    const std::size_t p0 = triangles[ar];
+    const std::size_t pr = triangles[a];
+    const std::size_t pl = triangles[al];
+    const std::size_t p1 = triangles[bl];
+
+    const bool illegal = in_circle(coords[2 * p0], coords[2 * p0 + 1], coords[2 * pr], coords[2 * pr + 1], coords[2 * pl], coords[2 * pl + 1], coords[2 * p1], coords[2 * p1 + 1]);
+
+    if(illegal)
+    {
+      triangles[a] = p1;
+      triangles[b] = p0;
+
+      auto hbl = halfedges[bl];
+
+      // Edge swapped on the other side of the hull (rare).
+      // Fix the halfedge reference
+      if(hbl == INVALID_INDEX)
+      {
+        std::size_t e = hull_start;
+        do
+        {
+          if(hull_tri[e] == bl)
+          {
+            hull_tri[e] = a;
+            break;
+          }
+          e = hull_prev[e];
+        } while(e != hull_start);
+      }
+      link(a, hbl);
+      link(b, halfedges[ar]);
+      link(ar, bl);
+      std::size_t br = b0 + (b + 1) % 3;
+
+      if(i < m_edge_stack.size())
+      {
+        m_edge_stack[i] = br;
+      }
+      else
+      {
+        m_edge_stack.push_back(br);
+      }
+      i++;
+    }
+    else
+    {
+      if(i > 0)
+      {
+        i--;
+        a = m_edge_stack[i];
+        continue;
+      }
+      else
+      {
+        break;
+      }
+    }
+  }
+  return ar;
+}
+
+std::size_t Delaunator::hash_key(const double x, const double y) const
+{
+  const double dx = x - m_center.x();
+  const double dy = y - m_center.y();
+  return fast_mod(static_cast<std::size_t>(std::llround(std::floor(pseudo_angle(dx, dy) * static_cast<double>(m_hash_size)))), m_hash_size);
+}
+
+std::size_t Delaunator::add_triangle(std::size_t i0, std::size_t i1, std::size_t i2, std::size_t a, std::size_t b, std::size_t c)
+{
+  std::size_t t = triangles.size();
+  triangles.push_back(i0);
+  triangles.push_back(i1);
+  triangles.push_back(i2);
+  link(t, a);
+  link(t + 1, b);
+  link(t + 2, c);
+  return t;
+}
+
+void Delaunator::link(const std::size_t a, const std::size_t b)
+{
+  std::size_t s = halfedges.size();
+  if(a == s)
+  {
+    halfedges.push_back(b);
+  }
+  else if(a < s)
+  {
+    halfedges[a] = b;
+  }
+  else
+  {
+    throw std::runtime_error("Cannot link edge");
+  }
+  if(b != INVALID_INDEX)
+  {
+    std::size_t s2 = halfedges.size();
+    if(b == s2)
+    {
+      halfedges.push_back(a);
+    }
+    else if(b < s2)
+    {
+      halfedges[b] = a;
+    }
+    else
+    {
+      throw std::runtime_error("Cannot link edge");
+    }
+  }
+}
+
+} // namespace nx::delaunator
diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.h b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.h
new file mode 100644
index 0000000000..5b7fbbddd9
--- /dev/null
+++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/delaunator.h
@@ -0,0 +1,188 @@
+/*
+* MIT License
+
+Copyright (c) 2018 Volodymyr Bilonenko
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+                                                             copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in all
+   copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+       SOFTWARE.
+           */
+
+#pragma once
+
+#ifdef DELAUNATOR_HEADER_ONLY
+#define INLINE inline
+#else
+#define INLINE
+#endif
+
+#include <limits>
+#include <ostream>
+#include <vector>
+
+namespace nx
+{
+
+namespace delaunator
+{
+
+constexpr std::size_t INVALID_INDEX = (std::numeric_limits<std::size_t>::max)();
+
+class Point
+{
+public:
+  Point(double x, double y)
+  : m_x(x)
+  , m_y(y)
+  {
+  }
+  Point()
+  : m_x(0)
+  , m_y(0)
+  {
+  }
+
+  double x() const
+  {
+    return m_x;
+  }
+
+  double y() const
+  {
+    return m_y;
+  }
+
+  double magnitude2() const
+  {
+    return m_x * m_x + m_y * m_y;
+  }
+
+  static double determinant(const Point& p1, const Point& p2)
+  {
+    return p1.m_x * p2.m_y - p1.m_y * p2.m_x;
+  }
+
+  static Point vector(const Point& p1, const Point& p2)
+  {
+    return Point(p2.m_x - p1.m_x, p2.m_y - p1.m_y);
+  }
+
+  static double dist2(const Point& p1, const Point& p2)
+  {
+    Point vec = vector(p1, p2);
+    return vec.m_x * vec.m_x + vec.m_y * vec.m_y;
+  }
+
+  static bool equal(const Point& p1, const Point& p2, double span)
+  {
+    double dist = dist2(p1, p2) / span;
+
+    // ABELL - This number should be examined to figure how how
+    // it correlates with the breakdown of calculating determinants.
+    return dist < 1e-20;
+  }
+
+private:
+  double m_x;
+  double m_y;
+};
+
+inline std::ostream& operator<<(std::ostream& out, const Point& p)
+{
+  out << p.x() << "/" << p.y();
+  return out;
+}
+
+class Points
+{
+public:
+  using const_iterator = Point const*;
+
+  Points(const std::vector<double>& coords)
+  : m_coords(coords)
+  {
+  }
+
+  const Point& operator[](size_t offset)
+  {
+    return reinterpret_cast<const Point&>(*(m_coords.data() + (offset * 2)));
+  };
+
+  Points::const_iterator begin() const
+  {
+    return reinterpret_cast<const Point*>(m_coords.data());
+  }
+  Points::const_iterator end() const
+  {
+    return reinterpret_cast<const Point*>(m_coords.data() + m_coords.size());
+  }
+  size_t size() const
+  {
+    return m_coords.size() / 2;
+  }
+
+private:
+  const std::vector<double>& m_coords;
+};
+
+class Delaunator
+{
+
+public:
+  std::vector<double> const& coords;
+  Points m_points;
+
+  // 'triangles' stores the indices to the 'X's of the input
+  // 'coords'.
+  std::vector<std::size_t> triangles;
+
+  // 'halfedges' store indices into 'triangles'.  If halfedges[X] = Y,
+  // It says that there's an edge from X to Y where a) X and Y are
+  // both indices into triangles and b) X and Y are indices into different
+  // triangles in the array.  This allows you to get from a triangle to
+  // its adjacent triangle.  If the a triangle edge has no adjacent triangle,
+  // its half edge will be INVALID_INDEX.
+  std::vector<std::size_t> halfedges;
+
+  std::vector<std::size_t> hull_prev;
+  std::vector<std::size_t> hull_next;
+
+  // This contains indexes into the triangles array.
+  std::vector<std::size_t> hull_tri;
+  std::size_t hull_start;
+
+  INLINE Delaunator(std::vector<double> const& in_coords);
+  INLINE double get_hull_area();
+  INLINE double get_triangle_area();
+
+private:
+  std::vector<std::size_t> m_hash;
+  Point m_center;
+  std::size_t m_hash_size;
+  std::vector<std::size_t> m_edge_stack;
+
+  INLINE std::size_t legalize(std::size_t a);
+  INLINE std::size_t hash_key(double x, double y) const;
+  INLINE std::size_t add_triangle(std::size_t i0, std::size_t i1, std::size_t i2, std::size_t a, std::size_t b, std::size_t c);
+  INLINE void link(std::size_t a, std::size_t b);
+};
+
+} // namespace delaunator
+
+} // namespace nx
+
+#undef INLINE
diff --git a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt
index ae5d959ff8..e7bb780248 100644
--- a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt
+++ b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt
@@ -117,28 +117,28 @@ if(EXISTS "${DREAM3D_DATA_DIR}" AND SIMPLNX_DOWNLOAD_TEST_FILES)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 12_IN625_GBCD.tar.gz SHA512 f696a8af181505947e6fecfdb1a11fda6c762bba5e85fea8d484b1af00bf18643e1d930d48f092ee238d1c19c9ce7c4fb5a8092d17774bda867961a1400e9cea)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_5_align_sections_mutual_information.tar.gz SHA512 363c52bf3b7677e06508fdb92350457d37bbd58ab2d26259e61ec9317f9146b85ec50c58151f69a0b38fc487f79af2896c8fef4e98d1e87c29af1cf19cbf884b)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_5_h5ebsd_exemplar.tar.gz SHA512 9e5936b6e8e6fd51ecc7c17521ff29898bb0e8bc3d516741c47ae65577b4f2e0f3df8f485f19ecf93f1848b70aeeb98b91c2c5742965107acabeaa40182d898d)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_ImportH5Data.tar.gz SHA512 68cf620d28a515d04040d66c4f8fa9d8d46707c65138b366b47f40d1d56a36c40e2e2cd0c9d35168b68c2cb5b2ce95650fa2a0d4e6ffce980d0dd0654908f40d)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_Small_IN100_GBCD.tar.gz SHA512 543e3bdcee24ff9e5cd80dfdedc39ceef1529e7172cebd01d8e5518292ffdf4e0eb2e79d75854cb3eaca5c60e19c861ca67f369e21b81c306edb66327f47a1e3)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_align_sections_misorientation.tar.gz SHA512 2343654a8bcb344fcc613b6715c1e0b2c780fedbdf06cc8e5306b23f9d7908d5eef8faff7e3f0dd6f7ac734a1c6e2b376832bed38548288cd5e9e0af1b5602a8)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME caxis_data.tar.gz SHA512 56468d3f248661c0d739d9acd5a1554abc700bf136586f698a313804536916850b731603d42a0b93aae47faf2f7ee49d4181b1c3e833f054df6f5c70b5e041dc)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_ebsd_segment_features.tar.gz SHA512 acbb493a0668e0115ac49d4fedbbf7600759b9a66deb5d1004c2749a61d2bad2fcc60344bf72b2aeda5c8c098f458949dd1f8d58cb21682fa1393dfb7d0a1b84)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_find_gbcd_metric_based.tar.gz SHA512 46032c758acc982eccaff38c3252a4f063b6ff5dc9ba3f33ed7f123c6f10771e1a2bdcbff8aab0fae7c91c03bb6025f49f0edbb085b3946b93b87980d31e4543)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_ImportH5Data.tar.gz SHA512 68cf620d28a515d04040d66c4f8fa9d8d46707c65138b366b47f40d1d56a36c40e2e2cd0c9d35168b68c2cb5b2ce95650fa2a0d4e6ffce980d0dd0654908f40d)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_read_ang_data.tar.gz SHA512 1777431a1623e2fffdc2daad9be51a72c3bdf6a83a33893827892c98a811991e21f1cf636e036604d0bbc523d8ca0b9d655c28be3d0e89ccafc1486dfa0bd0c7)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_read_ctf_data_2.tar.gz SHA512 f397fa3bf457615a90a4b48eaafded2aa4952b41ccb28d9da6a83adc38aea9c22f2bb5a955f251edeca9ef8265b6bf1d74e829b1340f45cf52620a237aad1707)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_Small_IN100_GBCD.tar.gz SHA512 543e3bdcee24ff9e5cd80dfdedc39ceef1529e7172cebd01d8e5518292ffdf4e0eb2e79d75854cb3eaca5c60e19c861ca67f369e21b81c306edb66327f47a1e3)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_stats_test.tar.gz SHA512 cb1d1c004ab5a3cb29cc10f7b6c291dd5819e57e303242b8162cd0b268dea24b1d3e5e3811ec4f5ee216179b6eb6b81c76ee7d37b7c49e83c6ab336147b4b14e)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME PoleFigure_Exemplars.tar.gz SHA512 7c90c69020a4fd4de10da81c1bc87499f5aededca57f3fc1aae3c60aad344cab1bdad2927d957a8d5848aca599028a9ade532dabc784f099bce2de39f61060e9)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_Ang_Files.tar.gz SHA512 79e9f6948d4e8e06187e11216a67596fa786ffd2700e51f594ad014090383eb8bcc003e14de2e88082aa9ae512cc4fc9cee22c80066fc54f38c3ebc75267eb5b)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_dream3d.tar.gz SHA512 6dd8a3412532bdc7481f7781c7087b4477c6a1efbe6b214f997dad30c53c59714a751be522f084b98065fe75100c74df901bb8af2f512ef47344d8f7941575cf)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_h5ebsd.tar.gz SHA512 31e606285ea9e8235dcb5f608fd2b252a5ab1492abd975e5ec33a21d083aa9720fe16fb8f752742c140f40e963d692f1a46256b9d36e96b1b09796c1e4ea3db9)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME align_sections.tar.gz SHA512 b6892e437df86bd79bd2f1d8f48e44d05bfe38b3453058744320bfaf1b1dc461a484edc9e593f6b9de4ad4d04c41b5dbd0a30e6fc605341d046aec4c3062e33e)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME bad_data_neighbor_orientation_check.tar.gz SHA512 6452cfb1f2394c10050082256f60a2068cfad78ef742e9e35b1d6e63b3fb7c35c9fe7bbe093bed4dbb4e758c49ec6da7b1f7e2473838a0421f39fbdd9f4a2f76)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME convert_hex_grid_to_square_grid_test.tar.gz SHA512 bb672ebbe2540ba493ad95bea95dac1f85b5634ac3311b5aa774ce3d2177103d1b45a13225221993dd40f0cbe02daf20ccd209d4ae0cab0bf034d97c5b234ba4)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME feature_boundary_neighbor_slip_transmission_1.tar.gz SHA512 5b81c1e42dfd88f1a629ba1d7b378391ed74895dcb21ee69e65de879e9c40f98dfba2959d6d0d491508dbcce0b1e6ee480f7ad9e9ed570a7419eab83113acef3)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME H5Oina_Test_Data.tar.gz SHA512 346573ac6b96983680078e8b0a401aa25bd9302dff382ca86ae4e503ded6db3947c4c5611ee603db519d8a8dc6ed35b044a7bfea9880fade5ab54479d140ea03 )
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME INL_writer.tar.gz SHA512 7d723351e51e84540abfbc38e69a6014852ba34808f9d216a27063a616bcfbd5eb708405305fd83334e48c9ca133d3d0be797c05040e4a115cc612e385d9ada6)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME neighbor_orientation_correlation.tar.gz SHA512 122367452174ade2f24dde7a4610bddc4f147a223722d9b30c1df9eaa2cd2bf25e1c7957aba83f3f9de79b4eadd79339b848f9530d1ebf44c69244ea5442cf85)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME PoleFigure_Exemplars_v2.tar.gz SHA512 1c7169e01e66cff074bd58fc188ec69d26df10c4eb61feeb2b123cab945ea44c0419b5d76ab422fc6c585e0a4078bf796f785e7e6ccccc84cc1c518631afb589)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_Ang_Files.tar.gz SHA512 79e9f6948d4e8e06187e11216a67596fa786ffd2700e51f594ad014090383eb8bcc003e14de2e88082aa9ae512cc4fc9cee22c80066fc54f38c3ebc75267eb5b)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_dream3d.tar.gz SHA512 6dd8a3412532bdc7481f7781c7087b4477c6a1efbe6b214f997dad30c53c59714a751be522f084b98065fe75100c74df901bb8af2f512ef47344d8f7941575cf)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_h5ebsd.tar.gz SHA512 31e606285ea9e8235dcb5f608fd2b252a5ab1492abd975e5ec33a21d083aa9720fe16fb8f752742c140f40e963d692f1a46256b9d36e96b1b09796c1e4ea3db9)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME so3_cubic_high_ipf_001.tar.gz SHA512 dfe4598cd4406e8b83f244302dc4fe0d4367527835c5ddd6567fe8d8ab3484d5b10ba24a8bb31db269256ec0b5272daa4340eedb5a8b397755541b32dd616b85)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME feature_boundary_neighbor_slip_transmission_1.tar.gz SHA512 5b81c1e42dfd88f1a629ba1d7b378391ed74895dcb21ee69e65de879e9c40f98dfba2959d6d0d491508dbcce0b1e6ee480f7ad9e9ed570a7419eab83113acef3)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_read_ang_data.tar.gz SHA512 1777431a1623e2fffdc2daad9be51a72c3bdf6a83a33893827892c98a811991e21f1cf636e036604d0bbc523d8ca0b9d655c28be3d0e89ccafc1486dfa0bd0c7)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_read_ctf_data_2.tar.gz SHA512 f397fa3bf457615a90a4b48eaafded2aa4952b41ccb28d9da6a83adc38aea9c22f2bb5a955f251edeca9ef8265b6bf1d74e829b1340f45cf52620a237aad1707)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME write_stats_gen_odf_angle_file.tar.gz SHA512 be3f663aae1f78e5b789200421534ed9fe293187ec3514796ac8177128b34ded18bb9a98b8e838bb283f9818ac30dc4b19ec379bdd581b1a98eb36d967cdd319)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME INL_writer.tar.gz SHA512 7d723351e51e84540abfbc38e69a6014852ba34808f9d216a27063a616bcfbd5eb708405305fd83334e48c9ca133d3d0be797c05040e4a115cc612e385d9ada6)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME convert_hex_grid_to_square_grid_test.tar.gz SHA512 bb672ebbe2540ba493ad95bea95dac1f85b5634ac3311b5aa774ce3d2177103d1b45a13225221993dd40f0cbe02daf20ccd209d4ae0cab0bf034d97c5b234ba4)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME H5Oina_Test_Data.tar.gz SHA512 346573ac6b96983680078e8b0a401aa25bd9302dff382ca86ae4e503ded6db3947c4c5611ee603db519d8a8dc6ed35b044a7bfea9880fade5ab54479d140ea03 )
 
 endif()
 
diff --git a/src/Plugins/OrientationAnalysis/test/WritePoleFigureTest.cpp b/src/Plugins/OrientationAnalysis/test/WritePoleFigureTest.cpp
index 43509483e7..ee4ee01c8a 100644
--- a/src/Plugins/OrientationAnalysis/test/WritePoleFigureTest.cpp
+++ b/src/Plugins/OrientationAnalysis/test/WritePoleFigureTest.cpp
@@ -36,7 +36,7 @@ void CompareComponentsOfArrays(const DataStructure& dataStructure, const DataPat
   usize exemplaryNumComp = exemplaryDataArray.getNumberOfComponents();
   usize generatedNumComp = generatedDataArray.getNumberOfComponents();
 
-  REQUIRE(exemplaryNumComp == 4);
+  REQUIRE(exemplaryNumComp == 3);
   REQUIRE(generatedNumComp == 3);
 
   REQUIRE(compIndex < exemplaryNumComp);
@@ -62,10 +62,10 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-1", "[OrientationAnalysis]
 {
   Application::GetOrCreateInstance()->loadPlugins(unit_test::k_BuildDir.view(), true);
 
-  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars.tar.gz", "PoleFigure_Exemplars");
+  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars_v2.tar.gz", "PoleFigure_Exemplars_v2");
 
   // Read the Small IN100 Data set
-  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
+  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars_v2/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
   DataStructure dataStructure = UnitTest::LoadDataStructure(baseDataFilePath);
 
   // Instantiate the filter, a DataStructure object and an Arguments Object
@@ -86,8 +86,8 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-1", "[OrientationAnalysis]
   args.insertOrAssign(WritePoleFigureFilter::k_UseMask_Key, std::make_any<bool>(false));
   args.insertOrAssign(WritePoleFigureFilter::k_ImageGeometryPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr Discrete Pole Figure [CALCULATED]"})));
 
-  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure [CALCULATED]", "Cell Data", fmt::format("{}{}", k_ImagePrefix, 1)});
-  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure", "CellData", "Image"});
+  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure [CALCULATED]", "Cell Data", fmt::format("Phase_{}", 1)});
+  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure", "Cell Data", "Phase_1"});
 
   args.insertOrAssign(WritePoleFigureFilter::k_CellEulerAnglesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "EulerAngles"})));
   args.insertOrAssign(WritePoleFigureFilter::k_CellPhasesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "Phases"})));
@@ -115,10 +115,10 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-2", "[OrientationAnalysis]
 {
   Application::GetOrCreateInstance()->loadPlugins(unit_test::k_BuildDir.view(), true);
 
-  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars.tar.gz", "PoleFigure_Exemplars");
+  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars_v2.tar.gz", "PoleFigure_Exemplars_v2");
 
   // Read the Small IN100 Data set
-  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
+  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars_v2/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
   DataStructure dataStructure = UnitTest::LoadDataStructure(baseDataFilePath);
 
   // Instantiate the filter, a DataStructure object and an Arguments Object
@@ -130,7 +130,7 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-2", "[OrientationAnalysis]
   args.insertOrAssign(WritePoleFigureFilter::k_LambertSize_Key, std::make_any<int32>(64));
   args.insertOrAssign(WritePoleFigureFilter::k_NumColors_Key, std::make_any<int32>(32));
   args.insertOrAssign(WritePoleFigureFilter::k_GenerationAlgorithm_Key, std::make_any<ChoicesParameter::ValueType>(1));
-  args.insertOrAssign(WritePoleFigureFilter::k_ImageLayout_Key, std::make_any<ChoicesParameter::ValueType>(1));
+  args.insertOrAssign(WritePoleFigureFilter::k_ImageLayout_Key, std::make_any<ChoicesParameter::ValueType>(0));
   args.insertOrAssign(WritePoleFigureFilter::k_OutputPath_Key, std::make_any<FileSystemPathParameter::ValueType>(fs::path(fmt::format("{}/Dir1/Dir2", unit_test::k_BinaryTestOutputDir))));
   args.insertOrAssign(WritePoleFigureFilter::k_ImagePrefix_Key, std::make_any<StringParameter::ValueType>(k_ImagePrefix));
   args.insertOrAssign(WritePoleFigureFilter::k_ImageSize_Key, std::make_any<int32>(1024));
@@ -139,12 +139,12 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-2", "[OrientationAnalysis]
   args.insertOrAssign(WritePoleFigureFilter::k_UseMask_Key, std::make_any<bool>(true));
   args.insertOrAssign(WritePoleFigureFilter::k_ImageGeometryPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked [CALCULATED]"})));
 
-  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked [CALCULATED]", "Cell Data", fmt::format("{}{}", k_ImagePrefix, 1)});
-  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked", "CellData", "Image"});
+  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked [CALCULATED]", "Cell Data", fmt::format("Phase_{}", 1)});
+  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked", "Cell Data", "Phase_1"});
 
   args.insertOrAssign(WritePoleFigureFilter::k_CellEulerAnglesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "EulerAngles"})));
   args.insertOrAssign(WritePoleFigureFilter::k_CellPhasesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "Phases"})));
-  args.insertOrAssign(WritePoleFigureFilter::k_MaskArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "ThresholdArray"})));
+  args.insertOrAssign(WritePoleFigureFilter::k_MaskArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "Mask"})));
   args.insertOrAssign(WritePoleFigureFilter::k_CrystalStructuresArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "CellEnsembleData", "CrystalStructures"})));
   args.insertOrAssign(WritePoleFigureFilter::k_MaterialNameArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "CellEnsembleData", "MaterialName"})));
 
@@ -169,10 +169,10 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-3", "[OrientationAnalysis]
 {
   Application::GetOrCreateInstance()->loadPlugins(unit_test::k_BuildDir.view(), true);
 
-  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars.tar.gz", "PoleFigure_Exemplars");
+  const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "PoleFigure_Exemplars_v2.tar.gz", "PoleFigure_Exemplars_v2");
 
   // Read the Small IN100 Data set
-  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
+  auto baseDataFilePath = fs::path(fmt::format("{}/PoleFigure_Exemplars_v2/fw-ar-IF1-aptr12-corr.dream3d", unit_test::k_TestFilesDir));
   DataStructure dataStructure = UnitTest::LoadDataStructure(baseDataFilePath);
 
   // Instantiate the filter, a DataStructure object and an Arguments Object
@@ -180,25 +180,25 @@ TEST_CASE("OrientationAnalysis::WritePoleFigureFilter-3", "[OrientationAnalysis]
   Arguments args;
 
   // Create default Parameters for the filter.
-  args.insertOrAssign(WritePoleFigureFilter::k_Title_Key, std::make_any<StringParameter::ValueType>("fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked Color"));
+  args.insertOrAssign(WritePoleFigureFilter::k_Title_Key, std::make_any<StringParameter::ValueType>("fw-ar-IF1-aptr12-corr Pole Figure Masked Color"));
   args.insertOrAssign(WritePoleFigureFilter::k_LambertSize_Key, std::make_any<int32>(64));
   args.insertOrAssign(WritePoleFigureFilter::k_NumColors_Key, std::make_any<int32>(32));
   args.insertOrAssign(WritePoleFigureFilter::k_GenerationAlgorithm_Key, std::make_any<ChoicesParameter::ValueType>(0));
-  args.insertOrAssign(WritePoleFigureFilter::k_ImageLayout_Key, std::make_any<ChoicesParameter::ValueType>(2));
+  args.insertOrAssign(WritePoleFigureFilter::k_ImageLayout_Key, std::make_any<ChoicesParameter::ValueType>(0));
   args.insertOrAssign(WritePoleFigureFilter::k_OutputPath_Key, std::make_any<FileSystemPathParameter::ValueType>(fs::path(fmt::format("{}/Dir1/Dir2", unit_test::k_BinaryTestOutputDir))));
   args.insertOrAssign(WritePoleFigureFilter::k_ImagePrefix_Key, std::make_any<StringParameter::ValueType>(k_ImagePrefix));
   args.insertOrAssign(WritePoleFigureFilter::k_ImageSize_Key, std::make_any<int32>(1024));
   args.insertOrAssign(WritePoleFigureFilter::k_SaveAsImageGeometry_Key, std::make_any<bool>(true));
   args.insertOrAssign(WritePoleFigureFilter::k_WriteImageToDisk, std::make_any<bool>(true));
   args.insertOrAssign(WritePoleFigureFilter::k_UseMask_Key, std::make_any<bool>(true));
-  args.insertOrAssign(WritePoleFigureFilter::k_ImageGeometryPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked Color [CALCULATED]"})));
+  args.insertOrAssign(WritePoleFigureFilter::k_ImageGeometryPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr Pole Figure Masked Color [CALCULATED]"})));
 
-  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked Color [CALCULATED]", "Cell Data", fmt::format("{}{}", k_ImagePrefix, 1)});
-  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Discrete Pole Figure Masked Color", "CellData", "Image"});
+  DataPath calculatedImageData({"fw-ar-IF1-aptr12-corr Pole Figure Masked Color [CALCULATED]", "Cell Data", fmt::format("Phase_{}", 1)});
+  DataPath exemplarImageData({"fw-ar-IF1-aptr12-corr Pole Figure Masked Color", "Cell Data", "Phase_1"});
 
   args.insertOrAssign(WritePoleFigureFilter::k_CellEulerAnglesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "EulerAngles"})));
   args.insertOrAssign(WritePoleFigureFilter::k_CellPhasesArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "Phases"})));
-  args.insertOrAssign(WritePoleFigureFilter::k_MaskArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "ThresholdArray"})));
+  args.insertOrAssign(WritePoleFigureFilter::k_MaskArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "Cell Data", "Mask"})));
   args.insertOrAssign(WritePoleFigureFilter::k_CrystalStructuresArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "CellEnsembleData", "CrystalStructures"})));
   args.insertOrAssign(WritePoleFigureFilter::k_MaterialNameArrayPath_Key, std::make_any<DataPath>(DataPath({"fw-ar-IF1-aptr12-corr", "CellEnsembleData", "MaterialName"})));
 
diff --git a/src/simplnx/Utilities/RTree.hpp b/src/simplnx/Utilities/RTree.hpp
index 85d9e93678..1678b9cc4f 100644
--- a/src/simplnx/Utilities/RTree.hpp
+++ b/src/simplnx/Utilities/RTree.hpp
@@ -1,6 +1,11 @@
 #ifndef RTREE_H
 #define RTREE_H
 
+/* ****************************************************************************
+ * This file comes from the https://github.com/nushoin/RTree repository.
+ * The file is under the MIT license as indicated within that repository
+ **************************************************************************** */
+
 // NOTE This file compiles under MSVC 6 SP5 and MSVC .Net 2003 it may not work on other compilers without modification.
 
 // NOTE These next few lines may be win32 specific, you may need to modify them to compile on other platform