diff --git a/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx new file mode 100644 index 0000000000..21ae159511 --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx @@ -0,0 +1,452 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkImageStencilData.cxx + + 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 "vtkImageStencilRaster.h" + +#include "vtkImageStencilSource.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkDemandDrivenPipeline.h" +#include "vtkDataSetAttributes.h" +#include "vtkDataArray.h" +#include "vtkObjectFactory.h" +#include "vtkMath.h" + +#include +#include + +//---------------------------------------------------------------------------- +// tolerance for float-to-int conversion in stencil operations + +#define VTK_STENCIL_TOL 7.62939453125e-06 + +//---------------------------------------------------------------------------- +vtkImageStencilRaster::vtkImageStencilRaster(const int extent[2]) +{ + int rsize = extent[1] - extent[0] + 1; + + // the "raster" is a sequence of pointer-pairs, where the first pointer + // points to the first value in the raster line, and the second pointer + // points to one location past the last vale in the raster line. The + // difference is the number of x values stored in the raster line. + this->Raster = new double*[2*static_cast(rsize)]; + + // the extent is the range of y values (one line per y) + this->Extent[0] = extent[0]; + this->Extent[1] = extent[1]; + + // tolerance should be larger than expected roundoff errors + this->Tolerance = VTK_STENCIL_TOL; + + // no extent is used initially + this->UsedExtent[0] = 0; + this->UsedExtent[1] = -1; +} + +//---------------------------------------------------------------------------- +vtkImageStencilRaster::~vtkImageStencilRaster() +{ + if (this->UsedExtent[1] >= this->UsedExtent[0]) + { + size_t i = 2*static_cast(this->UsedExtent[0] - this->Extent[0]); + size_t imax = 2*static_cast(this->UsedExtent[1] - this->Extent[0]); + + do + { + if (this->Raster[i]) + { + delete [] this->Raster[i]; + } + i += 2; + } + while (i <= imax); + } + delete [] this->Raster; +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::PrepareForNewData(const int allocateExtent[2]) +{ + if (this->UsedExtent[1] >= this->UsedExtent[0]) + { + // reset and re-use the allocated raster lines + size_t i = 2*static_cast(this->UsedExtent[0]-this->Extent[0]); + size_t imax=2*static_cast(this->UsedExtent[1]-this->Extent[0]); + do + { + this->Raster[i+1] = this->Raster[i]; + i += 2; + } + while (i <= imax); + } + + if (allocateExtent && allocateExtent[1] >= allocateExtent[1]) + { + this->PrepareExtent(allocateExtent[0], allocateExtent[1]); + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::PrepareExtent(int ymin, int ymax) +{ + // this does not do any allocation, it just initializes any + // raster lines are not already part of the UsedExtent, and + // then expands the UsedExtent to include [ymin, ymax] + + if (this->UsedExtent[1] < this->UsedExtent[0]) + { + size_t i = 2*static_cast(ymin - this->Extent[0]); + size_t imax = 2*static_cast(ymax - this->Extent[0]); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[0] = ymin; + this->UsedExtent[1] = ymax; + + return; + } + + if (ymin < this->UsedExtent[0]) + { + size_t i = 2*static_cast(ymin - this->Extent[0]); + size_t imax = 2*static_cast(this->UsedExtent[0]-this->Extent[0]-1); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[0] = ymin; + } + + if (ymax > this->UsedExtent[1]) + { + size_t i = 2*static_cast(this->UsedExtent[1]+1 - this->Extent[0]); + size_t imax = 2*static_cast(ymax - this->Extent[0]); + + do + { + this->Raster[i] = 0; + } + while (++i <= imax); + + this->UsedExtent[1] = ymax; + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::InsertPoint(int y, double x) +{ + size_t pos = 2*static_cast(y - this->Extent[0]); + double* &rhead = this->Raster[pos]; + double* &rtail = this->Raster[pos+1]; + + // current size is the diff between the tail and the head + size_t n = rtail - rhead; + + // no allocation on this raster line yet + if (rhead == 0) + { + rhead = new double[2]; + rtail = rhead; + } + // grow whenever size reaches a power of two + else if (n > 1 && (n & (n-1)) == 0) + { + double *ptr = new double[2*n]; + for (size_t i = 0; i < n; i++) + { + ptr[i] = rhead[i]; + } + delete [] rhead; + rhead = ptr; + rtail = ptr + n; + } + + // insert the value + *rtail++ = x; +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::InsertLine( + const double pt1[2], const double pt2[2], + bool inflection1, bool inflection2) +{ + double x1 = pt1[0]; + double x2 = pt2[0]; + double y1 = pt1[1]; + double y2 = pt2[1]; + + // swap end points if necessary + if (y1 > y2) + { + x1 = pt2[0]; + x2 = pt1[0]; + y1 = pt2[1]; + y2 = pt1[1]; + bool tmp = inflection1; + inflection1 = inflection2; + inflection2 = tmp; + } + + // find min and max of x values + double xmin = x1; + double xmax = x2; + if (x1 > x2) + { + xmin = x2; + xmax = x1; + } + + // check for parallel to the x-axis + if (y1 == y2) + { + return; + } + + double ymin = y1; + double ymax = y2; + + // if an end is an inflection point, include a tolerance + ymin -= inflection1*this->Tolerance; + ymax += inflection2*this->Tolerance; + + // Integer y values for start and end of line + int iy1, iy2; + iy1 = this->Extent[0]; + iy2 = this->Extent[1]; + + // Check for out of bounds + if (ymax < iy1 || ymin >= iy2) + { + return; + } + + // Guard against extentY + if (ymin >= iy1) + { + iy1 = vtkMath::Floor(ymin) + 1; + } + if (ymax < iy2) + { + iy2 = vtkMath::Floor(ymax); + } + + // Expand allocated extent if necessary + if (iy1 < this->UsedExtent[0] || + iy2 > this->UsedExtent[1]) + { + this->PrepareExtent(iy1, iy2); + } + + // Precompute values for a Bresenham-like line algorithm + double grad = (x2 - x1)/(y2 - y1); + double delta = (iy1 - y1)*grad; + + // Go along y and place each x in the proper raster line + for (int y = iy1; y <= iy2; y++) + { + double x = x1 + delta; + // incrementing delta has less roundoff error than incrementing x, + // since delta will typically be smaller than x + delta += grad; + + // clamp x (because of tolerance, it might not be in range) + if (x < xmin) + { + x = xmin; + } + else if (x > xmax) + { + x = xmax; + } + + this->InsertPoint(y, x); + } +} + +//---------------------------------------------------------------------------- +void vtkImageStencilRaster::FillStencilData( + vtkImageStencilData *data, const int extent[6], int xj, int yj) +{ + if (xj != 0) + { + // slices are stacked in the x direction + int xmin = extent[2*xj]; + int xmax = extent[2*xj+1]; + int ymin = this->UsedExtent[0]; + int ymax = this->UsedExtent[1]; + int zmin = extent[0]; + int zmax = extent[1]; + + for (int idY = ymin; idY <= ymax; idY++) + { + size_t pos = 2*static_cast(idY - this->Extent[0]); + double *rline = this->Raster[pos]; + double *rlineEnd = this->Raster[pos+1]; + + if (rline == 0) + { + continue; + } + + vtkstd::sort(rline, rlineEnd); + + int xy[2]; + xy[2-xj] = idY; + + int lastr = VTK_INT_MIN; + + size_t l = rlineEnd - rline; + l = l - (l & 1); // force l to be an even number + for (size_t k = 0; k < l; k += 2) + { + double x1 = rline[k] - this->Tolerance; + double x2 = rline[k+1] + this->Tolerance; + + // make sure one of the ends is in bounds + if (x2 < xmin || x1 >= xmax) + { + continue; + } + + // clip the line segment with the bounds + int r1 = xmin; + int r2 = xmax; + + if (x1 >= xmin) + { + r1 = vtkMath::Floor(x1) + 1; + } + if (x2 < xmax) + { + r2 = vtkMath::Floor(x2); + } + + // ensure no overlap occurs with previous + if (r1 <= lastr) + { + r1 = lastr + 1; + } + lastr = r2; + + for (int idX = r1; idX <= r2; idX++) + { + xy[xj-1] = idX; + data->InsertNextExtent(zmin, zmax, xy[0], xy[1]); + } + } + } + } + else + { + // slices stacked in the y or z direction + int zj = 3 - yj; + int xmin = extent[0]; + int xmax = extent[1]; + int ymin = this->UsedExtent[0]; + int ymax = this->UsedExtent[1]; + int zmin = extent[2*zj]; + int zmax = extent[2*zj+1]; + + // convert each raster line into extents for the stencil + for (int idY = ymin; idY <= ymax; idY++) + { + size_t pos = 2*static_cast(idY - this->Extent[0]); + double *rline = this->Raster[pos]; + double *rlineEnd = this->Raster[pos+1]; + + if (rline == 0) + { + continue; + } + + vtkstd::sort(rline, rlineEnd); + + int lastr = VTK_INT_MIN; + + // go through each raster line and fill the stencil + size_t l = rlineEnd - rline; + l = l - (l & 1); // force l to be an even number + for (size_t k = 0; k < l; k += 2) + { + int yz[2]; + + yz[yj-1] = idY; + yz[2-yj] = zmin; + + double x1 = rline[k] - this->Tolerance; + double x2 = rline[k+1] + this->Tolerance; + + if (x2 < xmin || x1 >= xmax) + { + continue; + } + + int r1 = xmin; + int r2 = xmax; + + if (x1 >= xmin) + { + r1 = vtkMath::Floor(x1) + 1; + } + if (x2 < xmax) + { + r2 = vtkMath::Floor(x2); + } + + // ensure no overlap occurs between extents + if (r1 <= lastr) + { + r1 = lastr + 1; + } + lastr = r2; + + if (r2 >= r1) + { + data->InsertNextExtent(r1, r2, yz[0], yz[1]); + } + } + } + + // copy the result to all other slices + if (zmin < zmax) + { + for (int idY = ymin; idY <= ymax; idY++) + { + int r1, r2; + int yz[2]; + + yz[yj-1] = idY; + yz[2-yj] = zmin; + + int iter = 0; + while (data->GetNextExtent(r1, r2, xmin, xmax, yz[0], yz[1], iter)) + { + for (int idZ = zmin + 1; idZ <= zmax; idZ++) + { + yz[2-yj] = idZ; + data->InsertNextExtent(r1, r2, yz[0], yz[1]); + } + yz[2-yj] = zmin; + } + } + } + } +} diff --git a/Modules/MitkExt/Algorithms/vtkImageStencilRaster.h b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.h new file mode 100644 index 0000000000..1815d054fa --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.h @@ -0,0 +1,98 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkImageStencilData.h + + 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 vtkImageStencilData - efficient description of an image stencil +// .SECTION Description +// vtkImageStencilData describes an image stencil in a manner which is +// efficient both in terms of speed and storage space. The stencil extents +// are stored for each x-row across the image (multiple extents per row if +// necessary) and can be retrieved via the GetNextExtent() method. +// .SECTION see also +// vtkImageStencilSource vtkImageStencil + +#ifndef __vtkImageStencilRaster_h +#define __vtkImageStencilRaster_h + +#include "vtkImageStencilData.h" +#include "vtkDataObject.h" +#include "MitkExtExports.h" + + +//BTX +// Description: +// This is a helper class for stencil creation. It is a raster with +// infinite resolution in the X direction (approximately, since it uses +// double precision). Lines that represent polygon edges can be drawn +// into this raster, and then filled given a tolerance. +class MitkExt_EXPORT vtkImageStencilRaster +{ +public: + // Description: + // Create a raster with the specified whole y extent. + vtkImageStencilRaster(const int wholeExtent[2]); + + // Description: + // Destructor. + ~vtkImageStencilRaster(); + + // Description: + // Reset the raster to its original state, but keep the same whole + // extent. Pre-allocate the specified 1D allocateExtent, which must be + // within the whole extent. + void PrepareForNewData(const int allocateExtent[2] = 0); + + // Description: + // Insert a line into the raster, given the two end points. + // The "inflect1" and "inflect2" should be set if you want + // to add a small vertical tolerance to either endpoints. + void InsertLine(const double p1[2], const double p2[2], + bool inflect1, bool inflect2); + + // Description: + // Fill the specified extent of a vtkImageStencilData with the raster, + // after permuting the raster according to xj and yj. + void FillStencilData(vtkImageStencilData *data, const int extent[6], + int xj = 0, int yj = 1); + + // Description: + // The tolerance for float-to-int conversions. + void SetTolerance(double tol) { this->Tolerance = tol; } + double GetTolerance() { return this->Tolerance; } + +protected: + // Description: + // Ensure that the raster is initialized for the specified range + // of y values, which must be within the Extent. + void PrepareExtent(int ymin, int ymax); + + // Description: + // Insert an x point into the raster. If the y value is larger + // than the y extent, the extent will grow automatically. + void InsertPoint(int y, double x); + + int Extent[2]; + int UsedExtent[2]; + double **Raster; + double Tolerance; + +private: + vtkImageStencilRaster(const vtkImageStencilRaster&); // Not implemented. + void operator=(const vtkImageStencilRaster&); // Not implemented. +}; +//ETX + +#endif + + + diff --git a/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx new file mode 100644 index 0000000000..bb172686ab --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx @@ -0,0 +1,609 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkLassoStencilSource.cxx + + 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 "vtkLassoStencilSource.h" + +#include "vtkMath.h" +#include "vtkPoints.h" +#include "vtkCardinalSpline.h" +#include "vtkDataArray.h" +#include "vtkImageData.h" +#include "vtkImageStencilData.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkObjectFactory.h" +#include "vtkStreamingDemandDrivenPipeline.h" +#include "vtkSmartPointer.h" + +#include "vtkImageStencilRaster.h" + +#include +#include + +vtkStandardNewMacro(vtkLassoStencilSource); +vtkCxxSetObjectMacro(vtkLassoStencilSource, Points, vtkPoints); + +//---------------------------------------------------------------------------- +class vtkLSSPointMap : public vtkstd::map > +{ +}; + +//---------------------------------------------------------------------------- +vtkLassoStencilSource::vtkLassoStencilSource() +{ + this->SetNumberOfInputPorts(0); + + this->Shape = vtkLassoStencilSource::POLYGON; + this->SliceOrientation = 2; + this->Points = NULL; + this->SplineX = vtkCardinalSpline::New(); + this->SplineY = vtkCardinalSpline::New(); + + this->PointMap = new vtkLSSPointMap(); +} + +//---------------------------------------------------------------------------- +vtkLassoStencilSource::~vtkLassoStencilSource() +{ + this->SetPoints(NULL); + if (this->SplineX) + { + this->SplineX->Delete(); + this->SplineX = NULL; + } + if (this->SplineY) + { + this->SplineY->Delete(); + this->SplineY = NULL; + } + if (this->PointMap) + { + delete this->PointMap; + this->PointMap = NULL; + } +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::PrintSelf(ostream& os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os,indent); + + os << indent << "Shape: " << this->GetShapeAsString() << "\n"; + os << indent << "Points: " << this->Points << "\n"; + os << indent << "SliceOrientation: " << this->GetSliceOrientation() << "\n"; + os << indent << "SlicePoints: " << this->PointMap->size() << "\n"; +} + +//---------------------------------------------------------------------------- +const char *vtkLassoStencilSource::GetShapeAsString() +{ + switch (this->Shape) + { + case vtkLassoStencilSource::POLYGON: + return "Polygon"; + case vtkLassoStencilSource::SPLINE: + return "Spline"; + } + return ""; +} + +//---------------------------------------------------------------------------- +unsigned long int vtkLassoStencilSource::GetMTime() +{ + unsigned long mTime = this->vtkImageStencilSource::GetMTime(); + + if ( this->Points != NULL ) + { + unsigned long t = this->Points->GetMTime(); + if (t > mTime) + { + mTime = t; + } + } + + if ( !this->PointMap->empty() ) + { + vtkLSSPointMap::iterator iter = this->PointMap->begin(); + while ( iter != this->PointMap->end() ) + { + unsigned long t = iter->second->GetMTime(); + if (t > mTime) + { + mTime = t; + } + iter++; + } + } + + return mTime; +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::SetSlicePoints(int i, vtkPoints *points) +{ + vtkLSSPointMap::iterator iter = this->PointMap->find(i); + if (iter != this->PointMap->end()) + { + if (iter->second == points) + { + return; + } + else if (points == 0) + { + this->PointMap->erase(iter); + } + else + { + iter->second = points; + } + } + else + { + if (points == NULL) + { + return; + } + else + { + this->PointMap->insert(iter, vtkLSSPointMap::value_type(i, points)); + } + } + + this->Modified(); +} + +//---------------------------------------------------------------------------- +void vtkLassoStencilSource::RemoveAllSlicePoints() +{ + this->PointMap->clear(); +} + +//---------------------------------------------------------------------------- +vtkPoints *vtkLassoStencilSource::GetSlicePoints(int i) +{ + vtkLSSPointMap::iterator iter = this->PointMap->find(i); + if (iter != this->PointMap->end()) + { + return iter->second; + } + return NULL; +} + +//---------------------------------------------------------------------------- +// tolerance for stencil operations + +#define VTK_STENCIL_TOL 7.62939453125e-06 + +//---------------------------------------------------------------------------- +// Compute a reduced extent based on the bounds of the shape. +static void vtkLassoStencilSourceSubExtent( + vtkPoints *points, + const double origin[3], const double spacing[3], + const int extent[6], int subextent[6]) +{ + double bounds[6]; + points->GetBounds(bounds); + + for (int i = 0; i < 3; i++) + { + double emin = (bounds[2*i] - origin[i])/spacing[i] - VTK_STENCIL_TOL; + double emax = (bounds[2*i+1] - origin[i])/spacing[i] + VTK_STENCIL_TOL; + + subextent[2*i] = extent[2*i]; + subextent[2*i+1] = extent[2*i+1]; + + if (extent[2*i] < emin) + { + subextent[2*i] = VTK_INT_MAX; + if (extent[2*i+1] >= emin) + { + subextent[2*i] = vtkMath::Floor(emin) + 1; + } + } + + if (extent[2*i+1] > emax) + { + subextent[2*i+1] = VTK_INT_MIN; + if (extent[2*i] <= emax) + { + subextent[2*i+1] = vtkMath::Floor(emax); + } + } + + } +} + +//---------------------------------------------------------------------------- +// Rasterize a polygon into the stencil +static int vtkLassoStencilSourcePolygon( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + const int extent[6], const double origin[3], const double spacing[3], + int xj, int yj) +{ + // get the bounds of the polygon + int subextent[6]; + vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent); + + // allocate the raster + raster->PrepareForNewData(&subextent[2*yj]); + + // rasterize each line + vtkIdType n = points->GetNumberOfPoints(); + double p[3]; + double p0[2], p1[2], p2[2], p3[2]; + + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + points->GetPoint(0, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + double dx = p1[0] - p0[0]; + double dy = p1[1] - p0[1]; + if (dx*dx + dy*dy <= VTK_STENCIL_TOL*VTK_STENCIL_TOL) + { + n -= 1; + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + } + + points->GetPoint(1, p); + p2[0] = (p[xj] - origin[xj])/spacing[xj]; + p2[1] = (p[yj] - origin[yj])/spacing[yj]; + + // inflection means the line changes vertical direction + bool inflection1, inflection2; + inflection1 = ( (p1[1] - p0[1])*(p2[1] - p1[1]) <= 0 ); + + for (vtkIdType i = 0; i < n; i++) + { + points->GetPoint((i+2)%n, p); + p3[0] = (p[xj] - origin[xj])/spacing[xj]; + p3[1] = (p[yj] - origin[yj])/spacing[yj]; + + inflection2 = ( (p2[1] - p1[1])*(p3[1] - p2[1]) <= 0 ); + + raster->InsertLine(p1, p2, inflection1, inflection2); + + p0[0] = p1[0]; p0[1] = p1[1]; + p1[0] = p2[0]; p1[1] = p2[1]; + p2[0] = p3[0]; p2[1] = p3[1]; + inflection1 = inflection2; + } + + raster->FillStencilData(data, extent, xj, yj); + + return 1; +} + + +//---------------------------------------------------------------------------- +// Generate the splines for the given set of points. The splines +// will be closed if the final point is equal to the first point. +// The parametric value for the resulting spline will be valid over +// the range [0, tmax] where the tmax value is returned by reference. +static void vtkLassoStencilSourceCreateSpline(vtkPoints *points, + const double origin[3], const double spacing[3], + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline, + double &tmax, double &dmax) +{ + // initialize the spline + xspline->RemoveAllPoints(); + yspline->RemoveAllPoints(); + xspline->ClosedOff(); + yspline->ClosedOff(); + + // get the number of points and the first/last point + vtkIdType n = points->GetNumberOfPoints(); + double p[3]; + double p0[2], p1[2]; + + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + points->GetPoint(0, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + // factor between real distance and parametric distance + double f = 1.0; + // the length of the implicit segment for closed loops + double lastd = 0; + + // aspect ratio + double xf = 1.0; + double yf = 1.0; + if (spacing[xj] > spacing[yj]) + { + xf = spacing[xj]/spacing[yj]; + } + else + { + yf = spacing[yj]/spacing[xj]; + } + + // if first and last point are same, spline is closed + double dx = (p1[0] - p0[0])*xf; + double dy = (p1[1] - p0[1])*yf; + double d2 = dx*dx + dy*dy; + while (d2 <= VTK_STENCIL_TOL*VTK_STENCIL_TOL && n > 1) + { + n -= 1; + points->GetPoint(n-1, p); + p0[0] = (p[xj] - origin[xj])/spacing[xj]; + p0[1] = (p[yj] - origin[yj])/spacing[yj]; + + xspline->ClosedOn(); + yspline->ClosedOn(); + + // vtkSpline considers the parametric length of the implicit + // segment of closed loops to be unity, so set "f" so that + // multiplying the real length of that segment by "f" gives unity. + dx = (p1[0] - p0[0])*xf; + dy = (p1[1] - p0[1])*yf; + d2 = dx*dx + dy*dy; + lastd = sqrt(d2); + if (lastd > 0) + { + f = 1.0/lastd; + } + } + + // Add all the points to the spline. + double d = 0.0; + for (vtkIdType i = 0; i < n; i++) + { + p0[0] = p1[0]; p0[1] = p1[1]; + + points->GetPoint(i, p); + p1[0] = (p[xj] - origin[xj])/spacing[xj]; + p1[1] = (p[yj] - origin[yj])/spacing[yj]; + + dx = (p1[0] - p0[0])*xf; + dy = (p1[1] - p0[1])*yf; + + d += sqrt(dx*dx + dy*dy); + + double t = f*d; + + xspline->AddPoint(t, p1[0]); + yspline->AddPoint(t, p1[1]); + } + + // Do the spline precomputations + xspline->Compute(); + yspline->Compute(); + + // The spline is valid over t = [0, tmax] + d += lastd; + tmax = f*d; + dmax = d; +} + +//---------------------------------------------------------------------------- +// Rasterize a spline contour into the stencil +static int vtkLassoStencilSourceSpline( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + const int extent[6], const double origin[3], const double spacing[3], + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline) +{ + // create the splines + double tmax, dmax; + vtkLassoStencilSourceCreateSpline( + points, origin, spacing, xj, yj, xspline, yspline, tmax, dmax); + + if (dmax <= VTK_STENCIL_TOL) + { + return 1; + } + + // get the bounds of the polygon as a first guess of the spline bounds + int subextent[6]; + vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent); + + // allocate the raster + raster->PrepareForNewData(&subextent[2*yj]); + + // go around the spline + vtkIdType n = vtkMath::Floor(dmax)+1; + double delta = tmax/n; + + double p0[2], p1[2], p2[2], p3[2]; + + double t = tmax; + if (xspline->GetClosed()) + { + t = (n-1)*tmax/n; + } + else + { + n = n + 1; + } + + p0[0] = xspline->Evaluate(t); + p0[1] = yspline->Evaluate(t); + + t = 0; + p1[0] = xspline->Evaluate(t); + p1[1] = yspline->Evaluate(t); + + t = delta; + p2[0] = xspline->Evaluate(t); + p2[1] = yspline->Evaluate(t); + + // inflection means the line changes vertical direction + bool inflection1, inflection2; + inflection1 = ( (p1[1] - p0[1])*(p2[1] - p1[1]) <= 0 ); + + for (vtkIdType i = 0; i < n; i++) + { + t += delta; + if (i == n-2) + { + t = 0; + } + + p3[0] = xspline->Evaluate(t); + p3[1] = yspline->Evaluate(t); + + inflection2 = ( (p2[1] - p1[1])*(p3[1] - p2[1]) <= 0 ); + + raster->InsertLine(p1, p2, inflection1, inflection2); + + p0[0] = p1[0]; p0[1] = p1[1]; + p1[0] = p2[0]; p1[1] = p2[1]; + p2[0] = p3[0]; p2[1] = p3[1]; + inflection1 = inflection2; + } + + raster->FillStencilData(data, extent, xj, yj); + + return 1; +} + +//---------------------------------------------------------------------------- +static int vtkLassoStencilSourceExecute( + vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster, + int extent[6], double origin[3], double spacing[3], int shape, + int xj, int yj, vtkSpline *xspline, vtkSpline *yspline) +{ + int result = 1; + + if (points == 0 || points->GetNumberOfPoints() < 3) + { + return 1; + } + + switch (shape) + { + case vtkLassoStencilSource::POLYGON: + result = vtkLassoStencilSourcePolygon( + points, data, raster, extent, origin, spacing, xj, yj); + break; + case vtkLassoStencilSource::SPLINE: + result = vtkLassoStencilSourceSpline( + points, data, raster, extent, origin, spacing, xj, yj, + xspline, yspline); + break; + } + + return result; +} + + +//---------------------------------------------------------------------------- +int vtkLassoStencilSource::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + int extent[6]; + double origin[3]; + double spacing[3]; + int result = 1; + + this->Superclass::RequestData(request, inputVector, outputVector); + + vtkInformation *outInfo = outputVector->GetInformationObject(0); + vtkImageStencilData *data = vtkImageStencilData::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent); + outInfo->Get(vtkDataObject::ORIGIN(), origin); + outInfo->Get(vtkDataObject::SPACING(), spacing); + + int slabExtent[6]; + slabExtent[0] = extent[0]; slabExtent[1] = extent[1]; + slabExtent[2] = extent[2]; slabExtent[3] = extent[3]; + slabExtent[4] = extent[4]; slabExtent[5] = extent[5]; + + int xj = 0; + int yj = 1; + int zj = 2; + + if (this->SliceOrientation == 0) + { + xj = 1; + yj = 2; + zj = 0; + } + else if (this->SliceOrientation == 1) + { + xj = 0; + yj = 2; + zj = 1; + } + + vtkImageStencilRaster raster(&extent[2*yj]); + raster.SetTolerance(VTK_STENCIL_TOL); + + int zmin = extent[2*zj]; + int zmax = extent[2*zj+1]; + + vtkLSSPointMap::iterator iter = this->PointMap->lower_bound(zmin); + vtkLSSPointMap::iterator maxiter = this->PointMap->upper_bound(zmax); + + while (iter != maxiter && result != 0) + { + this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1)); + + int i = iter->first; + vtkPoints *points = iter->second; + + // fill in the slices with no SlicePoints + if (this->Points && i > slabExtent[2*zj]) + { + slabExtent[2*zj+1] = i-1; + + result = vtkLassoStencilSourceExecute( + this->Points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + } + + // do the slice with its SlicePoints + if (result) + { + slabExtent[2*zj] = i; + slabExtent[2*zj+1] = i; + + result = vtkLassoStencilSourceExecute( + points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + + slabExtent[2*zj] = slabExtent[2*zj+1] + 1; + } + + ++iter; + } + + this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1)); + + // fill in the rest + if (result && slabExtent[2*zj] <= zmax) + { + slabExtent[2*zj+1] = zmax; + + result = vtkLassoStencilSourceExecute( + this->Points, data, &raster, slabExtent, origin, spacing, + this->Shape, xj, yj, this->SplineX, this->SplineY); + + this->SetProgress(1.0); + } + + return result; +} diff --git a/Modules/MitkExt/Algorithms/vtkLassoStencilSource.h b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.h new file mode 100644 index 0000000000..86ada7a3b7 --- /dev/null +++ b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.h @@ -0,0 +1,107 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkLassoStencilSource.h + + 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 vtkLassoStencilSource - Create a stencil from a contour +// .SECTION Description +// vtkLassoStencilSource will create an image stencil from a +// set of points that define a contour. Its output can be +// used with vtkImageStecil or other vtk classes that apply +// a stencil to an image. +// .SECTION See Also +// vtkROIStencilSource vtkPolyDataToImageStencil +// .SECTION Thanks +// Thanks to David Gobbi for contributing this class to VTK. + +#ifndef __vtkLassoStencilSource_h +#define __vtkLassoStencilSource_h + + +#include "vtkImageStencilSource.h" +#include "MitkExtExports.h" + +class vtkPoints; +class vtkSpline; +class vtkLSSPointMap; + +class MitkExt_EXPORT vtkLassoStencilSource : public vtkImageStencilSource +{ +public: + static vtkLassoStencilSource *New(); + vtkTypeMacro(vtkLassoStencilSource, vtkImageStencilSource); + void PrintSelf(ostream& os, vtkIndent indent); + +//BTX + enum { + POLYGON = 0, + SPLINE = 1, + }; +//ETX + + // Description: + // The shape to use, default is "Polygon". The spline is a + // cardinal spline. Bezier splines are not yet supported. + vtkGetMacro(Shape, int); + vtkSetClampMacro(Shape, int, POLYGON, SPLINE); + void SetShapeToPolygon() { this->SetShape(POLYGON); }; + void SetShapeToSpline() { this->SetShape(SPLINE); }; + virtual const char *GetShapeAsString(); + + // Description: + // The points that make up the lassoo. The loop does not + // have to be closed, the last point will automatically be + // connected to the first point by a straight line segment. + virtual void SetPoints(vtkPoints *points); + vtkGetObjectMacro(Points, vtkPoints); + + // Description: + // The slice orientation. The default is 2, which is XY. + // Other values are 0, which is YZ, and 1, which is XZ. + vtkGetMacro(SliceOrientation, int); + vtkSetClampMacro(SliceOrientation, int, 0, 2); + + // Description: + // The points for a particular slice. This will override the + // points that were set by calling SetPoints() for the slice. + // To clear the setting, call SetSlicePoints(slice, NULL). + virtual void SetSlicePoints(int i, vtkPoints *points); + virtual vtkPoints *GetSlicePoints(int i); + + // Description: + // Remove points from all slices. + virtual void RemoveAllSlicePoints(); + + // Description: + // Overload GetMTime() to include the timestamp on the points. + unsigned long GetMTime(); + +protected: + vtkLassoStencilSource(); + ~vtkLassoStencilSource(); + + virtual int RequestData(vtkInformation *, vtkInformationVector **, + vtkInformationVector *); + + int Shape; + int SliceOrientation; + vtkPoints *Points; + vtkSpline *SplineX; + vtkSpline *SplineY; + vtkLSSPointMap *PointMap; + +private: + vtkLassoStencilSource(const vtkLassoStencilSource&); // Not implemented. + void operator=(const vtkLassoStencilSource&); // Not implemented. +}; + +#endif diff --git a/Modules/MitkExt/files.cmake b/Modules/MitkExt/files.cmake index 8316b38ec7..6458616ce0 100644 --- a/Modules/MitkExt/files.cmake +++ b/Modules/MitkExt/files.cmake @@ -1,196 +1,198 @@ SET(CPP_FILES Algorithms/mitkMaskAndCutRoiImageFilter.cpp Algorithms/mitkBoundingObjectToSegmentationFilter.cpp Algorithms/vtkPointSetSlicer.cxx Algorithms/mitkCoreExtObjectFactory.cpp Algorithms/mitkImageStatisticsCalculator.cpp Algorithms/mitkAngleCorrectByPointFilter.cpp Algorithms/mitkAutoCropImageFilter.cpp Algorithms/mitkBoundingObjectCutter.cpp Algorithms/mitkCalculateSegmentationVolume.cpp Algorithms/mitkContourSetToPointSetFilter.cpp Algorithms/mitkContourUtils.cpp Algorithms/mitkCorrectorAlgorithm.cpp Algorithms/mitkCylindricToCartesianFilter.cpp Algorithms/mitkDiffImageApplier.cpp Algorithms/mitkDopplerToStrainRateFilter.cpp Algorithms/mitkExtractDirectedPlaneImageFilter.cpp Algorithms/mitkExtractDirectedPlaneImageFilterNew.cpp Algorithms/mitkExtractImageFilter.cpp Algorithms/mitkGeometryClipImageFilter.cpp Algorithms/mitkGeometryDataSource.cpp Algorithms/mitkHeightFieldSurfaceClipImageFilter.cpp Algorithms/mitkImageToLookupTableFilter.cpp Algorithms/mitkImageToSurfaceFilter.cpp Algorithms/mitkInterpolateLinesFilter.cpp Algorithms/mitkLabeledImageToSurfaceFilter.cpp Algorithms/mitkLabeledImageVolumeCalculator.cpp Algorithms/mitkLookupTableSource.cpp Algorithms/mitkManualSegmentationToSurfaceFilter.cpp Algorithms/mitkMaskImageFilter.cpp Algorithms/mitkMeshSource.cpp Algorithms/mitkNonBlockingAlgorithm.cpp Algorithms/mitkOverwriteSliceImageFilter.cpp Algorithms/mitkOverwriteDirectedPlaneImageFilter.cpp Algorithms/mitkPadImageFilter.cpp Algorithms/mitkPlaneCutFilter.cpp Algorithms/mitkPlaneFit.cpp Algorithms/mitkPlanesPerpendicularToLinesFilter.cpp Algorithms/mitkPointLocator.cpp Algorithms/mitkPointSetToCurvedGeometryFilter.cpp Algorithms/mitkPointSetToGeometryDataFilter.cpp Algorithms/mitkPointSetIndexToWorldTransformFilter.cpp Algorithms/mitkSurfaceIndexToWorldTransformFilter.cpp Algorithms/mitkPolygonToRingFilter.cpp Algorithms/mitkProbeFilter.cpp Algorithms/mitkSegmentationSink.cpp Algorithms/mitkShapeBasedInterpolationAlgorithm.cpp Algorithms/mitkShowSegmentationAsSurface.cpp Algorithms/mitkSimpleHistogram.cpp Algorithms/mitkSimpleUnstructuredGridHistogram.cpp Algorithms/mitkSurfaceToImageFilter.cpp Algorithms/mitkUnstructuredGridHistogram.cpp Algorithms/mitkUnstructuredGridSource.cpp + Algorithms/vtkImageStencilRaster.cxx + Algorithms/vtkLassoStencilSource.cxx Algorithms/mitkVolumeVisualizationImagePreprocessor.cpp Controllers/mitkMovieGenerator.cpp Controllers/mitkMultiStepper.cpp Controllers/mitkSegmentationInterpolationController.cpp Controllers/mitkToolManager.cpp DataManagement/mitkAffineTransformationOperation.cpp DataManagement/mitkApplyDiffImageOperation.cpp DataManagement/mitkBoundingObject.cpp DataManagement/mitkBoundingObjectGroup.cpp DataManagement/mitkCellOperation.cpp DataManagement/mitkColorConversions.cpp DataManagement/mitkColorSequence.cpp DataManagement/mitkColorSequenceCycleH.cpp DataManagement/mitkColorSequenceHalfTones.cpp DataManagement/mitkColorSequenceRainbow.cpp DataManagement/mitkCompressedImageContainer.cpp DataManagement/mitkCone.cpp DataManagement/mitkContour.cpp DataManagement/mitkContourSet.cpp DataManagement/mitkCuboid.cpp DataManagement/mitkCylinder.cpp DataManagement/mitkDataStorageSelection.cpp DataManagement/mitkDelegateManager.cpp DataManagement/mitkDrawOperation.cpp DataManagement/mitkEllipsoid.cpp DataManagement/mitkExternAbstractTransformGeometry.cpp DataManagement/mitkExtrudedContour.cpp DataManagement/mitkFrameOfReferenceUIDManager.cpp DataManagement/mitkGridRepresentationProperty.cpp DataManagement/mitkGridVolumeMapperProperty.cpp DataManagement/mitkItkBaseDataAdapter.cpp DataManagement/mitkLabeledImageLookupTable.cpp DataManagement/mitkLineOperation.cpp DataManagement/mitkMesh.cpp DataManagement/mitkObjectSet.cpp DataManagement/mitkOrganTypeProperty.cpp DataManagement/mitkPlaneLandmarkProjector.cpp DataManagement/mitkPlane.cpp DataManagement/mitkPropertyManager.cpp DataManagement/mitkPropertyObserver.cpp DataManagement/mitkSeedsImage.cpp DataManagement/mitkSeedsImageLookupTableSource.cpp DataManagement/mitkSphereLandmarkProjector.cpp # DataManagement/mitkUSLookupTableSource.cpp DataManagement/mitkUnstructuredGrid.cpp DataManagement/mitkVideoSource.cpp DataManagement/vtkObjectSet.cpp IO/mitkObjFileIOFactory.cpp IO/mitkObjFileReader.cpp IO/mitkPACSPlugin.cpp IO/mitkParRecFileIOFactory.cpp IO/mitkParRecFileReader.cpp IO/mitkStlVolumeTimeSeriesIOFactory.cpp IO/mitkStlVolumeTimeSeriesReader.cpp IO/mitkUnstructuredGridVtkWriter.cpp IO/mitkUnstructuredGridVtkWriterFactory.cpp IO/mitkVtkUnstructuredGridIOFactory.cpp IO/mitkVtkUnstructuredGridReader.cpp IO/mitkVtkVolumeTimeSeriesIOFactory.cpp IO/mitkVtkVolumeTimeSeriesReader.cpp Interactions/mitkAutoSegmentationTool.cpp Interactions/mitkConferenceEventMapper.cpp Interactions/mitkConnectPointsInteractor.cpp Interactions/mitkContourInteractor.cpp Interactions/mitkContourTool.cpp #Interactions/mitkCoordinateSupplier.cpp #Interactions/mitkDisplayCoordinateOperation.cpp #Interactions/mitkDisplayInteractor.cpp Interactions/mitkAffineInteractor3D.cpp Interactions/mitkDisplayPointSetInteractor.cpp #Interactions/mitkDisplayVectorInteractor.cpp Interactions/mitkExtrudedContourInteractor.cpp Interactions/mitkFeedbackContourTool.cpp Interactions/mitkInteractionDebug.cpp Interactions/mitkInteractionDebugger.cpp Interactions/mitkPaintbrushTool.cpp Interactions/mitkPointInteractor.cpp Interactions/mitkPointSelectorInteractor.cpp #Interactions/mitkPositionTracker.cpp Interactions/mitkSeedsInteractor.cpp Interactions/mitkSegTool2D.cpp Interactions/mitkSegmentationsProcessingTool.cpp Interactions/mitkSetRegionTool.cpp Interactions/mitkSocketClient.cpp Interactions/mitkSurfaceDeformationInteractor3D.cpp Interactions/mitkSurfaceInteractor.cpp Interactions/mitkTool.cpp Interactions/mitkAddContourTool.cpp Interactions/mitkAutoCropTool.cpp Interactions/mitkBinaryThresholdTool.cpp Interactions/mitkCalculateGrayValueStatisticsTool.cpp Interactions/mitkCalculateVolumetryTool.cpp Interactions/mitkCorrectorTool2D.cpp Interactions/mitkCreateSurfaceTool.cpp Interactions/mitkEraseRegionTool.cpp Interactions/mitkFillRegionTool.cpp Interactions/mitkRegionGrowingTool.cpp Interactions/mitkSubtractContourTool.cpp Interactions/mitkDrawPaintbrushTool.cpp Interactions/mitkErasePaintbrushTool.cpp Interactions/mitkMorphologicTool.cpp Interactions/mitkErodeTool.cpp Interactions/mitkDilateTool.cpp Interactions/mitkOpeningTool.cpp Interactions/mitkClosingTool.cpp Interactions/mitkBinaryThresholdULTool.cpp Interactions/mitkPixelManipulationTool.cpp Interactions/mitkRegionGrow3DTool.cpp Rendering/mitkContourMapper2D.cpp Rendering/mitkContourSetMapper2D.cpp Rendering/mitkContourSetVtkMapper3D.cpp Rendering/mitkContourVtkMapper3D.cpp Rendering/mitkEnhancedPointSetVtkMapper3D.cpp Rendering/mitkImageBackground2D.cpp Rendering/mitkLineMapper2D.cpp # Rendering/mitkLineVtkMapper3D.cpp Rendering/mitkMeshMapper2D.cpp Rendering/mitkMeshVtkMapper3D.cpp Rendering/mitkNativeRenderWindowInteractor.cpp Rendering/mitkSplineMapper2D.cpp Rendering/mitkSplineVtkMapper3D.cpp Rendering/mitkUnstructuredGridMapper2D.cpp Rendering/mitkUnstructuredGridVtkMapper3D.cpp Rendering/mitkVectorImageMapper2D.cpp Rendering/vtkUnstructuredGridMapper.cpp Rendering/vtkMaskedGlyph2D.cpp Rendering/vtkMaskedGlyph3D.cpp Rendering/vtkMitkVolumeTextureMapper3D.cpp Rendering/vtkMitkOpenGLVolumeTextureMapper3D.cpp Rendering/mitkGPUVolumeMapper3D.cpp Rendering/vtkMitkGPUVolumeRayCastMapper.cpp Rendering/vtkMitkOpenGLGPUVolumeRayCastMapper.cpp Rendering/vtkMitkOpenGLGPUVolumeRayCastMapperShaders.cpp ) IF(WIN32 AND NOT MINGW) SET(CPP_FILES Controllers/mitkMovieGeneratorWin32.cpp ${CPP_FILES} ) ENDIF(WIN32 AND NOT MINGW)