VectorField.cpp

Go to the documentation of this file.

/*=========================================================================

Program:   vtkINRIA3D
Module:    $Id: VectorField.cpp,v 1.6 2011/01/25 23:48:18 hawle Exp $
Language:  C++
Author:    $Author: hawle $
Date:      $Date: 2011/01/25 23:48:18 $
Version:   $Revision: 1.6 $

Copyright (c) 2007 INRIA - Asclepios Project. All rights reserved.
See Copyright.txt for details.

This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#include "VectorField.h"

#include "vtkInteractorObserver.h"
#include "vtkInteractorStyleSwitch.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRendererCollection.h"
#include "vtkProp.h"
#include "vtkTextActor.h"
#include "vtkCoordinate.h"
#include "vtkProperty.h"
#include "vtkProperty2D.h"
#include "vtkTextProperty.h"
#include <vtkCornerAnnotation.h>
#include <vtkLightCollection.h>
#include <vtkLight.h>
#include "vtkOrientationAnnotation.h"

#include "assert.h"
#include <iostream>
#include <sstream>
#include <cmath>

#include "NIREPvtkViewImage2DCommand.h"
#include "vtkViewImage2DFullCommand.h"
#include "NIREPvtkInteractorStyleImage2D.h"
#include "vtkCamera.h"
#include "vtkImageReslice.h"
#include "vtkRenderWindow.h"
#include "vtkTransform.h"
#include "vtkImageMapToWindowLevelColors.h"
#include "vtkImageMapToColors.h"
#include "vtkImageActor.h"
#include "vtkImageData.h"
#include "vtkActor.h"
#include "vtkPolyDataMapper.h"
#include "vtkLineSource.h"
#include "vtkLookupTable.h"
#include "vtkImageBlendWithMask.h"
#include "vtkImageBlend.h"
#include <vtkPolyDataMapper2D.h>
#include <vtkActor2D.h>

#include <vtkColorTransferFunction.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkExtractGeometry.h>
#include <vtkDataSet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkPlane.h>
#include <vtkPlaneCollection.h>
#include <vtkPolyDataMapper2D.h>
#include <vtkActor2D.h>
#include <vtkClipDataSet.h>
#include <vtkCutter.h>
#include <vtkBox.h>
#include <vtkPolyDataWriter.h>
#include <vtkPolyDataNormals.h>
#include <vtkTransform.h>
#include <vtkMath.h>
#include <vtkRulerWidget.h>

//Jeffrey Hawley added this
#include <vtkArrowSource.h>
#include <vtkGlyph3D.h>
#include <vtkVertex.h>
#include <vtkPoints.h>
#include <vtkStructuredGridReader.h>
#include <vtkVoxel.h>
#include <vtkCellArray.h>
#include <vtkDataArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkStructuredGrid.h>
#include <vtkGlyph2D.h>
#include <vtkLookupTableManager.h>
//end of addition

vtkCxxRevisionMacro(VectorField, "$Revision: 1.6 $");
vtkStandardNewMacro(VectorField);


VectorField::VectorField():NIREPvtkViewImage2D()
{
  //Jeffrey Hawley added this
  m_vectorActor = vtkActor::New();
//  InitializeArrays();
  //end of addition

  this->FirstRender = 1;
  this->ShowSliceNumber  = 1;
  this->InteractionStyle = SELECT_INTERACTION;
  this->LeftButtonInteractionStyle = SELECT_INTERACTION;
  this->MiddleButtonInteractionStyle = SELECT_INTERACTION;
  this->RightButtonInteractionStyle = SELECT_INTERACTION;
  this->WheelInteractionStyle = SELECT_INTERACTION;
  this->RulerWidgetVisibility = 0;
  this->Show2DAxis = 1;

  this->Conventions = RADIOLOGIC;
  
  this->InitialParallelScale = 1.0;


  this->ImageReslice     = vtkImageReslice::New();
  this->ImageActor       = NULL;

  vtkRenderer *rend = this->GetRenderer();

  this->WindowLevelForCorner = vtkImageMapToWindowLevelColors::New();
  this->WindowLevel      = vtkImageMapToColors::New();
  this->MaskFilter       = vtkImageBlendWithMask::New();
  this->Blender          = vtkImageBlend::New();

  this->OrientationMatrix = vtkMatrix4x4::New();
  this->OrientationMatrix->Identity();
  this->ScreenToRealWorldMatrix = vtkMatrix4x4::New();
  this->ScreenToRealWorldMatrix->Zero();
  this->ScreenToRealWorldMatrix->SetElement (2, 2, 1.0);

  this->BGActor          = 0;
  this->BGWindowLevel    = 0;
  this->BGImage = 0;

  
  this->HorizontalLineSource = vtkLineSource::New();
  this->VerticalLineSource   = vtkLineSource::New();
  this->HorizontalLineActor  = vtkActor::New();
  this->VerticalLineActor    = vtkActor::New();
  
  this->DataSetCutPlane      = vtkPlane::New();
  this->DataSetCutBox        = vtkBox::New();
  

  this->DataSetCutPlane->SetOrigin (0,0,0);
  this->DataSetCutPlane->SetNormal (0,0,1);
  this->DataSetCutBox->SetBounds (0,0,0,0,0,0);
  this->BoxThickness = 2;

  // set the filters properties
  this->Blender->SetBlendModeToNormal();
  //this->Blender->SetOpacity (0,0.25 ); GEC 7/10/2010 changed to 50% opacity
  //this->Blender->SetOpacity (1, 0.75);
  this->Blender->SetOpacity (0,0.5);
  this->Blender->SetOpacity (1,0.5);
 
  
  
  // set up the vtk pipeline  
  this->ImageReslice->SetOutputDimensionality(2);
  this->ImageReslice->InterpolateOff();
  this->ImageReslice->SetInputConnection( this->WindowLevel->GetOutputPort() );
  

  this->AuxInput     = this->WindowLevel->GetOutput();
  this->ResliceInput = this->WindowLevel->GetOutput();

 
  // Interactor Style
  //this->InitInteractorStyle(VectorField::SELECT_INTERACTION);
  this->SetInteractionStyle( NIREPvtkViewImage2D::SELECT_INTERACTION );

  
  // Initialize cursor lines
  vtkPolyDataMapper* mapper =  vtkPolyDataMapper::New();
  mapper->SetInputConnection( this->HorizontalLineSource->GetOutputPort() );

  this->HorizontalLineActor->SetMapper(mapper);
  this->HorizontalLineActor->GetProperty()->SetColor (1.0,0.0,0.0);
  mapper->Delete();
  this->HorizontalLineActor->SetVisibility (0);
  this->HorizontalLineActor->PickableOff();
  this->HorizontalLineActor->DragableOff();
  

  
  vtkPolyDataMapper* mapper2 =  vtkPolyDataMapper::New();
  mapper2->SetInputConnection( this->VerticalLineSource->GetOutputPort() );
  this->VerticalLineActor->SetMapper(mapper2);
  this->VerticalLineActor->GetProperty()->SetColor (1.0,0.0,0.0);
  mapper2->Delete();
  this->VerticalLineActor->SetVisibility (0);
  this->VerticalLineActor->PickableOff();
  this->VerticalLineActor->DragableOff();

  this->GetCornerAnnotation()->SetWindowLevel ( this->WindowLevelForCorner );
  
  this->SetOrientation (NIREPvtkViewImage::AXIAL_ID);

  if( NIREPvtkViewImage2D::GetViewImage2DDisplayConventions()==0 )
  {
    this->SetConventionsToRadiological();
  }
  else
  {
    this->SetConventionsToNeurological();
  }

  this->RulerWidget = vtkRulerWidget::New();
  this->RulerWidget->KeyPressActivationOff();
  
}


VectorField::~VectorField()
{
}

//void VectorField::InternalUpdate()
//{
//  this->UpdateImageActor();
//  //this->UpdatePosition();
//  this->InitializeImagePositionAndSize();
//  
//  //  this->UpdatePosition();
//  
//}


//void VectorField::Initialize()
//{
//  this->Superclass::Initialize();
//
//  if (this->GetRenderer())
//  {
//    this->GetRenderer()->TwoSidedLightingOff();
//  }
//  
//
//}


void VectorField::UpdateImageActor()
{
  if( !this->GetImage()  )
  {
    return;
  }
  //vtkPolyData *pd = vtkPolyData::New();
  //int * dimensions = this->ImageReslice->GetOutput()->GetDimensions();
  //this->ImageReslice->GetOutput()->SetScalarTypeToDouble();
  //double *scalarPointer = (double*)this->ImageReslice->GetOutput()->GetScalarPointer();
  //vtkPoints *points = vtkPoints::New();
  //for(int i=0; i<dimensions[0]*dimensions[1]*dimensions[2]; i=i+3) {
  //    points->InsertNextPoint(scalarPointer[i],scalarPointer[i+1],scalarPointer[i+2]);
  //}

  //pd->SetPoints(points);
  //
  //vtkArrowSource* as = vtkArrowSource::New();
  //vtkGlyph3D* gl = vtkGlyph3D::New();
  //gl->SetInput(pd);
  //gl->SetSource(as->GetOutput());
  //gl->SetScaleModeToScaleByVector();
  //gl->SetScaleFactor(0.1);
  //vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
  //mapper->SetInputConnection(gl->GetOutputPort());
  //    
  //m_vectorActor->SetMapper(mapper);
  this->FirstRender = 1;
}



//Changes:
//-Will need to do two ImagesReslices
void VectorField::UpdatePosition ()
{
  if( !this->GetImage() )
  {
    return;
  }
  
  double x=0;
  double y=0;
  double max_x=0;
  double max_y=0;
  double min_x=0;
  double min_y=0;
  double pos[4];

  this->GetCurrentPoint(pos);
  pos[3] = 0.0;

  double* spacing = this->GetImage()->GetSpacing();
  double* origin  = this->GetImage()->GetOrigin();

  // check if pos lies inside image bounds
  //if( pos[0]<imBounds[0] || pos[0]>imBounds[1] ||
  //    pos[1]<imBounds[2] || pos[1]>imBounds[3] ||
  //    pos[2]<imBounds[4] || pos[2]>imBounds[5])
  //{
  //  // we are outside image bounds
  //  return;
  //}

// round to the nearest integer slice
  pos[0] = vtkMath::Round ((pos[0]-origin[0])/spacing[0] )*spacing[0]+origin[0];
  pos[1] = vtkMath::Round ((pos[1]-origin[1])/spacing[1] )*spacing[1]+origin[1];
  pos[2] = vtkMath::Round ((pos[2]-origin[2])/spacing[2] )*spacing[2]+origin[2];

  /*
    double *bounds = this->ImageActor->GetBounds();
    min_x = bounds[0];
    max_x = bounds[1];
    min_y = bounds[2];
    max_y = bounds[3];
  */

  double min_bounds[4] = { this->GetWholeMinPosition(0),
                           this->GetWholeMinPosition(1),
                           this->GetWholeMinPosition(2),
                           0.0 };

  double max_bounds[4] = { this->GetWholeMaxPosition(0),
                           this->GetWholeMaxPosition(1),
                           this->GetWholeMaxPosition(2),
                           0.0 };
  

  // new efficient way to get the reslice origin and the axis coordinates.
  double reslice[4], pos2[4];
  this->ScreenToRealWorldMatrix->MultiplyPoint (pos, reslice);
  this->ImageReslice->SetResliceAxesOrigin(reslice[0], reslice[1], reslice[2]);


  vtkMatrix4x4* mat = vtkMatrix4x4::New();
  vtkMatrix4x4::Transpose (this->OrientationMatrix, mat);
  mat->MultiplyPoint (pos, pos2);
  mat->MultiplyPoint (min_bounds, min_bounds);
  mat->MultiplyPoint (max_bounds, max_bounds);
  mat->Delete();

  x = pos2[0];
  y = pos2[1];

  max_x = max_bounds[0];
  max_y = max_bounds[1];
  
  min_x = min_bounds[0];
  min_y = min_bounds[1];

  this->HorizontalLineSource->SetPoint1(min_x, y, 0.001);
  this->HorizontalLineSource->SetPoint2(max_x, y, 0.001);   
  this->VerticalLineSource->SetPoint1(x, min_y, 0.001);   
  this->VerticalLineSource->SetPoint2(x, max_y, 0.001);

  this->ImageReslice->Update(); // needed to update input Extent
 
  int imCoor[3];
  this->GetCurrentVoxelCoordinates(imCoor);
  int dims[3];
  this->GetImage()->GetDimensions (dims);
      
  std::ostringstream os2;
  std::ostringstream os;
  //os << "Slice: ";

  os2 << "Zoom: " << this->GetZoom()*100.0 << " \%\n";
    
  switch( this->Orientation )
  {
      case NIREPvtkViewImage::AXIAL_ID :
        os << "Image size: " << dims[0] << " x " << dims[1] << "\n";
        os << "Voxel size: " << spacing[0] << " x " << spacing[1] << "\n";
        //os << imCoor[2] << " / " << dims[2]-1 << std::endl;
        os << "X: " << imCoor[0] << " px Y: " << imCoor[1] << " px Value: "
           << this->GetCurrentPointDoubleValue() << std::endl;
        os << "X: " << pos[0] << " mm Y: " << pos[1] << " mm\n";
        if ( this->ShowSliceNumber )
          os2 << "Slice: " << imCoor[2] << " / " << dims[2]-1 << std::endl;
        os2 << "Location: " << pos[2] << " mm";
        break;
          
        
      case NIREPvtkViewImage::CORONAL_ID :
        os << "Image size: " << dims[0] << " x " << dims[2] << "\n";
        os << "Voxel size: " << spacing[0] << " x " << spacing[2] << "\n";
        //os << imCoor[1] << " / " << dims[1]-1 << std::endl;
        os << "X: " << imCoor[0] << " px Z: " << imCoor[2] << " px Value: "
           << this->GetCurrentPointDoubleValue() << std::endl;
        os << "X: " << pos[0] << " mm Z: " << pos[2] << " mm\n";
        if ( this->ShowSliceNumber )
          os2 << "Slice: " << imCoor[1] << " / " << dims[1]-1 << std::endl;
        os2 << "Location: " << pos[1] << " mm";
        break;
        
        
      case NIREPvtkViewImage::SAGITTAL_ID :
        os << "Image size: " << dims[1] << " x " << dims[2] << "\n";
        os << "Voxel size: " << spacing[1] << " x " << spacing[2] << "\n";
        //os << imCoor[0] << " / " << dims[0]-1 << std::endl;
        os << "Y: " << imCoor[1] << " px Z: " << imCoor[2] << " px Value: "
           << this->GetCurrentPointDoubleValue() << std::endl;
        os << "Y: " << pos[1] << " mm Z: " << pos[2] << " mm\n";
        if ( this->ShowSliceNumber )
          os2 << "Slice: " << imCoor[0] << " / " << dims[0]-1 << std::endl;
        os2 << "Location: " << pos[0] << " mm";
        break;
  }
  os << "<window_level>";
  
  this->SetUpLeftAnnotation( os.str().c_str() );
  this->SetDownLeftAnnotation( os2.str().c_str() );
  //  }
  //  else
  //    this->SetUpRightAnnotation("");
  
  unsigned int direction = this->GetOrthogonalAxis (this->GetOrientation());
  this->DataSetCutPlane->SetOrigin (reslice[0], reslice[1], reslice[2]);
  
  switch(direction)
  {
      case X_ID :
        //this->DataSetCutPlane->SetNormal (1,0,0);
        this->DataSetCutBox->SetBounds (this->DataSetCutPlane->GetOrigin()[0],this->DataSetCutPlane->GetOrigin()[0]+this->BoxThickness,
                                        this->GetWholeMinPosition(1),this->GetWholeMaxPosition(1),
                                        this->GetWholeMinPosition(2),this->GetWholeMaxPosition(2));

        break;
      case Y_ID :
        //this->DataSetCutPlane->SetNormal (0,1,0);
        this->DataSetCutBox->SetBounds (this->GetWholeMinPosition(0),this->GetWholeMaxPosition(0),
                                        this->DataSetCutPlane->GetOrigin()[1],this->DataSetCutPlane->GetOrigin()[1]+this->BoxThickness,
                                        this->GetWholeMinPosition(2),this->GetWholeMaxPosition(2));
        break;
      case Z_ID :
        //this->DataSetCutPlane->SetNormal (0,0,1);
        this->DataSetCutBox->SetBounds (this->GetWholeMinPosition(0),this->GetWholeMaxPosition(0),
                                        this->GetWholeMinPosition(1),this->GetWholeMaxPosition(1),
                                        this->DataSetCutPlane->GetOrigin()[2],this->DataSetCutPlane->GetOrigin()[2]+this->BoxThickness);
        break;
  }


  if( this->DataSetList.size() )
  {

    this->ResetAndRestablishZoomAndCamera();

    /*
      We need to correct for the origin of the actor. Indeed, the ImageActor
      has always position 0 in Z in axial view, in X in sagittal view and
      in Y in coronal view. The projected dataset have an origin that depends
      on the required slice and can be negative. In that case, the projected
      data are behind the image actor and thus not visible. Here, we correct
      this by translating the actor so that it becomes visible.
     */
    for( unsigned int i=0; i<this->DataSetActorList.size(); i++)
    {
      double Pos[3];
      this->DataSetActorList[i]->GetPosition (Pos);

      switch(direction)
      {
          case X_ID :
            Pos[0] = -1.0*reslice[0] + 1.0;
            break;

          case Y_ID:
            Pos[1] = -1.0*reslice[1] + 1.0;
            break;

          case Z_ID:
            Pos[2] = -1.0*reslice[2] + 1.0;
            break;
      }

      this->DataSetActorList[i]->SetPosition (Pos);
    }

  }

  this->Modified();
}



void VectorField::SetWindow (double w)
{

  if( w<0.0 )
  {
    w = 0.0;
  }

  double shiftScaleWindow = this->GetShift() + w*this->GetScale();
  
  NIREPvtkViewImage::SetWindow ( shiftScaleWindow );
  this->WindowLevelForCorner->SetWindow( shiftScaleWindow );
  
  double v_min = this->GetLevel() - 0.5*this->GetWindow();
  double v_max = this->GetLevel() + 0.5*this->GetWindow();

  if( this->GetLookupTable() && this->WindowLevel->GetLookupTable())
  {
  
    this->GetLookupTable()->SetRange ( (v_min-0.5*this->GetShift())/this->GetScale(),
                                       (v_max-1.5*this->GetShift())/this->GetScale());
    this->WindowLevel->GetLookupTable()->SetRange (v_min, v_max);
  }

  if(this->BGWindowLevel && this->BGWindowLevel->GetLookupTable())
    this->BGWindowLevel->GetLookupTable()->SetRange( v_min, v_max);
}



void VectorField::SetLevel (double l)
{

  double shiftScaleLevel = this->GetShift() + l*this->GetScale();

  NIREPvtkViewImage::SetLevel ( shiftScaleLevel );
  this->WindowLevelForCorner->SetLevel( shiftScaleLevel );
    
  double v_min = this->GetLevel() - 0.5*this->GetWindow();
  double v_max = this->GetLevel() + 0.5*this->GetWindow();
  
  if( this->GetLookupTable() )
  {
    this->GetLookupTable()->SetRange ( (v_min-0.5*this->GetShift())/this->GetScale(),
                                       (v_max-1.5*this->GetShift())/this->GetScale());
    this->WindowLevel->GetLookupTable()->SetRange (v_min, v_max);
  }
}



//Changes:
//-Either need to create another SetImage, or have this SetImage place the images into an array
//-Look into RegisterImage, I might have to do something else.
//-I think I should create another SetImage
void VectorField::SetImage(vtkImageData* image)
{
  //vtkStructuredGridReader *reader = vtkStructuredGridReader::New();
  //reader->SetFileName("subset.vtk");
  //reader->Update();

  //vtkArrowSource* as = vtkArrowSource::New();
  //vtkGlyph3D* gl = vtkGlyph3D::New();
  //gl->SetInputConnection(reader->GetOutputPort());
  //gl->SetSource(as->GetOutput());
  //gl->SetScaleModeToScaleByVector();
  //gl->SetScaleFactor(0.1);

  //vtkLookupTable* lu = vtkLookupTable::New();
  //lu->SetNumberOfTableValues(256);
  //lu->SetHueRange(0.6667, 0);
  //lu->SetTableRange(0, 5);
  //lu->SetVectorMode(vtkScalarsToColors::MAGNITUDE);
  //lu->Build();

  //vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
  //mapper->SetInputConnection(gl->GetOutputPort());
  //mapper->SetLookupTable(lu);
  //gl->Update();
  //m_vectorActor->SetMapper(mapper);

  //this->AddActor(m_vectorActor);

  //NIREPvtkViewImage::SetImage( image );
  if(!image)
  {
    return;
  }

  int* extent = image->GetExtent();
  if( extent[1]<extent[0] || extent[3]<extent[2] || extent[5]<extent[4] )
  {
    vtkErrorMacro ( << "Image extent is not valid: " << extent[0] << " "
                    << extent[1] << " "
                    << extent[2] << " "
                    << extent[3] << " "
                    << extent[4] << " "
                    << extent[5]);
    return;
  }

  

  NIREPvtkViewImage::SetImage( image );

  // check if there is a mask image. If yes, then we check
  // if the new image size and spacing agrees with the mask image.
  // If not, we remove the mask image
  if( this->GetMaskImage() )
  {
    int*    dims        = image->GetDimensions();
    double* spacing     = image->GetSpacing();
    int*    maskDims    = this->GetMaskImage()->GetDimensions();
    double* maskSpacing = this->GetMaskImage()->GetSpacing();

    if( dims[0]!=maskDims[0] || dims[1]!=maskDims[1] || dims[2]!=maskDims[2] ||
        spacing[0]!=maskSpacing[0] || spacing[1]!=maskSpacing[1] || spacing[2]!=maskSpacing[2] )
    {
      this->RemoveMaskImage();
    }
  }

  // should check also the overlapping image
  int type = image->GetScalarType();
  int components = image->GetNumberOfScalarComponents();
  if( image->GetScalarType() == VTK_UNSIGNED_CHAR  && (image->GetNumberOfScalarComponents()==3 || image->GetNumberOfScalarComponents()==4) )
  {
    this->AuxInput = image;
  }
  else if(image->GetScalarType() == VTK_FLOAT && image->GetNumberOfScalarComponents()==3) {
    this->AuxInput = image;
  }
  else
  {
    this->AuxInput = this->WindowLevel->GetOutput();

    this->WindowLevel->SetInput (image);
    double range[2];
    image->GetScalarRange (range);
    if ( this->WindowLevel->GetLookupTable() )
    {
      this->WindowLevel->GetLookupTable()->SetRange (range);
    }
  }


  if( this->GetOverlappingImage() )
  {
    this->Blender->SetInput (0, this->AuxInput );
  }
  else
  {
    if( this->GetMaskImage() )
    {
      this->MaskFilter->SetImageInput ( this->AuxInput );
    }
    else
    {
      this->ImageReslice->SetInput ( this->AuxInput );
      this->ResliceInput = this->AuxInput;
    }
  }

  //this->ImageActor->SetInput( this->ImageReslice->GetOutput() );

  if(this->GetImage())
  {
      this->AddActor(this->HorizontalLineActor);
      this->AddActor(this->VerticalLineActor);

      //this->AddActor(this->ImageActor);
      vtkStructuredGrid *pd = vtkStructuredGrid::New();
      //vtkPolyData *pd = vtkPolyData::New();
      this->ImageReslice->Update();
      int * dimensions = this->ImageReslice->GetOutput()->GetDimensions();

      int numberOfTuples = dimensions[0]*dimensions[1]*dimensions[2] /3;
      this->ImageReslice->GetOutput()->SetScalarTypeToDouble();
      double *scalarPointer = (double*)this->ImageReslice->GetOutput()->GetScalarPointer();
      vtkPoints *vtkpoints = vtkPoints::New();
      vtkpoints->SetNumberOfPoints(numberOfTuples);
      //vtkpoints->SetNumberOfPoints(7);
      vtkCellArray *vertices = vtkCellArray::New();
      vtkIdType pid;


      vtkDoubleArray *dataArray = vtkDoubleArray::New();
      dataArray->SetNumberOfComponents(3);
      dataArray->Allocate(numberOfTuples);
      dataArray->SetName("Momentum");

      vtkDoubleArray *dataArrayTwo = vtkDoubleArray::New();
      dataArrayTwo->SetNumberOfComponents(1);
      dataArrayTwo->Allocate(numberOfTuples);
      dataArrayTwo->SetName("default");

      vtkDoubleArray *dataArrayThree = vtkDoubleArray::New();
      dataArrayThree->SetNumberOfComponents(3);
      dataArrayThree->Allocate(numberOfTuples);

      
      int numberOfComponents = 3;
      double *ptrThree = dataArrayThree->WritePointer(0,numberOfTuples*numberOfComponents);

      int dim0,dim1,dim2;
      dim0=dim1=dim2=0;



      for (int i=0; i<numberOfTuples; i = i++) {
        for (int j=0; j<numberOfComponents; j++) {
          dim0+=50;
          if(dim0 > dimensions[0]*50) {
            dim0 = 0;
            dim1+=50;
          }
          if(dim1 >dimensions[1]*50) {
            dim1 = 0;
            dim2+=50;
          }
          if(dim2 > dimensions[2]*50) {
            dim2 = 0;
          }
          if(j == 0) {
            *ptrThree = dim0;
          }
          else if(j == 1) {
            *ptrThree = dim1;
          }
          else {
            *ptrThree = dim2;
          }
            ptrThree++;
        }
      }


      vtkpoints->SetData(dataArrayThree);

      
      numberOfComponents = 3;
      double *ptr = dataArray->WritePointer(0,numberOfTuples*numberOfComponents);
      for (int i=0; i<numberOfTuples; i = i++) {
        for (int j=0; j<numberOfComponents; j++) {
          *ptr = scalarPointer[(i*numberOfComponents)+j];
          ptr++;
        }
      }

      numberOfComponents = 1;
      double *ptrTwo = dataArrayTwo->WritePointer(0,numberOfTuples*numberOfComponents);
      for (int i=0; i<numberOfTuples; i = i++) {
        for (int j=0; j<numberOfComponents; j++) {
          //*ptrTwo = scalarPointer[(i*numberOfComponents)+j];
          //*ptrTwo = cellScalars[(i*numberOfComponents)+j];
          *ptrTwo = 1;
          ptrTwo++;
        }
      }

      pd->SetPoints(vtkpoints);
      //pd->InsertNextCell(vertex->GetCellType(),vertex->GetPointIds());
      vtkDataSetAttributes *a = pd->GetPointData();
      a->SetVectors(dataArray);
      vtkDataSetAttributes *a2=pd->GetCellData();
      a2->SetScalars(dataArrayTwo);
      //pd->SetVerts(vertices);
  
      vtkArrowSource* as = vtkArrowSource::New();
      //vtkGlyph3D* gl = vtkGlyph3D::New();
      vtkGlyph2D* gl = vtkGlyph2D::New();
      gl->SetInput(pd);
      gl->SetSource(as->GetOutput());
      gl->SetScaleModeToScaleByVector();
      //gl->SetColorModeToColorByVector();
      gl->SetScaleFactor(0.1);

      gl->Update();

      //double range[] = {0,0};

      //gl->GetOutput()->GetPointData()->GetScalars()->GetRange(range,2);

      vtkLookupTable* lu = vtkLookupTable::New();
      lu->SetNumberOfTableValues(256);
      lu->SetHueRange(0.6667, 0);
      //lu->SetTableRange(range);
      lu->SetVectorMode(vtkScalarsToColors::COMPONENT);
      lu->Build();
      this->GetWindowLevel()->SetLookupTable(lu);

      //vtkLookupTable *temp = vtkLookupTableManager::GetLookupTable(2);
      //temp->SetVectorModeToMagnitude();
      //temp->SetRange(range);
      //temp->Build();
      //this->SetLookupTable(temp);

      vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
      mapper->SetInputConnection(gl->GetOutputPort());
      mapper->SetLookupTable(lu);
      //mapper->SetLookupTable(temp);
      //mapper->SelectColorArray("Momentum");
      //mapper->SetUseLookupTableScalarRange(1);
      mapper->SetScalarMode(VTK_SCALAR_MODE_USE_POINT_FIELD_DATA);
      gl->Update();
      m_vectorActor->SetMapper(mapper);

      this->AddActor(m_vectorActor);

      // Save the camera focal and position, and zoom, before calling InternalUpdate()
      double focal[3], pos[3], zoom;
      this->GetCameraFocalAndPosition (focal, pos);
      zoom = this->GetZoom();
      

      this->SetupAnnotations();
      this->InternalUpdate();

      
      // Check that the current displayed point lies within the bounds of the image. If not,
      // we camp it to the nearest acceptable position.
      double *wextent = image->GetBounds();
      double position[3];
      this->GetCurrentPoint ( position );
      position[0]<wextent[0]?position[0]=wextent[0]:position[0];
      position[0]>wextent[1]?position[0]=wextent[1]:position[0];
      position[1]<wextent[2]?position[1]=wextent[2]:position[1];
      position[1]>wextent[3]?position[1]=wextent[3]:position[1];
      position[2]<wextent[4]?position[2]=wextent[4]:position[2];
      position[2]>wextent[5]?position[2]=wextent[5]:position[2];
      this->SetCurrentPoint ( position );


      // restore the zoom and focal in case this is not the first time we render the image.
      if( !this->GetFirstImage() )
      {
        this->SetZoom ( zoom );
        this->SetCameraFocalAndPosition (focal, pos);
        this->SetFirstImage (0);
      }

      if( this->RulerWidgetVisibility )
        this->RulerWidget->On();

  }
}

void VectorField::SetLookupTable (vtkScalarsToColors* lut)
{
    
  if( !lut )
  {
    return;
  }
    
  NIREPvtkViewImage::SetLookupTable (lut);

  double v_min = this->GetLevel() - 0.5*this->GetWindow();
  double v_max = this->GetLevel() + 0.5*this->GetWindow();

  lut->SetRange ( (v_min-0.5*this->GetShift())/this->GetScale(),
                  (v_max-1.5*this->GetShift())/this->GetScale() );

  
  vtkLookupTable* realLut = vtkLookupTable::SafeDownCast (lut);

  if( !realLut )
  {
    std::cerr << "Error: Cannot cast vtkScalarsToColors to vtkLookupTable." << std::endl;
    return;
  }
  
  vtkLookupTable* newLut = vtkLookupTable::New();
  newLut->DeepCopy (realLut);
  newLut->SetRange (v_min, v_max);
  this->WindowLevel->SetLookupTable (newLut);
  newLut->Delete();
  
}


void VectorField::SetOrientation(unsigned int p_orientation)
{
  if( (p_orientation > vtkViewImage::NB_PLAN_IDS - 1) || (this->Orientation == p_orientation) )
  {
    return;
  }

  this->Orientation = p_orientation;

  
  // setup the OrientationMatrix and the ScreenToRealWorldMatrix. Then, we have no
  // need to check the view's orientation to determine what slice to display: a
  // simple matrix multiplication will tell it for us.
  this->OrientationMatrix->Zero();
  this->ScreenToRealWorldMatrix->Zero();

  
  if( this->Orientation == vtkViewImage::SAGITTAL_ID )
  {
    this->ImageReslice->SetResliceAxesDirectionCosines( 0, -1, 0,
                                                        0, 0, 1,
                                                        1, 0, 0);
    this->OrientationMatrix->SetElement (0, 2, 0.0);
    this->OrientationMatrix->SetElement (1, 0, -1.0);
    this->OrientationMatrix->SetElement (2, 1, 1.0);

    this->ScreenToRealWorldMatrix->SetElement (0, 0, 1.0);
  }

  
  if( this->Orientation == vtkViewImage::CORONAL_ID )
  {
    if( this->Conventions == RADIOLOGIC )
    {
      this->OrientationMatrix->SetElement(0, 0, 1.0);
      this->ImageReslice->SetResliceAxesDirectionCosines(1,  0, 0,
                                                         0,  0, 1,
                                                         0,  1, 0); 
    }
    else
    {  
      this->OrientationMatrix->SetElement(0, 0, -1.0);
      this->ImageReslice->SetResliceAxesDirectionCosines(-1,  0, 0,
                                                         0,  0, 1,
                                                         0,  1, 0); 
    }
    
    this->OrientationMatrix->SetElement(2, 1, 1.0);
    this->ScreenToRealWorldMatrix->SetElement (1, 1, 1.0);
  }

  
  if( this->Orientation == vtkViewImage::AXIAL_ID )
  {

    
    if( this->Conventions == RADIOLOGIC )
    {
      this->OrientationMatrix->SetElement(0, 0, 1.0 );
      this->ImageReslice->SetResliceAxesDirectionCosines(1, 0, 0,
                                                         0, 1, 0,
                                                         0, 0, 1);
    }
    else
    {  
      this->OrientationMatrix->SetElement(0, 0, -1.0 );
      this->ImageReslice->SetResliceAxesDirectionCosines(-1, 0, 0,
                                                         0, 1, 0,
                                                         0, 0, 1);
    }

    this->OrientationMatrix->SetElement(1, 1, 1.0 );
    this->ScreenToRealWorldMatrix->SetElement (2, 2, 1.0);
  }

  unsigned int direction = this->GetOrthogonalAxis (this->GetOrientation());
  switch(direction)
  {
      case X_ID :
        this->DataSetCutPlane->SetNormal (1,0,0);
        break;
      case Y_ID :
        this->DataSetCutPlane->SetNormal (0,1,0);
        break;
      case Z_ID :
        this->DataSetCutPlane->SetNormal (0,0,1);
        break;
  }

  //this->ImageReslice->Modified();

  this->UpdateImageActor(); // make sure update extent are up to date
  this->UpdatePosition(); // make sure reslicer origin are ok
  this->InitializeImagePositionAndSize(); // reset and center the view
  this->SetupAnnotations(); // make sure annotations are ok
  
  this->Modified();
  
}



vtkActor* VectorField::AddDataSet (vtkDataSet* dataset,  vtkProperty* property)
{

  bool doit = true;

  if (!dataset)
    doit = false;
  
  if( this->HasDataSet (dataset) )
  {
    doit = false;
  }
  vtkImageData* imagedata = NULL;
  imagedata = vtkImageData::SafeDownCast(dataset);
  
  if (imagedata)
  {
    this->SetImage(imagedata);
  }
  else
  {
    
    if ( !this->GetImage() )
    {
      doit = false;
    }
    // don't constrain the memory of input datasets anymore.
    if (doit)
    {
      
      vtkMatrix4x4* matrix = vtkMatrix4x4::New();
      for (unsigned int i=0; i<3; i++)
      {
        for (unsigned int j=0; j<3; j++)
        {
          matrix->SetElement(i,j,this->ImageReslice->GetResliceAxes()->GetElement(j,i));
        }
        matrix->SetElement(i,3,0);
      }
      matrix->SetElement(3,3,1);
      
      vtkCutter* cutter = vtkCutter::New();
      cutter->SetCutFunction (this->DataSetCutPlane);

      // Very strangely in some cases (ex : landmarks)
      // the cutter increments the RefCount of the input dataset by 2
      // making some memory leek...
      // I could not manage to know what is wrong here
      
      cutter->SetInput (dataset);
      cutter->Update();
      
      if (!cutter->GetOutput())
      {
        vtkWarningMacro(<< "Unable to cut this dataset...");
        matrix->Delete();
        cutter->Delete();
        return NULL;
      }
      
      
      vtkPolyDataMapper* mapper = vtkPolyDataMapper::New();
      mapper->SetInput (cutter->GetOutput());
      
      vtkActor* actor = vtkActor::New();
      actor->SetUserMatrix (matrix);
      actor->SetMapper (mapper);
      if (property)
      {
        actor->SetProperty (property);
      }

      actor->PickableOff();
      
      this->AddActor (actor);
      this->DataSetList.push_back (dataset);
      this->DataSetActorList.push_back (actor);

      this->ResetAndRestablishZoomAndCamera();

      actor->Delete();
      mapper->Delete();
      matrix->Delete();
      cutter->Delete();
    }
  }

  
  return this->GetDataSetActor(dataset);
  
}


vtkActor* VectorField::SyncAddPolyData (vtkPolyData* polydata,  vtkProperty* property, double thickness)
{

  if( this->IsLocked() )
  {
    return NULL;
  }
  

  vtkActor* actor = this->AddPolyData (polydata, property, thickness);
  
  this->Lock();
  for( unsigned int i=0; i<this->Children.size(); i++)
  {
//     VectorField* view = dynamic_cast<VectorField*> (this->Children[i]);
    VectorField* view = VectorField::SafeDownCast (this->Children[i]);
    if( view )
    {
      view->SyncAddPolyData (polydata, property, thickness);
    }
  }
  this->UnLock();

  return actor;
  
  
}