diff --git a/Core/Code/Algorithms/mitkVtkImageIdxReslice.cpp b/Core/Code/Algorithms/mitkVtkImageIdxReslice.cpp index 496fc38bdb..2cb24a5013 100644 --- a/Core/Code/Algorithms/mitkVtkImageIdxReslice.cpp +++ b/Core/Code/Algorithms/mitkVtkImageIdxReslice.cpp @@ -1,933 +1,901 @@ /*========================================================================= Program: Visualization Toolkit Module: mitkVtkImageIdxReslice.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 "mitkVtkImageIdxReslice.h" #include "vtkImageData.h" #include "vtkImageStencilData.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMath.h" #include "vtkObjectFactory.h" #include "vtkStreamingDemandDrivenPipeline.h" #include "vtkTransform.h" #include "vtkDataSetAttributes.h" #include "vtkGarbageCollector.h" #include "vtkTemplateAliasMacro.h" // turn off 64-bit ints when templating over all types # undef VTK_USE_INT64 # define VTK_USE_INT64 0 # undef VTK_USE_UINT64 # define VTK_USE_UINT64 0 #include #include #include vtkStandardNewMacro(mitkVtkImageIdxReslice); bool Overwrite_Mode = false; //mitkVtkImageIdxReslice* self_obj; //-------------------------------------------------------------------------- // The 'floor' function on x86 and mips is many times slower than these // and is used a lot in this code, optimize for different CPU architectures template inline int vtkResliceFloor(double x, F &f) { #if defined mips || defined sparc || defined __ppc__ x += 2147483648.0; unsigned int i = static_cast(x); f = x - i; return static_cast(i - 2147483648U); #elif defined i386 || defined _M_IX86 union { double d; unsigned short s[4]; unsigned int i[2]; } dual; dual.d = x + 103079215104.0; // (2**(52-16))*1.5 f = dual.s[0]*0.0000152587890625; // 2**(-16) return static_cast((dual.i[1]<<16)|((dual.i[0])>>16)); #elif defined ia64 || defined __ia64__ || defined IA64 x += 103079215104.0; long long i = static_cast(x); f = x - i; return static_cast(i - 103079215104LL); #else double y = floor(x); f = x - y; return static_cast(y); #endif } inline int vtkResliceRound(double x) { #if defined mips || defined sparc || defined __ppc__ return static_cast(static_cast(x + 2147483648.5) - 2147483648U); #elif defined i386 || defined _M_IX86 union { double d; unsigned int i[2]; } dual; dual.d = x + 103079215104.5; // (2**(52-16))*1.5 return static_cast((dual.i[1]<<16)|((dual.i[0])>>16)); #elif defined ia64 || defined __ia64__ || defined IA64 x += 103079215104.5; long long i = static_cast(x); return static_cast(i - 103079215104LL); #else return static_cast(floor(x+0.5)); #endif } //---------------------------------------------------------------------------- mitkVtkImageIdxReslice::mitkVtkImageIdxReslice() { this->GetOutput()->SetScalarTypeToUnsignedInt(); } //---------------------------------------------------------------------------- mitkVtkImageIdxReslice::~mitkVtkImageIdxReslice() { Overwrite_Mode = false; } //---------------------------------------------------------------------------- // constants for different boundary-handling modes #define VTK_RESLICE_BACKGROUND 0 // use background if out-of-bounds #define VTK_RESLICE_WRAP 1 // wrap to opposite side of image #define VTK_RESLICE_MIRROR 2 // mirror off of the boundary #define VTK_RESLICE_BORDER 3 // use a half-voxel border #define VTK_RESLICE_NULL 4 // do nothing to *outPtr if out-of-bounds //---------------------------------------------------------------------------- // rounding functions for each type, where 'F' is a floating-point type #if (VTK_USE_INT8 != 0) template inline void vtkResliceRound(F val, vtkTypeInt8& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_UINT8 != 0) template inline void vtkResliceRound(F val, vtkTypeUInt8& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_INT16 != 0) template inline void vtkResliceRound(F val, vtkTypeInt16& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_UINT16 != 0) template inline void vtkResliceRound(F val, vtkTypeUInt16& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_INT32 != 0) template inline void vtkResliceRound(F val, vtkTypeInt32& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_UINT32 != 0) template inline void vtkResliceRound(F val, vtkTypeUInt32& rnd) { rnd = vtkResliceRound(val); } #endif #if (VTK_USE_FLOAT32 != 0) template inline void vtkResliceRound(F val, vtkTypeFloat32& rnd) { rnd = val; } #endif #if (VTK_USE_FLOAT64 != 0) template inline void vtkResliceRound(F val, vtkTypeFloat64& rnd) { rnd = val; } #endif //---------------------------------------------------------------------------- // clamping functions for each type #if (VTK_USE_INT8 != 0) template inline void vtkResliceClamp(F val, vtkTypeInt8& clamp) { if (val < -128.0) { val = -128.0; } if (val > 127.0) { val = 127.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_UINT8 != 0) template inline void vtkResliceClamp(F val, vtkTypeUInt8& clamp) { if (val < 0) { val = 0; } if (val > 255.0) { val = 255.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_INT16 != 0) template inline void vtkResliceClamp(F val, vtkTypeInt16& clamp) { if (val < -32768.0) { val = -32768.0; } if (val > 32767.0) { val = 32767.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_UINT16 != 0) template inline void vtkResliceClamp(F val, vtkTypeUInt16& clamp) { if (val < 0) { val = 0; } if (val > 65535.0) { val = 65535.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_INT32 != 0) template inline void vtkResliceClamp(F val, vtkTypeInt32& clamp) { if (val < -2147483648.0) { val = -2147483648.0; } if (val > 2147483647.0) { val = 2147483647.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_UINT32 != 0) template inline void vtkResliceClamp(F val, vtkTypeUInt32& clamp) { if (val < 0) { val = 0; } if (val > 4294967295.0) { val = 4294967295.0; } vtkResliceRound(val,clamp); } #endif #if (VTK_USE_FLOAT32 != 0) template inline void vtkResliceClamp(F val, vtkTypeFloat32& clamp) { clamp = val; } #endif #if (VTK_USE_FLOAT64 != 0) template inline void vtkResliceClamp(F val, vtkTypeFloat64& clamp) { clamp = val; } #endif //---------------------------------------------------------------------------- // Perform a wrap to limit an index to [0,range). // Ensures correct behaviour when the index is negative. inline int vtkInterpolateWrap(int num, int range) { if ((num %= range) < 0) { num += range; // required for some % implementations } return num; } //---------------------------------------------------------------------------- // Perform a mirror to limit an index to [0,range). inline int vtkInterpolateMirror(int num, int range) { if (num < 0) { num = -num - 1; } int count = num/range; num %= range; if (count & 0x1) { num = range - num - 1; } return num; } //---------------------------------------------------------------------------- // If the value is within one half voxel of the range [0,inExtX), then // set it to "0" or "inExtX-1" as appropriate. inline int vtkInterpolateBorder(int &inIdX0, int &inIdX1, int inExtX, double fx) { if (inIdX0 >= 0 && inIdX1 < inExtX) { return 0; } if (inIdX0 == -1 && fx >= 0.5) { inIdX1 = inIdX0 = 0; return 0; } if (inIdX0 == inExtX - 1 && fx < 0.5) { inIdX1 = inIdX0; return 0; } return 1; } inline int vtkInterpolateBorderCheck(int inIdX0, int inIdX1, int inExtX, double fx) { if ((inIdX0 >= 0 && inIdX1 < inExtX) || (inIdX0 == -1 && fx >= 0.5) || (inIdX0 == inExtX - 1 && fx < 0.5)) { return 0; } return 1; } //---------------------------------------------------------------------------- // Do nearest-neighbor interpolation of the input data 'inPtr' of extent // 'inExt' at the 'point'. The result is placed at 'outPtr'. // If the lookup data is beyond the extent 'inExt', set 'outPtr' to // the background color 'background'. // The number of scalar components in the data is 'numscalars' template static int vtkNearestNeighborInterpolation(T *&outPtr, const T *inPtr, const int inExt[6], const vtkIdType inInc[3], int numscalars, const F point[3], int mode, const T *background) { int inIdX0 = vtkResliceRound(point[0]) - inExt[0]; int inIdY0 = vtkResliceRound(point[1]) - inExt[2]; int inIdZ0 = vtkResliceRound(point[2]) - inExt[4]; int inExtX = inExt[1] - inExt[0] + 1; int inExtY = inExt[3] - inExt[2] + 1; int inExtZ = inExt[5] - inExt[4] + 1; if (inIdX0 < 0 || inIdX0 >= inExtX || inIdY0 < 0 || inIdY0 >= inExtY || inIdZ0 < 0 || inIdZ0 >= inExtZ) { if (mode == VTK_RESLICE_WRAP) { inIdX0 = vtkInterpolateWrap(inIdX0, inExtX); inIdY0 = vtkInterpolateWrap(inIdY0, inExtY); inIdZ0 = vtkInterpolateWrap(inIdZ0, inExtZ); } else if (mode == VTK_RESLICE_MIRROR) { inIdX0 = vtkInterpolateMirror(inIdX0, inExtX); inIdY0 = vtkInterpolateMirror(inIdY0, inExtY); inIdZ0 = vtkInterpolateMirror(inIdZ0, inExtZ); } else if (mode == VTK_RESLICE_BACKGROUND || mode == VTK_RESLICE_BORDER) { do { *outPtr++ = *background++; } while (--numscalars); return 0; } else { return 0; } } inPtr += inIdX0*inInc[0]+inIdY0*inInc[1]+inIdZ0*inInc[2]; do { if(!Overwrite_Mode) { *outPtr++ = *inPtr++; } else { *(const_cast(inPtr)) = *outPtr++; inPtr++; } } while (--numscalars); return 1; } //-------------------------------------------------------------------------- // get appropriate interpolation function according to interpolation mode // and scalar type template static void vtkGetResliceInterpFunc(mitkVtkImageIdxReslice *self, int (**interpolate)(void *&outPtr, const void *inPtr, const int inExt[6], const vtkIdType inInc[3], int numscalars, const F point[3], int mode, const void *background)) { int dataType = self->GetOutput()->GetScalarType(); switch (dataType) { vtkTemplateAliasMacro(*((int (**)(VTK_TT *&outPtr, const VTK_TT *inPtr, const int inExt[6], const vtkIdType inInc[3], int numscalars, const F point[3], int mode, const VTK_TT *background))interpolate) = \ &vtkNearestNeighborInterpolation); default: interpolate = 0; } } //---------------------------------------------------------------------------- // Some helper functions for 'RequestData' //---------------------------------------------------------------------------- //-------------------------------------------------------------------------- // pixel copy function, templated for different scalar types template struct vtkImageResliceSetPixels { static void Set(void *&outPtrV, const void *inPtrV, int numscalars, int n) { const T* inPtr = static_cast(inPtrV); T* outPtr = static_cast(outPtrV); for (int i = 0; i < n; i++) { const T *tmpPtr = inPtr; int m = numscalars; do { *outPtr++ = *tmpPtr++; } while (--m); } outPtrV = outPtr; } // optimized for 1 scalar components static void Set1(void *&outPtrV, const void *inPtrV, int vtkNotUsed(numscalars), int n) { const T* inPtr = static_cast(inPtrV); T* outPtr = static_cast(outPtrV); T val = *inPtr; for (int i = 0; i < n; i++) { *outPtr++ = val; } outPtrV = outPtr; } }; // get a pixel copy function that is appropriate for the data type static void vtkGetSetPixelsFunc(mitkVtkImageIdxReslice *self, void (**setpixels)(void *&out, const void *in, int numscalars, int n)) { int dataType = self->GetOutput()->GetScalarType(); int numscalars = self->GetOutput()->GetNumberOfScalarComponents(); switch (numscalars) { case 1: switch (dataType) { vtkTemplateAliasMacro( *setpixels = &vtkImageResliceSetPixels::Set1 ); default: setpixels = 0; } default: switch (dataType) { vtkTemplateAliasMacro( *setpixels = &vtkImageResliceSetPixels::Set ); default: setpixels = 0; } } } //---------------------------------------------------------------------------- // Convert background color from float to appropriate type template static void vtkAllocBackgroundPixelT(mitkVtkImageIdxReslice *self, T **background_ptr, int numComponents) { *background_ptr = new T[numComponents]; T *background = *background_ptr; for (int i = 0; i < numComponents; i++) { if (i < 4) { vtkResliceClamp(self->GetBackgroundColor()[i], background[i]); } else { background[i] = 0; } } } static void vtkAllocBackgroundPixel(mitkVtkImageIdxReslice *self, void **rval, int numComponents) { switch (self->GetOutput()->GetScalarType()) { vtkTemplateAliasMacro(vtkAllocBackgroundPixelT(self, (VTK_TT **)rval, numComponents)); } } static void vtkFreeBackgroundPixel(mitkVtkImageIdxReslice *self, void **rval) { switch (self->GetOutput()->GetScalarType()) { vtkTemplateAliasMacro(delete [] *((VTK_TT **)rval)); } *rval = 0; } //---------------------------------------------------------------------------- // helper function for clipping of the output with a stencil static int vtkResliceGetNextExtent(vtkImageStencilData *stencil, int &r1, int &r2, int rmin, int rmax, int yIdx, int zIdx, void *&outPtr, void *background, int numscalars, void (*setpixels)(void *&out, const void *in, int numscalars, int n), int &iter) { // trivial case if stencil is not set if (!stencil) { if (iter++ == 0) { r1 = rmin; r2 = rmax; return 1; } return 0; } // for clearing, start at last r2 plus 1 int clear1 = r2 + 1; if (iter == 0) { // if no 'last time', start at rmin clear1 = rmin; } int rval = stencil->GetNextExtent(r1, r2, rmin, rmax, yIdx, zIdx, iter); int clear2 = r1 - 1; if (rval == 0) { clear2 = rmax; } setpixels(outPtr, background, numscalars, clear2 - clear1 + 1); return rval; } //---------------------------------------------------------------------------- // This function simply clears the entire output to the background color, // for cases where the transformation places the output extent completely // outside of the input extent. static void vtkImageResliceClearExecute(mitkVtkImageIdxReslice *self, vtkImageData *, void *, vtkImageData *outData, void *outPtr, int outExt[6], int id) { int numscalars; int idY, idZ; vtkIdType outIncX, outIncY, outIncZ; int scalarSize; unsigned long count = 0; unsigned long target; void *background; void (*setpixels)(void *&out, const void *in, int numscalars, int n); // for the progress meter target = static_cast ((outExt[5]-outExt[4]+1)*(outExt[3]-outExt[2]+1)/50.0); target++; // Get Increments to march through data outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ); scalarSize = outData->GetScalarSize(); numscalars = outData->GetNumberOfScalarComponents(); // allocate a voxel to copy into the background (out-of-bounds) regions vtkAllocBackgroundPixel(self, &background, numscalars); // get the appropriate function for pixel copying vtkGetSetPixelsFunc(self, &setpixels); // Loop through output voxels for (idZ = outExt[4]; idZ <= outExt[5]; idZ++) { for (idY = outExt[2]; idY <= outExt[3]; idY++) { if (id == 0) { // update the progress if this is the main thread if (!(count%target)) { self->UpdateProgress(count/(50.0*target)); } count++; } // clear the pixels to background color and go to next row setpixels(outPtr, background, numscalars, outExt[1]-outExt[0]+1); outPtr = static_cast( static_cast(outPtr) + outIncY*scalarSize); } outPtr = static_cast( static_cast(outPtr) + outIncZ*scalarSize); } vtkFreeBackgroundPixel(self, &background); } //---------------------------------------------------------------------------- // This function executes the filter for any type of data. It is much simpler // in structure than vtkImageResliceOptimizedExecute. static void vtkImageResliceExecute(mitkVtkImageIdxReslice *self, vtkImageData *inData, void *inPtr, vtkImageData *outData, void *outPtr, int outExt[6], int id) { int numscalars; int idX, idY, idZ; int idXmin, idXmax, iter; vtkIdType outIncX, outIncY, outIncZ; int scalarSize; int inExt[6]; vtkIdType inInc[3]; unsigned long count = 0; unsigned long target; double point[4]; double f; double *inSpacing, *inOrigin, *outSpacing, *outOrigin, inInvSpacing[3]; void *background; int (*interpolate)(void *&outPtr, const void *inPtr, const int inExt[6], const vtkIdType inInc[3], int numscalars, const double point[3], int mode, const void *background); void (*setpixels)(void *&out, const void *in, int numscalars, int n); // the 'mode' species what to do with the 'pad' (out-of-bounds) area int mode = VTK_RESLICE_BACKGROUND; if (self->GetMirror()) { mode = VTK_RESLICE_MIRROR; } else if (self->GetWrap()) { mode = VTK_RESLICE_WRAP; } else if (self->GetBorder()) { mode = VTK_RESLICE_BORDER; } // the transformation to apply to the data vtkAbstractTransform *transform = self->GetResliceTransform(); vtkMatrix4x4 *matrix = self->GetResliceAxes(); // for conversion to data coordinates inOrigin = inData->GetOrigin(); inSpacing = inData->GetSpacing(); outOrigin = outData->GetOrigin(); outSpacing = outData->GetSpacing(); // save effor later: invert inSpacing inInvSpacing[0] = 1.0/inSpacing[0]; inInvSpacing[1] = 1.0/inSpacing[1]; inInvSpacing[2] = 1.0/inSpacing[2]; // find maximum input range inData->GetExtent(inExt); // for the progress meter target = static_cast ((outExt[5]-outExt[4]+1)*(outExt[3]-outExt[2]+1)/50.0); target++; // Get Increments to march through data inData->GetIncrements(inInc); outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ); scalarSize = outData->GetScalarSize(); numscalars = inData->GetNumberOfScalarComponents(); // allocate a voxel to copy into the background (out-of-bounds) regions vtkAllocBackgroundPixel(self, &background, numscalars); // get the appropriate functions for interpolation and pixel copying vtkGetResliceInterpFunc(self, &interpolate); vtkGetSetPixelsFunc(self, &setpixels); // get the stencil vtkImageStencilData *stencil = self->GetStencil(); // Loop through output voxels for (idZ = outExt[4]; idZ <= outExt[5]; idZ++) { for (idY = outExt[2]; idY <= outExt[3]; idY++) { if (id == 0) { // update the progress if this is the main thread if (!(count%target)) { self->UpdateProgress(count/(50.0*target)); } count++; } iter = 0; // if there is a stencil, it is applied here while (vtkResliceGetNextExtent(stencil, idXmin, idXmax, outExt[0], outExt[1], idY, idZ, outPtr, background, numscalars, setpixels, iter)) { for (idX = idXmin; idX <= idXmax; idX++) { // convert to data coordinates point[0] = idX*outSpacing[0] + outOrigin[0]; point[1] = idY*outSpacing[1] + outOrigin[1]; point[2] = idZ*outSpacing[2] + outOrigin[2]; // apply ResliceAxes matrix if (matrix) { point[3] = 1.0; matrix->MultiplyPoint(point, point); f = 1.0/point[3]; point[0] *= f; point[1] *= f; point[2] *= f; } // apply ResliceTransform if (transform) { transform->InternalTransformPoint(point, point); } // convert back to voxel indices point[0] = (point[0] - inOrigin[0])*inInvSpacing[0]; point[1] = (point[1] - inOrigin[1])*inInvSpacing[1]; point[2] = (point[2] - inOrigin[2])*inInvSpacing[2]; // interpolate output voxel from input data set interpolate(outPtr, inPtr, inExt, inInc, numscalars, point, mode, background); } } outPtr = static_cast( static_cast(outPtr) + outIncY*scalarSize); } outPtr = static_cast( static_cast(outPtr) + outIncZ*scalarSize); } vtkFreeBackgroundPixel(self, &background); } void mitkVtkImageIdxReslice::SetOverwriteMode(bool b){ Overwrite_Mode = b; } void mitkVtkImageIdxReslice::SetInputSlice(vtkImageData* slice){ this->SetOutput(slice); } -void mitkVtkImageIdxReslice::SetOutputExtent(int ex[]){SetOutputExtent(ex[0],ex[1],ex[2],ex[3],ex[4],ex[5]);} -void mitkVtkImageIdxReslice::SetOutputExtent(int xMin, int xMax, int yMin, int yMax, int zMin, int zMax){ - //if(Overwrite_Mode){ - /*if(xMin<0) - { - xMax-=xMin; - xMin=0; - } - if(yMin<0) - { - yMax-=yMin; - yMin=0; - }*/ - //} - - if(Overwrite_Mode){ - this->OutputOrigin[0] = -xMin; - this->OutputOrigin[1] = -yMin; - this->OutputOrigin[2] = 0.0; - }else{ - this->OutputOrigin[0] = xMin; - this->OutputOrigin[1] = yMin; - } - Superclass::SetOutputExtent( xMin, xMax, yMin, yMax, zMin, zMax); -} - -void mitkVtkImageIdxReslice::SetOutputOrigin(double or[]){this->SetOutputOrigin(or[0],or[1],or[2]);} -void mitkVtkImageIdxReslice::SetOutputOrigin(double x, double y, double z){ - //Superclass::SetOutputOrigin( x,y,z ); - -} - //---------------------------------------------------------------------------- // This method is passed a input and output region, and executes the filter // algorithm to fill the output from the input. // It just executes a switch statement to call the correct function for // the regions data types. void mitkVtkImageIdxReslice::ThreadedRequestData( vtkInformation *vtkNotUsed(request), vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *vtkNotUsed(outputVector), vtkImageData ***inData, vtkImageData **outData, int outExt[6], int id) { vtkDebugMacro(<< "Execute: inData = " << inData[0][0] << ", outData = " << outData[0]); // this filter expects that input is the same type as output. if (inData[0][0]->GetScalarType() != outData[0]->GetScalarType()) { vtkErrorMacro(<< "Execute: input ScalarType, " << inData[0][0]->GetScalarType() << ", must match out ScalarType " << outData[0]->GetScalarType()); return; } int inExt[6]; inData[0][0]->GetExtent(inExt); // check for empty input extent if (inExt[1] < inExt[0] || inExt[3] < inExt[2] || inExt[5] < inExt[4]) { return; } // Get the output pointer void *outPtr = outData[0]->GetScalarPointerForExtent(outExt); if (this->HitInputExtent == 0) { vtkImageResliceClearExecute(this, inData[0][0], 0, outData[0], outPtr, outExt, id); return; } // Now that we know that we need the input, get the input pointer void *inPtr = inData[0][0]->GetScalarPointerForExtent(inExt); vtkImageResliceExecute(this, inData[0][0], inPtr, outData[0], outPtr, outExt, id); } diff --git a/Core/Code/Algorithms/mitkVtkImageIdxReslice.h b/Core/Code/Algorithms/mitkVtkImageIdxReslice.h index 14a7059cb4..1036564857 100644 --- a/Core/Code/Algorithms/mitkVtkImageIdxReslice.h +++ b/Core/Code/Algorithms/mitkVtkImageIdxReslice.h @@ -1,57 +1,52 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision: $ 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. =========================================================================*/ #ifndef mitkVtkImageIdxReslice_h_Included #define mitkVtkImageIdxReslice_h_Included #include #include "MitkExports.h" /** \brief A Filter based on vtkImageReslice that extracts Coordinates in the 3D volume not pixel values. * */ class MITK_CORE_EXPORT mitkVtkImageIdxReslice : public vtkImageReslice { public: static mitkVtkImageIdxReslice *New(); vtkTypeMacro(mitkVtkImageIdxReslice, vtkImageReslice); void SetOverwriteMode(bool b); void SetInputSlice(vtkImageData* slice); - virtual void SetOutputExtent(int ex[]); - virtual void SetOutputExtent(int xMin, int xMax, int yMin, int yMax, int zMin, int zMax); - - virtual void SetOutputOrigin(double or[]); - virtual void SetOutputOrigin(double x, double y, double z); protected: mitkVtkImageIdxReslice(); virtual ~mitkVtkImageIdxReslice(); virtual void ThreadedRequestData(vtkInformation *vtkNotUsed(request), vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *vtkNotUsed(outputVector), vtkImageData ***inData, vtkImageData **outData, int outExt[6], int id); }; #endif //mitkVtkImageIdxReslice_h_Included diff --git a/Modules/MitkExt/Interactions/mitkSegTool2D.cpp b/Modules/MitkExt/Interactions/mitkSegTool2D.cpp index a512b9b550..a8df5fdfbb 100644 --- a/Modules/MitkExt/Interactions/mitkSegTool2D.cpp +++ b/Modules/MitkExt/Interactions/mitkSegTool2D.cpp @@ -1,365 +1,374 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ 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 "mitkSegTool2D.h" #include "mitkToolManager.h" #include "mitkDataStorage.h" #include "mitkBaseRenderer.h" #include "mitkPlaneGeometry.h" #include "mitkExtractImageFilter.h" #include "mitkExtractDirectedPlaneImageFilter.h" //Include of the new ImageExtractor #include "mitkExtractDirectedPlaneImageFilterNew.h" #include "mitkPlanarCircle.h" #include "mitkOverwriteSliceImageFilter.h" #include "mitkOverwriteDirectedPlaneImageFilter.h" #include "mitkGetModuleContext.h" //Includes for 3DSurfaceInterpolation #include "mitkImageToContourFilter.h" #include "mitkSurfaceInterpolationController.h" - +#include +#include #include #include #define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) mitk::SegTool2D::SegTool2D(const char* type) :Tool(type), m_LastEventSender(NULL), m_LastEventSlice(0), m_Contourmarkername ("Position"), m_ShowMarkerNodes (true), m_3DInterpolationEnabled (true) { // great magic numbers CONNECT_ACTION( 80, OnMousePressed ); CONNECT_ACTION( 90, OnMouseMoved ); CONNECT_ACTION( 42, OnMouseReleased ); CONNECT_ACTION( 49014, OnInvertLogic ); } mitk::SegTool2D::~SegTool2D() { - m_Map = NULL; - } bool mitk::SegTool2D::OnMousePressed (Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( positionEvent->GetSender()->GetMapperID() != BaseRenderer::Standard2D ) return false; // we don't want anything but 2D m_LastEventSender = positionEvent->GetSender(); m_LastEventSlice = m_LastEventSender->GetSlice(); return true; } bool mitk::SegTool2D::OnMouseMoved (Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::OnMouseReleased(Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::OnInvertLogic(Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::DetermineAffectedImageSlice( const Image* image, const PlaneGeometry* plane, int& affectedDimension, int& affectedSlice ) { assert(image); assert(plane); // compare normal of plane to the three axis vectors of the image Vector3D normal = plane->GetNormal(); Vector3D imageNormal0 = image->GetSlicedGeometry()->GetAxisVector(0); Vector3D imageNormal1 = image->GetSlicedGeometry()->GetAxisVector(1); Vector3D imageNormal2 = image->GetSlicedGeometry()->GetAxisVector(2); normal.Normalize(); imageNormal0.Normalize(); imageNormal1.Normalize(); imageNormal2.Normalize(); imageNormal0.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal0.Get_vnl_vector()) ); imageNormal1.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal1.Get_vnl_vector()) ); imageNormal2.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal2.Get_vnl_vector()) ); double eps( 0.00001 ); // transversal if ( imageNormal2.GetNorm() <= eps ) { affectedDimension = 2; } // sagittal else if ( imageNormal1.GetNorm() <= eps ) { affectedDimension = 1; } // frontal else if ( imageNormal0.GetNorm() <= eps ) { affectedDimension = 0; } else { affectedDimension = -1; // no idea return false; } // determine slice number in image Geometry3D* imageGeometry = image->GetGeometry(0); Point3D testPoint = imageGeometry->GetCenter(); Point3D projectedPoint; plane->Project( testPoint, projectedPoint ); Point3D indexPoint; imageGeometry->WorldToIndex( projectedPoint, indexPoint ); affectedSlice = ROUND( indexPoint[affectedDimension] ); MITK_DEBUG << "indexPoint " << indexPoint << " affectedDimension " << affectedDimension << " affectedSlice " << affectedSlice; // check if this index is still within the image if ( affectedSlice < 0 || affectedSlice >= static_cast(image->GetDimension(affectedDimension)) ) return false; return true; } mitk::Image::Pointer mitk::SegTool2D::GetAffectedImageSliceAs2DImage(const PositionEvent* positionEvent, const Image* image) { if (!positionEvent) return NULL; assert( positionEvent->GetSender() ); // sure, right? unsigned int timeStep = positionEvent->GetSender()->GetTimeStep( image ); // get the timestep of the visible part (time-wise) of the image // first, we determine, which slice is affected const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); if ( !image || !planeGeometry ) return NULL; - vtkSmartPointer mapReslice = vtkSmartPointer::New(); - mitk::ExtractSliceFilter::Pointer extractor = mitk::ExtractSliceFilter::New(mapReslice); + vtkSmartPointer reslice = vtkSmartPointer::New(); + reslice->SetOverwriteMode(false); + reslice->Modified(); + mitk::ExtractSliceFilter::Pointer extractor = mitk::ExtractSliceFilter::New(reslice); extractor->SetInput( image ); extractor->SetTimeStep( timeStep ); extractor->SetWorldGeometry( planeGeometry ); extractor->SetVtkOutputRequest(false); extractor->SetResliceTransformByGeometry( image->GetTimeSlicedGeometry()->GetGeometry3D( timeStep ) ); - mapReslice->Modified(); + extractor->Modified(); extractor->Update(); - Image::Pointer slice = extractor->GetOutput(); - m_Map = mapReslice->GetMap(); + Image::Pointer slice = extractor->GetOutput(); Vector3D axis0 = slice->GetGeometry()->GetAxisVector(0); Vector3D axis1 = slice->GetGeometry()->GetAxisVector(1); axis0.Normalize(); axis1.Normalize(); Point3D origin = slice->GetGeometry()->GetOrigin(); int offsetX = extractor->GetVtkOutput()->GetExtent()[0]; int offsetY = extractor->GetVtkOutput()->GetExtent()[2]; + + Vector3D spacing = slice->GetGeometry()->GetSpacing(); - origin += (axis0 * offsetX) + (axis1 * offsetY); + origin += (axis0 * (offsetX * spacing[0])) + (axis1 * (offsetY * spacing[1])); slice->GetGeometry()->SetOrigin(origin); return slice; } mitk::Image::Pointer mitk::SegTool2D::GetAffectedWorkingSlice(const PositionEvent* positionEvent) { DataNode* workingNode( m_ToolManager->GetWorkingData(0) ); if ( !workingNode ) return NULL; Image* workingImage = dynamic_cast(workingNode->GetData()); if ( !workingImage ) return NULL; return GetAffectedImageSliceAs2DImage( positionEvent, workingImage ); } mitk::Image::Pointer mitk::SegTool2D::GetAffectedReferenceSlice(const PositionEvent* positionEvent) { DataNode* referenceNode( m_ToolManager->GetReferenceData(0) ); if ( !referenceNode ) return NULL; Image* referenceImage = dynamic_cast(referenceNode->GetData()); if ( !referenceImage ) return NULL; return GetAffectedImageSliceAs2DImage( positionEvent, referenceImage ); } void mitk::SegTool2D::WriteBackSegmentationResult (const PositionEvent* positionEvent, Image* slice) { const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); DataNode* workingNode( m_ToolManager->GetWorkingData(0) ); Image* image = dynamic_cast(workingNode->GetData()); - - mitk::OverwriteSliceFilter::Pointer overwriter = mitk::OverwriteSliceFilter::New(); - overwriter->SetInput(image); - overwriter->SetInputSlice(slice); - overwriter->SetInputMap(m_Map); - overwriter->Modified(); - overwriter->Update(); + vtkSmartPointer reslice = vtkSmartPointer::New(); + reslice->SetInputSlice(slice->GetVtkImageData(this->m_TimeStep)); + reslice->SetOverwriteMode(true); + reslice->Modified(); + + mitk::ExtractSliceFilter::Pointer extractor = mitk::ExtractSliceFilter::New(reslice); + extractor->SetInput( image ); + extractor->SetTimeStep( this->m_TimeStep ); + extractor->SetWorldGeometry( planeGeometry ); + extractor->SetVtkOutputRequest(true); + extractor->SetResliceTransformByGeometry( image->GetTimeSlicedGeometry()->GetGeometry3D( this->m_TimeStep ) ); + extractor->Modified(); + extractor->Update(); + image->Modified(); if ( m_3DInterpolationEnabled ) { slice->DisconnectPipeline(); unsigned int pos = this->AddContourmarker(positionEvent); ImageToContourFilter::Pointer contourExtractor = ImageToContourFilter::New(); contourExtractor->SetInput(slice); contourExtractor->Update(); mitk::Surface::Pointer contour = contourExtractor->GetOutput(); mitk::ServiceReference serviceRef = mitk::GetModuleContext()->GetServiceReference(); PlanePositionManagerService* service = dynamic_cast(mitk::GetModuleContext()->GetService(serviceRef)); mitk::SurfaceInterpolationController::GetInstance()->AddNewContour( contour, service->GetPlanePosition(pos)); contour->DisconnectPipeline(); } } void mitk::SegTool2D::SetShowMarkerNodes(bool status) { m_ShowMarkerNodes = status; } void mitk::SegTool2D::Enable3DInterpolation(bool status) { m_3DInterpolationEnabled = status; } unsigned int mitk::SegTool2D::AddContourmarker ( const PositionEvent* positionEvent ) { const mitk::Geometry2D* plane = dynamic_cast (dynamic_cast< const mitk::SlicedGeometry3D*>( positionEvent->GetSender()->GetSliceNavigationController()->GetCurrentGeometry3D())->GetGeometry2D(0)); mitk::ServiceReference serviceRef = mitk::GetModuleContext()->GetServiceReference(); PlanePositionManagerService* service = dynamic_cast(mitk::GetModuleContext()->GetService(serviceRef)); unsigned int size = service->GetNumberOfPlanePositions(); unsigned int id = service->AddNewPlanePosition(plane, positionEvent->GetSender()->GetSliceNavigationController()->GetSlice()->GetPos()); mitk::PlanarCircle::Pointer contourMarker = mitk::PlanarCircle::New(); contourMarker->SetGeometry2D( const_cast(plane)); std::stringstream markerStream; mitk::DataNode* workingNode (m_ToolManager->GetWorkingData(0)); markerStream << m_Contourmarkername ; markerStream << " "; markerStream << id+1; DataNode::Pointer rotatedContourNode = DataNode::New(); rotatedContourNode->SetData(contourMarker); rotatedContourNode->SetProperty( "name", StringProperty::New(markerStream.str()) ); rotatedContourNode->SetProperty( "isContourMarker", BoolProperty::New(true)); rotatedContourNode->SetBoolProperty( "PlanarFigureInitializedWindow", true, positionEvent->GetSender() ); rotatedContourNode->SetProperty( "includeInBoundingBox", BoolProperty::New(false)); rotatedContourNode->SetProperty( "helper object", mitk::BoolProperty::New(!m_ShowMarkerNodes)); if (plane) { if ( id == size ) { m_ToolManager->GetDataStorage()->Add(rotatedContourNode, workingNode); } else { mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isContourMarker", mitk::BoolProperty::New(true)); mitk::DataStorage::SetOfObjects::ConstPointer markers = m_ToolManager->GetDataStorage()->GetDerivations(workingNode,isMarker); for ( mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) { std::string nodeName = (*iter)->GetName(); unsigned int t = nodeName.find_last_of(" "); unsigned int markerId = atof(nodeName.substr(t+1).c_str())-1; if(id == markerId) { return id; } } m_ToolManager->GetDataStorage()->Add(rotatedContourNode, workingNode); } } return id; } void mitk::SegTool2D::InteractiveSegmentationBugMessage( const std::string& message ) { MITK_ERROR << "********************************************************************************" << std::endl << " " << message << std::endl << "********************************************************************************" << std::endl << " " << std::endl << " If your image is rotated or the 2D views don't really contain the patient image, try to press the button next to the image selection. " << std::endl << " " << std::endl << " Please file a BUG REPORT: " << std::endl << " http://bugs.mitk.org" << std::endl << " Contain the following information:" << std::endl << " - What image were you working on?" << std::endl << " - Which region of the image?" << std::endl << " - Which tool did you use?" << std::endl << " - What did you do?" << std::endl << " - What happened (not)? What did you expect?" << std::endl; } diff --git a/Modules/MitkExt/Interactions/mitkSegTool2D.h b/Modules/MitkExt/Interactions/mitkSegTool2D.h index 6d702c6b8a..75c59fd53b 100644 --- a/Modules/MitkExt/Interactions/mitkSegTool2D.h +++ b/Modules/MitkExt/Interactions/mitkSegTool2D.h @@ -1,144 +1,143 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ 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. =========================================================================*/ #ifndef mitkSegTool2D_h_Included #define mitkSegTool2D_h_Included #include "mitkCommon.h" #include "MitkExtExports.h" #include "mitkTool.h" #include "mitkImage.h" #include "mitkStateEvent.h" #include "mitkPositionEvent.h" #include "mitkPlanePositionManager.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkInteractionConst.h" #include #include -#include -#include + namespace mitk { class BaseRenderer; /** \brief Abstract base class for segmentation tools. \sa Tool \ingroup Interaction \ingroup ToolManagerEtAl Implements 2D segmentation specific helper methods, that might be of use to all kind of 2D segmentation tools. At the moment these are: - Determination of the slice where the user paints upon (DetermineAffectedImageSlice) - Projection of a 3D contour onto a 2D plane/slice SegTool2D tries to structure the interaction a bit. If you pass "PressMoveRelease" as the interaction type of your derived tool, you might implement the methods OnMousePressed, OnMouseMoved, and OnMouseReleased. Yes, your guess about when they are called is correct. \warning Only to be instantiated by mitk::ToolManager. $Author$ */ class MitkExt_EXPORT SegTool2D : public Tool { public: mitkClassMacro(SegTool2D, Tool); /** \brief Calculates for a given Image and PlaneGeometry, which slice of the image (in index corrdinates) is meant by the plane. \return false, if no slice direction seems right (e.g. rotated planes) \param affectedDimension The image dimension, which is constant for all points in the plane, e.g. Transversal --> 2 \param affectedSlice The index of the image slice */ static bool DetermineAffectedImageSlice( const Image* image, const PlaneGeometry* plane, int& affectedDimension, int& affectedSlice ); void SetShowMarkerNodes(bool); void Enable3DInterpolation(bool); protected: SegTool2D(); // purposely hidden SegTool2D(const char*); // purposely hidden virtual ~SegTool2D(); virtual bool OnMousePressed (Action*, const StateEvent*); virtual bool OnMouseMoved (Action*, const StateEvent*); virtual bool OnMouseReleased(Action*, const StateEvent*); virtual bool OnInvertLogic (Action*, const StateEvent*); /** \brief Extract the slice of an image that the user just scribbles on. \return NULL if SegTool2D is either unable to determine which slice was affected, or if there was some problem getting the image data at that position. */ Image::Pointer GetAffectedImageSliceAs2DImage(const PositionEvent*, const Image* image); /** \brief Extract the slice of the currently selected working image that the user just scribbles on. \return NULL if SegTool2D is either unable to determine which slice was affected, or if there was some problem getting the image data at that position, or just no working image is selected. */ Image::Pointer GetAffectedWorkingSlice(const PositionEvent*); /** \brief Extract the slice of the currently selected reference image that the user just scribbles on. \return NULL if SegTool2D is either unable to determine which slice was affected, or if there was some problem getting the image data at that position, or just no reference image is selected. */ Image::Pointer GetAffectedReferenceSlice(const PositionEvent*); void WriteBackSegmentationResult (const PositionEvent*, Image*); /** \brief Adds a new node called Contourmarker to the datastorage which holds a mitk::PlanarFigure. By selecting this node the slicestack will be reoriented according to the PlanarFigure's Geometry */ unsigned int AddContourmarker ( const PositionEvent* ); void InteractiveSegmentationBugMessage( const std::string& message ); unsigned int* m_Map; private: BaseRenderer* m_LastEventSender; unsigned int m_LastEventSlice; //The prefix of the contourmarkername. Suffix is a consecutive number const std::string m_Contourmarkername; bool m_ShowMarkerNodes; bool m_3DInterpolationEnabled; }; } // namespace #endif