Merge topic '0016736-vtkMPIImage-2ToPwer31-elements'

27eb09b1 VTK Bug #16736: Add large data support&test for vtkMPIImageReader

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !1683

IO/MPIImage/Testing/Cxx/CMakeLists.txt

@@ -3,5 +3,6 @@ include(vtkMPI)
 vtk_add_test_mpi(${vtk-module}CxxTests-MPI tests
   TESTING_DATA
   ParallelIso.cxx
+  ParallelIso2.cxx
   )
 vtk_test_mpi_executable(${vtk-module}CxxTests-MPI tests)

IO/MPIImage/Testing/Cxx/ParallelIso.cxx

@@ -44,6 +44,8 @@
 
 #include "vtkDebugLeaks.h"
 
+namespace {
+
 static const float ISO_START=4250.0;
 static const float ISO_STEP=-1250.0;
 static const int ISO_NUM=3;
@@ -231,6 +233,8 @@ void MyMain( vtkMultiProcessController *controller, void *arg )
 }
 
 
+} // end anon namespace
+
 int ParallelIso( int argc, char* argv[] )
 {
   // This is here to avoid false leak messages from vtkDebugLeaks when

IO/MPIImage/Testing/Cxx/ParallelIso2.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    ParallelIso2.cxx
+
+  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+  All rights reserved.
+  See Copyright.txt or http://www.kitware.com/Copyright.htm 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 notice for more information.
+
+=========================================================================*/
+// This example demonstrates the use of data parallelism in VTK. The
+// pipeline ( vtkMPIImageReader -> vtkContourFilter -> vtkElevationFilter )
+// is created in parallel and each process is assigned 1 piece to process.
+// All satellite processes send the result to the first process which
+// collects and renders them.
+
+#include <mpi.h>
+
+#include "vtkActor.h"
+#include "vtkAppendPolyData.h"
+#include "vtkCamera.h"
+#include "vtkConeSource.h"
+#include "vtkContourFilter.h"
+#include "vtkDataSet.h"
+#include "vtkElevationFilter.h"
+#include "vtkMath.h"
+#include "vtkMPIController.h"
+#include "vtkMPIImageReader.h"
+#include "vtkPolyData.h"
+#include "vtkPolyDataMapper.h"
+#include "vtkTestUtilities.h"
+#include "vtkRegressionTestImage.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderer.h"
+#include "vtkWindowToImageFilter.h"
+#include "vtkImageData.h"
+#include "vtkStreamingDemandDrivenPipeline.h"
+#include "vtkInformation.h"
+
+#include "vtkDebugLeaks.h"
+
+namespace {
+
+static const float ISO_START=4250.0;
+static const float ISO_STEP=-1250.0;
+static const int ISO_NUM=3;
+// Just pick a tag which is available
+static const int ISO_VALUE_RMI_TAG=300;
+static const int ISO_OUTPUT_TAG=301;
+
+struct ParallelIsoArgs_tmp
+{
+  int* retVal;
+  int argc;
+  char** argv;
+};
+
+struct ParallelIsoRMIArgs_tmp
+{
+  vtkContourFilter* ContourFilter;
+  vtkMultiProcessController* Controller;
+  vtkElevationFilter* Elevation;
+};
+
+// call back to set the iso surface value.
+void SetIsoValueRMI(void *localArg, void* vtkNotUsed(remoteArg),
+                    int vtkNotUsed(remoteArgLen), int vtkNotUsed(id))
+{
+  ParallelIsoRMIArgs_tmp* args = (ParallelIsoRMIArgs_tmp*)localArg;
+
+  float val;
+
+  vtkMultiProcessController* contrl = args->Controller;
+  int myid = contrl->GetLocalProcessId();
+  int numProcs = contrl->GetNumberOfProcesses();
+
+  vtkContourFilter *iso = args->ContourFilter;
+  val = iso->GetValue(0);
+  iso->SetValue(0, val + ISO_STEP);
+  args->Elevation->UpdatePiece(myid, numProcs, 0);
+
+  contrl->Send(args->Elevation->GetOutput(), 0, ISO_OUTPUT_TAG);
+}
+
+
+// This will be called by all processes
+void MyMain( vtkMultiProcessController *controller, void *arg )
+{
+  vtkMPIImageReader *reader;
+  vtkContourFilter *iso;
+  vtkElevationFilter *elev;
+  int myid, numProcs;
+  float val;
+  ParallelIsoArgs_tmp* args = reinterpret_cast<ParallelIsoArgs_tmp*>(arg);
+
+  // Obtain the id of the running process and the total
+  // number of processes
+  myid = controller->GetLocalProcessId();
+  numProcs = controller->GetNumberOfProcesses();
+
+  // Create the reader, the data file name might have
+  // to be changed depending on where the data files are.
+  char* fname = vtkTestUtilities::ExpandDataFileName(args->argc, args->argv,
+                                                   "Data/headsq/quarter");
+  reader = vtkMPIImageReader::New();
+  reader->SetDataByteOrderToLittleEndian();
+  reader->SetDataExtent(0, 63, 0, 63, 1, 93);
+  reader->SetFilePrefix(fname);
+  reader->SetDataSpacing(3.2, 3.2, 1.5);
+  delete[] fname;
+
+
+  // Iso-surface.
+  iso = vtkContourFilter::New();
+  iso->SetInputConnection(reader->GetOutputPort());
+  iso->SetValue(0, ISO_START);
+  iso->ComputeScalarsOff();
+  iso->ComputeGradientsOff();
+
+  // Compute a different color for each process.
+  elev = vtkElevationFilter::New();
+  elev->SetInputConnection(iso->GetOutputPort());
+  val = (myid+1) / static_cast<float>(numProcs);
+  elev->SetScalarRange(val, val+0.001);
+
+  // Make sure all processes update at the same time.
+  elev->UpdatePiece(myid, numProcs, 0);
+
+  if (myid != 0)
+    {
+    // If I am not the root process
+    ParallelIsoRMIArgs_tmp args2;
+    args2.ContourFilter = iso;
+    args2.Controller = controller;
+    args2.Elevation = elev;
+
+    // Last, set up a RMI call back to change the iso surface value.
+    // This is done so that the root process can let this process
+    // know that it wants the contour value to change.
+    controller->AddRMI(SetIsoValueRMI, (void *)&args2, ISO_VALUE_RMI_TAG);
+    controller->ProcessRMIs();
+    }
+  else
+    {
+    // Create the rendering part of the pipeline
+    vtkAppendPolyData *app = vtkAppendPolyData::New();
+    vtkRenderer *ren = vtkRenderer::New();
+    vtkRenderWindow *renWindow = vtkRenderWindow::New();
+    vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+    vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
+    vtkActor *actor = vtkActor::New();
+    vtkCamera *cam = vtkCamera::New();
+    renWindow->AddRenderer(ren);
+    iren->SetRenderWindow(renWindow);
+    ren->SetBackground(0.9, 0.9, 0.9);
+    renWindow->SetSize( 400, 400);
+    mapper->SetInputConnection(app->GetOutputPort());
+    actor->SetMapper(mapper);
+    ren->AddActor(actor);
+    cam->SetFocalPoint(100, 100, 65);
+    cam->SetPosition(100, 450, 65);
+    cam->SetViewUp(0, 0, -1);
+    cam->SetViewAngle(30);
+    cam->SetClippingRange(177.0, 536.0);
+    ren->SetActiveCamera(cam);
+
+    // loop through some iso surface values.
+    for (int j = 0; j < ISO_NUM; ++j)
+      {
+      for (int i = 1; i < numProcs; ++i)
+        {
+        // trigger the RMI to change the iso surface value.
+        controller->TriggerRMI(i, ISO_VALUE_RMI_TAG);
+        }
+
+      // set the local value
+      iso->SetValue(0, iso->GetValue(0) + ISO_STEP);
+      elev->UpdatePiece(myid, numProcs, 0);
+
+      for (int i = 1; i < numProcs; ++i)
+        {
+        vtkPolyData* pd = vtkPolyData::New();
+        controller->Receive(pd, i, ISO_OUTPUT_TAG);
+        if (j == ISO_NUM - 1)
+          {
+          app->AddInputData(pd);
+          }
+        pd->Delete();
+        }
+      }
+
+    // Tell the other processors to stop processing RMIs.
+    for (int i = 1; i < numProcs; ++i)
+      {
+      controller->TriggerRMI(i, vtkMultiProcessController::BREAK_RMI_TAG);
+      }
+
+    vtkPolyData* outputCopy = vtkPolyData::New();
+    outputCopy->ShallowCopy(elev->GetOutput());
+    app->AddInputData(outputCopy);
+    outputCopy->Delete();
+    app->Update();
+    outputCopy->RemoveGhostCells();
+    renWindow->Render();
+
+    *(args->retVal) =
+      vtkRegressionTester::Test(args->argc, args->argv, renWindow, 10);
+
+    if ( *(args->retVal) == vtkRegressionTester::DO_INTERACTOR)
+      {
+      iren->Start();
+      }
+
+    // Clean up
+    app->Delete();
+    ren->Delete();
+    renWindow->Delete();
+    iren->Delete();
+    mapper->Delete();
+    actor->Delete();
+    cam->Delete();
+    }
+
+  // clean up objects in all processes.
+  reader->Delete();
+  iso->Delete();
+  elev->Delete();
+}
+
+
+} // end anon namespace
+
+int ParallelIso2( int argc, char* argv[] )
+{
+  // This is here to avoid false leak messages from vtkDebugLeaks when
+  // using mpich. It appears that the root process which spawns all the
+  // main processes waits in MPI_Init() and calls exit() when
+  // the others are done, causing apparent memory leaks for any objects
+  // created before MPI_Init().
+  MPI_Init(&argc, &argv);
+
+  // Note that this will create a vtkMPIController if MPI
+  // is configured, vtkThreadedController otherwise.
+  vtkMPIController* controller = vtkMPIController::New();
+
+  controller->Initialize(&argc, &argv, 1);
+
+  // Added for regression test.
+  // ----------------------------------------------
+  int retVal = 1;
+  ParallelIsoArgs_tmp args;
+  args.retVal = &retVal;
+  args.argc = argc;
+  args.argv = argv;
+  // ----------------------------------------------
+
+  controller->SetSingleMethod(MyMain, &args);
+  controller->SingleMethodExecute();
+
+  controller->Finalize();
+  controller->Delete();
+
+  return !retVal;
+}

IO/MPIImage/Testing/Data/Baseline/ParallelIso2.png.md5

+4ceb35879beb33c50c5c58fc02453dc9

IO/MPIImage/vtkMPIImageReader.cxx

@@ -52,6 +52,7 @@
 #include "vtkPointData.h"
 #include "vtkTransform.h"
 #include "vtkType.h"
+#include <algorithm>
 
 // This macro can be wrapped around MPI function calls to easily report errors.
 // Reporting errors is more important with file I/O because, unlike network I/O,
@@ -295,16 +296,25 @@ void vtkMPIImageReader::ReadSlice(int slice, const int extent[6], void *buffer)
   length *= extent[3]-extent[2]+1;
   if (this->GetFileDimensionality() == 3) length *= extent[5]-extent[4]+1;
 
-  if (length > VTK_INT_MAX)
+  vtkIdType pos = 0;
+  while (length > pos)
     {
-    vtkErrorMacro(<< "Cannot read more than " << VTK_INT_MAX
-                  << " bytes at a time.");
+    MPI_Status stat;
+    // we know this will fit in an int because it can't exceed VTK_INT_MAX.
+   const int remaining = static_cast<int>(std::min(length - pos,
+                                  static_cast<vtkIdType>(VTK_INT_MAX)));
+   MPICall(MPI_File_read(file.Handle, (static_cast<char*>(buffer)) + pos, remaining,
+                         MPI_BYTE, &stat));
+   int rd = 0;
+   MPICall(MPI_Get_elements(&stat, MPI_BYTE, &rd));
+   if (MPI_UNDEFINED == rd)
+      {
+      vtkErrorMacro("Error obtaining number of values read in " << remaining <<
+      "-byte read.");
+      }
+      pos += static_cast<vtkIdType>(rd);
     }
 
-  // Do the read.  This is a coordinated parallel operation for efficiency.
-  MPICall(MPI_File_read_all(file.Handle, buffer, static_cast<int>(length),
-                            MPI_BYTE, MPI_STATUS_IGNORE));
-
   MPICall(MPI_File_close(&file.Handle));
 }
 #else // VTK_USE_MPI_IO