diff --git a/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx index 6b00e1f26f..43d62a862c 100644 --- a/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx +++ b/Modules/MitkExt/Algorithms/vtkImageStencilRaster.cxx @@ -1,452 +1,453 @@ /*========================================================================= 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 "mitkvtkImageStencilRaster.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/vtkLassoStencilSource.cxx b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx index dbe66d903d..e86d4ad535 100644 --- a/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx +++ b/Modules/MitkExt/Algorithms/vtkLassoStencilSource.cxx @@ -1,609 +1,610 @@ /*========================================================================= 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 "mitkvtkLassoStencilSource.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 "mitkvtkImageStencilRaster.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/Testing/Data/testimage.dcm b/Modules/MitkExt/Testing/Data/testimage.dcm new file mode 100644 index 0000000000..ecf7b12d33 Binary files /dev/null and b/Modules/MitkExt/Testing/Data/testimage.dcm differ diff --git a/Modules/MitkExt/Testing/files.cmake b/Modules/MitkExt/Testing/files.cmake index 7062436d27..1256295f5e 100644 --- a/Modules/MitkExt/Testing/files.cmake +++ b/Modules/MitkExt/Testing/files.cmake @@ -1,39 +1,40 @@ SET(MODULE_TESTS mitkAutoCropImageFilterTest.cpp mitkBoundingObjectCutterTest.cpp mitkContourMapper2DTest.cpp mitkContourTest.cpp mitkCoreExtObjectFactoryTest mitkDataNodeExtTest.cpp mitkExternalToolsTest.cpp + mitkImageStatisticsCalculatorTest.cpp mitkMeshTest.cpp mitkMultiStepperTest.cpp mitkOrganTypePropertyTest.cpp mitkPipelineSmartPointerCorrectnessTest.cpp mitkPlaneFitTest.cpp mitkPointLocatorTest.cpp # mitkSegmentationInterpolationTest.cpp # mitkTestTemplate.cpp mitkToolManagerTest.cpp ) SET(MODULE_IMAGE_TESTS mitkUnstructuredGridVtkWriterTest.cpp mitkCompressedImageContainerTest.cpp mitkCylindricToCartesianFilterTest.cpp mitkExtractImageFilterTest.cpp mitkManualSegmentationToSurfaceFilterTest.cpp mitkOverwriteSliceImageFilterTest.cpp mitkSurfaceToImageFilterTest.cpp ) SET(MODULE_CUSTOM_TESTS mitkLabeledImageToSurfaceFilterTest.cpp ) SET(MODULE_TESTIMAGES US4DCyl.pic.gz Pic3D.pic.gz Pic2DplusT.pic.gz BallBinary30x30x30.pic.gz Png2D-bw.png binary.stl ball.stl ) diff --git a/Modules/MitkExt/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/MitkExt/Testing/mitkImageStatisticsCalculatorTest.cpp new file mode 100644 index 0000000000..49f8125849 --- /dev/null +++ b/Modules/MitkExt/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -0,0 +1,471 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2008-02-25 17:27:17 +0100 (Mo, 25 Feb 2008) $ +Version: $Revision: 7837 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html 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 notices for more information. + +=========================================================================*/ + +#include "mitkStandardFileLocations.h" +#include "mitkDicomSeriesReader.h" +#include "mitkTestingMacros.h" +#include "mitkImageStatisticsCalculator.h" +#include "mitkPlanarPolygon.h" + +#include "mitkDicomSeriesReader.h" +#include + +#include "vtkStreamingDemandDrivenPipeline.h" + + +//#include + +/** + * \brief Test class for mitkImageStatisticsCalculator + * + * This test covers: + * - instantiation of an ImageStatisticsCalculator class + * - correctness of statistics when using PlanarFigures for masking + */ +class mitkImageStatisticsCalculatorTestClass +{ + +public: + + struct testCase + { + int id; + mitk::PlanarFigure::Pointer figure; + double mean; + double sd; + }; + + +// calculate statistics for the given image and planarpolygon +static const mitk::ImageStatisticsCalculator::Statistics TestStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) +{ + mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); + statisticsCalculator->SetImage( image ); + statisticsCalculator->SetMaskingModeToPlanarFigure(); + statisticsCalculator->SetPlanarFigure( polygon ); + + statisticsCalculator->ComputeStatistics(); + + return statisticsCalculator->GetStatistics(); +} + + +// returns a vector of defined test-cases +static std::vector InitializeTestCases( mitk::Geometry2D::Pointer geom ) +{ + std::vector testCases; + + { + /***************************** + * one whole white pixel + * -> mean of 255 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 255.0; + test.sd = 0.0; + + testCases.push_back( test ); + } + + { + /***************************** + * half pixel in x-direction (white) + * -> mean of 255 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 10.0 ; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 255.0; + test.sd = 0.0; + + testCases.push_back( test ); + } + + { + /***************************** + * half pixel in diagonal-direction (white) + * -> mean of 255 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; + figure1->SetControlPoint( 2, pnt3, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 255.0; + test.sd = 0.0; + + testCases.push_back( test ); + } + + + { + /***************************** + * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) + * -> mean of 191.25 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 191.25; + test.sd = 127.5; + + testCases.push_back( test ); + } + + + { + /***************************** + * whole pixel (white) + half pixel (gray) in x-direction + * -> mean of 191.5 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 191.50; + test.sd = 89.80; + + testCases.push_back( test ); + } + + { + /***************************** + * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction + * -> mean of 191.5 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 191.5; + test.sd = 89.80; + + testCases.push_back( test ); + } + + { + /***************************** + * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction + * -> mean of 127.66 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); + figure1->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; + figure1->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; + figure1->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; + figure1->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; + figure1->SetControlPoint( 3, pnt4, true ); + figure1->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure1; + test.mean = 127.66; + test.sd = 127.5; + + testCases.push_back( test ); + } + + + { + /***************************** + * whole pixel (gray) + * -> mean of 128 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); + figure2->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; + figure2->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; + figure2->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; + figure2->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; + figure2->SetControlPoint( 3, pnt4, true ); + figure2->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure2; + test.mean = 128.0; + test.sd = 0.0; + + testCases.push_back( test ); + } + + { + /***************************** + * whole pixel (gray) + half pixel (white) in y-direction + * -> mean of 191.5 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); + figure2->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; + figure2->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; + figure2->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; + figure2->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; + figure2->SetControlPoint( 3, pnt4, true ); + figure2->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure2; + test.mean = 191.5; + test.sd = 89.80; + + testCases.push_back( test ); + } + + { + /***************************** + * 2 whole pixel (white) + 2 whole pixel (black) in y-direction + * -> mean of 127.66 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); + figure2->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; + figure2->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; + figure2->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; + figure2->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; + figure2->SetControlPoint( 3, pnt4, true ); + figure2->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure2; + test.mean = 127.66; + test.sd = 127.5; + + testCases.push_back( test ); + } + + + { + /***************************** + * 9 whole pixels (white) + 3 half pixels (white) + * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)] + * -> mean of 204.0 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); + figure2->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; + figure2->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; + figure2->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; + figure2->SetControlPoint( 2, pnt3, true ); + mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; + figure2->SetControlPoint( 3, pnt4, true ); + figure2->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure2; + test.mean = 204.0; + test.sd = 105.58; + + testCases.push_back( test ); + } + + { + /***************************** + * half pixel (white) + whole pixel (white) + half pixel (black) + * -> mean of 212.66 expected + ******************************/ + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); + figure2->SetGeometry2D( geom ); + mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; + figure2->PlaceFigure( pnt1 ); + + mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; + figure2->SetControlPoint( 1, pnt2, true ); + mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; + figure2->SetControlPoint( 2, pnt3, true ); + figure2->GetPolyLine(0); + + testCase test; + test.id = testCases.size(); + test.figure = figure2; + test.mean = 212.66; + test.sd = 73.32; + + testCases.push_back( test ); + } + + + + return testCases; +} + +// loads the test image +static mitk::Image::Pointer GetTestImage() +{ + mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); + + const std::string filename = locator->FindFile("testimage.dcm", "Modules/MitkExt/Testing/Data"); + if (filename.empty()) + { + MITK_ERROR << "Could not find test file"; + return NULL; + } + else + { + MITK_INFO << "Found testimage.dcm"; + } + + + itk::FilenamesContainer file; + file.push_back( filename ); + + mitk::DicomSeriesReader* reader = new mitk::DicomSeriesReader; + + mitk::DataNode::Pointer node = reader->LoadDicomSeries( file, false, false ); + mitk::Image::Pointer image = dynamic_cast( node->GetData() ); + + return image; +} + +}; + +int mitkImageStatisticsCalculatorTest(int argc, char* argv[]) +{ + // always start with this! + MITK_TEST_BEGIN("mitkImageStatisticsCalculatorTest") + + //QCoreApplication app(argc, argv); + + mitk::Image::Pointer image = mitkImageStatisticsCalculatorTestClass::GetTestImage(); + MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Loading test image" ); + + mitk::Geometry2D::Pointer geom = image->GetSlicedGeometry()->GetGeometry2D(0); + + std::vector allTestCases = + mitkImageStatisticsCalculatorTestClass::InitializeTestCases( geom ); + + + for ( int i=0; i