vtkThresholdingLineFactory.cpp

#include "MantidVatesAPI/vtkThresholdingLineFactory.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkFloatArray.h"
#include "vtkSmartPointer.h"
#include "vtkLine.h"
#include "MantidGeometry/MDGeometry/Coordinate.h"
#include "MDDataObjects/MDIndexCalculator.h"
#include <vector>
#include <boost/math/special_functions/fpclassify.hpp> 
#include "MantidAPI/IMDWorkspace.h"

using Mantid::API::IMDWorkspace;

namespace Mantid
{

  namespace VATES
  {

    vtkThresholdingLineFactory::vtkThresholdingLineFactory(ThresholdRange_scptr thresholdRange, const std::string& scalarName) : m_scalarName(scalarName),
      m_thresholdRange(thresholdRange)
    {
    }

      /**
  Assigment operator
  @param other : vtkThresholdingLineFactory to assign to this instance from.
  @return ref to assigned current instance.
  */
  vtkThresholdingLineFactory& vtkThresholdingLineFactory::operator=(const vtkThresholdingLineFactory& other)
  {
    if(this != &other)
    {
      this->m_scalarName = other.m_scalarName;
      this->m_thresholdRange = other.m_thresholdRange;
      this->m_workspace = other.m_workspace;
    }
    return *this;
  }

  /**
  Copy Constructor
  @param other : instance to copy from.
  */
  vtkThresholdingLineFactory::vtkThresholdingLineFactory(const vtkThresholdingLineFactory& other)
  {
   this->m_scalarName = other.m_scalarName;
   this->m_thresholdRange = other.m_thresholdRange;
   this->m_workspace = other.m_workspace;
  }

    vtkDataSet* vtkThresholdingLineFactory::create() const
    {
      validate();
      //use the successor factory's creation method if this type cannot handle the dimensionality of the workspace.
      size_t nonIntegratedSize = m_workspace->getNonIntegratedDimensions().size();
      if(nonIntegratedSize != vtkDataSetFactory::OneDimensional)
      {
        return m_successor->create();
      }
      else
      {
        const int nBinsX = static_cast<int>( m_workspace->getXDimension()->getNBins() );

        const double maxX = m_workspace-> getXDimension()->getMaximum();
        const double minX = m_workspace-> getXDimension()->getMinimum();

        double incrementX = (maxX - minX) / (nBinsX-1);

        const int imageSize = nBinsX;
        vtkPoints *points = vtkPoints::New();
        points->Allocate(static_cast<int>(imageSize));

        vtkFloatArray * signal = vtkFloatArray::New();
        signal->Allocate(imageSize);
        signal->SetName(m_scalarName.c_str());
        signal->SetNumberOfComponents(1);

        //The following represent actual calculated positions.
        double posX;

        UnstructuredPoint unstructPoint;
        float signalScalar;
        const int nPointsX = nBinsX;
        Column column(nPointsX);

        signal_t maxThreshold = m_thresholdRange->getMaximum();
        signal_t minThreshold = m_thresholdRange->getMinimum();

        //Loop through dimensions
        for (int i = 0; i < nPointsX; i++)
        {
          posX = minX + (i * incrementX); //Calculate increment in x;
          
          signalScalar = static_cast<float>(m_workspace->getSignalNormalizedAt(i));

            if (boost::math::isnan( signalScalar ) || (signalScalar < minThreshold) || (signalScalar > maxThreshold))
            {
              //Flagged so that topological and scalar data is not applied.
              unstructPoint.isSparse = true;
            }
            else
            {
              if (i < (nBinsX -1))
              {
                signal->InsertNextValue(static_cast<float>(signalScalar));
              }
              unstructPoint.isSparse = false;
            }
            unstructPoint.pointId = points->InsertNextPoint(posX, 0, 0);
            column[i] = unstructPoint;

          }
        
        points->Squeeze();
        signal->Squeeze();

        vtkUnstructuredGrid *visualDataSet = vtkUnstructuredGrid::New();
        visualDataSet->Allocate(imageSize);
        visualDataSet->SetPoints(points);
        visualDataSet->GetCellData()->SetScalars(signal);

        for (int i = 0; i < nBinsX - 1; i++)
        {
            //Only create topologies for those cells which are not sparse.
            if (!column[i].isSparse)
            {
              vtkLine* line = vtkLine::New();
              line->GetPointIds()->SetId(0, column[i].pointId);
              line->GetPointIds()->SetId(1, column[i + 1].pointId);
              visualDataSet->InsertNextCell(VTK_LINE, line->GetPointIds());
            }
        }

        points->Delete();
        signal->Delete();
        visualDataSet->Squeeze();
        return visualDataSet;
      }
    }

    vtkDataSet* vtkThresholdingLineFactory::createMeshOnly() const
    {
      throw std::runtime_error("::createMeshOnly() does not apply for this type of factory.");
    }

    vtkFloatArray* vtkThresholdingLineFactory::createScalarArray() const
    {
      throw std::runtime_error("::createScalarArray() does not apply for this type of factory.");
    }

    void vtkThresholdingLineFactory::initialize(Mantid::API::Workspace_sptr wspace_sptr)
    {
      m_workspace = boost::dynamic_pointer_cast<IMDWorkspace>(wspace_sptr);
      validate();
      // When the workspace can not be handled by this type, take action in the form of delegation.
      size_t nonIntegratedSize = m_workspace->getNonIntegratedDimensions().size();
      if(nonIntegratedSize != vtkDataSetFactory::OneDimensional)
      {
        if(this->hasSuccessor())
        {
          m_successor->initialize(m_workspace);
          return;
        }
        else
        {
          throw std::runtime_error("There is no successor factory set for this vtkThresholdingLineFactory type");
        }
      }

      //Setup range values according to whatever strategy object has been injected.
      m_thresholdRange->setWorkspace(m_workspace);
      m_thresholdRange->calculate();
    }

    void vtkThresholdingLineFactory::validate() const
    {
      if(NULL == m_workspace.get())
      {
        throw std::runtime_error("IMDWorkspace is null");
      }
    }

    vtkThresholdingLineFactory::~vtkThresholdingLineFactory()
    {

    }
  }
}