[vtkusers] Re: MarchingCubes

Wed Jul 28 18:15:33 EDT 2004

Sean McInerney wrote:
> 
> Jonathan,
> 
>    That's a cool idea!
> 
>    Foreach desired slice plane for the ImageData volume that you will
> ultimately produce, do something like this:
> 
> 1.) Use vtkClipPolyData to get the 2D contour of the given slice.
> 
> 2.) Use vtkPolyDataToImageStencil to produce ImageStencilData
> 
> 3.) Use vtkImageStencil to produce the ImageData
>      (need to do a little filtering to produce the desired binary image).
> 
> 4.) Use vtkAppendFilter to fuse the 2D binary images into a 3D volume.
>      Hopefully, the ImageStencilData retains a 3D extent.
> 
> 5.) Voila! Binary volume from PolyData ... maybe.
> 
> -Sean

Thanks for the help!

If someone is inspirated to actually implement this (or has already
completed this), please notify me!
Anyway, here is a stupid but working solution. Code is sometimes ugly,
but it works (only voxels at border of each slice are "colored"). Output
is AIR analyse format, recognized on ITK and usable on VTK if you know
the parameters and consider the .img file as the raw volume.

All personal includes in:
http://www.greyc.ismra.fr/~bailleul/Ressources/summary.html (see gpl
code section)
for those who really want to end up testing the program.

I really should have been more investigating on VTK abilities before
doing all this by hand, but I was VTK newbie when I coded this and had
most of the necessary code at hand on my personal libraries.

#include "vtkDataArray.h"
#include "vtkFloatArray.h"
#include "vtkPointData.h"
#include "vtkSphereSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkActor2D.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPointPicker.h"
#include "vtkCellPicker.h"
#include "vtkPicker.h"
#include "vtkPolyDataReader.h"
#include "vtkPolyData.h"
#include "vtkProperty.h"
#include "vtkCamera.h"
#include "vtkCommand.h"
#include "vtkCallbackCommand.h"
#include "vtkIdList.h"
#include "vtkConeSource.h"
#include "vtkSphereSource.h"
#include "vtkGlyph3D.h"
#include "vtkCellArray.h"
#include "vtkTextSource.h"
#include "vtkTextProperty.h"
#include "vtkTextMapper.h"
#include "vtkVectorText.h"
#include "vtkCoordinate.h"
#include "vtkPolyDataReader.h"
#include "vtkPolyDataMapper.h"
#include "vtkPolyDataNormals.h"

// these are "personal" includes
#define DATATYPE unsigned short int 
#include "macros.h"
#include "AIR_Volume.h"

#include "Point.h"
#include "Triangle.h"
#include "Vector.h"
#include "Plane.h"
#include "Buffer_vector.h"

#include "Tools.h"

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <string>

#define DEBUG_MODE false
#define DEBUG2_MODE false
#define SLOTTYPE Segment<double>

#define LANDK_COLOR_START 150
#define OUTL_COLOR 200

/* GLOS (Graphic Library in Open Source), an ANSI Common Lisp OpenGL
subset.
   Copyright (C) 2000 the GLOS development team
(http://glos.sourceforge.net) */

/* global variables - due to function pointers */
AIR_Volume<DATATYPE> g_out_volume;
int g_color, g_current_z;

static void
usage(int argc, char **argv) 
{
  if ((argc - 1) < 3)
    {
      printf("Traces the outline of a polygonal .vtk object file into
given AIR volume\n");
      printf("Usage: %s <vtk object filename> <output AIR volume>
<outline color / input AIR volume> <opt: normals half-length> <normals
tag>\n", argv[0]);
      printf("AIR volume: if only .hdr exists, builds .img ; otherwise
loads .img and write onto it.\n");
      printf("normals half-length: if provided, traces recomputed
normals centered at vertices and of length 1+2*HL \n");
      printf("if AIR volume provided instead of color, traced voxels
take their values from input volume \n");
      printf("normals tag: 0 (default): landmarks = 0, normals of color
the ldk index. 1: landks = landk idx, normals of bright different
colors. 2: only normals are displayed, no polygon! \n");
      exit(1);
    }  
}

void 
PrintPoint(double* p, ostream &stream)
{
  stream << " [" << p[0] << "; " << p[1] << "; " << p[2] << "] "; 
}

/* local function called by bresenham2d. Uses globals for dynamic */
void 
writePoint(int x, int y)
{
  g_out_volume.write_Point_std(Point<int>(x, y, g_current_z),
(DATATYPE)g_color);
}

/*
--------------------------------------------------------------------------
*/
/* Main program */
/*
--------------------------------------------------------------------------
*/

int 
main(int argc, char* argv[])
{
  usage(argc, argv);

  // Air volume 
  int dim_X, dim_Y, dim_Z;
  AIR_Volume<DATATYPE> in_volume;
  BOOL fixed_color_mode = TRUE;
  BOOL only_normals_mode = FALSE;
  BOOL normal_ldkidx_mode = FALSE;
  BOOL normal_bright_mode = FALSE;

  BOOL go_on;
  vtkIdType nb_cell_pts;
  vtkIdType* cell_pts_ids; // array of points identifiers

  //vtkPoints *poly_points; // points of polydata object
  vtkPolyData* polydata;
  int nb_points = 0, triangles_number = 0;
  int slice;
  int normal_length = 0;

  // vtk vertices
  //  double *A, *B, *C;
  Point<double> A, B, C;
  A.set_read_flip_xy(true);
  B.set_read_flip_xy(true);
  C.set_read_flip_xy(true);

  Triangle<double> triangle;
  Plane<double> plane;
  Vector<double> N(0.0, 0.0, 1.0);

  // arguments parsing
  char* vtk_filename = argv[1];
  char* air_filename = argv[2];
  in_volume.change_files(argv[3]); 
  fixed_color_mode = (! in_volume.get_VOL_files().exists_hdr_p());
  g_color = (fixed_color_mode) ? atoi(argv[3]) : OUTL_COLOR;
  if ((argc - 1) >= 4) {normal_length = atoi(argv[4]);}
  if ((argc - 1) == 5) 
    switch (atoi(argv[5]))
      {
      case 1: normal_bright_mode = TRUE; break;
      case 2: only_normals_mode = TRUE; break;
      default: normal_ldkidx_mode = TRUE;
      }
  else
    normal_ldkidx_mode = TRUE;

  // Air volume operations
  g_out_volume.change_files(air_filename);
  assert(g_out_volume.get_VOL_files().exists_hdr_p());
  g_out_volume.update_definitions();

  dim_X = g_out_volume.get_x_dim(); dim_Y = g_out_volume.get_y_dim(); 
  dim_Z = g_out_volume.get_z_dim();

  // array of Buffer_vector of intersection segments
  Buffer_vector<SLOTTYPE>* seg_vectors = new
Buffer_vector<SLOTTYPE>[dim_Z](10);

  g_out_volume.allocate(); 
  // if an img file exists, load it and write into. Otherwise, build it.
  if (g_out_volume.get_VOL_files().exists_hdr_p())
    g_out_volume.update_data();

  // vtk file reading  
  vtkPolyDataReader *reader = vtkPolyDataReader::New();
  reader -> SetFileName(vtk_filename);
  reader -> Update();

  polydata = reader -> GetOutput();  
  //poly_points = polydata -> GetPoints();
  nb_points = polydata -> GetNumberOfPoints();

  vtkCellArray* triangles = polydata -> GetPolys();
  triangles_number = polydata -> GetNumberOfPolys();

  //vtkIdList *ptsIdList = vtkIdList::New();

  /* for each volume plane, store (without order) segments resulting
     from intersection with any triangle */
  if (triangles_number && (! only_normals_mode))
    {

      vtkPolyDataMapper* poly_map = vtkPolyDataMapper::New();
      poly_map -> SetInput(polydata);
      double bounds[6]; poly_map -> GetBounds(bounds); 

      polydata -> BuildCells();

      // Process in turn every z-plane between Z-bounds
      for (slice = int(bounds[4] - 3); slice < int(bounds[5] + 3);
slice++) // for (int slice = 44; slice <= 44; slice++)
	{
	  //cout << " * " << bounds[4] << " : " << bounds[5];
	  // current plane cuts through middle of current voxel Z-slice
	  plane.reset(Point<double>(0.0, 0.0, double(slice + 0.5)), N);      

	  // Process in turn all triangles of the vtk shape
	  if (DEBUG2_MODE) {cout << endl << "plane: " << plane;}

	  for (vtkIdType t = 0; t < triangles_number; t++)
	    {	 
	      /** changed for vtk 4.4. warning, if data not hard
		  stored, pointers A/B/C point value of last read C! */
	      polydata -> GetCellPoints(t, nb_cell_pts, cell_pts_ids);
	      assert(nb_cell_pts == 3); // only triangles expected
	      A.set(polydata -> GetPoint(cell_pts_ids[0]));
	      B.set(polydata -> GetPoint(cell_pts_ids[1]));
	      C.set(polydata -> GetPoint(cell_pts_ids[2])); 

	      // current triangle is rebuilt from its vertices
	      // triangle.reset(Point<double>(A, FALSE), Point<double>(B,
FALSE), Point<double>(C, FALSE));
	      triangle.reset(A, B, C);

	      // for current plane, process on this triangle
	      if (triangle.intersects_plane_p(plane))
		{
		  if (DEBUG_MODE)
		    cout << endl << "slice: " << slice << "inter! triangle:" <<
triangle 
			 << "segt: " << triangle.intersection_plane(plane);

		  seg_vectors[slice].add_element(triangle.intersection_plane(plane));
		}
	    }
	}

      /* Run through saved segments and trace lines into outut volumes
*/
      // NB: devised BEFORE adjunction of bresenham inside AIR_VOL
      {
	g_out_volume.set_pen_color(g_color);

	int count = 0;
	int x1, y1, x2, y2;
	SLOTTYPE cur_segt;
	for (g_current_z = 0; g_current_z < dim_Z; g_current_z++)
	  {
	    count = seg_vectors[g_current_z].get_elements();	    
	    for (int segt = 0; segt < count; segt++)
	      {
		cur_segt = seg_vectors[g_current_z][segt];
		g_out_volume.write_Segment(cur_segt.round_coords());
	      }
	  }
      }
    }

  // if we also want to trace normals
  if (normal_length > 0)
    {
      Point<double> readpt; 
      Point<int> roundedpt; 
      Vector<double> readnormal;

      readpt.set_read_flip_xy(true);
      readnormal.set_read_flip_xy(true);

      /* (re)compute normals and make them available in an array */
      vtkPolyDataNormals* normals = vtkPolyDataNormals::New();
      normals -> SetInput(polydata);
      normals -> ConsistencyOn();
      normals -> SplittingOff();
      normals -> Update();

      vtkDataArray* norm1; // = vtkDataArray::New();
      norm1 = normals -> GetOutput() -> GetPointData() -> GetNormals();
      norm1 -> SetNumberOfTuples(nb_points);
      norm1 -> SetNumberOfComponents(3);

      // run through every (landmark) vertices
      for (vtkIdType p = 0; p < nb_points; p++)
	{
	  readpt.set(polydata -> GetPoint(p)); //set_vtk
	  cout << endl << p << ":" << readpt;
	  roundedpt = readpt.round_coords();

	  readnormal.set(norm1 -> GetTuple(p));  //set_vtk
	  if (DEBUG_MODE) {cout << endl << p << ":" << roundedpt << readnormal
<< endl;}
	  // write normal segment
	  g_out_volume.set_pen_color( normal_ldkidx_mode ? p :
LANDK_COLOR_START + p );
	  g_out_volume.set_voxels_around(roundedpt, readnormal, normal_length);
	  // (over)write landmark
	  g_out_volume.write_Point(roundedpt, normal_ldkidx_mode ? 0 : p);
	}
    }

  /* if non-fixed color mode, all str outline voxels and normal voxels
     must take values from given input volume */
  if (! fixed_color_mode)
    {
      in_volume.update_definitions();
      in_volume.allocate(); 
      in_volume.update_data();

      int z, x, y; 
      Voxel<DATATYPE> outvox(0);

      for (z = 0; z < dim_Z; z++)
	for (x = 0; x < dim_X; x++)
	  for (y = 0; y < dim_Y; y++)
	    {
	      outvox = g_out_volume.read_Voxel(x, y, z);
	      if (outvox.get_value() == OUTL_COLOR)
		{
		  outvox.set_value(in_volume.read_Value(outvox.toPoint()));
		  g_out_volume.write_Voxel(outvox);
		}
	    }
    }

  // creates/updates .img file on disk
  g_out_volume.write_data();

  return EXIT_SUCCESS;
}

> Jonathan Bailleul wrote:
> > Goodwin Lawlor wrote:
> >
> >>Something like vtkPolyDataToImageStencil?
> >
> >
> > I do not understand what it does!
> >
> > Anyway, I did a simple program writing (slice by slice) the contours of
> > a polydata into an RAW (air/analyse) image some time ago, if that can
> > help...
> >
> >
> >
> >
> >>"Leila Baghdadi" <baghdadi at sickkids.ca> wrote in message
> >>news:1090948771.980.6.camel at mouse18.phenogenomics.ca...
> >>
> >>>Dear Bill,
> >>>
> >>>Hope you are doing well,
> >>>
> >>>I was wondering since you are smart enough to have come up with
> >>>MarchingCubes,
> >>>
> >>>have you come up with something that does the opposite,
> >>>
> >>>i.e Input (mesh)  , output (binary image)
> >>>
> >>>
> >>>Yours Sincerely,
> >>>
> >>>
> >>>Leila
> >>>
> _______________________________________________
> This is the private VTK discussion list.
> Please keep messages on-topic. Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers

-- 
-----------------------------------
Jonathan BAILLEUL, Doctorant
GREYC Image - Université de Caen 
http://www.greyc.ismra.fr/~bailleul