[vtk-developers] potential speedup for the vtkpolydatatoimagestencil
Mark Roden
mmroden at gmail.com
Mon Feb 27 15:24:35 EST 2012
I'm working with both attempts at speedups at the moment. Our final
goal is a binary image of the data that is perfectly in register with
the image; I was getting too many problems confused at once on Friday.
Once we have the binary data, then we have a variety of display
mechanisms, including contours. If we see no contours, it's because
no binary data was made.
There are two methods suggested.
The first is to use the extruder, polydata to stencil, then the
stencil class to produce a binary image. As discussed, this approach
is slow because the extruder produces polygons, and when the polygons
are fed into the polydata to stencil class, a slower code path is
used. I just don't see how the stencil approach works without the
extruder. I don't see any vtk sample code that uses the stencil
without the extruder, and it just seems that the extruder is both
necessary and incredibly slow for this application. David, if you
could please shed some light on how to use the faster code path--
right now, it produces blank binary images.
The second is to iterate through the points in the image and, for each
point, check to see if it's valid using vtkPolygon::PointInPolygon.
This approach is producing a very garbled set of data, ie, is very
wrong. I'm attaching my code, in case I've made some obvious mistake,
but I believe I'm following the basic idea you laid out, Karthik. The
problem is the asserts that are fairly constantly firing, so I believe
that something's not quite right here.
Honestly, I don't care which approach gets used, I just need this
method to be faster.
vtkImageData* vtkPipelines::createBinaryImage(
vtkImageData * const binaryOrgan,
vtkPolyData * const data,
const EOrientation &eOrientation) {
if (data->GetNumberOfPoints() == 0) {
vtkImageData* returned = vtkImageData::New();
returned->DeepCopy(binaryOrgan);
return returned;
}
double spacing[3];
binaryOrgan->GetSpacing(spacing);
double origin[3];
binaryOrgan->GetOrigin(origin);
int extent[6];
binaryOrgan->GetExtent(extent);
double imageBounds[6];
binaryOrgan->GetBounds(imageBounds);
double ROIBounds[6];
data->GetBounds(ROIBounds);
double pt[3];
int displayPt[3];
double n[3];
int numPolys = data->GetNumberOfPolys();
int zmin, zmax;
vtkImageData *drawImage = vtkImageData::New();
drawImage->DeepCopy(binaryOrgan);
QVector<double> thePoint(3,0);
for (int i=0, index = 0; i<numPolys; ++i, ++index) {
int cellType = data->GetCellType(index);
while (cellType != VTK_POLYGON &&
cellType != VTK_QUAD &&
cellType != VTK_TRIANGLE)
++index;
vtkIdList *pointList = vtkIdList::New();
data->GetCellPoints(index, pointList);
double* points =
static_cast<double*>(data->GetPoints()->GetData()->GetVoidPointer(0));
int numPoints = pointList->GetNumberOfIds();
vtkPolygon::ComputeNormal(numPoints, points, n);
thePoint[2] = data->GetPoint(pointList->GetId(0))[2];
pt[2] = thePoint[2];
for (double y=ROIBounds[2]; y<=ROIBounds[3]; y += spacing[1]) {
for (double x=ROIBounds[0]; x<=ROIBounds[1]; x += spacing[0]) {
thePoint[0] = x;
thePoint[1] = y;
pt[0] = x;
pt[1] = y;
int inside = vtkPolygon::PointInPolygon(
pt,
numPoints,
points,
ROIBounds,
n);//constantly fires assertions, always zero.
}
}
pointList->Delete();
}
vtkImageData *returnImage = vtkImageData::New();
returnImage->DeepCopy(drawImage);
return returnImage;
}
On Fri, Feb 24, 2012 at 9:51 PM, Karthik Krishnan
<karthik.krishnan at kitware.com> wrote:
>
> On Sat, Feb 25, 2012 at 9:06 AM, Mark Roden <mmroden at gmail.com> wrote:
>>
>> That isn't working. I get an image with no contours if I do that.
>>
>> The polydata come from an rtstruct generated by Eclipse (ie, I know
>> it's 'good'). If I just use the polydata directly, with no extrusion,
>> then nothing shows up. If I use the extrusion, the data show up.
>
>
> Its likely that the polydata that needs to be shown is slightly occluded by
> the image actor. That's why extrusion works since it creates a wall out of
> the polydata with some height.
>
> Assuming that the imageActor and the polydata are on the same plane, this
> may be fixed by setting the mapper's
> ResolveCoincidentTopologyToPolygonOffset; although this is still susceptible
> to precision issues.
>
> Better yet, you may displace the polydata's points by a tiny amount towards
> the camera.
>
>>
>>
>> There are plenty of examples of this throughout the vtk examples;
>>
>> http://vtk.org/gitweb?p=VTK.git;a=blob;f=Examples/GUI/Python/ImageTracerWidget.py
>> is one.
>>
>> If I'm misunderstanding here how to display a contour on data-- ie, if
>> I can do this without the extrusion or the stencil-- I'd love to know.
>>
>> On Fri, Feb 24, 2012 at 7:24 PM, Karthik Krishnan
>> <karthik.krishnan at kitware.com> wrote:
>> >
>> >
>> > On Sat, Feb 25, 2012 at 3:13 AM, Mark Roden <mmroden at gmail.com> wrote:
>> >>
>> >> We need to display the contour itself, not voxels that are inside the
>> >> contour.
>> >
>> >
>> > Then why would you be generating the stencil ? Just cut the polydata
>> > with
>> > the slice plane and display it.
>> >
>> >
>> >>
>> >> It seems that Eclipse has trained everyone to expect this
>> >> appearance, and the rtstruct itself is based around this premise.
>> >>
>> >> Unless we can get a contour through that approach?
>> >>
>> >> On Fri, Feb 24, 2012 at 1:40 PM, Karthik Krishnan
>> >> <karthik.krishnan at kitware.com> wrote:
>> >> > But isn't vtkPolyDataToImageStencil, extrusion etc overkill for your
>> >> > specific problem. Why not compute the mask by looping over all voxels
>> >> > within
>> >> > the bounds of the contour on a slice and using
>> >> > vtkPolygon::PointInPolygon
>> >> > for every such voxel. It evaluates inside-outside by using rayfiring
>> >> > and
>> >> > that should run much faster; in the order of milliseconds
>> >> >
>> >> > On Sat, Feb 25, 2012 at 2:21 AM, Mark Roden <mmroden at gmail.com>
>> >> > wrote:
>> >> >>
>> >> >> Hi David,
>> >> >>
>> >> >> Sorry for being extremely tardy on continuing this thread, but the
>> >> >> issue took a back burner to more pressing problems that arose pretty
>> >> >> suddenly.
>> >> >>
>> >> >> Anyway, I think I see what's going on here, and I have some
>> >> >> hypotheses.
>> >> >>
>> >> >> We have a vtkPolyData read in by GDCM. This structure is read in
>> >> >> via
>> >> >> double precision (a change I've made recently that may or may not
>> >> >> have
>> >> >> any bearing).
>> >> >>
>> >> >> If we use the following series of vtk classes, the data can be
>> >> >> visualized on a CT Dicom image (whose z coordinates are in float
>> >> >> space, as best as I can tell):
>> >> >>
>> >> >> vtkLinearExtrusionFilter
>> >> >> vtkPolyDataToImageStencil
>> >> >> vtkImageStencil
>> >> >>
>> >> >> This series of operations, on one test data set, requires 25 seconds
>> >> >> on a fairly fast machine to produce contours for all the organs in
>> >> >> the
>> >> >> structure set. On our testers' machines, that number is roughly
>> >> >> triple, meaning that they are going nuts with waiting times.
>> >> >>
>> >> >> If we remove the extruder, then the time drops to 10 seconds, a
>> >> >> pretty
>> >> >> significant speedup (using vtk git source from 10 jan 2012).
>> >> >> Problem
>> >> >> is, now the data are not visible on the CT image at all. I suspect
>> >> >> that they are not visible because the coordinates do not match
>> >> >> _exactly_ in z, but I really have no way to back up that hypothesis.
>> >> >>
>> >> >> It seems that the extruder, in this case, is acting as an error
>> >> >> range
>> >> >> in z; ie, allowing the stencil a bit more leeway to finding the
>> >> >> proper
>> >> >> plane in Z. Trouble is, that code path is triggering the slow path
>> >> >> you talked about.
>> >> >>
>> >> >> Is there a way to set a kind of z tolerance to the stencil to
>> >> >> achieve
>> >> >> the same effect? Or am I off-base as to why the contours would not
>> >> >> appear on the data?
>> >> >>
>> >> >> Thanks,
>> >> >> Mark
>> >> >>
>> >> >>
>> >> >> On Wed, Jan 25, 2012 at 4:20 PM, David Gobbi <david.gobbi at gmail.com>
>> >> >> wrote:
>> >> >> > Hi Mark,
>> >> >> >
>> >> >> > In all my previous replies to your questions, I assumed that your
>> >> >> > data
>> >> >> > consisted of a surface (i.e. with polygonal faces). If your data
>> >> >> > is
>> >> >> > a
>> >> >> > series of polylines, then that changes things... particularly with
>> >> >> > respect to what version of VTK you are using.
>> >> >> >
>> >> >> > If your data consists of a series 2D polyline contours in 3D
>> >> >> > space,
>> >> >> > prior to VTK 5.8, vtkPolyDataToImageStencil could not use this
>> >> >> > kind
>> >> >> > of
>> >> >> > data, it required a 3D surface and it would cut that surface to
>> >> >> > generate contours. In VTK 5.8, the ability to directly process a
>> >> >> > series of 2D contours was added, but very inefficiently. In VTK
>> >> >> > git,
>> >> >> > the code for handling a series of 2D contours was modified to
>> >> >> > become
>> >> >> > much more efficient.
>> >> >> >
>> >> >> > For the vtkPolyDataToImageStencil in VTK 5.8, the PolyDataCutter
>> >> >> > function does this: if the input contains any only polylines, then
>> >> >> > the
>> >> >> > "if (cell->GetCellDimension() == 1)" branch is used. If the input
>> >> >> > contains any polygons, however, then the "if
>> >> >> > (cell->GetCellDimension()
>> >> >> > == 2)" branch is used. Only this latter branch actually does any
>> >> >> > cutting/contouring of the data.
>> >> >> >
>> >> >> > In VTK git, the code for polylines uses its own efficient
>> >> >> > subroutine
>> >> >> > called PolyDataSelector. So for polyline contours, this is the
>> >> >> > best
>> >> >> > version to use. Also, for polyline contours, the "loose end" code
>> >> >> > will
>> >> >> > not find any loose ends (and therefore not eat up any CPU) if the
>> >> >> > polylines that make up your contours are already closed. You can
>> >> >> > always run your contours through vtkCleanPolyData and vtkStripper
>> >> >> > to
>> >> >> > ensure that this is the case.
>> >> >> >
>> >> >> > - David
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> >
>> >> >> > On Wed, Jan 25, 2012 at 4:08 PM, Mark Roden <mmroden at gmail.com>
>> >> >> > wrote:
>> >> >> >> OK, so I've looked further into this problem, and it seems that I
>> >> >> >> was
>> >> >> >> optimizing the wrong thing.
>> >> >> >>
>> >> >> >> If I comment out everything other than the PolyDataCutter method
>> >> >> >> (specifically, the call to 'contour' in there), then contours
>> >> >> >> take
>> >> >> >> the
>> >> >> >> same amount of time to load as they do with everything turned on.
>> >> >> >> All
>> >> >> >> of the time is being spent in the Contour routine.
>> >> >> >>
>> >> >> >> However, rtstructs are already organized into contours. Each set
>> >> >> >> of
>> >> >> >> three numbers corresponds to a point in 3D space, and then each
>> >> >> >> succeeding group corresponds to the next point in the contour.
>> >> >> >> From
>> >> >> >> what I can determine, there is no need for this PolyDataCutter
>> >> >> >> routine.
>> >> >> >>
>> >> >> >> How can I tell the vtkPolyDataToImageStencil class that the data
>> >> >> >> are
>> >> >> >> already properly organized, that each point that follows is
>> >> >> >> connected
>> >> >> >> to the one previous, and in the case of a CLOSED contour, the
>> >> >> >> last
>> >> >> >> connects to the first?
>> >> >> >>
>> >> >> >> On Mon, Jan 23, 2012 at 11:27 AM, David Gobbi
>> >> >> >> <david.gobbi at gmail.com>
>> >> >> >> wrote:
>> >> >> >>> On Mon, Jan 23, 2012 at 12:09 PM, Mark Roden <mmroden at gmail.com>
>> >> >> >>> wrote:
>> >> >> >>>> Hi David,
>> >> >> >>>>
>> >> >> >>>> Thanks for the response, I'll take a look at the
>> >> >> >>>> clipclosedsurface
>> >> >> >>>> class. That may be what I want anyway (does it allow for any
>> >> >> >>>> one
>> >> >> >>>> plane to have a hollow interior?), but if you've already done
>> >> >> >>>> the
>> >> >> >>>> hashing there, then I'll take a look at porting it to this
>> >> >> >>>> class
>> >> >> >>>> as
>> >> >> >>>> well. Or should there be a shared routine that both classes
>> >> >> >>>> use,
>> >> >> >>>> if
>> >> >> >>>> the work done is very similar?
>> >> >> >>>
>> >> >> >>> A shared routine would be ideal, I had hoped to eventually write
>> >> >> >>> a
>> >> >> >>> vtkCutClosedSurface filter for that purpose.
>> >> >> >>>
>> >> >> >>> The vtkClipClosedSurface class does allow holes. It is
>> >> >> >>> described
>> >> >> >>> here: http://vtk.org/Wiki/VTK/Closed_Surface_Clipping
>> >> >> >>>
>> >> >> >>> - David
>> >> >> >>>
>> >> >> >>>
>> >> >> >>>
>> >> >> >>>> On Mon, Jan 23, 2012 at 11:04 AM, David Gobbi
>> >> >> >>>> <david.gobbi at gmail.com>
>> >> >> >>>> wrote:
>> >> >> >>>>> On Mon, Jan 23, 2012 at 11:19 AM, Mark Roden
>> >> >> >>>>> <mmroden at gmail.com>
>> >> >> >>>>> wrote:
>> >> >> >>>>>> Hi David,
>> >> >> >>>>>>
>> >> >> >>>>>> I'm looking to make the vtkpolydatatoimagestencil class
>> >> >> >>>>>> faster.
>> >> >> >>>>>> Right
>> >> >> >>>>>> now, on a core i7 machine and in using release-compiled code,
>> >> >> >>>>>> translating a body mask in an rtstruct to pixels using this
>> >> >> >>>>>> filter
>> >> >> >>>>>> is
>> >> >> >>>>>> prohibitively expensive in time (upwards of a minute for the
>> >> >> >>>>>> mask).
>> >> >> >>>>>>
>> >> >> >>>>>> There are three possible speedups to do here, in my mind. I
>> >> >> >>>>>> wanted
>> >> >> >>>>>> to
>> >> >> >>>>>> clear them by you, because I don't want to fork vtk to solve
>> >> >> >>>>>> this
>> >> >> >>>>>> problem, but the problem has become a serious problem for us.
>> >> >> >>>>>>
>> >> >> >>>>>> First, treat each plane independently of one another, and
>> >> >> >>>>>> then
>> >> >> >>>>>> have
>> >> >> >>>>>> each plane processed by different threads. This change is
>> >> >> >>>>>> fairly
>> >> >> >>>>>> straightforward theoretically, and something you mentioned
>> >> >> >>>>>> you
>> >> >> >>>>>> were
>> >> >> >>>>>> looking into a while back
>> >> >> >>>>>>
>> >> >> >>>>>>
>> >> >> >>>>>> (http://www.vtk.org/pipermail/vtkusers/2011-January/114538.html).
>> >> >> >>>>>> Is
>> >> >> >>>>>> this something you're still investigating?
>> >> >> >>>>>
>> >> >> >>>>> There were three things that I was considering:
>> >> >> >>>>> 1) multi-threading so that each CPU gets N slices to work on
>> >> >> >>>>> 2) increasing the efficiency of the polygon cutting code, I
>> >> >> >>>>> added
>> >> >> >>>>> some
>> >> >> >>>>> nice cutting code to vtkClipClosedSurface and planned to
>> >> >> >>>>> eventually
>> >> >> >>>>> also use it in vtkPolyDataToImageStencil
>> >> >> >>>>> 3) improving the efficiency or extent insertion for the
>> >> >> >>>>> stencils
>> >> >> >>>>>
>> >> >> >>>>> So far I've only done #3, which was the least important but
>> >> >> >>>>> was
>> >> >> >>>>> the
>> >> >> >>>>> easiest. Right now my own apps are bottlenecking on my
>> >> >> >>>>> segmentation
>> >> >> >>>>> algorithms, rather than on vtkPolyDataToImageStencil, so
>> >> >> >>>>> improving
>> >> >> >>>>> vtkPolyDataToImageStencil hasn't been a high priority.
>> >> >> >>>>>
>> >> >> >>>>>> For my work on other projects, I've found that the intel tbb
>> >> >> >>>>>> (http://threadingbuildingblocks.org/) makes multithreading
>> >> >> >>>>>> this
>> >> >> >>>>>> kind
>> >> >> >>>>>> of work very easy; the library is free and works with any c++
>> >> >> >>>>>> project
>> >> >> >>>>>> that uses the C++0x standard and can use lambda expressions.
>> >> >> >>>>>
>> >> >> >>>>> VTK will have to continue to support pre-C++0x compilers for a
>> >> >> >>>>> long
>> >> >> >>>>> time, several more years at least. So if threading is to be
>> >> >> >>>>> done,
>> >> >> >>>>> it
>> >> >> >>>>> should be done with VTK's threading classes.
>> >> >> >>>>>
>> >> >> >>>>>> Second, change the interior while/for loop collision
>> >> >> >>>>>> detection
>> >> >> >>>>>> to
>> >> >> >>>>>> be a
>> >> >> >>>>>> hash table. Right now, in the ThreadedExecute function,
>> >> >> >>>>>> there
>> >> >> >>>>>> is
>> >> >> >>>>>> this
>> >> >> >>>>>> code:
>> >> >> >>>>>>
>> >> >> >>>>>> for (vtkIdType i = 0; i < numberOfPoints; ++i)
>> >> >> >>>>>> {
>> >> >> >>>>>> ...
>> >> >> >>>>>> while( lines->GetNextCell(npts, pointIds) )
>> >> >> >>>>>> {
>> >> >> >>>>>> for (vtkIdType j = 0; j < npts; ++j)
>> >> >> >>>>>> if ( pointIds[j] == i )
>> >> >> >>>>>> {
>> >> >> >>>>>>
>> >> >> >>>>>> But what if collision detection was changed to uses a hash
>> >> >> >>>>>> table
>> >> >> >>>>>> where
>> >> >> >>>>>> the hashing function automatically detected point collisions
>> >> >> >>>>>> through a
>> >> >> >>>>>> single pass through the data?
>> >> >> >>>>>
>> >> >> >>>>> I use a hash vtkClipClosedSurface to accelerate the clipping
>> >> >> >>>>> (i.e.
>> >> >> >>>>> #2
>> >> >> >>>>> on my to-do list above). Take a look at the
>> >> >> >>>>> vtkClipClosedSurface
>> >> >> >>>>> code, specifically the vtkCCSEdgeLocator class.
>> >> >> >>>>>
>> >> >> >>>>>> Third, there does not appear to be an iterator over the
>> >> >> >>>>>> points
>> >> >> >>>>>> vector,
>> >> >> >>>>>> but a Get and Set function. These functions appear to be
>> >> >> >>>>>> pretty
>> >> >> >>>>>> slow,
>> >> >> >>>>>> and go through several thunking layers before data can
>> >> >> >>>>>> actually
>> >> >> >>>>>> be
>> >> >> >>>>>> set
>> >> >> >>>>>> or not. Is there an iterator class for points as there is
>> >> >> >>>>>> for
>> >> >> >>>>>> lines?
>> >> >> >>>>>> If not, how hard would it be to create such a class?
>> >> >> >>>>>
>> >> >> >>>>> The vtkPoints::GetData() method returns an array (either a
>> >> >> >>>>> vtkFloatArray or a vtkDoubleArray) that contains the points.
>> >> >> >>>>> Once
>> >> >> >>>>> you
>> >> >> >>>>> have this array, the GetTupleValue() method is a purely
>> >> >> >>>>> inlined
>> >> >> >>>>> method
>> >> >> >>>>> that can be used to efficiently get the points. If you know
>> >> >> >>>>> ahead
>> >> >> >>>>> of
>> >> >> >>>>> time whether the points are double or float, then this is
>> >> >> >>>>> probably
>> >> >> >>>>> the
>> >> >> >>>>> most efficient way of accessing them, apart from getting the
>> >> >> >>>>> raw
>> >> >> >>>>> "float *" or "double *".
>> >> >> >>>>>
>> >> >> >>>>> - David
>> >> >> _______________________________________________
>> >> >> Powered by www.kitware.com
>> >> >>
>> >> >> Visit other Kitware open-source projects at
>> >> >> http://www.kitware.com/opensource/opensource.html
>> >> >>
>> >> >> Follow this link to subscribe/unsubscribe:
>> >> >> http://www.vtk.org/mailman/listinfo/vtk-developers
>> >> >>
>> >> >
>> >> >
>> >> >
>> >> > --
>> >> > --
>> >> > karthik
>> >> >
>> >
>> >
>> >
>> >
>> > --
>> > --
>> > karthik
>> >
>
>
>
>
> --
> --
> karthik
>
More information about the vtk-developers
mailing list