5 #include <itkImageRegionIterator.h>
7 #include <itkImageDuplicator.h>
8 #include <itkFFTConvolutionImageFilter.h>
14 m_HotspotRadiusinMM(6.2035049089940),
15 m_HotspotMustBeCompletelyInsideImage(true),
20 m_InternalMaskUpdateTime = 0;
28 m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension());
29 m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension());
49 if(radiusInMillimeter != m_HotspotRadiusinMM)
51 m_HotspotRadiusinMM = radiusInMillimeter;
58 return m_HotspotRadiusinMM;
63 return m_HotspotMustBeCompletelyInsideImage;
68 if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage)
70 m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage;
78 if (IsUpdateRequired())
82 throw std::runtime_error(
"Error: image empty!" );
87 throw std::runtime_error(
"Error: invalid time step!" );
93 imageTimeSelector->UpdateLargestPossibleRegion();
96 m_internalImage = timeSliceImage;
97 m_internalMask2D =
nullptr;
98 m_internalMask3D =
nullptr;
100 if ( m_Mask !=
nullptr )
105 if ( m_internalImage->GetDimension() == 3 )
110 else if ( m_internalImage->GetDimension() == 2 )
117 throw std::runtime_error(
"Error: invalid image dimension" );
123 if ( m_internalImage->GetDimension() == 3 )
127 else if ( m_internalImage->GetDimension() == 2 )
133 throw std::runtime_error(
"Error: invalid image dimension" );
153 if (label != m_Label)
163 return m_ConvolutionImageMinIndex;
169 return m_ConvolutionImageMaxIndex;
172 template <
typename TPixel,
unsigned int VImageDimension >
174 HotspotMaskGenerator::CalculateExtremaWorld(
const itk::Image<TPixel, VImageDimension>* inputImage,
176 double neccessaryDistanceToImageBorderInMM,
179 typedef itk::Image< TPixel, VImageDimension >
ImageType;
180 typedef itk::Image< unsigned short, VImageDimension >
MaskImageType;
182 typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MaskImageIteratorType;
183 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> InputImageIndexIteratorType;
185 typename ImageType::SpacingType spacing = inputImage->GetSpacing();
188 minMax.Defined =
false;
189 minMax.MaxIndex.set_size(VImageDimension);
190 minMax.MaxIndex.set_size(VImageDimension);
192 typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion();
194 bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 );
195 if (keepDistanceToImageBorders)
197 long distanceInPixels[VImageDimension];
198 for(
unsigned short dimension = 0; dimension < VImageDimension; ++dimension)
204 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5);
207 allowedExtremaRegion.ShrinkByRadius(distanceInPixels);
210 InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion);
215 typename ImageType::IndexType maxIndex;
216 typename ImageType::IndexType minIndex;
218 for(
unsigned short i = 0; i < VImageDimension; ++i)
224 if (maskImage !=
nullptr)
226 MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());
227 typename ImageType::IndexType imageIndex;
229 typename ImageType::IndexType maskIndex;
231 for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
233 imageIndex = maskIndex = maskIt.GetIndex();
235 if(maskIt.Get() == label)
237 if( allowedExtremaRegion.IsInside(imageIndex) )
239 imageIndexIt.SetIndex( imageIndex );
240 double value = imageIndexIt.Get();
241 minMax.Defined =
true;
244 if( value > maxValue )
246 maxIndex = imageIndexIt.GetIndex();
250 if(value < minValue )
252 minIndex = imageIndexIt.GetIndex();
261 for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt)
263 double value = imageIndexIt.Get();
264 minMax.Defined =
true;
267 if( value > maxValue )
269 maxIndex = imageIndexIt.GetIndex();
273 if(value < minValue )
275 minIndex = imageIndexIt.GetIndex();
281 minMax.MaxIndex.set_size(VImageDimension);
282 minMax.MinIndex.set_size(VImageDimension);
284 for(
unsigned int i = 0; i < minMax.MaxIndex.size(); ++i)
286 minMax.MaxIndex[i] = maxIndex[i];
289 for(
unsigned int i = 0; i < minMax.MinIndex.size(); ++i)
291 minMax.MinIndex[i] = minIndex[i];
294 minMax.Max = maxValue;
295 minMax.Min = minValue;
300 template <
unsigned int VImageDimension>
301 itk::Size<VImageDimension>
302 HotspotMaskGenerator::CalculateConvolutionKernelSize(
double spacing[VImageDimension],
305 typedef itk::Image< float, VImageDimension > KernelImageType;
306 typedef typename KernelImageType::SizeType SizeType;
309 for(
unsigned int i = 0; i < VImageDimension; ++i)
311 maskSize[i] =
static_cast<int>( 2 * radiusInMM / spacing[i]);
314 if(maskSize[i] % 2 == 0 )
322 template <
unsigned int VImageDimension>
324 HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(
double mmPerPixel[VImageDimension],
327 std::stringstream ss;
328 for (
unsigned int i = 0; i < VImageDimension; ++i)
331 if (i < VImageDimension -1)
334 MITK_DEBUG <<
"Update convolution kernel for spacing (" << ss.str() <<
") and radius " << radiusInMM <<
"mm";
337 double radiusInMMSquared = radiusInMM * radiusInMM;
338 typedef itk::Image< float, VImageDimension > KernelImageType;
342 typedef typename KernelImageType::SizeType SizeType;
343 SizeType maskSize = this->CalculateConvolutionKernelSize<VImageDimension>(mmPerPixel, radiusInMM);
346 convolutionMaskCenterIndex.Fill(0.0);
347 for(
unsigned int i = 0; i < VImageDimension; ++i)
349 convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1);
352 typedef typename KernelImageType::IndexType IndexType;
356 typedef typename KernelImageType::RegionType RegionType;
357 RegionType maskRegion;
358 maskRegion.SetSize(maskSize);
359 maskRegion.SetIndex(maskIndex);
361 convolutionKernel->SetRegions(maskRegion);
362 convolutionKernel->SetSpacing(mmPerPixel);
363 convolutionKernel->Allocate();
366 typedef itk::ImageRegionIteratorWithIndex<KernelImageType> MaskIteratorType;
367 MaskIteratorType maskIt(convolutionKernel,maskRegion);
369 int numberOfSubVoxelsPerDimension = 2;
370 int numberOfSubVoxels = ::pow( static_cast<float>(numberOfSubVoxelsPerDimension), static_cast<float>(VImageDimension) );
371 double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension;
372 double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels;
373 double maskValue = 0.0;
375 double distanceSquared = 0.0;
377 typedef itk::ContinuousIndex<double, VImageDimension> ContinuousIndexType;
378 for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
380 ContinuousIndexType indexPoint(maskIt.GetIndex());
382 for (
unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
384 voxelPosition[dimension] = indexPoint[dimension];
390 for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0;
391 subVoxelOffset[0] < +0.5;
392 subVoxelOffset[0] += subVoxelSizeInPixels)
394 for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0;
395 subVoxelOffset[1] < +0.5;
396 subVoxelOffset[1] += subVoxelSizeInPixels)
398 for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0;
399 subVoxelOffset[2] < +0.5;
400 subVoxelOffset[2] += subVoxelSizeInPixels)
402 subVoxelIndexPosition = voxelPosition + subVoxelOffset;
404 (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0]
405 + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1]
406 + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2];
408 if (distanceSquared <= radiusInMMSquared)
410 maskValue += valueOfOneSubVoxel;
415 maskIt.Set( maskValue );
418 return convolutionKernel;
421 template <
typename TPixel,
unsigned int VImageDimension>
423 HotspotMaskGenerator::GenerateConvolutionImage(
const itk::Image<TPixel, VImageDimension>* inputImage )
425 double mmPerPixel[VImageDimension];
426 for (
unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
428 mmPerPixel[dimension] = inputImage->GetSpacing()[dimension];
432 typedef itk::Image< float, VImageDimension > KernelImageType;
433 typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel<VImageDimension>(mmPerPixel, m_HotspotRadiusinMM);
437 typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
440 ConvolutionImageType> ConvolutionFilterType;
443 typedef itk::ConstantBoundaryCondition<InputImageType, InputImageType> BoundaryConditionType;
444 BoundaryConditionType boundaryCondition;
445 boundaryCondition.SetConstant(0.0);
447 if (m_HotspotMustBeCompletelyInsideImage)
450 convolutionFilter->SetBoundaryCondition(&boundaryCondition);
453 convolutionFilter->SetInput(inputImage);
454 convolutionFilter->SetKernelImage(convolutionKernel);
455 convolutionFilter->SetNormalize(
true);
456 MITK_DEBUG <<
"Update Convolution image for hotspot search";
457 convolutionFilter->UpdateLargestPossibleRegion();
460 convolutionImage->SetSpacing( inputImage->GetSpacing() );
462 return convolutionImage;
465 template <
typename TPixel,
unsigned int VImageDimension>
467 HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image<TPixel, VImageDimension>* maskImage,
468 itk::Point<double, VImageDimension> sphereCenter,
469 double sphereRadiusInMM )
472 typedef itk::ImageRegionIteratorWithIndex<MaskImageType> MaskImageIteratorType;
474 MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());
476 typename MaskImageType::IndexType maskIndex;
480 for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
482 maskIndex = maskIt.GetIndex();
483 maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition);
484 maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 );
488 template <
typename TPixel,
unsigned int VImageDimension>
490 HotspotMaskGenerator::CalculateHotspotMask(itk::Image<TPixel, VImageDimension>* inputImage,
495 typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
496 typedef itk::Image< float, VImageDimension > KernelImageType;
497 typedef itk::Image< unsigned short, VImageDimension >
MaskImageType;
501 if (convolutionImage.IsNull())
503 MITK_ERROR <<
"Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error).";
504 throw std::logic_error(
"Empty convolution image in CalculateHotspotStatistics()");
509 if (maskImage ==
nullptr)
512 typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion();
513 typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing();
515 typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection();
516 maskImage->SetRegions(maskRegion);
517 maskImage->Allocate();
518 maskImage->SetOrigin(maskOrigin);
519 maskImage->SetSpacing(maskSpacing);
520 maskImage->SetDirection(maskDirection);
522 maskImage->FillBuffer(1);
528 double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0;
529 ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label);
531 bool isHotspotDefined = convolutionImageInformation.
Defined;
533 if (!isHotspotDefined)
535 MITK_ERROR <<
"No origin of hotspot-sphere was calculated!";
549 hotspotMaskITK->SetOrigin(inputImage->GetOrigin());
550 hotspotMaskITK->SetSpacing(inputImage->GetSpacing());
551 hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion());
552 hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion());
553 hotspotMaskITK->SetDirection(inputImage->GetDirection());
554 hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel());
555 hotspotMaskITK->Allocate();
556 hotspotMaskITK->FillBuffer(1);
558 typedef typename InputImageType::IndexType IndexType;
559 IndexType maskCenterIndex;
560 for (
unsigned int d =0; d< VImageDimension;++d)
562 maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d];
566 inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter);
568 FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM);
574 m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex;
575 m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex;
579 bool HotspotMaskGenerator::IsUpdateRequired()
const
581 unsigned long thisClassTimeStamp = this->GetMTime();
583 unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime();
584 unsigned long inputImageTimeStamp =
m_inputImage->GetMTime();
586 if (thisClassTimeStamp > m_InternalMaskUpdateTime)
591 if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp)
596 if (internalMaskTimeStamp > m_InternalMaskUpdateTime)
mitk::Image::Pointer m_inputImage
vnl_vector< int > GetHotspotIndex()
Returns the image index where the hotspot is located.
void SetHotspotRadiusInMM(double radiusInMillimeter)
Set the radius of the hotspot (in MM)
bool GetHotspotMustBeCompletelyInsideImage() const
itk::SmartPointer< Self > Pointer
const double & GetHotspotRadiusinMM() const
DataCollection - Class to facilitate loading/accessing structured data.
void SetLabel(unsigned short label)
If a maskGenerator is set, this detemines which mask value is used.
void SetInputImage(mitk::Image::Pointer inputImage)
Set the input image. Required for this class.
Image::Pointer GrabItkImageMemory(itk::SmartPointer< ItkOutputImageType > &itkimage, mitk::Image *mitkImage=nullptr, const BaseGeometry *geometry=nullptr, bool update=true)
Grabs the memory of an itk::Image (with a specific type) and puts it into an mitk::Image.The memory is managed by the mitk::Image after calling this function. The itk::Image remains valid until the mitk::Image decides to free the memory.
mitk::Image::Pointer GetMask()
Computes and returns the hotspot mask. The hotspot mask has the same size as the input image...
void SetMask(MaskGenerator::Pointer mask)
Set a mask (can be nullptr if no mask is desired)
map::core::discrete::Elements< 3 >::InternalImageType ImageType
itk::Image< double, 3 > InputImageType
#define AccessFixedDimensionByItk_2(mitkImage, itkImageTypeFunction, dimension, arg1, arg2)
itk::Image< unsigned char, 3 > MaskImageType
void SetTimeStep(unsigned int timeStep)
SetTimeStep is used to set the time step for which the mask is to be generated.
mitk::Image::Pointer m_InternalMask
vnl_vector< int > GetConvolutionImageMinIndex()
Returns the index where the convolution image is minimal (darkest spot in image)
void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage)
Define whether the hotspot must be completely inside the image. Default is true.
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.