[vtkusers] Python memory buildup

[Randy Heiland]

vtk-Pythoneers,

While I begin to look thru vtkFieldData (and others), I thought I'd throw this
question out to the (Python) group...

My situation is similar to what's recently been discussed on the list -- i.e.,
I have a time series of data and want to simply update fields.  My problem is
that, in Python, memory is being alloc'd and not released - eventually causing
memory to be exhausted.

The crux of the problem is found in the following snippet of my script.  For
completion, I've attached the full script below, as well as a C pgm to generate
a dummy struct grid file and two field files.

The memory buildup only occurs in a Python script.  Doing the equivalent
pipeline in C++ has no such problem - memory is apparently released properly.
 I haven't tried Tcl.

thanks,
--Randy

-----------------------------------
fieldReader = vtkDataObjectReader()
fieldReader.DebugOn()
distScalars = vtkScalars()

# loop and watch the memory buildup for this process...
for count in range(1,500,1):
  print '---------  ' + str(count) + '   ----------------'
  if (count % 2) == 1:
    fname = 'dist1.vtk'
  else:
    fname = 'dist2.vtk'
  print 'fname = ',fname

  fieldReader.SetFileName(fname)
  fieldReader.Update()

#  foo = fieldReader.GetOutput()                 # doesn't buildup memory
#  foo = fieldReader.GetOutput().GetFieldData()  # doesn't buildup memory
  foo = fieldReader.GetOutput().GetFieldData().GetArray(0)  # culprit

#  distScalars.SetData(fieldReader.GetOutput().GetFieldData().GetArray(0))
#  sgReader.GetOutput().GetPointData().SetScalars(distScalars)

#  renWin.Render()
-------------------------------------

>>>>>>>>>>>>> full Python script <<<<<<<<<<<<<<<<<<<
#!/usr/local/bin/python
#
#
import os
try:
  VTK_DATA = os.environ['VTK_DATA']
except KeyError:
  VTK_DATA = '../../../vtkdata/'

from libVTKCommonPython import *
from libVTKGraphicsPython import *
from libVTKLocalPython import *

#from colors import *

ren = vtkRenderer()
renWin = vtkRenderWindow()
renWin.AddRenderer(ren)
renWin.SetSize(300,300)
iren = vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

sgReader = vtkStructuredGridReader()
#sgReader.DebugOn()
sgReader.SetFileName('distSG.vtk')
sgReader.Update()
distSG = sgReader.GetOutput()

sgBox = vtkStructuredGridOutlineFilter()
sgBox.SetInput(distSG)

sgBoxMapper = vtkPolyDataMapper()
sgBoxMapper.ImmediateModeRenderingOn()
sgBoxMapper.SetInput(sgBox.GetOutput())

sgBoxActor = vtkActor()
sgBoxActor.SetMapper(sgBoxMapper)
sgBoxActor.GetProperty().SetColor(0.5,0.5,0.5)

_range = distSG.GetPointData().GetScalars().GetRange()
print 'dist scalar range= ',_range
ren.RemoveActor(sgBoxActor)

#-----------------------
contours1 = vtkContourFilter()
contours1.SetInput(distSG)

contours1Mapper = vtkDataSetMapper()
contours1Mapper.ImmediateModeRenderingOn()
contours1Mapper.SetInput(contours1.GetOutput())
#contours1Mapper.ScalarVisibilityOn()
contours1Mapper.ScalarVisibilityOff()
contours1Actor = vtkActor()
contours1Actor.SetMapper(contours1Mapper)

contours1.GenerateValues(6,_range[0], _range[1])
#contours1.GenerateValues(3,100, 200)
ren.AddActor(contours1Actor)

#=====================================
ren.AddActor(sgBoxActor)
ren.AddActor(contours1Actor)

renWin.Render()

fieldReader = vtkDataObjectReader()
fieldReader.DebugOn()
distScalars = vtkScalars()
#distScalars.DebugOn()

# loop and watch the memory buildup for this process...
for count in range(1,500,1):
  print '---------  ' + str(count) + '   ----------------'
  if (count % 2) == 1:
    fname = 'dist1.vtk'
  else:
    fname = 'dist2.vtk'
  print 'fname = ',fname

  fieldReader.SetFileName(fname)
  fieldReader.Update()

#  foo = fieldReader.GetOutput()                 # doesn't buildup memory
#  foo = fieldReader.GetOutput().GetFieldData()  # doesn't buildup memory
  foo = fieldReader.GetOutput().GetFieldData().GetArray(0)  # culprit

  del foo

#  distScalars.SetData(fieldReader.GetOutput().GetFieldData().GetArray(0))
#  sgReader.GetOutput().GetPointData().SetScalars(distScalars)

#  renWin.Render()

print '----------> done !!'
iren.Initialize()
iren.Start()


>>>>>>>>>>>>> C pgm to generate data files <<<<<<<<<<<<<<<<<<<

/*-----------
  genDist.c - Generate a VTK struct grid file (containing a single scalar
              "radial distance" field)
  and also generate a couple of field files that will be used to update the
  scalar field.


cc -O -o genDist genDist.c -lm

----------- */
#include <stdio.h>
#include <math.h>

#define SGFILE "distSG.vtk"
#define FLDFILE1 "dist1.vtk"
#define FLDFILE2 "dist2.vtk"
#define EPS 0.2

int appxEqual(float val1, float val2)
{
/*  if ( (val1 < (val2+EPS)) && (val1 > (val2-EPS)) ) */
  if ( (val1 <= (val2+EPS)) && (val1 >= (val2)) )
    return 1;
  else
    return 0;
}

void main(int argc, char **argv)
{
  FILE *fpout;
  int nx,ny,idx,idy;
  float x0,y0,xval,yval, pt[3], sval;
  float theta, perturb, radians;

  fpout = fopen(SGFILE,"w");

  nx = 500;
  ny = nx;
  x0 = -(float)nx / 2.0;
  y0 = -(float)ny / 2.0;

  /* ---- header ---- */
  fprintf(fpout,"# vtk DataFile Version 2.1\n");
  fprintf(fpout,"my dummy data\n");
  fprintf(fpout,"BINARY\n");
/*  fprintf(fpout,"ASCII\n"); */
  fprintf(fpout,"DATASET STRUCTURED_GRID\n");
  fprintf(fpout,"DIMENSIONS %d %d 1\n", nx,ny);
  fprintf(fpout,"POINTS %d float\n", nx*ny);

  /* ---- points ---- */
  yval = y0;
  pt[2] = 0.0;
  for (idy=0; idy<ny; idy++) {
    xval = x0;
    for (idx=0; idx<nx; idx++) {
      pt[0] = xval;
      pt[1] = yval;
      fwrite(pt,sizeof(float),3,fpout);
/*      fprintf(fpout,"%f %f %f\n",pt ); */
      xval += 1.0;
    }
    yval += 1.0;
  }

  fprintf(fpout,"POINT_DATA %d\n", nx*ny);
  fprintf(fpout,"SCALARS distance float\n");
  fprintf(fpout,"LOOKUP_TABLE default\n");

  yval = y0;
  for (idy=0; idy<ny; idy++) {
    xval = x0;
    for (idx=0; idx<nx; idx++) {
      sval = sqrt(xval*xval + yval*yval);    /* distance function */
      fwrite( &sval,sizeof(float),1,fpout);
      xval += 1.0;
    }
    yval += 1.0;
  }
  printf("---> %s\n",SGFILE);
  fclose(fpout);

  /*------------------------------------------------*/
  /* Create 1st field field */
  fpout = fopen(FLDFILE1,"w");
  /* ---- header ---- */
  fprintf(fpout,"# vtk DataFile Version 2.1\n");
  fprintf(fpout,"Distance field data: perturbation at theta = PI/4\n");
  fprintf(fpout,"BINARY\n");
  fprintf(fpout,"FIELD all 1\n");
  fprintf(fpout,"distance 1 %d float\n", nx*ny);

  radians = M_PI/4.0;
  printf("radians %f\n", radians);
  perturb = nx/20.0;
  yval = y0;
  for (idy=0; idy<ny; idy++) {
    xval = x0;
    for (idx=0; idx<nx; idx++) {
      sval = sqrt(xval*xval + yval*yval);
      theta = atan(yval/xval);
      if (appxEqual(theta, radians))
      {
        sval += perturb;
      }
      fwrite( &sval,sizeof(float),1,fpout);
      xval += 1.0;
    }
    yval += 1.0;
  }
  printf("---> %s\n",FLDFILE1);
  fclose(fpout);

  /*------------------------------------------------*/
  /* Create 2nd field field */
  fpout = fopen(FLDFILE2,"w");
  /* ---- header ---- */
  fprintf(fpout,"# vtk DataFile Version 2.1\n");
  fprintf(fpout,"Distance field data: perturbation at theta = -PI/4\n");
  fprintf(fpout,"BINARY\n");
  fprintf(fpout,"FIELD all 1\n");
  fprintf(fpout,"distance 1 %d float\n", nx*ny);

  radians = -M_PI/4.0;
  printf("radians %f\n", radians);
  yval = y0;
  for (idy=0; idy<ny; idy++) {
    xval = x0;
    for (idx=0; idx<nx; idx++) {
      sval = sqrt(xval*xval + yval*yval);
      theta = atan(yval/xval);
      if (appxEqual(theta, radians))
      {
        sval += perturb;
      }
      fwrite( &sval,sizeof(float),1,fpout);
      xval += 1.0;
    }
    yval += 1.0;
  }
  printf("---> %s\n",FLDFILE2);
  fclose(fpout);

}