[vtkusers] contouring point data

[Randall Hand]

Attached to this email find my vtkMetaballs.cxx/h and a simple example
program using it.  This is a class I made about ayear ago to let you specify
points in space, and then sample it (with a metaballs algorithm) into a
rectilinear space.

You should be able to use this but replace my vtkXMLImageDataWriter with a
simple vtkContourFilter to get the result you want.  Or better yet, let it
write the VTI to disk and load that in ParaView to play with.

-- 
----------------------------------------
Randall Hand
Visualization Scientist
ERDC MSRC-ITL

On 7/26/07, Erik Anderson <eranders at sci.utah.edu> wrote:
>
> Hi all,
>      I have been trying to get iso-contours generated for
> point data for a while now with no luck.  Any help anyone
> can give would be much appreciated.  Below is a more
> complete description of the problem:
>
> I have a vtkUnstructuredGrid that is storing point data.
> In other words, every vertex in the dataset is a discrete
> point to be rendered.  Each point has several data arrays
> associated with it to be used for coloring, scaling, etc
> operations during glyph rendering.  I now want to generate
> isosurfaces for these data, but I cannot figure out how to
> do something like this from within VTK.
>
> One method that I can loosely translate into VTK language
> (but which I've had no success getting a pipeline working
> with) would be as follows:
>
> embed the point set in a structured grid
> for each vertex in the grid:
>     evaluate an implicit function (meta-ball like
> function)
>
> apply marching cubes or other isosurface extraction
> technique
> render surface
>
> The problem with the above approach is that I have been
> unable to construct a pipeline that would evaluate a
> function taking input from each point in the dataset.
>
> Any help would be greatly appreciated!
>
> Thank you,
> Erik Anderson
> _______________________________________________
> This is the private VTK discussion list.
> Please keep messages on-topic. Check the FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20070726/2d8c4d69/attachment.htm>
-------------- next part --------------
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkMetaballs.cxx,v $

  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.

=========================================================================*/
#include "vtkMetaballs.h"
#include "vtkObjectFactory.h"

vtkCxxRevisionMacro(vtkMetaballs, "$Revision: 1.2 $");
vtkStandardNewMacro(vtkMetaballs);

// Construct sphere with center at (0,0,0) and radius=0.5.
vtkMetaballs::vtkMetaballs()
{
  mb = NULL;
  mbSize = 0;
  mbCount = 0;
  mbExt = 1000;

  Resize(mbExt);
}

vtkMetaballs::~vtkMetaballs() 
{
    Reset();
}

void vtkMetaballs::Allocate(const vtkIdType sz, const vtkIdType ext)
{
    mbExt = ext;
    if (mb) {
        delete mb;
        mb = NULL;
    }
    mbCount = 0;
    mbSize = 0;
    Resize(sz);
}


void vtkMetaballs::Resize(const vtkIdType ext) 
{
    PointInfo *ptr;

    ptr = new PointInfo[mbSize + ext];
    if (mb) {
        // copy old into new
        memcpy(ptr, mb, sizeof(PointInfo) * mbSize);
        delete mb;
    }

    mb = ptr;
    mbSize += ext;
}

void vtkMetaballs::Reset(void)
{
    Allocate(1000, 1000);
}

void vtkMetaballs::AddMetaball(double x[3], double g)
{
    if (mbCount == mbSize-1) {
        Resize(mbExt);
    }
    this->mb[this->mbCount].p[0] = x[0];
    this->mb[this->mbCount].p[1] = x[1];
    this->mb[this->mbCount].p[2] = x[2];
    this->mb[this->mbCount].g = g;
    this->mbCount++;
}

// Evaluate metaball equation 
double vtkMetaballs::EvaluateFunction(double x[3])
{
  int i;
  PointInfo *point;
  double value = 0.0;
  double pValue, pDist;
  if (this->mbCount == 0)
      return 0;

  point = this->mb;
  for(i=0; i<this->mbCount; i++, point++) {
      pValue = point->g;
      pDist = (x[0] - point->p[0])*(x[0] - point->p[0]) +
              (x[1] - point->p[1])*(x[1] - point->p[1]) +
              (x[2] - point->p[2])*(x[2] - point->p[2]);
      value += (pValue / pDist);
  }
  //printf("EvalFun: %lf, %lf, %lf = %lf\n", x[0], x[1], x[2], value);
  return value;
}

// Evaluate sphere gradient.
void vtkMetaballs::EvaluateGradient(double x[3], double n[3])
{
  int i;
  PointInfo *point;
  double pValue;
  double pPoint[3];
  n[1] = n[2] = n[0] = 0;
  
  if (this->mbCount == 0) 
      return;

  point = this->mb;
  for(i=0; i<this->mbCount; i++, point++) {
      pValue = point->g;
      n[0] += (pValue / ((x[0] - point->p[0])*(x[0] - point->p[0])));
      n[1] += (pValue / ((x[1] - point->p[1])*(x[1] - point->p[1])));
      n[2] += (pValue / ((x[2] - point->p[2])*(x[2] - point->p[2])));
  }
}

void vtkMetaballs::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
}
-------------- next part --------------
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkMetaballs.h,v $

  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.

=========================================================================*/
// .NAME vtkMetaballs - implicit function for a sphere
// .SECTION Description
// vtkMetaballs computes the implicit function and/or gradient for a sphere.
// vtkMetaballs is a concrete implementation of vtkImplicitFunction.

#ifndef __vtkMetaballs_h
#define __vtkMetaballs_h

#include "vtkImplicitFunction.h"
#include <vtkPoints.h>
#include <vtkDoubleArray.h>

typedef struct {
    double p[3];
    double g;
} PointInfo;

class VTK_FILTERING_EXPORT vtkMetaballs : public vtkImplicitFunction
{
public:
  vtkTypeRevisionMacro(vtkMetaballs,vtkImplicitFunction);
  void PrintSelf(ostream& os, vtkIndent indent);

  // Description
  // Construct sphere with center at (0,0,0) and radius=0.5.
  static vtkMetaballs *New();

  // Description
  // Evaluate sphere equation ((x-x0)^2 + (y-y0)^2 + (z-z0)^2) - R^2.
  double EvaluateFunction(double x[3]);
  double EvaluateFunction(double x, double y, double z)
    {return this->vtkImplicitFunction::EvaluateFunction(x, y, z); } ;

  // Description
  // Evaluate sphere gradient.
  void EvaluateGradient(double x[3], double n[3]);

  void AddMetaball(double x[3], double g);
  void Reset();

  void Allocate(const vtkIdType sz, const vtkIdType ext=1000);

protected:
  vtkMetaballs();
  ~vtkMetaballs();

  void Resize(const vtkIdType ext);
  PointInfo *mb;
  unsigned int mbCount, mbSize;
  unsigned int mbExt;
private:
  vtkMetaballs(const vtkMetaballs&);  // Not implemented.
  void operator=(const vtkMetaballs&);  // Not implemented.
};

#endif


-------------- next part --------------
#include <vtkDataSet.h>
#include <vtkDataArray.h>
#include <vtkDataSetReader.h>
#include <vtkSampleFunction.h>
#include <vtkXMLImageDataWriter.h>
#include "vtkMetaballs.h"
#include <vtkPointData.h>
#include <vtkImageData.h>
#include <vtkObject.h>
#include <vtkCallbackCommand.h>
#include <vtkCommand.h>
#include "stat.h"

#define SCALE 1.0
static void ShowProgress( vtkObject *caller,
        unsigned long vtkNotUsed(eventId), void *, void *av) 
{
    double amount;
    amount = *((double*)av);
    printf("DEBUG: %s at %.1f%% Completion\n", caller->GetClassName(),
            (amount)*100.0);
}

int main (int argc, char *argv[]) {
    int i;
    double point[3];
    vtkCallbackCommand *vtkProgress = vtkCallbackCommand::New();
    vtkProgress->SetCallback(ShowProgress);

    vtkDataSetReader *reader = vtkDataSetReader::New();
    printf("Reading droplet file \"%s\"\n", argv[1]);
    reader->SetFileName(argv[1]);
            reader->ReadAllFieldsOn();
            reader->ReadAllScalarsOn();
            reader->ReadAllVectorsOn();
            reader->ReadAllTensorsOn();
            reader->ReadAllNormalsOn();
            reader->ReadAllColorScalarsOn();
            reader->ReadAllTCoordsOn();
    reader->Update();

    vtkDataSet *pdata = reader->GetOutput();
    PrintStatistics(pdata);
    vtkDataArray *radius = reader->GetOutput()->GetPointData()->GetArray("Radius");
    vtkMetaballs *mb = vtkMetaballs::New();
    printf("-> %i points...\n", pdata->GetNumberOfPoints());
    for(i=0; i<pdata->GetNumberOfPoints(); i++) {
        pdata->GetPoint(i, point);
        point[0] *= SCALE;
        point[1] *= SCALE;
        point[2] *= SCALE;
        mb->AddMetaball(point, radius->GetTuple1(i) * SCALE);
    }
    
    vtkSampleFunction *sampler = vtkSampleFunction::New();
    sampler->AddObserver(vtkCommand::ProgressEvent, vtkProgress);
    sampler->SetImplicitFunction(mb);
    sampler->SetOutputScalarTypeToDouble();
    sampler->SetSampleDimensions(50,50,50);
    {
        double *bounds;
        bounds = pdata->GetBounds();
        sampler->SetModelBounds(bounds[0]*SCALE, bounds[1]*SCALE,
                                bounds[2]*SCALE, bounds[3]*SCALE,
                                bounds[4]*SCALE, bounds[5]*SCALE);
    }
                            
    sampler->ComputeNormalsOff();
    sampler->CappingOff();
    PrintStatistics(sampler->GetOutput());

    printf("Sampling...\n");
    sampler->Update();
    vtkXMLImageDataWriter *writer = vtkXMLImageDataWriter::New();
    printf("Writing output \"%s\"...\n", argv[2]);
    writer->AddObserver(vtkCommand::ProgressEvent, vtkProgress);
    writer->SetFileName(argv[2]);
    writer->SetInput(sampler->GetOutput());
    writer->Write();
}