diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp index cbbe73198c..242690fbac 100644 --- a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -1,547 +1,546 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "itkMultiGaussianImageSource.h" #include #include #include #include #include #include #include #include // Commandline:(for exmaple) mitkMultiGaussianTest C:/temp/output C:/temp/inputFile.xml // //For Example: inputFile.xml // // // // // // // // // // //For Example: output.xml // // // // // // // // // // bool IsInOtherROI(int, itk::MultiGaussianImageSource>::VectorType, itk::MultiGaussianImageSource>::VectorType, itk::MultiGaussianImageSource>::VectorType, itk::MultiGaussianImageSource>::VectorType, itk::MultiGaussianImageSource>::VectorType, itk::MultiGaussianImageSource>::VectorType); int mitkMultiGaussianTest(int argc, char* argv[]) { // Read the parmaeter from a .xml File. // In the inputFile.xml we find the characteristics of the Gaussian and the ROI's. Hier we can have more then one ROI -> we find the hot spot for each of the ROI's; we can set the entire HotSpot to be in the image or just its midpoint, but not necessary the whole HotSpot. const unsigned int Dimension = 3; typedef double PixelType; typedef itk::DOMNode::Pointer DOMNodeType; typedef itk::Image ImageType; typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; std::string outputFilename = argv[1], name; int numberOfImages; double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude, hotSpotRadiusInMM; unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, numberOfLabels; itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec, ROImaxSizeX, ROIminSizeX, ROImaxSizeY, ROIminSizeY, ROImaxSizeZ, ROIminSizeZ, label; itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; MultiGaussianImageSource::Pointer gaussianGenerator; itk::DOMNodeXMLWriter::Pointer xmlWriter; itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; DOMNodeType domTestCase, domTestImage, domGaussian, domSegmentation, domStatistics, domROI; ImageType::SizeValueType size[3]; std::stringstream ss; double radius = pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10; char * fileNamePointer; std::string attributeValue; double value; bool entireHotSpotInImage; int maxSize, minSize; std::string filename = argv[2]; itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); xmlReader->SetFileName( filename ); xmlReader->Update(); itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); typedef std::vector NodeList; // read test image parameters, fill result structure NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; attributeValue = testimage->GetAttribute("image-rows"); std::stringstream(attributeValue) >> size[0]; attributeValue = testimage->GetAttribute("image-columns"); std::stringstream(attributeValue) >> size[1]; attributeValue = testimage->GetAttribute("image-slices"); std::stringstream(attributeValue) >> size[2]; attributeValue = testimage->GetAttribute( "numberOfGaussians" ); std::stringstream(attributeValue) >> numberOfGaussians; attributeValue = testimage->GetAttribute( "spacingX" ); std::stringstream(attributeValue) >> spacing[0]; attributeValue = testimage->GetAttribute( "spacingY" ); std::stringstream(attributeValue) >> spacing[1]; attributeValue = testimage->GetAttribute( "spacingZ" ); std::stringstream(attributeValue) >> spacing[2]; attributeValue = testimage->GetAttribute( "entireHotSpotInImage" ); std::stringstream(attributeValue) >> entireHotSpotInImage; std::cout << "Read size parameters (x,y,z): " << size[0] << ", " << size[1] << ", " << size[2] << "\n" << std::endl; std::cout << "Read spacing parameters (x,y,z): " << spacing[0] << ", " << spacing[1] << ", " << spacing[2]<< "\n" << std::endl; NodeList gaussians; testimage->GetChildren("gaussian", gaussians); itk::DOMNode* gaussian; for(int i = 0; i < numberOfGaussians ; ++i) { gaussian = gaussians[i]; //TODO attributeValue = gaussian->GetAttribute( "centerIndexX" ); std::stringstream(attributeValue) >> value; centerXVec.push_back(value * spacing[0]); attributeValue = gaussian->GetAttribute( "centerIndexY" ); std::stringstream(attributeValue) >> value; centerYVec.push_back(value * spacing[1]); attributeValue = gaussian->GetAttribute( "centerIndexZ" ); std::stringstream(attributeValue) >> value; centerZVec.push_back(value * spacing[2]); std::cout << "Read center of Gaussian (x,y,z) in mm: " << centerXVec[i] << ", " << centerYVec[i] << ", " << centerZVec[i] << "\n" << std::endl; attributeValue = gaussian->GetAttribute( "deviationX" ); std::stringstream(attributeValue) >> value; sigmaXVec.push_back(value * spacing[0]); attributeValue = gaussian->GetAttribute( "deviationY" ); std::stringstream(attributeValue) >> value; sigmaYVec.push_back(value * spacing[1]); attributeValue = gaussian->GetAttribute( "deviationZ" ); std::stringstream(attributeValue) >> value; sigmaZVec.push_back(value * spacing[2]); std::cout << "Read deviation of Gaussian (x,y,z) in mm: " << sigmaXVec[i] << ", " << sigmaYVec[i] << ", " << sigmaZVec[i] << "\n" << std::endl; attributeValue = gaussian->GetAttribute( "altitude" ); std::stringstream(attributeValue) >> value; altitudeVec.push_back(value); std::cout << "Read altitude: " << altitudeVec[i] << "\n" << std::endl; } // read ROI's parameter NodeList segmentations; domRoot->GetChildren("segmentation", segmentations); MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One ROI image defined" ) itk::DOMNode* segmentation = segmentations[0]; attributeValue = segmentation->GetAttribute("numberOfLabels"); std::stringstream(attributeValue) >> numberOfLabels; attributeValue = segmentation->GetAttribute("hotspotRadiusInMM"); std::stringstream(attributeValue) >> hotSpotRadiusInMM; std::cout << "Read number of labels: " << numberOfLabels << std::endl; std::cout << "Read radius in mm : " << hotSpotRadiusInMM << std::endl; NodeList rois; segmentation->GetChildren("roi", rois); itk::DOMNode* roi; // for each label i take the ROI and set it to be the i'th element of ROImaxSize* and ROIminSize* ( * = X, Y, Z) for(int i = 0; i < numberOfLabels ; ++i) { roi = rois[i]; attributeValue = roi->GetAttribute( "label" ); std::stringstream(attributeValue) >> value; label.push_back(value); attributeValue = roi->GetAttribute( "maximumSizeX" ); std::stringstream(attributeValue) >> value; ROImaxSizeX.push_back(value); attributeValue = roi->GetAttribute( "minimumSizeX" ); std::stringstream(attributeValue) >> value; ROIminSizeX.push_back(value); attributeValue = roi->GetAttribute( "maximumSizeY" ); std::stringstream(attributeValue) >> value; ROImaxSizeY.push_back(value); attributeValue = roi->GetAttribute( "minimumSizeY" ); std::stringstream(attributeValue) >> value; ROIminSizeY.push_back(value); attributeValue = roi->GetAttribute( "maximumSizeZ" ); std::stringstream(attributeValue) >> value; ROImaxSizeZ.push_back(value); attributeValue = roi->GetAttribute( "minimumSizeZ" ); std::stringstream(attributeValue) >> value; ROIminSizeZ.push_back(value); std::cout << "Read ROI with label number: " << label[i] << " with min and max values in the x-, y-, z-Achse: [" << ROIminSizeX[i] << " " << ROImaxSizeX[i] <<"], [" << ROIminSizeY[i] << " " << ROImaxSizeY[i] <<"], [" << ROIminSizeZ[i] << " " << ROImaxSizeZ[i] <<"]\n" << std::endl; } // Check whether the ROI's are correct defined, i.e. whether the ROI's are disjoint for(int i = 1; i < numberOfLabels ; ++i) { // check whether the edges of the i'th ROI is in another ROI included (when yes -> ERROR) bool isInOtherROI = IsInOtherROI( i, ROIminSizeX, ROImaxSizeX, ROIminSizeY, ROImaxSizeY, ROIminSizeZ, ROImaxSizeZ ); if( isInOtherROI) { std::cout << "The ROI's in the different labels should be disjoint! Please define it correct. " << std::endl; return 0; } } //write test image parameter xmlWriter = itk::DOMNodeXMLWriter::New(); domTestCase = itk::DOMNode::New(); domTestCase->SetName("testcase"); domTestImage = itk::DOMNode::New(); domTestImage->SetName("testimage"); ss.str(""); ss << size[0]; domTestImage->SetAttribute("image-rows", ss.str()); ss.str(""); ss << size[1]; domTestImage->SetAttribute("image-columns", ss.str()); ss.str(""); ss << size[2]; domTestImage->SetAttribute("image-slices", ss.str()); ss.str(""); ss << numberOfGaussians; domTestImage->SetAttribute("numberOfGaussians", ss.str()); ss.str(""); ss << spacing[0]; domTestImage->SetAttribute("spacingX", ss.str()); ss.str(""); ss << spacing[1]; domTestImage->SetAttribute("spacingY", ss.str()); ss.str(""); ss << spacing[2]; domTestImage->SetAttribute("spacingZ", ss.str()); ss.str(""); ss << entireHotSpotInImage; domTestImage->SetAttribute("entireHotSpotInImage", ss.str()); domTestCase->AddChildAtBegin(domTestImage); for( unsigned int i = 0; i < numberOfGaussians; ++i) { domGaussian = itk::DOMNode::New() ; domGaussian->SetName("gaussian"); domTestImage->AddChildAtEnd(domGaussian); // write the midpoint and the daviation in pixel units centerX = centerXVec[i] / spacing[0]; ss.str(""); ss << centerX; // * spacing[0]; //static_cast( static_cast( centerX / spacing[0] + 0.9999 ) ); domGaussian->SetAttribute("centerIndexX", ss.str()); centerY = centerYVec[i] / spacing[1]; ss.str(""); ss << centerY; // * spacing[1]; //static_cast( static_cast( centerY / spacing[1] + 0.9999 ) ); domGaussian->SetAttribute("centerIndexY", ss.str()); centerZ = centerZVec[i] / spacing[2]; ss.str(""); ss << centerZ; // * spacing[2]; //static_cast( static_cast( centerZ / spacing[2] + 0.9999 ) ); domGaussian->SetAttribute("centerIndexZ", ss.str()); sigmaX = sigmaXVec[i] / spacing[0]; ss.str(""); ss << sigmaX; // * spacing[0]; // static_cast( static_cast( sigmaX / spacing[0] + 0.9999 ) ); domGaussian->SetAttribute("deviationX", ss.str()); sigmaY = sigmaYVec[i] / spacing[1]; ss.str(""); ss << sigmaY; // * spacing[1]; //static_cast( static_cast( sigmaY / spacing[1] + 0.9999 ) ); domGaussian->SetAttribute("deviationY", ss.str()); sigmaZ = sigmaZVec[i] / spacing[2]; ss.str(""); ss << sigmaZ; // * spacing[2]; //static_cast( static_cast( sigmaZ / spacing[2] + 0.9999 ) ); domGaussian->SetAttribute("deviationZ", ss.str()); altitude = altitudeVec[i]; ss.str(""); ss << altitude; domGaussian->SetAttribute("altitude", ss.str()); } radius = hotSpotRadiusInMM; // set the parameter for the gaussianGenerator gaussianGenerator = MultiGaussianImageSource::New(); gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); - gaussianGenerator->SetRadiusStepNumber(8); gaussianGenerator->SetRadius(radius); gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); domSegmentation = itk::DOMNode::New(); domSegmentation->SetName("segmentation"); domTestCase->AddChildAtEnd(domSegmentation); ss.str(""); ss << numberOfLabels; domSegmentation->SetAttribute("numberOfLabels", ss.str()); ss.str(""); ss << hotSpotRadiusInMM; domSegmentation->SetAttribute("hotspotRadiusInMM", ss.str()); // set the region of interest for each label i for (unsigned int i = 0; i < numberOfLabels; ++i) { // Set region of interest in index values. The entire HotSpot is in the image. if(entireHotSpotInImage) { // x axis region of interest------------------------------------------------------ maxSize = size[0] - static_cast((radius)/spacing[0] + 0.9999) + 1; minSize = 0.0 + static_cast((radius)/spacing[0]+ 0.9999) - 1; if( minSize >= maxSize ) { std::cout << "The sphere is larger then the image in the x axis!" << std::endl; } // the maximum in the x-Axis regionOfInterestMax.SetElement( 0, ( ROImaxSizeX[i] < maxSize ) ? ROImaxSizeX[i] : maxSize ); // the minimum in the x-Axis regionOfInterestMin.SetElement( 0, ( ROIminSizeX[i] > minSize ) ? ROIminSizeX[i] : minSize ); // y axis region of interest------------------------------------------------------ maxSize = size[1] - static_cast((radius)/spacing[1] + 0.9999) + 1; minSize = 0.0 + static_cast((radius)/spacing[1]+ 0.9999) - 1; if( minSize >= maxSize ) { std::cout << "The sphere is larger then the image in the y axis!" << std::endl; } // the maximum in the y-Axis regionOfInterestMax.SetElement( 1, ( ROImaxSizeY[i] < maxSize ) ? ROImaxSizeY[i] : maxSize ); // the minimum in the y-Axis regionOfInterestMin.SetElement( 1, ( ROIminSizeY[i] > minSize ) ? ROIminSizeY[i] : minSize ); // z axis region of interest------------------------------------------------------ maxSize = size[2] - static_cast((radius)/spacing[2] + 0.9999) + 1; minSize = 0.0 + static_cast((radius)/spacing[2]+ 0.9999) - 1; if( minSize >= maxSize ) { std::cout << "The sphere is larger then the image in the z axis!" << std::endl; } // the maximum in the z-Axis regionOfInterestMax.SetElement( 2, ( ROImaxSizeZ[i] < maxSize ) ? ROImaxSizeZ[i] : maxSize ); // the minimum in the z-Axis regionOfInterestMin.SetElement( 2, ( ROIminSizeZ[i] > minSize ) ? ROIminSizeZ[i] : minSize ); } // Set region of interest in index values. The midpoint of the HotSpot is in the image, but not necessary the whole HotSpot else { // x axis region of interest------------------------------------------------------ regionOfInterestMax.SetElement( 0, ROImaxSizeX[i] ); regionOfInterestMin.SetElement( 0, ROIminSizeX[i] ); // y axis region of interest------------------------------------------------------ regionOfInterestMax.SetElement( 1, ROImaxSizeY[i] ); regionOfInterestMin.SetElement( 1, ROIminSizeY[i] ); // z axis region of interest------------------------------------------------------ regionOfInterestMax.SetElement( 2, ROImaxSizeZ[i] ); regionOfInterestMin.SetElement( 2, ROIminSizeZ[i] ); } gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); gaussianGenerator->Update(); //write region of interest for the .xml file domROI = itk::DOMNode::New(); domROI->SetName("roi"); domSegmentation->AddChildAtEnd(domROI); ss.str(""); ss << label[i]; domROI->SetAttribute("label", ss.str()); ss.str(""); ss << ROImaxSizeX[i]; domROI->SetAttribute("maximumSizeX", ss.str()); ss.str(""); ss << ROIminSizeX[i]; domROI->SetAttribute("minimumSizeX", ss.str()); ss.str(""); ss << ROImaxSizeY[i]; domROI->SetAttribute("maximumSizeY", ss.str()); ss.str(""); ss << ROIminSizeY[i]; domROI->SetAttribute("minimumSizeY", ss.str()); ss.str(""); ss << ROImaxSizeZ[i]; domROI->SetAttribute("maximumSizeZ", ss.str()); ss.str(""); ss << ROIminSizeZ[i]; domROI->SetAttribute("minimumSizeZ", ss.str()); // Calculate the mean value and the midpoint of the wanted sphere. gaussianGenerator->CalculateTheMidpointAndTheMeanValueWithOctree(); //peak and peak coordinate domStatistics = itk::DOMNode::New(); domStatistics->SetName("statistic"); domTestCase->AddChildAtEnd(domStatistics); sphereMidpt = gaussianGenerator->GetSphereMidpoint(); ss.str(""); ss << sphereMidpt[0]; domStatistics->SetAttribute("hotspotIndexX", ss.str()); ss.str(""); ss << sphereMidpt[1]; domStatistics->SetAttribute("hotspotIndexY", ss.str()); ss.str(""); ss << sphereMidpt[2]; domStatistics->SetAttribute("hotspotIndexZ", ss.str()); ss.str(""); ss << gaussianGenerator->GetMaxMeanValue(); domStatistics->SetAttribute("mean", ss.str()); //maximum and maximum coordinate gaussianGenerator->CalculateMaxAndMinInSphere(); maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); ss.str(""); ss << maxValueIndexInSphere[0]; domStatistics->SetAttribute("maximumIndexX", ss.str()); ss.str(""); ss << maxValueIndexInSphere[1]; domStatistics->SetAttribute("maximumIndexY", ss.str()); ss.str(""); ss << maxValueIndexInSphere[2]; domStatistics->SetAttribute("maximumIndexZ", ss.str()); ss.str(""); ss << gaussianGenerator->GetMaxValueInSphere(); domStatistics->SetAttribute("maximum", ss.str()); //minimum and minimum coordinate minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); ss.str(""); ss << minValueIndexInSphere[0]; domStatistics->SetAttribute("minimumIndexX", ss.str()); ss.str(""); ss << minValueIndexInSphere[1]; domStatistics->SetAttribute("minimumIndexY", ss.str()); ss.str(""); ss << minValueIndexInSphere[2]; domStatistics->SetAttribute("minimumIndexZ", ss.str()); ss.str(""); ss << gaussianGenerator->GetMinValueInSphere(); domStatistics->SetAttribute("minimum", ss.str()); } // .xml (Data) ss.str(""); ss << outputFilename << ".xml"; name = ss.str(); fileNamePointer = (char*) name.c_str(); xmlWriter->SetFileName( fileNamePointer); xmlWriter->SetInput( domTestCase ); xmlWriter->Update(); ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); //.nrrd (Image) typedef itk::ImageFileWriter< ImageType > WriterType; WriterType::Pointer writer = WriterType::New(); ss.str(""); ss << outputFilename << ".nrrd"; name = ss.str(); fileNamePointer = (char*) name.c_str(); writer->SetFileName( fileNamePointer); writer->SetInput( gaussianImage ); writer->Update(); } // check whether the edges of the i'th ROI is in another ROI included bool IsInOtherROI(int i, itk::MultiGaussianImageSource>::VectorType ROIminSizeX, itk::MultiGaussianImageSource>::VectorType ROImaxSizeX, itk::MultiGaussianImageSource>::VectorType ROIminSizeY, itk::MultiGaussianImageSource>::VectorType ROImaxSizeY, itk::MultiGaussianImageSource>::VectorType ROIminSizeZ, itk::MultiGaussianImageSource>::VectorType ROImaxSizeZ ) { bool error = 0; std::vector xBound, yBound, zBound; xBound.push_back( ROIminSizeX[i] ); xBound.push_back( ROImaxSizeX[i] ); yBound.push_back( ROIminSizeY[i] ); yBound.push_back( ROImaxSizeY[i] ); zBound.push_back( ROIminSizeZ[i] ); zBound.push_back( ROImaxSizeZ[i] ); //for each ROI for( unsigned int j = 0; j < i; ++j ) { for( unsigned int x = 0; x < 2; ++x) { for( unsigned int y = 0; y < 2; ++y) { for( unsigned int z = 0; z < 2; ++z) { double edgeXCoord = xBound[x]; double edgeYCoord = yBound[y]; double edgeZCoord = zBound[z]; // check if the edge with coordinate [edgeXCoord; edgeYCoord; edgeZCoord] is inside the j'th ROI if ( ROIminSizeX[j] < edgeXCoord && edgeXCoord < ROImaxSizeX[j] && ROIminSizeY[j] < edgeYCoord && edgeYCoord < ROImaxSizeY[j] && ROIminSizeZ[j] < edgeZCoord && edgeZCoord < ROImaxSizeZ[j]) { error = 1; return error; } } } } } return error; } \ No newline at end of file diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/itkMultiGaussianImageSource.h index 8ce8ff63f7..ba4d250dd2 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.h +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.h @@ -1,394 +1,375 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_h #define __itkMultiGaussianImageSource_h #include "itkImageSource.h" #include "itkNumericTraits.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkImageFileWriter.h" #include namespace itk { /** \class MultiGaussianImageSource \brief 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 build as a sum of N gaussian function. This is defined by the following parameters: +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 build as a sum of N gaussian function. This is defined by the following parameters: -1. 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)]. +1. 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)]. -2. 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-the gaussian bell curve deviation in the x-axis is SigmaX(i), in the y-axis is SigmaY(i) and in the z-axis is SigmaZ(i). +2. 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). 3. 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 at the boundary of the image. Furthermore it can calculate the maximal und minimal pixel intensities and whose indices in the founded sphere. +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. \section Multigaussian-image-source-generator Multigaussian image source generator \subsection Abstract Abstract The class name is MultiGaussianImageSource and can be found in the directory MINT\Mint_Bin\CMakeExternals\Source\MITK\Modules\ImageStatistics. This class generates an image, which pixel intensities are the values of a multigauss function. The class allows calculating a spherical region of defined size which has a maximal mean value and returns the position and the mean value. Furthermore it can calculate the maximal und minimal pixel intensities and whose indices in the founded sphere. \subsection Generation-of-a-multigauss-image Generation of a multigauss image A multigauss function consists of the sum of \f$ N \f$ gauss function. The \f$ i \f$-th \linebreak (\f$0 \leq i \leq N \f$) gaussian is described with the following seven parameters: - \f$ x_0^{(i)} \f$ is the expectancy value in the \f$ x \f$-Axis - \f$ y_0^{(i)} \f$ is the expectancy value in the \f$ y \f$-Axis - \f$ z_0^{(i)} \f$ is the expectancy value in the \f$ z \f$-Axis - \f$ \sigma_x^{(i)} \f$ is the deviation in the \f$ x \f$-Axis - \f$ \sigma_y^{(i)} \f$ is the deviation in the \f$ y \f$-Axis - \f$ \sigma_z^{(i)} \f$ is the deviation in the \f$ z \f$-Axis - \f$ a^{(i)} \f$ is the altitude of the gaussian. A gauss function has the following form: \f{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]. \f} A multigauss function has then the form: \f{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]. \f} The multigauss function \f$f_N\f$ will be evaluated at each voxel coordinate to become the voxel intensity. \subsection Algorithm-for-calculating-statistic-in-a-sphere 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 \f$ \Phi(x) \f$ for \f$ x \f$ between \f$ -3.99 \f$ and \f$ 3.99 \f$ with step size \f$ 0.01 \f$. 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 \f$ \Phi \f$. 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 \f$ Q = [a_1, b_1]\times[a_2, b_2]\times[a_3, b_3] \f$ we apply Fubini's theorem for integrable functions and become for the integral the following: \f{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. \f} So we calculate three one dimensional integrals: \f{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]. \f} and become for the integral over \f$ Q \f$: \f{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]. \f} For the integral over the prism we take the half of the integral over the corresponding cuboid. Altogether we become for the mean value in the sphere: \f{align*}{ \left( \sum_{Q_k \text{ Cuboid}} m_k + \sum_{P_l \text{ Prism}} 0.5 m_l \right )/Volume(B), \f} 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. \subsection Input-and-output 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: +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 --------------------------------------------------- +\verbatim - - - - - - - - +\endverbatim ---------------------- OUTPUT --------------------------------------------------- +--------------------- OUTPUT -------------------------------------------------- +\verbatim - - - - - - - - - +\endverbatim \subsection Parameter-for-the-input Parameter for the input In the tag \a testimage we describe the image that we generate. Image rows/columns/slices gives the number of rows/columns/slices of the image; \a numberOfGaussians is the number of gauss functions (\f$ N \f$); spacing defines the extend of one voxel for each direction. The parameter \a entireHotSpotInImage determines whether the whole sphere is in the image included (\f$ = 1 \f$) or only the midpoint of the sphere is inside the image. NOTE: When the \a entireHotSpotInImage \f$ = 0 \f$ 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 \a gaussian we describe each gauss function as mentioned in the second section. -In the tag \a segmentation we define the radius of the wanted sphere in mm (\a hotspotRadiusInMM ). We can also set the number of labels (\a 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 \a roi we set label number and the index coordinates for the borders of the roi. +In the tag \a segmentation we define the radius of the wanted sphere in mm (\a hotspotRadiusInMM ). We can also set the number of labels (\a 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 \a roi we set label number and the index coordinates for the borders of the roi. */ template< typename TOutputImage > class ITK_EXPORT MultiGaussianImageSource:public ImageSource< TOutputImage > { public: /** Standard class typedefs. */ typedef MultiGaussianImageSource Self; typedef ImageSource< TOutputImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; /** Typedef for the output image PixelType. */ typedef typename TOutputImage::PixelType OutputImagePixelType; /** Typedef to describe the output image region type. */ typedef typename TOutputImage::RegionType OutputImageRegionType; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Basic types from the OutputImageType */ typedef typename TOutputImage::SizeType SizeType; typedef typename TOutputImage::IndexType IndexType; typedef typename TOutputImage::SpacingType SpacingType; typedef typename TOutputImage::PointType PointType; typedef typename SizeType::SizeValueType SizeValueType; typedef SizeValueType SizeValueArrayType[TOutputImage::ImageDimension]; typedef typename TOutputImage::SpacingValueType SpacingValueType; typedef SpacingValueType SpacingValueArrayType[TOutputImage::ImageDimension]; typedef typename TOutputImage::PointValueType PointValueType; typedef PointValueType PointValueArrayType[TOutputImage::ImageDimension]; typedef typename itk::ImageRegion<3> ::SizeValueType SizeRegionType; /** Typedef to describe the sphere radius type. */ typedef double RadiusType; /** Typedef to describe the standard vector type. */ typedef std::vector VectorType; /** Typedef to describe the itk vector type. */ typedef Vector ItkVectorType; /** Typedef to describe the ImageRegionIteratorWithIndex type. */ typedef ImageRegionIteratorWithIndex IteratorType; /** Typedef to describe the Poiner type at the output image. */ typedef typename TOutputImage::Pointer ImageType; typedef MapContainer MapContainerPoints; typedef MapContainer MapContainerRadius; /** Set/Get size of the output image. */ itkSetMacro(Size, SizeType); virtual void SetSize(SizeValueArrayType sizeArray); virtual const SizeValueType * GetSize() const; /** Set/Get spacing of the output image. */ itkSetMacro(Spacing, SpacingType); virtual void SetSpacing(SpacingValueArrayType spacingArray); virtual const SpacingValueType * GetSpacing() const; - /** Set/Get origin of the output image. */ + /** Set/Get origin of the output image. This programm works proper only with origin [0.0, 0.0, 0.0] */ itkSetMacro(Origin, PointType); virtual void SetOrigin(PointValueArrayType originArray); virtual const PointValueType * GetOrigin() const; /** Get the number of gaussian functions in the output image. */ virtual unsigned int GetNumberOfGaussians() const; /** Set the number of gaussian function. */ virtual void SetNumberOfGausssians( unsigned int ); /** Set/Get the radius of the sphere. */ virtual const RadiusType GetRadius() const; virtual void SetRadius( RadiusType radius ); - /** Set/Get the number of steps to traverse the radius of the sphere. */ - virtual const int GetRadiusStepNumber() const; - /** Set the number of steps to traverse the radius. */ - virtual void SetRadiusStepNumber( unsigned int stepNumber ); /** Get the maximal mean value in a sphere over all possible spheres with midpoint in the image. */ virtual const OutputImagePixelType GetMaxMeanValue() const; /** Get the index of the midpoint of a sphere with the maximal mean value.*/ virtual const IndexType GetSphereMidpoint() const; /** Calculates the value of the multigaussian function at a Point given by its coordinates [x, y, z]. */ virtual const double MultiGaussianFunctionValueAtPoint(double , double, double); /** Adds a multigaussian defined by the parameter: CenterX, CenterY, CenterZ, SigmaX, SigmaY, SigmaZ, Altitude. All parameters should have the same size, which determinates the number of the gaussian added. */ virtual void AddGaussian( VectorType, VectorType, VectorType, VectorType, VectorType, VectorType, VectorType); /** Calculates and set the index of the midpoint of the sphere with the maximal mean value as well as the mean value. */ virtual void CalculateTheMidpointAndTheMeanValueWithOctree(); /** Calculates and set the index an the value of maximulm and minimum in the wanted sphere. */ virtual void CalculateMaxAndMinInSphere(); /** Get the index in the sphere with maximal value. */ virtual const IndexType GetMaxValueIndexInSphere() const; /** Get the maximal value in the sphere. */ virtual const OutputImagePixelType GetMaxValueInSphere() const; /** Get the index in the sphere with minimal value. */ virtual const IndexType GetMinValueIndexInSphere() const; /** Get the minimal value in the sphere. */ virtual const OutputImagePixelType GetMinValueInSphere() const; /** Set the region of interest. */ virtual void SetRegionOfInterest(ItkVectorType, ItkVectorType); /** Write a .mps file to visualise the point in the sphere. */ virtual void WriteXMLToTestTheCuboidInsideTheSphere(); /**This recursive method realise the octree method. It subdivide a cuboid in eight cuboids, when this cuboid crosses the boundary of sphere. If the cuboid is inside the sphere, it calculates the integral. */ virtual void CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, int level); /**Calculate and return value of the integral of the gaussian in a cuboid region with the dimension 3: in the x-axis between xMin and xMax and in the y-axis between yMin and yMax and in the z-axis also between zMin and zMax. */ virtual double MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax); /** Inseret the midpoints of cuboid in a vector m_Midpoints, so that we can visualise it. */ virtual void InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius); /** Start the octree recursion in eigth directions for the sphere with midpoint globalCoordinateMidpointSphere. */ virtual void GenerateCuboidSegmentationInSphere( PointType globalCoordinateMidpointSphere ); /** Get the the values of the cumulative distribution function of the normal distribution. */ virtual double FunctionPhi(double value); /** Check if a cuboid with midpoint globalCoordinateMidpointCuboid and side length sideLength intersect the sphere with midpoint globalCoordinateMidpointSphere boundary. */ virtual unsigned int IntersectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double sideLength); /** Set the tabel values of the distribution function of the normal distribution. */ void SetNormalDistributionValues(); /** Set the minimum possible pixel value. By default, it is * NumericTraits::min(). */ itkSetClampMacro( Min, OutputImagePixelType, NumericTraits< OutputImagePixelType >::NonpositiveMin(), NumericTraits< OutputImagePixelType >::max() ); /** Check if a index is inside the image*/ bool IsInImage(IndexType index); /** Get the minimum possible pixel value. */ itkGetConstMacro(Min, OutputImagePixelType); /** Set the maximum possible pixel value. By default, it is * NumericTraits::max(). */ itkSetClampMacro( Max, OutputImagePixelType, NumericTraits< OutputImagePixelType >::NonpositiveMin(), NumericTraits< OutputImagePixelType >::max() ); /** Get the maximum possible pixel value. */ itkGetConstMacro(Max, OutputImagePixelType); protected: MultiGaussianImageSource(); ~MultiGaussianImageSource(); void PrintSelf(std::ostream & os, Indent indent) const; virtual void GenerateData(); virtual void GenerateOutputInformation(); private: MultiGaussianImageSource(const MultiGaussianImageSource &); //purposely not implemented void operator=(const MultiGaussianImageSource &); //purposely not implemented SizeType m_Size; //size of the output image SpacingType m_Spacing; //spacing PointType m_Origin; //origin OutputImagePixelType m_MaxValueInSphere; //maximal value in the wanted sphere IndexType m_MaxValueIndexInSphere; //index of the maximal value in the wanted sphere OutputImagePixelType m_MinValueInSphere; //minimal value in the wanted sphere IndexType m_MinValueIndexInSphere; //index of the minimal value in the wanted sphere unsigned int m_NumberOfGaussians; //number of Gaussians RadiusType m_Radius; //radius of the sphere unsigned int m_RadiusStepNumber; //number of steps to traverse the sphere radius OutputImagePixelType m_MeanValue; //mean value in the wanted sphere OutputImagePixelType m_ValueAtMidpoint; //value at the midpoint of the wanted sphere IndexType m_SphereMidpoint; //midpoint of the wanted sphere VectorType m_SigmaX; //deviation in the x-axis VectorType m_SigmaY; //deviation in the y-axis VectorType m_SigmaZ; //deviation in the z-axis VectorType m_CenterX; //x-coordinate of the mean value of Gaussians VectorType m_CenterY; //y-coordinate of the mean value of Gaussians VectorType m_CenterZ; //z-coordinate of the mean value of Gaussians VectorType m_Altitude; //amplitude ItkVectorType m_RegionOfInterestMax; //maximal values for the coordinates in the region of interest ItkVectorType m_RegionOfInterestMin; //minimal values for the coordinates in the region of interest typename TOutputImage::PixelType m_Min; //minimum possible value typename TOutputImage::PixelType m_Max; //maximum possible value PointType m_GlobalCoordinate; //physical coordiante of the sphere midpoint bool m_WriteMPS; //1 = write a MPS File to visualise the cuboid midpoints of one approximation of the sphere MapContainerPoints m_Midpoints; //the midpoints of the cuboids MapContainerRadius m_RadiusCuboid; //the radius ( = 0.5 * side length) of the cuboids (in the same order as the midpoints in m_Midpoints) double m_Volume; //the volume of the body, that approximize the sphere double m_NormalDistValues [410];//normal distribution values double m_meanValueTemp; //= m_Volume * meanValue in each sphere // The following variables are deprecated, and provided here just for // backward compatibility. It use is discouraged. mutable PointValueArrayType m_OriginArray; mutable SpacingValueArrayType m_SpacingArray; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiGaussianImageSource.hxx" #endif #endif diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx index b43cae0c1d..f43bfbcf91 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx @@ -1,1001 +1,972 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_hxx #define __itkMultiGaussianImageSource_hxx #include #include #include #include "itkMultiGaussianImageSource.h" #include "itkImageRegionIterator.h" #include "itkObjectFactory.h" #include "itkProgressReporter.h" #include "stdlib.h" namespace itk { /** * */ template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::MultiGaussianImageSource() { //Initial image is 100 wide in each direction. for ( unsigned int i = 0; i < TOutputImage::GetImageDimension(); i++ ) { m_Size[i] = 100; m_Spacing[i] = 1.0; m_Origin[i] = 0.0; m_SphereMidpoint[i] = 0; } m_NumberOfGaussians = 0; m_Radius = 1; m_RadiusStepNumber = 5; m_MeanValue = 0; m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); m_Max = NumericTraits< OutputImagePixelType >::max(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::~MultiGaussianImageSource() {} //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSize(SizeValueArrayType sizeArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( sizeArray[i] != this->m_Size[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Size[i] = sizeArray[i]; } } } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SizeValueType * MultiGaussianImageSource< TOutputImage > ::GetSize() const { return this->m_Size.GetSize(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSpacing(SpacingValueArrayType spacingArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( spacingArray[i] != this->m_Spacing[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Spacing[i] = spacingArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetOrigin(PointValueArrayType originArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( originArray[i] != this->m_Origin[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Origin[i] = originArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::PointValueType * MultiGaussianImageSource< TOutputImage > ::GetOrigin() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_OriginArray[i] = this->m_Origin[i]; } return this->m_OriginArray; } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SpacingValueType * MultiGaussianImageSource< TOutputImage > ::GetSpacing() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_SpacingArray[i] = this->m_Spacing[i]; } return this->m_SpacingArray; } //----------------------------------------------------------------------------------------------------------------------- /** * */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Max: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Max ) << std::endl; os << indent << "Min: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Min ) << std::endl; os << indent << "Origin: ["; unsigned int ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Origin[ii] << ", "; ++ii; } os << m_Origin[ii] << "]" << std::endl; os << indent << "Spacing: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Spacing[ii] << ", "; ++ii; } os << m_Spacing[ii] << "]" << std::endl; os << indent << "Size: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Size[ii] << ", "; ++ii; } os << m_Size[ii] << "]" << std::endl; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::GetNumberOfGaussians() const { return this->m_NumberOfGaussians; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::RadiusType MultiGaussianImageSource< TOutputImage > ::GetRadius() const { return this->m_Radius; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRadius( RadiusType radius ) { this->m_Radius = radius; } - //----------------------------------------------------------------------------------------------------------------------- - template< class TOutputImage > - const int - MultiGaussianImageSource< TOutputImage > - ::GetRadiusStepNumber() const - { - return this->m_RadiusStepNumber; - } - //----------------------------------------------------------------------------------------------------------------------- - template< class TOutputImage > - void - MultiGaussianImageSource< TOutputImage > - ::SetRadiusStepNumber( unsigned int stepNumber ) - { - this->m_RadiusStepNumber = stepNumber; - } + //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNumberOfGausssians( unsigned int n ) { this->m_NumberOfGaussians = n; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) { m_RegionOfInterestMax.operator=(roiMax); m_RegionOfInterestMin.operator=(roiMin); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxMeanValue() const { return m_MeanValue; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxValueInSphere() const { return m_MaxValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMaxValueIndexInSphere() const { return m_MaxValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMinValueInSphere() const { return m_MinValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMinValueIndexInSphere() const { return m_MinValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetSphereMidpoint() const { return m_SphereMidpoint; } //----------------------------------------------------------------------------------------------------------------------- /* Calculate and return value of the integral of the gaussian in a cuboid region with the dimension 3: in the x-axis between xMin and xMax and in the y-axis between yMin and yMax and in the z-axis also between zMin and zMax. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax) { double mean = 0; double summand0, summand1, summand2, value, factor; for(unsigned int n = 0; n < m_NumberOfGaussians; ++n) { //TODO //summand0 = FunctionPhi((xMax - m_CenterX[n] * m_Spacing[0]) / (m_SigmaX[n] * m_Spacing[0]) ) - FunctionPhi((xMin - m_CenterX[n] * m_Spacing[0]) / (m_SigmaX[n] * m_Spacing[0]) ); //summand1 = FunctionPhi((yMax - m_CenterY[n] * m_Spacing[1]) / (m_SigmaY[n] * m_Spacing[1]) ) - FunctionPhi((yMin - m_CenterY[n] * m_Spacing[1]) / (m_SigmaY[n] * m_Spacing[1]) ); //summand2 = FunctionPhi((zMax - m_CenterZ[n] * m_Spacing[2]) / (m_SigmaZ[n] * m_Spacing[2]) ) - FunctionPhi((zMin - m_CenterZ[n] * m_Spacing[2]) / (m_SigmaZ[n] * m_Spacing[2]) ); summand0 = FunctionPhi((xMax - m_CenterX[n]) / m_SigmaX[n] ) - FunctionPhi((xMin - m_CenterX[n]) / m_SigmaX[n] ); summand1 = FunctionPhi((yMax - m_CenterY[n]) / m_SigmaY[n] ) - FunctionPhi((yMin - m_CenterY[n]) / m_SigmaY[n] ); summand2 = FunctionPhi((zMax - m_CenterZ[n]) / m_SigmaZ[n] ) - FunctionPhi((zMin - m_CenterZ[n]) / m_SigmaZ[n] ); value = summand0 * summand1 * summand2; factor = (m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) * pow(2.0 * itk::Math::pi, 1.5 ); // * / // ( ((xMax - m_CenterX[n]) / m_SigmaX[n] - (xMin - m_CenterX[n]) / m_SigmaX[n]) // * ((yMax - m_CenterY[n]) / m_SigmaY[n] - (yMin - m_CenterY[n]) / m_SigmaY[n]) // * ((zMax - m_CenterZ[n]) / m_SigmaZ[n] - (zMin - m_CenterZ[n]) / m_SigmaZ[n])); // 1 / ((xMax - xMin) * (yMax - yMin) * (zMax - zMin) * pow(2.0 * itk::Math::pi, 1.5 ));// (m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) / ((xMax - xMin) * (yMax - yMin) * (zMax - zMin)); //* pow(2.0 * itk::Math::pi, 1.5 )); mean = mean + factor * value * m_Altitude[n]; // m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin); } return mean; } //--------------------------------------------------------------------------------------------------------------------- /* Returns the linear interpolation of the values of the standard normal distribution function. This values could be seen in the vector m_NormalDistValues. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::FunctionPhi(double value) { double phiAtValue; //linear interpolation between the values int indexValue = static_cast( 100 * value); if( indexValue > 409 ) { return phiAtValue = 1.0; } else if( indexValue < -409 ) { return phiAtValue = 0.0; } else if( indexValue == 409 ) { return phiAtValue = m_NormalDistValues[409]; } else if( indexValue == -409 ) { return phiAtValue = 1.0 - m_NormalDistValues[409]; } else { bool negative = false; if (indexValue < 0.0) { negative = true; value = -value; } int indexUp = static_cast( 100 * value) + 1; int indexDown = static_cast( 100 * value); double alpha = 100.0 * value - static_cast(indexDown); phiAtValue = (1.0 - alpha) * m_NormalDistValues[indexDown] + alpha * m_NormalDistValues[indexUp] ; if(negative) { phiAtValue = 1.0 - phiAtValue; } return phiAtValue; } } //---------------------------------------------------------------------------------------------------------------------- /* Set the midpoint of the cuboid in a vector m_Midpoints. This cuboids discretise the sphere with the octree method. Set the radius of the cuboid ( = 0.5 * length of the side of the cuboid ) in the vector m_RadiusCuboid. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius) { MapContainerPoints::ElementIdentifier id = m_Midpoints.Size(); m_Midpoints.InsertElement(id, globalCoordinateMidpointCuboid); m_RadiusCuboid.InsertElement(id, cuboidRadius); } //---------------------------------------------------------------------------------------------------------------------- /* This recursive method realise the octree method. It subdivide a cuboid in eight cuboids, when this cuboid crosses the boundary of sphere. If the cuboid is inside the sphere, it calculates the integral and, if uncommented, insert the midpoint of the cuboid in the m_Midpoints vector. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, int level ) { double xMin, xMax, yMin, yMax, zMin, zMax; double cuboidRadiusRecursion = cuboidRadius; PointType newMidpoint; int intersect = IntersectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadiusRecursion); if( (intersect == 1) ) { if (level < 4) { // cuboid intersect the sphere -> call CalculateEdgesInSphere (this method) for the subdivided cuboid cuboidRadiusRecursion = cuboidRadiusRecursion / 2.0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadiusRecursion; newMidpoint[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadiusRecursion; newMidpoint[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadiusRecursion; this->CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadiusRecursion, level + 1 ); } } } } // last step of recursion -> on the boundary else { // Calculate the integral and take the half of it (because we are on the boundary) xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; //check if the boundary of the integral is inside the image if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { //double temp = this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0) / 2.0; } } } else if(intersect == 2) { // cuboid in the sphere // To insert the midpoint and the radius of the cuboid in a vector, so that we can visualise the midpoints, uncomment the next line and the line WriteXMLToTestTheCuboidInsideTheSphere(); // InsertPoints(globalCoordinateMidpointCuboid, cuboidRadius); // Calculate the integral boundary xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane bool yzPlaneAtOriginCrossXSection = xMin <= m_Origin[0] && xMax >= m_Origin[0]; bool yzPlaneAtImageBorderCrossXSection = xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]; bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane bool xzPlaneAtOriginCrossYSection = yMin <= m_Origin[1] && yMax >= m_Origin[1]; bool xzPlaneAtImageBorderCrossYSection = yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]; bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane bool xyPlaneAtOriginCrossZSection = zMin <= m_Origin[2] && zMax >= m_Origin[2]; bool xyPlaneAtImageBorderCrossZSection = zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]; bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0); } // check if the boundary of the image intersect the cuboid and if yes, change the limits of the cuboid to be only inside the image; therefor we cut the sphere and neglect the part of it outside the image else if( // one plane crosses the cuboid ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) || // two plane cross the cuboid (on the image edges possible) ( (yzPlaneAtOriginCrossXSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection)) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) ) { // x-Axis if(xMin <= m_Origin[0] && xMax >= m_Origin[0]) { xMin = m_Origin[0]; }else if(xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]) { xMax = m_Size[0] * m_Spacing[0]; } // y-Axis if(yMin <= m_Origin[1] && yMax >= m_Origin[1]) { yMin = m_Origin[1]; }else if(yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]) { yMax = m_Size[1] * m_Spacing[1]; } // z-Axis if(zMin <= m_Origin[2] && zMax >= m_Origin[2]) { zMin = m_Origin[2]; }else if(zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]) { zMax = m_Size[2] * m_Spacing[2]; } m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin) ; } } } //----------------------------------------------------------------------------------------------------------------------- /* Start the octree recursion in eigth directions for the sphere with midpoint globalCoordinateMidpointSphere and, if uncommented, write this in a file, so that we can visualise it. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateCuboidSegmentationInSphere(PointType globalCoordinateMidpointSphere) { double cuboidRadius = m_Radius * 0.5; PointType newMidpoint; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointSphere[0] + static_cast(i) * cuboidRadius; newMidpoint[1] = globalCoordinateMidpointSphere[1] + static_cast(k) * cuboidRadius; newMidpoint[2] = globalCoordinateMidpointSphere[2] + static_cast(j) * cuboidRadius; CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadius, 0); } } } if(m_WriteMPS) { m_WriteMPS = 0; // uncomment to write an .mps file to visualise the midpoints // std::cout << "Wrote .xml to visualise the midpoints." << std::endl; // WriteXMLToTestTheCuboidInsideTheSphere(); } } //---------------------------------------------------------------------------------------------------------------------- /* This class allows by the method CalculateTheMidpointAndTheMeanValueWithOctree() to find a sphere with a specified radius that has a maximal mean value over all sphere with that radius with midpoint inside or at the boundary of the image. We approximaze the sphere with the octree recursiv method. CalculateTheMidpointAndTheMeanValueWithOctree works as follows: 1. The for-loops traverse the region of interest and assume the current point to be the wanted sphere midpoint. 2. Calculate the mean value for that sphere. 3. Compare with the until-now-found-maximum and take the bigger one. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateTheMidpointAndTheMeanValueWithOctree() { m_MeanValue = 0.0; double meanValueTemp; PointType midpoint; MapContainerPoints::ElementIdentifier cuboidNumber = m_Midpoints.Size(); SetNormalDistributionValues(); double radius; //double xMin, xMax, yMin, yMax, zMin, zMax; m_WriteMPS = 1; PointType globalCoordinateMidpointSphere; IndexType index; OutputImageRegionType regionOfInterest; IndexType indexR; indexR.SetElement( 0, m_RegionOfInterestMin[0] ); indexR.SetElement( 1, m_RegionOfInterestMin[1] ); indexR.SetElement( 2, m_RegionOfInterestMin[2] ); regionOfInterest.SetIndex(indexR); SizeType sizeROI; sizeROI.SetElement( 0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] ); sizeROI.SetElement( 1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] ); sizeROI.SetElement( 2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] ); regionOfInterest.SetSize(sizeROI); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType regionOfInterestIterator(image, regionOfInterest); for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) { index = regionOfInterestIterator.GetIndex(); image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); m_Volume = 0.0; m_meanValueTemp = 0.0; this->GenerateCuboidSegmentationInSphere(globalCoordinateMidpointSphere); meanValueTemp = m_meanValueTemp / m_Volume; // std::cout << "index: " << index <<" meanvalue: " << meanValueTemp << std::endl; if(meanValueTemp > m_MeanValue) { m_GlobalCoordinate = globalCoordinateMidpointSphere; m_MeanValue = meanValueTemp; m_SphereMidpoint = index; } } } //---------------------------------------------------------------------------------------------------------------------- /* Check if a cuboid intersect the sphere boundary. Returns 2, if the cuboid is inside the sphere; returns 1, if the cuboid intersects the sphere boundary and 0, if the cuboid is out of the sphere. */ template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::IntersectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius) { unsigned int intersect = 1; PointType cuboidEdge; int count = 0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { cuboidEdge[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadius; cuboidEdge[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadius; cuboidEdge[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadius; count = count + ( globalCoordinateMidpointSphere.SquaredEuclideanDistanceTo(cuboidEdge) <= m_Radius * m_Radius ); } } } if ( count == 0 ) { // cuboid not in the sphere intersect = 0; } if (count == 8 ) { // cuboid in the sphere intersect = 2; } return intersect; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtPoint(double x, double y, double z) { double summand0, summand1, summand2, power, value = 0.0; // the for-loop represent the sum of the gaussian function for(unsigned int n =0; n < m_NumberOfGaussians; ++n) - { //TODO * spacing - //summand0 = ( x - m_CenterX[n] * m_Spacing[0] ) / ( m_SigmaX[n] * m_Spacing[0] ); - //summand1 = ( y - m_CenterY[n] * m_Spacing[1] ) / ( m_SigmaY[n] * m_Spacing[1] ); - //summand2 = ( z - m_CenterZ[n] * m_Spacing[2] ) / ( m_SigmaZ[n] * m_Spacing[2] ); + { summand0 = ( x - m_CenterX[n] ) / m_SigmaX[n]; summand1 = ( y - m_CenterY[n] ) / m_SigmaY[n]; summand2 = ( z - m_CenterZ[n] ) / m_SigmaZ[n]; power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; - value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); //* (1 / (pow(2.0 * itk::Math::pi, 1.5 ) * m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n]) ) ; + value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); } - - /*for(unsigned int n =0; n < m_NumberOfGaussians; ++n) - { - summand0 = (x - m_CenterX[n]) / m_SigmaX[n]; - summand1 = (y - m_CenterY[n]) / m_SigmaY[n]; - summand2 = (z - m_CenterZ[n]) / m_SigmaZ[n]; - - power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; - value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); - }*/ - return value; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::AddGaussian( VectorType x, VectorType y, VectorType z, VectorType sx, VectorType sy, VectorType sz, VectorType altitude) { for(unsigned int i = 0; i < x.size(); ++i) { // NOT ?? the x,y,z are given in index coordinates -> multiply by the spacing aand save the global coordinate m_CenterX.push_back( x[i] ); m_CenterY.push_back( y[i] ); m_CenterZ.push_back( z[i] ); m_SigmaX.push_back( sx[i] ); m_SigmaY.push_back( sy[i] ); m_SigmaZ.push_back( sz[i] ); m_Altitude.push_back(altitude[i]); } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateOutputInformation() { TOutputImage *output; IndexType index; index.Fill(0); output = this->GetOutput(0); typename TOutputImage::RegionType largestPossibleRegion; largestPossibleRegion.SetSize(this->m_Size); largestPossibleRegion.SetIndex(index); output->SetLargestPossibleRegion(largestPossibleRegion); output->SetSpacing(m_Spacing); output->SetOrigin(m_Origin); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateData() { itkDebugMacro(<< "Generating a image of scalars "); double valueReal; IndexType index; typedef typename TOutputImage::PixelType scalarType; typename TOutputImage::Pointer image = this->GetOutput(0); image = this->GetOutput(0); image->SetBufferedRegion( image->GetRequestedRegion() ); image->Allocate(); IteratorType imageIt(image, image->GetLargestPossibleRegion()); PointType globalCoordinate; for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) { valueReal = 0.0; index = imageIt.GetIndex(); image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); valueReal = MultiGaussianFunctionValueAtPoint(globalCoordinate[0], globalCoordinate[1] ,globalCoordinate[2]); imageIt.Set(valueReal); } } //----------------------------------------------------------------------------------------------------------------------- /*This method is used to write a .mps file, so that we can visualize the midpoints of the approximated sphere as a scatterplot (for example with MITK Workbench). */ template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::WriteXMLToTestTheCuboidInsideTheSphere() { std::stringstream ss; int numberSummand = 1.0; //write an .mps test file itk::DOMNodeXMLWriter::Pointer xmlWriter; typedef itk::DOMNode::Pointer DOMNodeType; DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; xmlWriter = itk::DOMNodeXMLWriter::New(); domXML = itk::DOMNode::New(); domXML->SetName("?xml"); domPointSetFile = itk::DOMNode::New(); domPointSetFile->SetName("point_set_file"); //domFileVersion = itk::DOMNode::New(); //domFileVersion->SetName("file_version"); domPointSet = itk::DOMNode::New(); domPointSet->SetName("point_set"); ss.str(""); ss << 1.0; domXML->SetAttribute("version", ss.str()); domXML->AddChildAtBegin(domPointSetFile); //domPointSetFile -> AddChildAtBegin(domFileVersion); domPointSetFile -> AddChildAtBegin(domPointSet); int cap = m_Midpoints.Size(); for(unsigned int iter = 0 ; iter < cap; ++iter) { domPoint = itk::DOMNode::New(); domPoint->SetName("point"); domX = itk::DOMNode::New(); domX->SetName("x"); domY = itk::DOMNode::New(); domY->SetName("y"); domZ = itk::DOMNode::New(); domZ->SetName("z"); domId = itk::DOMNode::New(); domId->SetName("id"); ss.str(""); ss << numberSummand; domId->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domId); double scaleFactor = 10.0; PointType point = m_Midpoints.GetElement( numberSummand - 1 ); ss.str(""); ss << point[0] * scaleFactor; domX->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domX); ss.str(""); ss << point[1] * scaleFactor; domY->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domY); ss.str(""); ss << point[2] * scaleFactor; domZ->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domZ); domPointSet -> AddChildAtEnd(domPoint); numberSummand += 1.0; } // .mps (Data) ss.str(""); ss << "C:/temp/CuboidsInTheSphere.mps"; std::string name = ss.str(); char * fileNamePointer = (char*) name.c_str(); xmlWriter->SetFileName( fileNamePointer); xmlWriter->SetInput( domXML ); xmlWriter->Update(); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateMaxAndMinInSphere() { IndexType index; MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value; m_MaxValueInSphere = std::numeric_limits::min(); m_MinValueInSphere = std::numeric_limits::max(); int radInt; OutputImageRegionType cuboidRegion; IndexType indexR; SizeType sizeR; int indexRegion, sizeRegion; for( unsigned int i = 0; i < 3; ++i ) { radInt = static_cast(m_Radius/m_Spacing[i]); indexRegion = m_SphereMidpoint[i] - radInt; if( m_Origin[i] > indexRegion ) { indexR.SetElement(i, m_Origin[i] ); } else { indexR.SetElement(i, indexRegion ); } sizeRegion = 2 *radInt + 1; if(indexR[i] + sizeRegion > m_Size[i]) { std::cout << "Not the entire sphere is in the image!" << std::endl; sizeR.SetElement(i, m_Size[i] - indexRegion ); } else { sizeR.SetElement(i, sizeRegion ); } } cuboidRegion.SetIndex(indexR); cuboidRegion.SetSize(sizeR); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType cuboidRegionOfInterestIterator(image, cuboidRegion); PointType globalCoordinate; for(cuboidRegionOfInterestIterator.GoToBegin(); !cuboidRegionOfInterestIterator.IsAtEnd(); ++cuboidRegionOfInterestIterator) { index = cuboidRegionOfInterestIterator.GetIndex(); if(IsInImage(index)) { image->TransformIndexToPhysicalPoint(cuboidRegionOfInterestIterator.GetIndex(), globalCoordinate); if( m_GlobalCoordinate.EuclideanDistanceTo(globalCoordinate) <= m_Radius ) { value = cuboidRegionOfInterestIterator.Get(); if(m_MaxValueInSphere < value) { m_MaxValueInSphere = value; m_MaxValueIndexInSphere = index; } if(m_MinValueInSphere > value) { m_MinValueInSphere = value; m_MinValueIndexInSphere = index; } } } } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > bool MultiGaussianImageSource< TOutputImage > ::IsInImage(IndexType index) { bool isInImage = 0; if(index[0] >= m_Origin[0] && index[0] <= m_Size[0] && index[1] >= m_Origin[1] && index[1] <= m_Size[1] && index[2] >= m_Origin[2] && index[2] <= m_Size[2]) { isInImage = 1; } return isInImage; } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNormalDistributionValues() { double temp[410] = { 0.5 , 0.50399 , 0.50798, 0.51197, 0.51595, 0.51994, 0.52392, 0.5279, 0.53188, 0.53586, 0.53983, 0.5438, 0.54776, 0.55172, 0.55567, 0.55962, 0.56356, 0.56749, 0.57142, 0.57535, 0.57926, 0.58317, 0.58706, 0.59095, 0.59483 , 0.59871, 0.60257, 0.60642, 0.61026, 0.61409, 0.61791, 0.62172, 0.62552, 0.6293, 0.63307, 0.63683, 0.64058, 0.64431, 0.64803, 0.65173, 0.65542, 0.6591, 0.66276, 0.6664, 0.67003, 0.67364, 0.67724, 0.68082, 0.68439, 0.68793, 0.69146, 0.69497, 0.69847, 0.70194, 0.7054, 0.70884, 0.71226, 0.71566, 0.71904, 0.7224, 0.72575, 0.72907, 0.73237, 0.73565, 0.73891, 0.74215, 0.74537, 0.74857, 0.75175, 0.7549, 0.75804, 0.76115, 0.76424, 0.7673, 0.77035, 0.77337, 0.77637, 0.77935, 0.7823, 0.78524, 0.78814, 0.79103, 0.79389, 0.79673, 0.79955, 0.80234, 0.80511, 0.80785, 0.81057, 0.81327, 0.81594, 0.81859, 0.82121, 0.82381, 0.82639, 0.82894, 0.83147, 0.83398, 0.83646, 0.83891, 0.84134, 0.84375, 0.84614, 0.84849, 0.85083, 0.85314, 0.85543, 0.85769, 0.85993, 0.86214, 0.86433, 0.8665, 0.86864, 0.87076, 0.87286, 0.87493, 0.87698, 0.879, 0.881, 0.88298, 0.88493, 0.88686, 0.88877, 0.89065, 0.89251, 0.89435, 0.89617, 0.89796, 0.89973, 0.90147, 0.9032, 0.9049, 0.90658, 0.90824, 0.90988, 0.91149, 0.91309, 0.91466, 0.91621, 0.91774, 0.91924, 0.92073, 0.9222, 0.92364, 0.92507, 0.92647, 0.92785, 0.92922, 0.93056, 0.93189, 0.93319, 0.93448, 0.93574, 0.93699, 0.93822, 0.93943, 0.94062, 0.94179, 0.94295, 0.94408, 0.9452, 0.9463, 0.94738, 0.94845, 0.9495, 0.95053, 0.95154, 0.95254, 0.95352, 0.95449, 0.95543, 0.95637, 0.95728, 0.95818, 0.95907, 0.95994, 0.9608, 0.96164, 0.96246, 0.96327, 0.96407, 0.96485, 0.96562, 0.96638, 0.96712, 0.96784, 0.96856, 0.96926, 0.96995, 0.97062, 0.97128, 0.97193, 0.97257, 0.9732, 0.97381, 0.97441, 0.975, 0.97558, 0.97615, 0.9767, 0.97725, 0.97778, 0.97831, 0.97882, 0.97932, 0.97982, 0.9803, 0.98077, 0.98124, 0.98169, 0.98214, 0.98257, 0.983, 0.98341, 0.98382, 0.98422, 0.98461, 0.985, 0.98537, 0.98574, 0.9861, 0.98645, 0.98679, 0.98713, 0.98745, 0.98778, 0.98809, 0.9884, 0.9887, 0.98899, 0.98928, 0.98956, 0.98983, 0.9901, 0.99036, 0.99061, 0.99086, 0.99111, 0.99134, 0.99158, 0.9918, 0.99202, 0.99224, 0.99245, 0.99266, 0.99286, 0.99305, 0.99324, 0.99343, 0.99361, 0.99379, 0.99396, 0.99413, 0.9943, 0.99446, 0.99461, 0.99477, 0.99492, 0.99506, 0.9952, 0.99534, 0.99547, 0.9956, 0.99573, 0.99585, 0.99598, 0.99609, 0.99621, 0.99632, 0.99643, 0.99653, 0.99664, 0.99674, 0.99683, 0.99693, 0.99702, 0.99711, 0.9972, 0.99728, 0.99736, 0.99744, 0.99752, 0.9976, 0.99767, 0.99774, 0.99781, 0.99788, 0.99795, 0.99801, 0.99807, 0.99813, 0.99819, 0.99825, 0.99831, 0.99836, 0.99841, 0.99846, 0.99851, 0.99856, 0.99861, 0.99865, 0.99869, 0.99874, 0.99878, 0.99882, 0.99886, 0.99889, 0.99893, 0.99896, 0.999, 0.99903, 0.99906, 0.9991, 0.99913, 0.99916, 0.99918, 0.99921, 0.99924, 0.99926, 0.99929, 0.99931, 0.99934, 0.99936, 0.99938, 0.9994, 0.99942, 0.99944, 0.99946, 0.99948, 0.9995, 0.99952, 0.99953, 0.99955, 0.99957, 0.99958, 0.9996, 0.99961, 0.99962, 0.99964, 0.99965, 0.99966, 0.99968, 0.99969, 0.9997, 0.99971, 0.99972, 0.99973, 0.99974, 0.99975, 0.99976, 0.99977, 0.99978, 0.99978, 0.99979, 0.9998, 0.99981, 0.99981, 0.99982, 0.99983, 0.99983, 0.99984, 0.99985, 0.99985, 0.99986, 0.99986, 0.99987, 0.99987, 0.99988, 0.99988, 0.99989, 0.99989, 0.9999, 0.9999, 0.9999, 0.99991, 0.99991, 0.99992, 0.99992, 0.99992, 0.99992, 0.99993, 0.99993, 0.99993, 0.99994, 0.99994, 0.99994, 0.99994, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99998, 0.99998, 0.99998, 0.99998 }; for(int i=0; i < 410; i++) { m_NormalDistValues[i] = temp[i]; } } } // end namespace itk #endif \ No newline at end of file