[vtkusers] Points to 3D-Surface
George Kamucha
kamucha at hfm.e-technik.uni-kassel.de
Thu Jul 12 06:04:54 EDT 2001
Hi krs
I got involved in another project and didn't have enough time to go through the program you sent me
until today. I am really surprised that this program is doing what you said: "I wrote a little
utility 2 days back to map teh difference between an isosurface and another isosurface that
approximates it, to see the errors. I compute the difference by computing the point closest to each of
the vertices of the approximate surface (using PointLocator on the original surface),computed the
distances and produced a scalar array, I then color mapped this via a lookuptable onto the fitted
surface to show the errors (say blue being right on, and red being the largest error)..".
In the program (attached), the lookuptable is set to the original data, mc_orig, and so I am left
wondering how it can produce the colormapping of the difference data. The difference scalars, sc_impl,
don't play any role in the rendering effects?! Well, I may be wrong, but unless this LUT uses some
sort of magic to decide where to place the two colours......?
"K. R. Subramanian" wrote:
> George:
>
> Here is the program. It wont run, as it requires data files, etc.., but
> essentially it is working
> with 2 structured points datasets (mc_orig and mc_impl) - the
> differences are mapped onto
> mc_impl and displayed via a lookup table..
>
> -- krs
>
> --
> Dr. K.R.Subramanian Phone: (301) 402-0042
> National Institutes of Health FAX: (301) 402-4080
> National Library of Medicine Email: krs at cs.uncc.edu
> Bldg 38A, Rm. B1N30D WWW: www.cs.uncc.edu/~krs
> 8600 Rockville Pike, Bethesda, MD 20894
>
> ------------------------------------------------------------------------
> #include <stdio.h>
> #include <vtkSphereSource.h>
> #include <vtkRenderer.h>
> #include <vtkRenderWindow.h>
> #include <vtkRenderWindowInteractor.h>
> #include <vtkPolyData.h>
> #include <vtkPolyDataMapper.h>
> #include <vtkPolyDataNormals.h>
> #include <vtkActor.h>
> #include <vtkInteractorStyleTrackball.h>
> #include <vtkScalars.h>
> #include <vtkStructuredPoints.h>
> #include <vtkMarchingCubes.h>
> #include <vtkPointSet.h>
> #include <vtkPoints.h>
> #include <vtkPointLocator.h>
> #include <vtkNormals.h>
> #include <vtkSphereSource.h>
> #include <vtkLookupTable.h>
> #include <fstream>
>
> vtkStructuredPoints *ReadOriginal(char*),
> *ReadImplicit(char*);
>
> main(int argc, char **argv)
> {
> if (argc < 2)
> {
> cout << "usage: diff dset-file1 dset-file2" << endl;
> exit (-1);
> }
>
> float iso_val = 100.0;
>
> vtkRenderer *renderer = vtkRenderer::New();
> renderer->SetBackground(0.1,0.3,0.8);
>
> vtkRenderWindow *renWin = vtkRenderWindow::New();
> renWin->AddRenderer(renderer);
> renWin->SetSize(600,600);
>
> vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
> iren->SetRenderWindow(renWin);
>
> vtkInteractorStyleTrackball *trackball = vtkInteractorStyleTrackball::New();
> trackball->SetTrackballModeToTrackball();
>
> iren->SetInteractorStyle (trackball);
>
> vtkStructuredPoints *sp = ReadOriginal(argv[1]);
>
> vtkMarchingCubes *mc_orig = vtkMarchingCubes::New();
> mc_orig->SetInput (sp);
> mc_orig->SetValue (1, iso_val);
>
> vtkPolyDataNormals *poly_norm = vtkPolyDataNormals::New();
> poly_norm->SetInput (mc_orig->GetOutput());
>
> vtkLookupTable *lut = vtkLookupTable::New();
> lut->SetHueRange (0.667, 0.0);
> lut->SetSaturationRange (1.0, 1.0);
> lut->SetValueRange (1.0, 1.0);
> lut->SetNumberOfColors (1024);
> lut->Build();
>
> vtkPolyDataMapper *pm = vtkPolyDataMapper::New();
> pm->SetInput (poly_norm->GetOutput());
> pm->ScalarVisibilityOff();
> pm->SetLookupTable (lut);
>
> vtkStructuredPoints *sp2 = ReadImplicit(argv[2]);
>
> vtkMarchingCubes *mc_impl = vtkMarchingCubes::New();
> mc_impl->SetInput (sp2);
> mc_impl->SetValue (1, 0.);
>
> vtkPolyDataNormals *poly_norm2 = vtkPolyDataNormals::New();
> poly_norm2->SetInput (mc_impl->GetOutput());
>
> vtkPolyDataMapper *pm2 = vtkPolyDataMapper::New();
> pm2->SetInput (poly_norm2->GetOutput());
>
> pm->Update();
> pm2->Update();
>
> // mc_orig - original, mc_impl - implicit
>
> vtkPolyData *pd_orig, *pd_impl;
>
> pd_orig = mc_orig->GetOutput();
> pd_impl = mc_impl->GetOutput();
>
> vtkPointLocator *pl = vtkPointLocator::New();
> pl->SetDataSet (pd_orig);
>
> vtkPoints *pts_impl = pd_impl->GetPoints();
>
> float pt_orig[3], pt_impl[3], d, dx, dy, dz;
> int id;
>
> vtkScalars *sc_impl = vtkScalars::New(VTK_FLOAT, 1);
> sc_impl->SetNumberOfScalars (pts_impl->GetNumberOfPoints());
>
> float d_max = -1e10, d_min = 1e10, d_avg, d_sum = 0.;;
>
> for (int i = 0; i < pts_impl->GetNumberOfPoints(); i++)
> {
> pts_impl->GetPoint (i, pt_impl);
> id = pl->FindClosestPoint (pt_impl);
> pd_orig->GetPoint (id, pt_orig);
> dx = pt_orig[0] - pt_impl[0];
> dy = pt_orig[1] - pt_impl[1];
> dz = pt_orig[2] - pt_impl[2];
>
> d = sqrt((double) (dx*dx+dy*dy+dz*dz));
>
> if (d < d_min) d_min = d;
> if (d > d_max) d_max = d;
>
> d_sum += d;
>
> sc_impl->SetScalar (i, (float) d);
> }
>
> d_avg = d_sum/pts_impl->GetNumberOfPoints();
>
> cout << "Distance:(min, max, avg)" << d_min <<","<<d_max<<","<<d_avg<<endl;
>
> sc_impl->Modified();
>
> pd_impl->GetPointData()->SetScalars (sc_impl);
> pd_impl->Modified();
>
> vtkActor *mc_actor = vtkActor::New();
> mc_actor->SetMapper (pm2);
> mc_actor->GetProperty()->SetOpacity (1.0);
>
> renderer->AddActor (mc_actor2);
>
> renWin->Render();
> iren->Start();
>
> renderer->Delete();
> renWin->Delete();
> iren->Delete();
> }
> /////////////////////////////////////////////////////////////////////////
> vtkStructuredPoints *ReadImplicit (char *filename)
> {
> FILE *fp;
> float *slice_ptr, *slice_buf;
>
> vtkStructuredPoints *sp;
>
> if ((fp = fopen (filename, "r")) == NULL)
> {
> fprintf (stderr, "fopen(): couldnt open %s\n", filename);
> exit (-1);
> }
>
> float version;
> char binary[50], type[50], dummy[50], dummy2[50];
> int bytes_per_cell, res[3];
>
> fscanf (fp, "%f %s %s", &version, binary, type);
> if (version != 2.0)
> {
> fprintf (stderr, "Version %f not supported.\n", version);
> exit (-1);
> }
> fscanf ( fp, "%s %d %d %d %s %d\n", dummy, res, res+1, res+2,
> dummy2, &bytes_per_cell);
>
> int slice_size = res[0]*res[1]*4;
>
> if ((slice_buf = new float[slice_size]) == NULL)
> {
> fprintf (stderr, "slicebuf: out of memory\n");
> exit (-1);
> }
>
> vtkScalars *scalars = vtkScalars::New (VTK_FLOAT, 1);
> scalars->SetNumberOfScalars (res[0]*res[1]*res[2]);
>
> int j, k, n = 0,
> xy_res = res[0]*res[1];
>
> float val;
>
> for (k = 0; k < res[2]; k++)
> {
> slice_ptr = slice_buf;
> fread ((char *) slice_ptr, 1, slice_size, fp);
> for (j = 0; j < xy_res; j++, slice_ptr++)
> {
> val = *slice_ptr;
> scalars->SetScalar (n++, val);
> }
> }
>
> sp = vtkStructuredPoints::New();
>
> sp->SetDimensions (res);
> sp->GetPointData()->SetScalars (scalars);
>
> fprintf (stderr, "\t%d X %d X %d voxels.\n", res[0], res[1], res[2]);
>
> return sp;
>
> }
> //////////////////////////////////////////////////////////////////////////
> vtkStructuredPoints *ReadOriginal (char *filename)
> {
> FILE *fp;
> unsigned
> char *slice_ptr, *slice_buf;
>
> vtkStructuredPoints *sp;
>
> if ((fp = fopen (filename, "r")) == NULL)
> {
> fprintf (stderr, "fopen(): couldnt open %s\n", filename);
> exit (-1);
> }
>
> float version;
> char binary[50], type[50], dummy[50], dummy2[50];
> int bytes_per_cell, res[3];
>
> fscanf (fp, "%f %s %s", &version, binary, type);
> if (version != 2.0)
> {
> fprintf (stderr, "Version %f not supported.\n", version);
> exit (-1);
> }
> fscanf ( fp, "%s %d %d %d %s %d\n", dummy, res, res+1, res+2,
> dummy2, &bytes_per_cell);
>
> int slice_size = res[0]*res[1];
>
> if ((slice_buf = new unsigned char[slice_size]) == NULL)
> {
> fprintf (stderr, "slicebuf: out of memory\n");
> exit (-1);
> }
>
> vtkScalars *scalars = vtkScalars::New (VTK_FLOAT, 1);
> scalars->SetNumberOfScalars (res[0]*res[1]*res[2]);
>
> int j, k, n = 0,
> xy_res = res[0]*res[1];
>
> float val;
>
> for (k = 0; k < res[2]; k++)
> {
> slice_ptr = slice_buf;
> fread ((char *) slice_ptr, 1, slice_size, fp);
> for (j = 0; j < xy_res; j++, slice_ptr++)
> {
> val = *slice_ptr;
> scalars->SetScalar (n++, val);
> }
> }
>
> sp = vtkStructuredPoints::New();
>
> sp->SetDimensions (res);
> sp->GetPointData()->SetScalars (scalars);
>
> fprintf (stderr, "\t%d X %d X %d voxels.\n", res[0], res[1], res[2]);
>
> return sp;
>
> }
> //////////////////////////////////////////////////////////////////////////
More information about the vtkusers
mailing list