diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp index 8cbac0c47f..e9bfbfd5d2 100644 --- a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -1,923 +1,923 @@ /*=================================================================== 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 // // // // // // // // // // // //bin\Release\ImageStatisticsTestDriver.exe mitkMultiGaussianTest C:/temp/TestImage/image 1 12 12 10 1 1 5 20 200 "2.5" "2.5" 3 // // bin\Release\ImageStatisticsTestDriver.exe mitkMultiGaussianTest C:/temp/HotSpotTestImage/ImageNeu30 5 30 30 15 1 5 15 20 200 1 1 2 int mitkMultiGaussianTest(int argc, char* argv[]) { //if ( argc != 14 || argc != 3 ) //{ // std::cerr << " 14 arguments expected: [outputFilename, numberOfImages, imageSizeX, imageSizeY, imageSizeZ, numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, spacingX, spacingY, spacingZ ] \n OR \n 3 argument expected : [outputFilename, input.xml] \n " << std::endl; // return EXIT_FAILURE; //} //// always start with this! //MITK_TEST_BEGIN("mitkMultiGaussianTest"); //MITK_TEST_CONDITION_REQUIRED(argc == 14, "Test called with 14 parameters"); if( argc == 14) { 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; unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian; itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; MultiGaussianImageSource::Pointer gaussianGenerator; itk::DOMNodeXMLWriter::Pointer xmlWriter; itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; DOMNodeType domTestCase, domTestImage, domGaussian, 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; if ( ! (std::istringstream(argv[2]) >> numberOfImages) ) numberOfImages = 0; if ( ! (std::istringstream(argv[3]) >> size[0]) ) size[0] = 0; if ( ! (std::istringstream(argv[4]) >> size[1]) ) size[1] = 0; if ( ! (std::istringstream(argv[5]) >> size[2]) ) size[2] = 0; if ( ! (std::istringstream(argv[6]) >> numberOfGaussians) ) numberOfGaussians = 0; if ( ! (std::istringstream(argv[7]) >> minWidthOfGaussian) ) minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; if ( ! (std::istringstream(argv[8]) >> maxWidthOfGaussian) ) maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; if ( ! (std::istringstream(argv[9]) >> minAltitudeOfGaussian) ) minAltitudeOfGaussian = 5; if ( ! (std::istringstream(argv[10]) >> maxAltitudeOfGaussian) ) maxAltitudeOfGaussian = 200; if ( ! (std::istringstream(argv[11]) >> spacing[0]) ) spacing[0] = 1; if ( ! (std::istringstream(argv[12]) >> spacing[1]) ) spacing[1] = 1; if ( ! (std::istringstream(argv[13]) >> spacing[2]) ) spacing[2] = 1; // Set region of interest in pixels regionOfInterestMax.SetElement(0, size[0] - static_cast((radius)/spacing[0] + 0.9999) - 1); regionOfInterestMax.SetElement(1, size[1] - static_cast((radius)/spacing[1] + 0.9999) - 1); regionOfInterestMax.SetElement(2, size[2] - static_cast((radius)/spacing[2] + 0.9999) - 1); regionOfInterestMin.SetElement(0, 0 + static_cast((radius)/spacing[0]+ 0.9999)); regionOfInterestMin.SetElement(1, 0 + static_cast((radius)/spacing[1]+ 0.9999)); regionOfInterestMin.SetElement(2, 0 + static_cast((radius)/spacing[2]+ 0.9999)); srand (time(NULL)); int numberAddGaussian = numberOfGaussians; unsigned int k = 0; unsigned int count = 0; while(k != numberOfImages && count < 500) // for(unsigned int k = 1; k <= numberOfImages; ++k) { ++count; gaussianGenerator = MultiGaussianImageSource::New(); gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); gaussianGenerator->SetRadiusStepNumber(8); gaussianGenerator->SetRadius(radius); gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); // DOM Node Writer 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()); domTestCase->AddChildAtBegin(domTestImage); for( unsigned int i = 0; i < numberAddGaussian; ++i) { domGaussian = itk::DOMNode::New() ; domGaussian->SetName("gaussian"); domTestImage->AddChildAtEnd(domGaussian); // generate the midpoint and the daviation in mm centerX = rand() % static_cast(size[0] * spacing[0]); ss.str(""); ss << centerX; domGaussian->SetAttribute("centerIndexX", ss.str()); centerY = rand() % static_cast(size[1] * spacing[1]); ss.str(""); ss << centerY; domGaussian->SetAttribute("centerIndexY", ss.str()); centerZ = rand() % static_cast(size[2] * spacing[2]); ss.str(""); ss << centerZ; domGaussian->SetAttribute("centerIndexZ", ss.str()); sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); ss.str(""); ss << sigmaX; domGaussian->SetAttribute("deviationX", ss.str()); sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); ss.str(""); ss << sigmaY; domGaussian->SetAttribute("deviationY", ss.str()); sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); ss.str(""); ss << sigmaZ; domGaussian->SetAttribute("deviationZ", ss.str()); altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); ss.str(""); ss << altitude; domGaussian->SetAttribute("altitude", ss.str()); centerXVec.push_back(centerX); centerYVec.push_back(centerY); centerZVec.push_back(centerZ); sigmaXVec.push_back(sigmaX); sigmaYVec.push_back(sigmaY); sigmaZVec.push_back(sigmaZ); altitudeVec.push_back(altitude); } gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); centerXVec.clear(); centerYVec.clear(); centerZVec.clear(); sigmaXVec.clear(); sigmaYVec.clear(); sigmaZVec.clear(); altitudeVec.clear(); try { gaussianGenerator->Update(); gaussianGenerator->CalculateMidpointAndMeanValue(); } catch (std::exception& e) { std::cout << "Error: " << e.what() << std::endl; } //region of interest domROI = itk::DOMNode::New(); domROI->SetName("roi"); domTestCase->AddChildAtEnd(domROI); ss.str(""); ss << radius; domROI->SetAttribute("hotspotRadiusInMM", ss.str()); ss.str(""); ss << regionOfInterestMax[0]; domROI->SetAttribute("maximumSizeX", ss.str()); ss.str(""); ss << regionOfInterestMin[0]; domROI->SetAttribute("minimumSizeX", ss.str()); ss.str(""); ss << regionOfInterestMax[1]; domROI->SetAttribute("maximumSizeY", ss.str()); ss.str(""); ss << regionOfInterestMin[1]; domROI->SetAttribute("minimumSizeY", ss.str()); ss.str(""); ss << regionOfInterestMax[2]; domROI->SetAttribute("maximumSizeZ", ss.str()); ss.str(""); ss << regionOfInterestMin[2]; domROI->SetAttribute("minimumSizeZ", ss.str()); //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("peak", ss.str()); //optimize the mean value in the sphere // gaussianGenerator->OptimizeMeanValue(); 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()); // test if the midpoint of the sphere is not the same as the maximum's index and saves only such a case if( sphereMidpt[0]!= maxValueIndexInSphere[0] && sphereMidpt[1]!= maxValueIndexInSphere[1] && sphereMidpt[2]!= maxValueIndexInSphere[2]) { k++; // .xml (Data) ss.str(""); if(k < 10){ ss << outputFilename <<"00"<< k <<".xml"; }else if(k < 100){ ss << outputFilename <<"0"<< k <<".xml"; }else{ ss << outputFilename << k <<".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(""); if(k < 10) { ss << outputFilename <<"00"<< k <<".nrrd"; }else if(k < 100) { ss << outputFilename <<"0"<< k <<".nrrd"; }else { ss << outputFilename << k <<".nrrd"; } name = ss.str(); fileNamePointer = (char*) name.c_str(); writer->SetFileName( fileNamePointer); writer->SetInput( gaussianImage ); writer->Update(); } } } //-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- // 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 ROI or just its midpoint, but not necessary the whole HotSpot. else { 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; itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; 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]; attributeValue = gaussian->GetAttribute( "centerIndexX" ); std::stringstream(attributeValue) >> value; centerXVec.push_back(value); attributeValue = gaussian->GetAttribute( "centerIndexY" ); std::stringstream(attributeValue) >> value; centerYVec.push_back(value); attributeValue = gaussian->GetAttribute( "centerIndexZ" ); std::stringstream(attributeValue) >> value; centerZVec.push_back(value); std::cout << "Read center of Gaussian (x,y,z): " << centerXVec[i] << ", " << centerYVec[i] << ", " << centerZVec[i] << "\n" << std::endl; attributeValue = gaussian->GetAttribute( "deviationX" ); std::stringstream(attributeValue) >> value; sigmaXVec.push_back(value); attributeValue = gaussian->GetAttribute( "deviationY" ); std::stringstream(attributeValue) >> value; sigmaYVec.push_back(value); attributeValue = gaussian->GetAttribute( "deviationZ" ); std::stringstream(attributeValue) >> value; sigmaZVec.push_back(value); std::cout << "Read deviation of Gaussian (x,y,z): " << 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 << "\n" << std::endl; std::cout << "Read radius in mm : " << hotSpotRadiusInMM << "\n" << 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 nummber: " << 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; + 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; } //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); int numberAddGaussian = numberOfGaussians; for( unsigned int i = 0; i < numberAddGaussian; ++i) { domGaussian = itk::DOMNode::New() ; domGaussian->SetName("gaussian"); domTestImage->AddChildAtEnd(domGaussian); // generate the midpoint and the daviation in mm centerX = centerXVec[i]; ss.str(""); ss << centerX; domGaussian->SetAttribute("centerIndexX", ss.str()); centerY = centerYVec[i]; ss.str(""); ss << centerY; domGaussian->SetAttribute("centerIndexY", ss.str()); centerZ = centerZVec[i]; ss.str(""); ss << centerZ; domGaussian->SetAttribute("centerIndexZ", ss.str()); sigmaX = sigmaXVec[i]; ss.str(""); ss << sigmaX; domGaussian->SetAttribute("deviationX", ss.str()); sigmaY = sigmaYVec[i]; ss.str(""); ss << sigmaY; domGaussian->SetAttribute("deviationY", ss.str()); sigmaZ = sigmaZVec[i]; ss.str(""); ss << sigmaZ; 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); - minSize = 0.0 + static_cast((radius)/spacing[0]+ 0.9999); + 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); - minSize = 0.0 + static_cast((radius)/spacing[1]+ 0.9999); + 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[1] + 0.9999); - minSize = 0.0 + static_cast((radius)/spacing[1]+ 0.9999); + maxSize = size[2] - 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 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->CalculateTheMidPointMeanValueWithOctree(); //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("peak", ss.str()); //optimize the mean value in the sphere // gaussianGenerator->OptimizeMeanValue(); 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()); } /* //The midpoint of the HotSpot is in the image, but not necessary the whole HotSpot else { // set the region of interest for each label i for (unsigned int i = 0; i < numberOfLabels; ++i) { for (unsigned int k = 0 ; k < 3; ++k) { maxSize = ROImaxSizeX[i]; minSize = ROIminSizeX[i]; if(minSize > maxSize) { std::cout << "The sphere is larger then the region of interest! Set the roi to be the whole image " << std::endl; maxSize = size[k]; minSize = 0.0; } // the maximum in the k-Axis regionOfInterestMax.SetElement(k, (maxSize < size[k]) ? maxSize : size[k] - 1.0 ); // the minimum in the k-Axis regionOfInterestMin.SetElement(k, (minSize > 0.0) ? minSize : 0.0 ); } gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); gaussianGenerator->Update(); //TODO // gaussianGenerator->GenerateCuboidSegmentationInSphere(); //gaussianGenerator->CalculateMidpointAndMeanValue(); gaussianGenerator->CalculateTheMidPointMeanValueInCuboid(); //region of interest domROI = itk::DOMNode::New(); domROI->SetName("roi"); domSegmentation->AddChildAtEnd(domROI); ss.str(""); ss << label[i]; domROI->SetAttribute("label", ss.str()); ss.str(""); ss << regionOfInterestMax[0]; domROI->SetAttribute("maximumSizeX", ss.str()); ss.str(""); ss << regionOfInterestMin[0]; domROI->SetAttribute("minimumSizeX", ss.str()); ss.str(""); ss << regionOfInterestMax[1]; domROI->SetAttribute("maximumSizeY", ss.str()); ss.str(""); ss << regionOfInterestMin[1]; domROI->SetAttribute("minimumSizeY", ss.str()); ss.str(""); ss << regionOfInterestMax[2]; domROI->SetAttribute("maximumSizeZ", ss.str()); ss.str(""); ss << regionOfInterestMin[2]; domROI->SetAttribute("minimumSizeZ", ss.str()); //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("peak", ss.str()); //optimize the mean value in the sphere // gaussianGenerator->OptimizeMeanValue(); 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(); // gaussianGenerator -> WriteXMLToTestTheCuboid(); } // MITK_TEST_END() } diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx index 00bb8de463..d9f140db80 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx @@ -1,1663 +1,1809 @@ /*========================================================================= * * 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) { 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); } //for(unsigned int n =0; n < m_NumberOfGaussians; ++n) //{ // 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 = 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] ???? TODO //} + //TODO 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; // TODO qubische interpolation //linear interpolation 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); - - // std::cout << "midpoint: " << globalCoordinateMidpointCuboid << " radius: " << cuboidRadius << " id: " << id << std::endl; - } //---------------------------------------------------------------------------------------------------------------------- /* recursive; This 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, we insert its midpoint 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; - //check if the boundary of the integral is inside the image - if( xMax <= m_Size[0] && xMin >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) - { - double temp = this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; - if(temp > 10000000) - { - std::cout << "temp: " << temp << std::endl; - } + // yz Plane + bool yzPlaneAtOriginCrossXSection = xMin <= m_Origin[0] && xMax >= m_Origin[0]; + bool yzPlaneAtImageBorderCrossXSection = xMin <= m_Size[0] && xMax >= m_Size[0]; + bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0]; + // xz Plane + bool xzPlaneAtOriginCrossYSection = yMin <= m_Origin[1] && yMax >= m_Origin[1]; + bool xzPlaneAtImageBorderCrossYSection = yMin <= m_Size[1] && yMax >= m_Size[1]; + bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1]; + // xy Plane + bool xyPlaneAtOriginCrossZSection = zMin <= m_Origin[2] && zMax >= m_Origin[2]; + bool xyPlaneAtImageBorderCrossZSection = zMin <= m_Size[2] && zMax >= m_Size[2]; + bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2]; + //check if the boundary of the integral is inside the image + // if( xMax <= m_Size[0] && xMin >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) + 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; } + /* // check if the boundary of the image intersect the cuboid and if yes change the limits of the cuboid + else if( // one cross + ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || + (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || + + (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || + (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || + + (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || + (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) + || // two crosses + ( (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) ) + ) + /* if( + ( ( xMin <= m_Origin[0] && xMax >= m_Origin[0] && yMin <= m_Origin[1] && yMax >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) || + ( xMin <= m_Origin[0] && xMax >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMin <= m_Size[2] && zMax >= m_Size[2]) || + (xMin <= m_Size[0] && xMax >= m_Size[0] && ( ( yMax <= m_Size[1] && yMin >= m_Origin[1] ) || ( zMax <= m_Size[2] && zMin >= m_Origin[2] ) ) ) || + (yMin <= m_Origin[1] && yMax >= m_Origin[1] && ( ( xMax <= m_Size[0] && xMin >= m_Origin[0] ) || ( zMax <= m_Size[2] && zMin >= m_Origin[2] ) ) ) || + (yMin <= m_Size[1] && yMax >= m_Size[1] && ( ( xMax <= m_Size[0] && xMin >= m_Origin[0] ) || ( zMax <= m_Size[2] && zMin >= m_Origin[2] ) ) ) || + (zMin <= m_Origin[2] && zMax >= m_Origin[2] && ( ( xMax <= m_Size[0] && xMin >= m_Origin[0] ) || ( yMax <= m_Size[1] && yMin >= m_Origin[1] ) ) ) || + (zMin <= m_Size[2] && zMax >= m_Size[2] && ( (xMax <= m_Size[0] && xMin >= m_Origin[0] ) || ( yMax <= m_Size[1] && yMin >= m_Origin[1] ) ) ) ) + || /* one axis crosses the cuboid + ( ( xMin <= m_Origin[0] && xMax >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) || + (xMin <= m_Size[0] && xMax >= m_Size[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) || + (yMin <= m_Origin[1] && yMax >= m_Origin[1] && xMax <= m_Size[0] && xMin >= m_Origin[0] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) || + (yMin <= m_Size[1] && yMax >= m_Size[1] && xMax <= m_Size[0] && xMin >= m_Origin[0] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) || + (zMin <= m_Origin[2] && zMax >= m_Origin[2] && xMax <= m_Size[0] && xMin >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] ) || + (zMin <= m_Size[2] && zMax >= m_Size[2] && xMax <= m_Size[0] && xMin >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] ) ) ) + { + // x-Axis + if(xMin <= m_Origin[0] && xMax >= m_Origin[0]) + { + xMin = m_Origin[0]; + }else if(xMin <= m_Size[0] && xMax >= m_Size[0]) + { + xMax = m_Size[0]; + } + // y-Axis + if(yMin <= m_Origin[1] && yMax >= m_Origin[1]) + { + yMin = m_Origin[1]; + }else if(yMin <= m_Size[1] && yMax >= m_Size[1]) + { + yMax = m_Size[1]; + } + // z-Axis + if(zMin <= m_Origin[2] && zMax >= m_Origin[2]) + { + zMin = m_Origin[2]; + }else if(zMin <= m_Size[2] && zMax >= m_Size[2]) + { + zMax = m_Size[2]; + } + 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; + m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin) / 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; - if( xMax <= m_Size[0] && xMin >= m_Origin[0] && yMax <= m_Size[1] && yMin >= m_Origin[1] && zMax <= m_Size[2] && zMin >= m_Origin[2] ) + + // yz Plane + bool yzPlaneAtOriginCrossXSection = xMin <= m_Origin[0] && xMax >= m_Origin[0]; + bool yzPlaneAtImageBorderCrossXSection = xMin <= m_Size[0] && xMax >= m_Size[0]; + bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0]; + // xz Plane + bool xzPlaneAtOriginCrossYSection = yMin <= m_Origin[1] && yMax >= m_Origin[1]; + bool xzPlaneAtImageBorderCrossYSection = yMin <= m_Size[1] && yMax >= m_Size[1]; + bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1]; + // xy Plane + bool xyPlaneAtOriginCrossZSection = zMin <= m_Origin[2] && zMax >= m_Origin[2]; + bool xyPlaneAtImageBorderCrossZSection = zMin <= m_Size[2] && zMax >= m_Size[2]; + bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[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 + // same as above + else if( // one cross + ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || + (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || + + (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || + (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || + + (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || + (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) + || // two crosses + ( (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] && xMax >= m_Size[0]) + { + xMax = m_Size[0]; + } + // y-Axis + if(yMin <= m_Origin[1] && yMax >= m_Origin[1]) + { + yMin = m_Origin[1]; + }else if(yMin <= m_Size[1] && yMax >= m_Size[1]) + { + yMax = m_Size[1]; + } + // z-Axis + if(zMin <= m_Origin[2] && zMax >= m_Origin[2]) + { + zMin = m_Origin[2]; + }else if(zMin <= m_Size[2] && zMax >= m_Size[2]) + { + zMax = m_Size[2]; + } + m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); + m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin) ; + + } } } //----------------------------------------------------------------------------------------------------------------------- /* For each voxel calculate the octree (start the recursion) and write this in a file, so that we can visualise it. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateCuboidSegmentationInSphere(PointType globalCoordinateMidpointSphere) { double cuboidRadius; PointType globalCoordinateMidpointCuboid, newMidpoint; - //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); cuboidRadius = m_Radius / 2.0; - m_Volume = 0.0; + // m_Volume = 0.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] = 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_dispVol) { - std::cout << "m_Volume: " << m_Volume < void MultiGaussianImageSource< TOutputImage > ::CalculateTheMidPointMeanValueWithOctree() { 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_dispVol = 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); //for(unsigned int i = 0; i < cuboidNumber; ++i ) //{ // midpoint = m_Midpoints.ElementAt(i); // radius = m_RadiusCuboid.ElementAt(i); // midpoint[0] = midpoint[0] + globalCoordinateMidpointSphere[0]; // midpoint[1] = midpoint[1] + globalCoordinateMidpointSphere[1]; // midpoint[2] = midpoint[2] + globalCoordinateMidpointSphere[2]; // xMin = midpoint[0] - radius; // xMax = midpoint[0] + radius; // yMin = midpoint[1] - radius; // yMax = midpoint[1] + radius; // zMin = midpoint[2] - radius; // zMax = midpoint[2] + radius; // // std::cout << " 2 * radius: " << 2.0 * radius << " xMin: " << xMin << " xMax: " << xMax << " xMax - 2 * radius: " << xMax - 2.0 * radius<< std::endl; // //std::cout << " 2 * radius: " << 2.0 * radius << " yMin: " << yMin << " yMax: " << yMax << " yMax - 2 * radius: " << yMax - 2.0 * radius<< std::endl; // //std::cout << " 2 * radius: " << 2.0 * radius << " zMin: " << zMin << " zMax: " << zMax << " zMax - 2 * radius: " << zMax - 2.0 * radius<< std::endl; // //// std::cout << "( 2 * radius) ^3: " << pow(2.0 * radius, 3.0) << " Produkt: " << (xMax - xMin) * (yMax - yMin) * (zMax - zMin) << std::endl; // m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin); // meanValueTemp = meanValueTemp + MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); //} meanValueTemp = m_meanValueTemp / m_Volume; //((4.0 / 3.0) * itk::Math::pi * m_Radius * m_Radius * m_Radius); - // std::cout << "index: " << index << " meanValue: " << meanValueTemp << std::endl; + //std::cout << " meanValue/ vol(Sphere): " << m_meanValueTemp / ((4.0 / 3.0) * itk::Math::pi * m_Radius * m_Radius * m_Radius) << std::endl; // std::cout << "m_Volume: " << m_Volume << " ... " << (4.0 / 3.0) * itk::Math::pi * m_Radius * m_Radius * m_Radius << std::endl; if(meanValueTemp > m_MeanValue) { m_GlobalCoordinate = globalCoordinateMidpointSphere; m_MeanValue = meanValueTemp; m_SphereMidpoint = index; } } } /* ////----------------------------------------------------------------------------------------------------------------------- //// TODO zu loeschen (ueberpruefen!) //template< class TOutputImage > //void // MultiGaussianImageSource< TOutputImage > // ::CalculateMidpoint() //{ // double valueMean = 0; // m_MeanValue = 0.0; // double sideLength = m_Radius; // double cuboidRadius; // double sideLengthNew = m_Radius; // double meanValueTemp; // bool intersect; // PointType globalCoordinateMidpointCuboid; // 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] +1); // sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); // sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); // regionOfInterest.SetSize(sizeROI); // typename TOutputImage::Pointer image = this->GetOutput(0); // IteratorType regionOfInterestIterator(image, regionOfInterest); // SetNormalDistributionValues(); // for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) // { // cuboidRadius = sideLength / 2.0; // meanValueTemp = 0.0; // index = regionOfInterestIterator.GetIndex(); // image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, // cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, // cuboidRadius, meanValueTemp); // meanValueTemp = meanValueTemp / (m_Radius * itk::Math::pi * itk::Math::pi); // if(meanValueTemp > m_MeanValue) // { // m_GlobalCoordinate = globalCoordinateMidpointSphere; // m_MeanValue = meanValueTemp; // m_SphereMidpoint = index; // } // //to test only one step!!! TODO // regionOfInterestIterator.GoToReverseBegin(); // } //} ////---------------------------------------------------------------------------------------------------------------------- //// TODO zu loeschen (ueberpruefen!) //template< class TOutputImage > //double // MultiGaussianImageSource< TOutputImage > // ::Quadtrees( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, double meanValueTemp) //{ // double xMin, xMax, yMin, yMax, zMin, zMax; // double minDistance = m_Spacing[0] / 2.0; // unsigned int intersect; // while(cuboidRadius > minDistance) // { // intersect = IntersectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius); // if( intersect == 1 ) // { // cuboidRadius = cuboidRadius / 2.0; // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; // globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; // globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; // meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); // } // else if( intersect == 2 ) // { // xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; // xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; // yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; // yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; // zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; // zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // meanValueTemp = meanValueTemp + MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax); // cuboidRadius = minDistance; // } // else if( intersect == 0 ) // { // cuboidRadius = minDistance; // } // } // return meanValueTemp; //} */ //---------------------------------------------------------------------------------------------------------------------- /* - Check if a cuboid intersect the sphere boundary. Returns 0, if the cuboid is inside the sphere; returns 1, if the cuboid intersects the sphere boundary and 2, if the cuboid is out of the sphere. + 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) { 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]) ) ; } /*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) { 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 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. The parameter RadiusStepNumber controls the accuracy of that calculation (the higher the value the higher the exactness). The algorithm works as follows: 1. the first three 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 (use sphere coordinates): 2.1. traverse the radius of the sphere with step size Radius divided by RadiusStepNumber (the for-loop with index i) 2.2. define a variable dist, which gives a approximately distance between the points at the sphere surface (here we take such a distance, that on the smallest sphere are located 8 points) 2.3. calculate the angles so that the points are equally spaced on the surface (for-loops with indexes j and k) 2.3. for all radius length add the values at the points on the sphere and divide by the number of points added (the values at each point is calculate with the method MultiGaussianFunctionValueAtPoint()) 3. Compare with the until-now-found-maximum and take the bigger one */ template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateMidpointAndMeanValue() { itkDebugMacro(<< "Generating a image of scalars "); double numberSummand = 0.0; int fiStepNumber, psiStepNumber; double x, y, z, temp, abst; MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value, meanValueTemp, valueAdd; PointType globalCoordinate; IndexType index; double riStep, fijStep, psikStep, ri, fij, psik, rNew; double dist = itk::Math::pi * m_Radius / static_cast< double >(2 * m_RadiusStepNumber); m_MeanValue = 0.0; riStep = m_Radius / static_cast< double >( m_RadiusStepNumber ); 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] +1); sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); regionOfInterest.SetSize(sizeROI); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType regionOfInterestIterator(image, regionOfInterest); for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) { numberSummand = 1.0; value = regionOfInterestIterator.Get(); m_ValueAtMidpoint = value; index = regionOfInterestIterator.GetIndex(); image->TransformIndexToPhysicalPoint(index, globalCoordinate); ri = riStep; //for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) //{ // angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); // //fij = 0.0; // fijStep = itk::Math::pi / static_cast< double >( angleStepNumberOverTwo ); // fij = fijStep; // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; // for(unsigned int j = 0; j < angleStepNumberOverTwo -1 ; ++j) // from 0 to pi j<= angleStep... // { // z = ri * cos(fij); // psikStep = 2.0 * itk::Math::pi / (2.0 * static_cast< double >( angleStepNumberOverTwo ) ); // psik = -itk::Math::pi + psikStep; // temp = ri * sin(fij); // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; // for(unsigned int k = 0; k < 2 * angleStepNumberOverTwo; ++k) // from -pi to pi // { // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; // x = temp * cos(psik); // y = temp * sin(psik); // numberSummand += 1.0; // // std::cout << "x, y, z: " << x <<" " << y << " " << z << "\n" < xMax ) xMax = tempX; // //if ( tempY > yMax ) yMax = tempY; // //if ( tempZ > zMax ) zMax = tempZ; // // } // fij = fij + fijStep; // } // ri = ri + riStep; //} // second implementation ... corrected! // for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) //{ // angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); // //fij = 0.0; // fijStep = itk::Math::pi / static_cast< double >( angleStepNumberOverTwo ); // fij = fijStep; // // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; // for(unsigned int j = 0; j < angleStepNumberOverTwo -1 ; ++j) // from 0 to pi j<= angleStep... // { // z = ri * cos(fij); // /* double l = 2.0 * ri / static_cast< double >( angleStepNumberOverTwo - 1); // double tempVar = (j + 1) * l; // abst = (tempVar < ri ) ? tempVar : (2.0 * ri - tempVar); // double testVar = abst * (ri - abst); // rNew = std::sqrt(abst * (ri - abst)); */ // if(fij < itk::Math::pi / 2.0) // { // rNew = sin(fij) * ri; // } // else // { // rNew = sin(itk::Math::pi - fij) * ri; // } // psiStepNumber = static_cast( 2.0 * itk::Math::pi * rNew / dist); // psikStep = 2.0 * itk::Math::pi / psiStepNumber ; // psik = -itk::Math::pi + psikStep; // temp = ri * sin(fij); // // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; // for(unsigned int k = 0; k < psiStepNumber; ++k) // from -pi to pi // { // // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; // x = temp * cos(psik); // y = temp * sin(psik); // numberSummand += 1.0; // valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); // value = value + valueAdd; // psik = psik + psikStep; // } // fij = fij + fijStep; // } // ri = ri + riStep; //} for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) { fiStepNumber = static_cast( (itk::Math::pi * ri / dist)); //fij = 0.0; fijStep = itk::Math::pi / (static_cast< double >( fiStepNumber ) * 2.0); fij = fijStep; // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; for(unsigned int j = 0; j < fiStepNumber -2 ; ++j) // from 0 to pi/2 j<= angleStep... { z = ri * cos(fij); if(fij < itk::Math::pi / 2.0) { rNew = sin(fij) * ri; } else { rNew = sin(itk::Math::pi - fij) * ri; } psiStepNumber = static_cast(( 2.0 * itk::Math::pi * rNew / dist) / 4.0); psikStep = itk::Math::pi / (static_cast< double >( psiStepNumber ) * 2.0); psik = psikStep; temp = ri * sin(fij); // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; if(psiStepNumber > 2){ for(unsigned int k = 0; k < psiStepNumber -2; ++k) // from 0 to pi/2 { // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; x = temp * cos(psik); y = temp * sin(psik); numberSummand += 8.0; valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint( - x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint( - x + globalCoordinate[0], - y + globalCoordinate[1], - z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], - y + globalCoordinate[1], z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], - y + globalCoordinate[1], - z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], - z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint(- x + globalCoordinate[0], y + globalCoordinate[1], - z + globalCoordinate[2]); value = value + valueAdd; valueAdd = MultiGaussianFunctionValueAtPoint(- x + globalCoordinate[0], - y + globalCoordinate[1], z + globalCoordinate[2]); value = value + valueAdd; psik = psik + psikStep; } } fij = fij + fijStep; } ri = ri + riStep; } meanValueTemp = value / numberSummand; if(meanValueTemp > m_MeanValue) { m_GlobalCoordinate = globalCoordinate; m_MeanValue = meanValueTemp; m_SphereMidpoint = index; } } } //------------------------------------------------------------------------------------------------------------------------- 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 > ::WriteXMLToTestTheCuboid() { // to create a point set std::stringstream ss; int fiStepNumber, psiStepNumber; double x, y, z, temp; MultiGaussianImageSource< TOutputImage >::OutputImagePixelType meanValueTemp; PointType globalCoordinate; double riStep, fijStep, psikStep, ri, fij, psik, rNew; int angleStepNumber = 4; double dist = 2.0 * itk::Math::pi * m_Radius / static_cast< double >(angleStepNumber * m_RadiusStepNumber); m_MeanValue = 0.0; riStep = m_Radius / static_cast< double >( m_RadiusStepNumber); 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_XCoordToTest.size(); for(unsigned int iter = 0 ; iter < cap; ++iter) { numberSummand += 1.0; 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; ss.str(""); ss << m_XCoordToTest[iter] * scaleFactor; domX->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domX); ss.str(""); ss << m_YCoordToTest[iter] * scaleFactor; domY->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domY); ss.str(""); ss << m_ZCoordToTest[iter] * scaleFactor; domZ->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domZ); domPointSet -> AddChildAtEnd(domPoint); } // .mps (Data) ss.str(""); ss << "C:/temp/CuboidMidpointsSet.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 > ::WriteXMLToTest() { // to create a point set double numberSummand = 0.0; std::stringstream ss; int fiStepNumber, psiStepNumber; double x, y, z, temp; MultiGaussianImageSource< TOutputImage >::OutputImagePixelType meanValueTemp; PointType globalCoordinate; double riStep, fijStep, psikStep, ri, fij, psik, rNew; int angleStepNumber = 4; double dist = 2.0 * itk::Math::pi * m_Radius / static_cast< double >(angleStepNumber * m_RadiusStepNumber); m_MeanValue = 0.0; riStep = m_Radius / static_cast< double >( m_RadiusStepNumber); 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); typename TOutputImage::Pointer image = this->GetOutput(0); image->TransformIndexToPhysicalPoint(m_SphereMidpoint, globalCoordinate); ri = riStep; for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) { fiStepNumber = static_cast( (itk::Math::pi * ri / dist) / 2.0); //fij = 0.0; fijStep = itk::Math::pi / (static_cast< double >( fiStepNumber ) * 2.0); fij = fijStep; // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; for(unsigned int j = 0; j < fiStepNumber - 2 ; ++j) // from 0 to pi/2 j<= angleStep... { z = ri * cos(fij); if(fij < itk::Math::pi / 2.0) { rNew = sin(fij) * ri; } else { rNew = sin(itk::Math::pi - fij) * ri; } psiStepNumber = static_cast(( 2.0 * itk::Math::pi * rNew / dist) / 4.0); psikStep = itk::Math::pi / (static_cast< double >( psiStepNumber ) * 2.0); psik = psikStep; temp = ri * sin(fij); for(unsigned int k = 0; k < psiStepNumber - 2; ++k) // from 0 to pi/2 { x = temp * cos(psik); y = temp * sin(psik); numberSummand += 1.0; 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 = 40.0; ss.str(""); ss << (x + globalCoordinate[0]) * scaleFactor; domX->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domX); ss.str(""); ss << (y + globalCoordinate[1]) * scaleFactor; domY->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domY); ss.str(""); ss << (z + globalCoordinate[2]) * scaleFactor; domZ->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domZ); domPointSet -> AddChildAtEnd(domPoint); psik = psik + psikStep; } fij = fij + fijStep; } ri = ri + riStep; } // .mps (Data) ss.str(""); ss << "C:/temp/SpherePointSet.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 > ::OptimizeMeanValue() { int radiusStepNumberOptimized = m_RadiusStepNumber * 4; int numberSummand = 1, angleStepNumberOverTwo; double x, y, z, temp; double riStep, fijStep, psikStep, ri, fij, psik; double dist = itk::Math::pi * m_Radius / (2 * radiusStepNumberOptimized); MultiGaussianImageSource< TOutputImage >::OutputImagePixelType valueAdd, value, m_MeanValue = 0; riStep = m_Radius / radiusStepNumberOptimized; ri = riStep; value = m_ValueAtMidpoint; for(unsigned int i = 0; i < radiusStepNumberOptimized; ++i) { angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); fij = 0.0; fijStep = itk::Math::pi / angleStepNumberOverTwo; for(unsigned int j = 0; j <= angleStepNumberOverTwo; ++j) // from 0 to pi { z = ri * cos(fij); psikStep = 2.0 * itk::Math::pi / (2.0 * angleStepNumberOverTwo); psik = -itk::Math::pi + psikStep; temp = ri * sin(fij); for(unsigned int k = 0; k < 2 * angleStepNumberOverTwo; ++k) // from -pi to pi { x = temp * cos(psik); y = temp * sin(psik); numberSummand++; valueAdd = MultiGaussianFunctionValueAtPoint(x + m_GlobalCoordinate[0] , y + m_GlobalCoordinate[1], z + m_GlobalCoordinate[2]); value = value + valueAdd; psik = psik + psikStep; } fij = fij + fijStep; } ri = ri + riStep; } m_MeanValue = value / numberSummand; } //------------------------------------------------------------------------------------------------------------------------------------------- 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