tutorials/medical/viewer.cc

// RedGRID by Orel & Mick
// Parallel VTK Viewer for Medical Tutorial

#include "vtkImageData.h"
#include "vtkPointData.h"
#include "vtkDataArray.h"
#include "vtkUnsignedShortArray.h"
#include "vtkMath.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include "vtkCamera.h"
#include "vtkContourFilter.h"
#include "vtkDataSet.h"
#include "vtkMultiProcessController.h"
#include "vtkMPIController.h"
#include "vtkAppendPolyData.h"
#include "vtkTimerLog.h"
#include "vtkProperty.h"
#include "vtkLookupTable.h"

#include <mpi.h>
#include <cassert>
#include <unistd.h>
#include <iostream>
#include "RedCORBA.hh"
#include "RedSYM.hh"
#include "CorbaTools.hh"
using namespace std;

//===========================================================================//

RedCORBA::Proxy * proxy;;
RedCORBA::Node * node;
RedSYM::Grid * grid;
RedSYM::DistributionInt * distribution; 
int prank, psize;
string sourceID, dataID;
RedSYM::TypeCode rs_typecode;
int vtk_typecode, sizeoftype;
CORBA::ORB_ptr orb;
RedCORBA::ConnectionID connection;
RedCORBA::ProxyChannel * channel;
static const int NODE_TAG=999;
float dx, dy, dz;
float ISO_VALUE = 900.0; // default

//===========================================================================//

// This will be called by all processes
void MyMain( vtkMultiProcessController *controller, void *arg )
{
  vtkImageData * imagedata;
  
  try {
    
    // ==== create node & proxy ==== //    
    if(prank == 0) 
      proxy =  new RedCORBA::Proxy("viewer3D",psize,orb,RedCORBA::_non_blocking); 
    node = new RedCORBA::Node("viewer3D",prank,psize,orb,RedCORBA::_blocking);
    if(prank == 0) proxy->waitActivation();  
    
    // ==== connection  ==== //    
    if(prank == 0) 
      connection = proxy->connect(sourceID);
    node->waitConnection(sourceID); 
    controller->Barrier();
    
    // ==== get remote domain ==== //
    connection = node->getConnectionID(sourceID);
    RedSYM::ParallelGrid * remote_grid = 
      (RedSYM::ParallelGrid *) node->getRemoteParallelObjectFromProxy(connection,dataID,RedSYM::_grid);
    RedSYM::BlockInt domain = remote_grid->getDomain();
    assert(domain.getNumberOfDimensions() == 3);

    // get grid spacing
    double dx = remote_grid->getGridSpacing(0);
    double dy = remote_grid->getGridSpacing(1);
    double dz = remote_grid->getGridSpacing(2);
    cout << "dx = " << dx << "  dy = " << dy << "  dz = " << dz << endl;

    // ==== create distribution ==== //
    
    // define topology
    RedSYM::AxisProc axisproc[3] = {0, 0, 0};
    RedSYM::AxisType axistype[3] = {RedSYM::_collapsed, RedSYM::_collapsed, RedSYM::_collapsed} ;
    RedSYM::AxisInfo axisinfo[3] = {0, 0, 0};
    
    if ((psize%4 == 0)&&(psize > 4)){
      axisproc[0] = psize/4; 
      axisproc[1] = 2; 
      axisproc[2] = 2;
      axistype[0] = RedSYM::_block; 
      axistype[1] = RedSYM::_block; 
      axistype[2] = RedSYM::_block;
    }
    else if ((psize%3 == 0)&&(psize > 3)){
      axisproc[0] = psize/3; 
      axisproc[1] = 3;
      axistype[0] = RedSYM::_block; 
      axistype[1] = RedSYM::_block;
    }
    else if ((psize%2 == 0)&&(psize > 3)){
      axisproc[0] = psize/2; 
      axisproc[1] = 2;
      axistype[0] = RedSYM::_block; 
      axistype[1] = RedSYM::_block;
    }
    else {
      axisproc[0] = psize;
      axistype[0] = RedSYM::_block;
    }

    // compute grid spacing
    // dx = 200.0/(domain[0].upperbound() - domain[0].lowerbound() + 1);
    // dy = 200.0/(domain[1].upperbound() - domain[1].lowerbound() + 1);
    // dz = 200.0/(domain[2].upperbound() - domain[2].lowerbound() + 1);
    // cout << "dx = " << dx << "  dy = " << dy << "  dz = " << dz << endl;

    // the distribution
    distribution = RedSYM::DistributionFactoryInt::createBlockCyclic(domain, 
                                                                     prank, 
                                                                     axisproc, 
                                                                     axistype, 
                                                                     axisinfo);
    
    // overlap
    assert(distribution->size() == 1);
    if ((*distribution)[0][0].lowerbound() != 0) (*distribution)[0][0].lowerbound()--;
    if ((*distribution)[0][1].lowerbound() != 0) (*distribution)[0][1].lowerbound()--;
    if ((*distribution)[0][2].lowerbound() != 0) (*distribution)[0][2].lowerbound()--;
    int x0, x1, y0, y1, z0, z1;
    x0 = (*distribution)[0][0].lowerbound(); 
    x1 = (*distribution)[0][0].upperbound();
    y0 = (*distribution)[0][1].lowerbound(); 
    y1 = (*distribution)[0][1].upperbound();     
    z0 = (*distribution)[0][2].lowerbound(); 
    z1 = (*distribution)[0][2].upperbound(); 

    // the typecode
    rs_typecode = remote_grid->getTypeCode(); 
    
    sizeoftype =  RedSYM::SizeOfType[rs_typecode];
    switch(rs_typecode) {
    case RedSYM::_octet : 
      vtk_typecode = VTK_CHAR;
      break;
    case RedSYM::_short : 
      vtk_typecode = VTK_SHORT;
      break;
    case RedSYM::_long : 
      vtk_typecode = VTK_LONG;
      break;
    case RedSYM::_unsigned_short : 
      vtk_typecode = VTK_UNSIGNED_SHORT;
      break;
    case RedSYM::_unsigned_long :
      vtk_typecode = VTK_UNSIGNED_LONG;
      break;
    case RedSYM::_float : 
      vtk_typecode = VTK_FLOAT;
      break;
    case RedSYM::_double : 
      vtk_typecode = VTK_DOUBLE;
      break;
    default:
      cerr << "unknown type code !!!" << endl;
    }
    
    // ==== image data ==== //    
    imagedata = vtkImageData::New();
    imagedata->SetDimensions(x1-x0+1, y1-y0+1, z1-z0+1)  ;
    imagedata->SetSpacing(dx, dy, dz);
    imagedata->SetOrigin(x0*dx,y0*dy,z0*dz);
    imagedata->SetScalarType(vtk_typecode)  ;
    imagedata->AllocateScalars()  ;
    
    // ==== create grid 3D ==== //
    grid = new RedSYM::Grid(dataID,prank,psize,3); // 3D 
    int serie0 = grid->addSerie("serie0",rs_typecode);
    grid->setDistribution(*distribution);
    grid->wrap(serie0, 0, imagedata->GetScalarPointer());
    node->addData(grid); 
    
    //  ==== set ready ==== //
    if(prank == 0) 
      proxy->ready(RedCORBA::_non_blocking); 
    node->ready(RedCORBA::_blocking); 
    cout << "Node " << prank << " Ready" << endl;
    controller->Barrier();  
    
    // ==== create channel ==== //
    if(prank == 0) {      
      cout << "create channel for " << dataID << "... " << flush;
      channel = proxy->createChannel(connection,dataID);
      cout << "done" << endl;
    }    
    controller->Barrier();  
    
    // get request
    if(prank == 0) {
      cout << "get data from reader... " << flush;
      channel->get(0,RedCORBA::_blocking);
      cout << "done" << endl;
    }
    controller->Barrier(); 
    
  }
  catch(RedCORBA::Exception & ex) { 
    cout << "RedCORBA Exception: " << ex <<  endl; 
    exit(EXIT_FAILURE);
  }
  catch(ColCOWS::ColCOWSException & ex) { 
    cout << "ColCOWS Exception: " << ex <<  endl; 
    exit(EXIT_FAILURE);
  }
  catch(...) { 
    cout << "Unknown Exception!" << endl; exit(EXIT_FAILURE);
  }
  
  // Iso-surface Filter
  vtkContourFilter * iso = vtkContourFilter::New();
  iso->SetInput(imagedata);
  iso->SetValue(0, ISO_VALUE);
  iso->ComputeScalarsOff();
  iso->ComputeGradientsOff();
  
  // ==== slave processes ==== //
  if (prank != 0) {
    controller->Send(iso->GetOutput(), 0, NODE_TAG);
  }

  // ==== master process ==== //
  else {

    vtkActor * isoActors[psize];   
    vtkPolyDataMapper * isoMappers[psize];

    vtkRenderer *ren = vtkRenderer::New();
    vtkRenderWindow *renWindow = vtkRenderWindow::New();
    vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
    vtkTimerLog *timer = vtkTimerLog::New();
    vtkCamera *cam = vtkCamera::New();

    isoMappers[0] = vtkPolyDataMapper::New();
    isoMappers[0]->SetInput(iso->GetOutput());
    isoActors[0] = vtkActor::New();
    isoActors[0]->SetMapper(isoMappers[0]); 
    ren->AddActor(isoActors[0]);
    
    // Collect the data from slaves and render it.
    for (int i = 1; i < psize; ++i) {
      vtkPolyData *pd;
      controller->Receive(pd, i, NODE_TAG);
      isoMappers[i] = vtkPolyDataMapper::New();
      isoMappers[i]->SetInput(pd);
      isoActors[i] = vtkActor::New();
      isoActors[i]->SetMapper(isoMappers[i]);
      ren->AddActor(isoActors[i]);
    }

    // Create Lookup Table
    vtkLookupTable *lookupTable = vtkLookupTable::New();
    lookupTable->SetTableRange (0, psize-1);
    lookupTable->SetHueRange (.70, 0);
    lookupTable->SetSaturationRange (1, 1);
    lookupTable->SetValueRange (1, 1);
    lookupTable->Build();


    // set one color per processor
    for (int i = 0; i < psize; ++i) {
      double procColor[3];
      lookupTable->GetColor(i,procColor);
      isoActors[i]->GetProperty()->SetColor(procColor); 
    }

    
    // Create the rendering part of the pipeline
    iren->SetRenderWindow(renWindow);
    renWindow->AddRenderer(ren);
    ren->SetBackground(0.9, 0.9, 0.9);
    renWindow->SetSize( 400, 400);
    cam->SetFocalPoint(100, 100, 100);
    cam->SetPosition(100, 450, 100);
    cam->SetViewUp(0, 0, -1);
    cam->SetViewAngle(30);
    cam->SetClippingRange(177.0, 536.0);
    ren->SetActiveCamera(cam);

    // Initialize the event loop and then start it.
    iren->Initialize();
    cout << "start parallel iso-surface calculation, wait few seconds..." << endl;
    iren->Start(); 

    // nb triangles
    for (int i = 0; i < psize; ++i) {
      int nb_triangles = isoMappers[i]->GetInput()->GetNumberOfCells();
      cerr << nb_triangles << " triangles on process " << i << endl;
    }
    
    // Tell the other processors to stop processing RMIs.
    for (int i = 1; i < psize; ++i) {
      controller->TriggerRMI(i, vtkMultiProcessController::BREAK_RMI_TAG); 
    }
    
    // Clean up
    ren->Delete();
    renWindow->Delete();
    iren->Delete();
    cam->Delete();
    timer->Delete();
  }
  
  try {
    
    // ==== shutdown ==== //    
    if(prank == 0) {
      proxy->disconnect(connection);
      proxy->deactivate();
    }
    controller->Barrier();  

    // if(prank == 0) delete proxy;
    // delete node;    
    // delete grid;
    // delete distribution;    
  }
  catch(RedCORBA::Exception & ex) { 
    cout << "RedCORBA Exception: " << ex <<  endl; 
    exit(EXIT_FAILURE);
  }
  catch(ColCOWS::ColCOWSException & ex) { 
    cout << "ColCOWS Exception: " << ex <<  endl; 
  }
  catch(...) { 
    cout << "Unknown Exception!" << endl; 
    exit(EXIT_FAILURE);
  }
  
  controller->Barrier();  
}

//===========================================================================//

int main(int argc, char** argv) 
{
  sourceID = "reader3D";
  dataID = "mygrid";

  if (argc >= 2) { 
    ISO_VALUE = atof(argv[1]);
  }
  else {
    cout << "Usage: " << " [iso_value] [ORB options]" << endl;
    return 0;
  }
  
  vtkMultiProcessController *controller;
  controller = vtkMPIController::New();
  controller->Initialize(&argc, &argv);
  
  if (controller->IsA("vtkThreadedController"))
    controller->SetNumberOfProcesses(2);
  
  prank = controller->GetLocalProcessId();
  psize = controller->GetNumberOfProcesses();
  cout << "prank = " << prank << endl;
  cout << "psize = " << psize << endl;
  
  controller->SetSingleMethod(MyMain, reinterpret_cast<void*>(argv[1]));
  
  // ====  start ORB ==== //
  orb = CorbaTools::startORB(argc,argv);
  
  // ====  start Controller ==== //
  controller->SingleMethodExecute();
  controller->Barrier(); 
  
  // ==== goodbye ==== //    
  CorbaTools::killORB();
  controller->Finalize();
  controller->Delete();
  cout << "return OK" << endl;
  
  return 0;
}

//===========================================================================//

// EOF

---