Generate an 3-dimensional multigaussian image. More...

#include <itkMultiGaussianImageSource.h>

Inheritance diagram for itk::MultiGaussianImageSource< TOutputImage >:

Collaboration diagram for itk::MultiGaussianImageSource< TOutputImage >:

Public Types
typedef MultiGaussianImageSource	Self

typedef ImageSource< TOutputImage >	Superclass

typedef SmartPointer< Self >	Pointer

typedef SmartPointer< const Self >	ConstPointer

typedef TOutputImage::PixelType	OutputImagePixelType

typedef TOutputImage::RegionType	OutputImageRegionType

typedef TOutputImage::SizeType	SizeType

typedef TOutputImage::IndexType	IndexType

typedef TOutputImage::SpacingType	SpacingType

typedef TOutputImage::PointType	PointType

typedef SizeType::SizeValueType	SizeValueType

typedef SizeValueType	SizeValueArrayType[TOutputImage::ImageDimension]

typedef TOutputImage::SpacingValueType	SpacingValueType

typedef SpacingValueType	SpacingValueArrayType[TOutputImage::ImageDimension]

typedef TOutputImage::PointValueType	PointValueType

typedef PointValueType	PointValueArrayType[TOutputImage::ImageDimension]

typedef itk::ImageRegion< 3 >::SizeValueType	SizeRegionType

typedef double	RadiusType

typedef std::vector< double >	VectorType

typedef Vector< double, TOutputImage::ImageDimension >	ItkVectorType

typedef ImageRegionIteratorWithIndex< TOutputImage >	IteratorType

typedef TOutputImage::Pointer	ImageType

typedef MapContainer< unsigned int, PointType >	MapContainerPoints

typedef MapContainer< unsigned int, double >	MapContainerRadius

Public Member Functions
virtual void	SetSize (SizeType _arg)

virtual void	SetSize (SizeValueArrayType sizeArray)

virtual const SizeValueType *	GetSize () const

virtual void	SetSpacing (SpacingType _arg)

virtual void	SetSpacing (SpacingValueArrayType spacingArray)

virtual const SpacingValueType *	GetSpacing () const

virtual void	SetOrigin (PointType _arg)

virtual void	SetOrigin (PointValueArrayType originArray)

virtual const PointValueType *	GetOrigin () const

virtual unsigned int	GetNumberOfGaussians () const

virtual void	SetNumberOfGausssians (unsigned int)

virtual RadiusType	GetRadius () const

virtual void	SetRadius (RadiusType radius)

virtual const OutputImagePixelType	GetMaxMeanValue () const

virtual const IndexType	GetSphereMidpoint () const

virtual double	MultiGaussianFunctionValueAtPoint (double, double, double)

virtual void	AddGaussian (VectorType centerX, VectorType centerY, VectorType centerZ, VectorType sigmaX, VectorType sigmaY, VectorType sigmaZ, VectorType altitude)

virtual void	CalculateTheMidpointAndTheMeanValueWithOctree ()

virtual void	CalculateMaxAndMinInSphere ()

virtual const IndexType	GetMaxValueIndexInSphere () const

virtual const OutputImagePixelType	GetMaxValueInSphere () const

virtual const IndexType	GetMinValueIndexInSphere () const

virtual const OutputImagePixelType	GetMinValueInSphere () const

virtual void	SetRegionOfInterest (ItkVectorType, ItkVectorType)

virtual void	WriteXMLToTestTheCuboidInsideTheSphere ()

virtual void	CalculateEdgesInSphere (PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, int level)

virtual double	MultiGaussianFunctionValueAtCuboid (double xMin, double xMax, double yMin, double yMax, double zMin, double zMax)

virtual void	InsertPoints (PointType globalCoordinateMidpointCuboid, double cuboidRadius)

virtual void	GenerateCuboidSegmentationInSphere (PointType globalCoordinateMidpointSphere)

virtual double	FunctionPhi (double value)

virtual unsigned int	IntersectTheSphere (PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double sideLength)

void	SetNormalDistributionValues ()

virtual void	SetMin (OutputImagePixelType _arg)

bool	IsInImage (IndexType index)

virtual OutputImagePixelType	GetMin () const

virtual void	SetMax (OutputImagePixelType _arg)

virtual OutputImagePixelType	GetMax () const

Static Public Member Functions
static Pointer	New ()

Protected Member Functions
	MultiGaussianImageSource ()

	~MultiGaussianImageSource () override

void	PrintSelf (std::ostream &os, Indent indent) const override

void	GenerateData () override

void	GenerateOutputInformation () override

Detailed Description

template<typename TOutputImage>
class itk::MultiGaussianImageSource< TOutputImage >

Generate an 3-dimensional multigaussian image.

This class defines an 3-dimensional Image, in which the value at one voxel equals the value of a multigaussian function evaluated at the voxel's coordinates. The multigaussian function is built as a sum of N gaussian function and is defined by the following parameters (Generation of a multigauss image):

CenterX, CenterY, CenterZ - vectors of the size of N determining the expectancy value at the x-, y- and the z-axis. That means: The i-th gaussian bell curve takes its maximal value at the voxel with index [CenterX(i); CenterY(i); Centerz(i)].
SigmaX, SigmaY, SigmaZ - vectors of the size of N determining the deviation at the x-, y- and the z-axis. That means: The width of the i-th gaussian bell curve is determined by the deviation in the x-axis, which is SigmaX(i), in the y-axis is SigmaY(i) and in the z-axis is SigmaZ(i).
Altitude - vector of the size of N determining the altitude: the i-th gaussian bell curve has a height of Altitude(i). This class allows by the method CalculateMidpointAndMeanValue() to find a sphere with a specified radius that has a maximal mean value over all sphere with that radius with midpoint inside or on the boundary of the image. Furthermore it can calculate the maximal und minimal pixel intensities and whose indices in the founded sphere.

To serve as a test tool for ImageStatisticsCalculator, esp. the "hotspot search" feature of this class, MultiGaussianImageSource is also able to calculate the position of a sphere that maximizes the mean value of the voxels within the sphere (Algorithm for calculating statistic in a sphere).

Generation of a multigauss image

A multigauss function consists of the sum of \( N \) gauss function. The \( i \)-th
( \(0 \leq i \leq N \)) gaussian is described with the following seven parameters (see above):

\( x_0^{(i)} \) is the expectancy value in the \( x \)-Axis
\( y_0^{(i)} \) is the expectancy value in the \( y \)-Axis
\( z_0^{(i)} \) is the expectancy value in the \( z \)-Axis
\( \sigma_x^{(i)} \) is the deviation in the \( x \)-Axis
\( \sigma_y^{(i)} \) is the deviation in the \( y \)-Axis
\( \sigma_z^{(i)} \) is the deviation in the \( z \)-Axis
\( a^{(i)} \) is the altitude of the gaussian.

A gauss function has the following form:

\begin{eqnarray} \nonumber f^{(i)}(x,y,z) = a^{(i)} exp \left[ - \left( \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} + \frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} + \frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right) \right]. \end{eqnarray}

A multigauss function has then the form:

\begin{align*} f_N(x,y,z) =& \sum_{i=0}^{N}f^{(i)}(x,y,z)\\ =&\sum_{0}^{N} a^{(i)} exp \left[ - \left( \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} + \frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} + \frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right) \right]. \end{align*}

The multigauss function \(f_N\) will be evaluated at each voxel coordinate to become the voxel intensity.

Algorithm for calculating statistic in a sphere

This section explains how we can find a sphere region which has a maximal mean value over all sphere regions with a fixed radius. Furthermore we want to calculate the maximal and minimal value in the wanted sphere.

To calculate the mean value in a sphere we integrate the gaussians over the whole sphere. The antiderivative is unknown as an explicit function, but we have a value table for the distribution function of the normal distribution \( \Phi(x) \) for \( x \) between \( -3.99 \) and \( 3.99 \) with step size \( 0.01 \). The only problem is that we cannot integrate over a spherical region, because we have an 3-dim integral and therefore are the integral limits dependent from each other and we cannot evaluate \( \Phi \). So we approximate the sphere with cuboids inside the sphere and prisms on the boundary of the sphere. We calculate these cuboids with the octree recursive method: We start by subdividing the wrapping box of the sphere in eight cuboids. Further we subdivide each cuboid in eight cuboids and check for each of them, whether it is inside or outside the sphere or whether it intersects the sphere surface. We save those of them, which are inside the sphere and continue to subdivide the cuboids that intersect the sphere until the recursion breaks. In the last step we take the half of the cuboids on the boundary and this are the prisms. Now we can calculate and sum the integrals over the cuboids and divide through the volume of the body to obtain the mean value.

For each cuboid \( Q = [a_1, b_1]\times[a_2, b_2]\times[a_3, b_3] \) we apply Fubini's theorem for integrable functions and become for the integral the following:

\begin{align*} m_k =& \sum_{i=0}^{N} \int_{Q} f^{(i)}(x,y,z)dx\\ =&\sum_{i=0}^{N} a^{(i)} \int_{Q} exp \left[ - \left( \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} + \frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} + \frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right) \right] dx \\ =& \sum_{i=0}^{N} a^{(i)} \int_{a_1}^{b_1} exp \left[ - \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} \right] dx \int_{a_2}^{b_2}exp \left[ -\frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} \right]dx \int_{a_3}^{b_3}exp \left[ -\frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right]dx. \end{align*}

So we calculate three one dimensional integrals:

\begin{align*} \int_{a}^{b} & exp \left[ - \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} \right] dx \\ =&\int_{-\infty}^{b} exp \left[ - \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} \right] dx - \int_{-\infty}^{a} exp \left[ - \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} \right] dx \\ =& \sigma_x^{(i)} \left[\int_{-\infty}^{(a - x_0^{(i)})/ \sigma_x^{(i)}} e^{-\frac{t^2}{2}} dt - \int_{-\infty}^{(b - x_0^{(i)})/ \sigma_x^{(i)}} e^{-\frac{t^2}{2}}dt \right] \\ =&\sigma_x^{(i)} \sqrt{(\pi)} \left[ \Phi \left( \frac{(a - x_0^{(i)})^2}{\sigma_x^{(i)}} \right) - \Phi \left ( \frac{(b - x_0^{(i)})^2}{\sigma_x^{(i)}} \right) \right]. \end{align*}

and become for the integral over \( Q \):

\begin{align*} m_k =& \sum_{i=0}^{N} \sigma_x^{(i)} \sigma_y^{(i)} \sigma_z^{(i)} \pi^{1.5} \left[ \Phi \left( \frac{(a_1 - x_0^{(i)})^2}{\sigma_x^{(i)}} \right) - \Phi \left ( \frac{(b_1 - x_0^{(i)})^2}{\sigma_x^{(i)}} \right) \right]\times \\ &\left[ \Phi \left( \frac{(a_2 - y_0^{(i)})^2}{\sigma_y^{(i)}} \right) - \Phi \left ( \frac{(b_2 - y_0^{(i)})^2}{\sigma_y^{(i)}} \right) \right]\times \left[ \Phi \left( \frac{(a_3 - z_0^{(i)})^2}{\sigma_z^{(i)}} \right) - \Phi \left ( \frac{(b_3 - z_0^{(i)})^2}{\sigma_z^{(i)}} \right) \right]. \end{align*}

For the integral over the prism we take the half of the integral over the corresponding cuboid.

Altogether we find the mean value in the sphere as:

\begin{align*} \left( \sum_{Q_k \text{ Cuboid}} m_k + \sum_{P_l \text{ Prism}} 0.5 m_l \right )/Volume(B), \end{align*}

where Volume(B) is the volume of the body that approximate the sphere.

Now we know how to calculate the mean value in a sphere for given midpoint and radius. So we assume each voxel in the given image to be the sphere midpoint and we calculate the mean value as described above. If the new mean value is greater than the "maximum-until-now", we take the new value to be the "maximum-until-now". Then we go to the next voxel and make the same calculation and so on. At the same time we save the coordinates of the midpoint voxel.

After we found the midpoint and the maximal mean value, we can calculate the maximum and the minimum in the sphere: we just traverse all the voxels in the region and take the maximum and minimum value and the respective coordinates.

Input and output

An example for an input in the command-line is:

  mitkMultiGaussianTest C:/temp/outputFile C:/temp/inputFile.xml

Here is outputFile the name of the gaussian image with extension .nrrd and at the same time the name of the output file with extension .xml, which is the same as the inputFile, only added the calculated mean value, max and min and the corresponding indexes in the statistic tag. Here we see an example for the input and output .xml file:

INPUT:

<testcase>
  <testimage image-rows="20" image-columns="20" image-slices="20" numberOfGaussians="2" spacingX="1" spacingY="1" spacingZ="1" entireHotSpotInImage="1">
    <gaussian centerIndexX="4" centerIndexY="16" centerIndexZ="10" deviationX="7" deviationY="7" deviationZ="7" altitude="200"/>
    <gaussian centerIndexX="18" centerIndexY="2" centerIndexZ="10" deviationX="1" deviationY="1" deviationZ="1" altitude="210"/>
  </testimage>
  <segmentation numberOfLabels="1" hotspotRadiusInMM="6.2035">
    <roi label="1" maximumSizeX="20" minimumSizeX="0" maximumSizeY="20" minimumSizeY="0" maximumSizeZ="20" minimumSizeZ="0"/>
  </segmentation>
</testcase>

OUTPUT:

<testcase>
  <testimage image-rows="20" image-columns="20" image-slices="20" numberOfGaussians="2" spacingX="1" spacingY="1" spacingZ="1" entireHotSpotInImage="1">
    <gaussian centerIndexX="4" centerIndexY="16" centerIndexZ="10" deviationX="7" deviationY="7" deviationZ="7" altitude="200"/>
    <gaussian centerIndexX="18" centerIndexY="2" centerIndexZ="10" deviationX="1" deviationY="1" deviationZ="1" altitude="210"/>
  </testimage>
  <segmentation numberOfLabels="1" hotspotRadiusInMM="6.2035">
    <roi label="1" maximumSizeX="20" minimumSizeX="0" maximumSizeY="20" minimumSizeY="0" maximumSizeZ="20" minimumSizeZ="0"/>
  </segmentation>
  <statistic hotspotIndexX="6" hotspotIndexY="13" hotspotIndexZ="10" peak="141.544" mean="141.544" maximumIndexX="4" maximumIndexY="16" maximumIndexZ="10" maximum="200" minimumIndexX="9" minimumIndexY="8" minimumIndexZ="8"  minimum="77.4272"/>
</testcase>

Parameter for the input

In the tag testimage we describe the image that we generate. Image rows/columns/slices gives the number of rows/columns/slices of the image; numberOfGaussians is the number of gauss functions ( \( N \)); spacing defines the extend of one voxel for each direction. The parameter entireHotSpotInImage determines whether the whole sphere is in the image included ( \( = 1 \)) or only the midpoint of the sphere is inside the image.

NOTE: When the entireHotSpotInImage \( = 0 \) it is possible that we find the midpoint of the sphere on the border of the image. In this case we cut the approximation of the sphere, so that we become a body, which is completely inside the image, but not a "sphere" anymore. To that effect is the volume of the body decreased and that could lead to unexpected results.

In the subtag gaussian we describe each gauss function as mentioned in the second section.

In the tag segmentation we define the radius of the wanted sphere in mm (hotspotRadiusInMM ). We can also set the number of labels (numberOfLabels ) to be an positive number and this determines the number of regions of interest (ROI). In each ROI we find the sphere with the wanted properties and midpoint inside the ROI, but not necessarily the whole sphere. In the subtag roi we set label number and the index coordinates for the borders of the roi.

Definition at line 194 of file itkMultiGaussianImageSource.h.