From 3f83aa7c95c283b6c3d47cb2a6085b83b009dfd7 Mon Sep 17 00:00:00 2001 From: Kent Williams Date: Wed, 1 May 2013 13:41:34 -0500 Subject: [PATCH] COMP: Remove #if 0 clauses I also removed 3 header files whose entire contents were inside #if 0, and weren't being included anywhere --- BRAINSABC/brainseg/AtlasDefinition.cxx | 16 - .../brainseg/AtlasRegistrationMethod.hxx | 92 +--- BRAINSABC/brainseg/BRAINSABC.cxx | 161 ++----- BRAINSABC/brainseg/BRAINSABCUtilities.hxx | 56 --- BRAINSABC/brainseg/DenoiseFiltering.h | 45 -- BRAINSABC/brainseg/EMSegmentationFilter.hxx | 31 -- BRAINSABC/brainseg/EMSplitPriorMST.h | 392 ------------------ BRAINSABC/brainseg/RegistrationParameters.h | 64 --- BRAINSCommonLib/genericRegistrationHelper.h | 24 -- BRAINSCommonLib/genericRegistrationHelper.hxx | 21 - BRAINSCommonLib/itkIO.h | 2 +- .../src/BRAINSHoughEyeDetector.hxx | 22 - .../src/BRAINSTransformFromFiducials.cxx | 14 - .../src/landmarksConstellationCommon.cxx | 36 -- .../src/landmarksConstellationCommon.h | 15 - .../src/tempSphereFromAreaEstimates.cxx | 53 --- .../liblinear-1.8/linear.cpp | 15 - .../BRAINSCutConfiguration/LandmarkType.h | 33 -- .../ProcessObjectBase.h | 5 - .../ImageRegionPlotter/ImageRegionPlotter.cxx | 12 - BRAINSDemonWarp/VValidationInputParser.hxx | 2 +- BRAINSDemonWarp/commandIterationupdate.h | 14 - BRAINSFit/BRAINSFit.cxx | 49 --- BRAINSFit/BRAINSFitEZ.cxx | 49 --- BRAINSResample/BRAINSResample.cxx | 2 +- .../itkGridForwardWarpImageFilterNew.hxx | 13 - BRAINSTalairach/vtkTalairachGrid.cxx | 338 --------------- .../BRAINSTransformConvert.cxx | 13 - DebugImageViewer/QDebugImageViewerWindow.cxx | 20 - GTRACT/Cmdline/gtractCoregBvalues.cxx | 14 - .../gtractTransformToDisplacementField.cxx | 8 - GTRACT/Common/itkAnatomicalBSplineFilter.cxx | 17 - GTRACT/Common/itkAnatomicalBSplineFilter.h | 24 -- .../Common/itkAnatomicalVersorRigidFilter.h | 18 - .../Common/itkEigenVectorToColorImageFilter.h | 17 - .../itkTimeSeriesVersorScaleSkewFilter.h | 17 - .../itkVectorImageRegisterAffineFilter.h | 19 - .../itkVectorImageRegisterAffineFilter.hxx | 33 -- .../itkVectorImageRegisterVersorRigidFilter.h | 20 +- ...VectorImageRegisterVersorScaleSkewFilter.h | 19 - ICCDEF/ApplyWarp.cxx | 322 ++++++++++++++ ICCDEF/ICCDEFWarp.h | 20 - ICCDEF/IccdefPreprocessor.h | 19 - ICCDEF/IccdefPreprocessor.hxx | 56 --- ICCDEF/IccdefRegistrator.hxx | 25 -- ICCDEF/commandIterationupdate.h | 23 - ICCDEF/iccdefRegistration_New.cxx | 53 --- ICCDEF/itkBrains2LandmarkReader.hxx | 7 - ICCDEF/itkICCDeformableFunction.h | 34 -- ICCDEF/itkICCDeformableFunction.hxx | 159 +------ ICCDEF/itkICCDeformableRegistrationFilter.h | 11 - ICCDEF/itkICCDeformableRegistrationFilter.hxx | 69 --- ...ltiResolutionICCDeformableRegistration.hxx | 11 - 53 files changed, 367 insertions(+), 2257 deletions(-) delete mode 100644 BRAINSABC/brainseg/EMSplitPriorMST.h delete mode 100644 BRAINSABC/brainseg/RegistrationParameters.h delete mode 100644 BRAINSCut/BRAINSCutConfiguration/LandmarkType.h create mode 100644 ICCDEF/ApplyWarp.cxx diff --git a/BRAINSABC/brainseg/AtlasDefinition.cxx b/BRAINSABC/brainseg/AtlasDefinition.cxx index 1ac2e8beb7..023f0d4496 100644 --- a/BRAINSABC/brainseg/AtlasDefinition.cxx +++ b/BRAINSABC/brainseg/AtlasDefinition.cxx @@ -65,22 +65,6 @@ XMLcharhandler(void *data, const char *txt, int txtlen) } } } -#if 0 -const char *AtlasDefinition::tissueTypes[] = - { - "WM", - "GM", - "BGM", - "CSF", - "VB", - "NOTCSF", - "NOTGM", - "NOTVB", - "NOTWM", - "AIR", - 0 - }; -#endif AtlasDefinition::AtlasDefinition() : m_LastWeight(0.0), diff --git a/BRAINSABC/brainseg/AtlasRegistrationMethod.hxx b/BRAINSABC/brainseg/AtlasRegistrationMethod.hxx index 4732eed86a..718630aeb2 100644 --- a/BRAINSABC/brainseg/AtlasRegistrationMethod.hxx +++ b/BRAINSABC/brainseg/AtlasRegistrationMethod.hxx @@ -13,7 +13,6 @@ // MI registration module #include "AtlasRegistrationMethod.h" -#include "RegistrationParameters.h" #include "vnl/vnl_math.h" @@ -605,96 +604,7 @@ AtlasRegistrationMethod << " " << m_AtlasToSubjectTransform->GetParameters() << std::endl ); } } - // End generating the best initial transform for atlas T1 to subject T1 -#if 0 - // Update the registration with all the other non-primary images. - { - for( unsigned int atlasIter = 1; atlasIter < m_AtlasOriginalImageList.size(); atlasIter++ ) - { - // Save code snippets for future integration 2012-10-18 - { // Set the fixed image - InternalImageType::Pointer currentWarpedIntraSubject = NULL; - // NOTE: This is to save memory, so just resample the images as needed - // to avoid keeping an extra copy of them - if( atlasReferenceImageIndex == 0 ) - { // First subject image does not need to be resampled - currentWarpedIntraSubject = m_IntraSubjectOriginalImageList[0]; - } - else - { // All other subject images must be resampled - typedef itk::ResampleImageFilter ResampleType; - typedef ResampleType::Pointer ResamplePointer; - - ResamplePointer resampler = ResampleType::New(); - resampler->SetInput(m_IntraSubjectOriginalImageList[atlasReferenceImageIndex]); - resampler->SetTransform(m_IntraSubjectTransforms[atlasIter]); - resampler->SetOutputParametersFromImage(m_IntraSubjectOriginalImageList[0]); - resampler->SetDefaultPixelValue(0); - resampler->Update(); - currentWarpedIntraSubject = resampler->GetOutput(); - } - atlasToSubjectRegistrationHelper->SetFixedVolume(currentWarpedIntraSubject); - } - std::string preprocessMovingString(""); - { // Set the moving image - // TODO: Just turn histogram matching off at this point. - // histogram matching often help in the registration for difficult - // cases - // Need to find some way to quantify how much it helps/hurts. - // Changed so FLAIR doesn't get histogram matched to non-existing - // atlas - // const bool histogramMatch=true;//Setting histogram matching to true - const bool histogramMatch = false; - // const bool histogramMatch=(m_InputVolumeTypes[atlasIter]=="FLAIR")?false:true; - if( histogramMatch ) - { - typedef itk::HistogramMatchingImageFilter HistogramMatchingFilterType; - HistogramMatchingFilterType::Pointer histogramfilter - = HistogramMatchingFilterType::New(); - - histogramfilter->SetInput( m_AtlasOriginalImageList[atlasIter] ); - histogramfilter->SetReferenceImage( atlasToSubjectRegistrationHelper->GetFixedVolume() ); - - histogramfilter->SetNumberOfHistogramLevels( 128 ); - histogramfilter->SetNumberOfMatchPoints( 2 ); - histogramfilter->ThresholdAtMeanIntensityOn(); - histogramfilter->Update(); - InternalImageType::Pointer equalizedMovingImage = histogramfilter->GetOutput(); - if( equalizedMovingImage->GetLargestPossibleRegion().GetSize() != - m_AtlasOriginalImageList[atlasIter]->GetLargestPossibleRegion().GetSize() ) - { - itkExceptionMacro(<< "Histogram equalized image has wrong size." ); - } - preprocessMovingString = "histogram equalized "; - atlasToSubjectRegistrationHelper->SetMovingVolume(equalizedMovingImage); - if( this->m_DebugLevel > 7 ) - { - typedef itk::ImageFileWriter ShortWriterType; - ShortWriterType::Pointer writer = ShortWriterType::New(); - writer->UseCompressionOn(); - - std::string fn = - this->m_OutputDebugDir + std::string("AtlasPostHistogram_") + ".nii.gz"; - - writer->SetInput( equalizedMovingImage ); - writer->SetFileName(fn.c_str() ); - writer->Update(); - muLogMacro( << __FILE__ << " " << __LINE__ << " " << std::endl ); - itkExceptionMacro(<< "Histogram match"); - } - } - else - { - atlasToSubjectRegistrationHelper->SetMovingVolume(m_AtlasOriginalImageList[atlasIter]); - } - } - // muLogMacro( << "Refining transform with information from atlas " << atlasIter << " of " << - // m_AtlasOriginalImageList.size() << " to subject image " << atlasIter << "." << std::endl); - // Code for updateing and refining the registration for all non-zero index locations. Not implemented yet. - } - } -#endif + // End generating the best initial transform for atlas T1 to subject T1 { muLogMacro(<< "Writing " << this->m_AtlasToSubjectTransformFileName << "." << std::endl); WriteTransformToDisk(m_AtlasToSubjectTransform, this->m_AtlasToSubjectTransformFileName); diff --git a/BRAINSABC/brainseg/BRAINSABC.cxx b/BRAINSABC/brainseg/BRAINSABC.cxx index b9d12ed3c2..f663148044 100644 --- a/BRAINSABC/brainseg/BRAINSABC.cxx +++ b/BRAINSABC/brainseg/BRAINSABC.cxx @@ -804,34 +804,6 @@ int main(int argc, char * *argv) } // Initialize with file read in FloatImageType::Pointer typewiseEqualizedToFirstImage = imgreader->GetOutput(); -#if 0 // This needs more testing. - // Now go looking to see if this image type has already been found, - // and equalize to the first image of this type if found. - for( unsigned int prevImageIndex = 0; prevImageIndex < i; prevImageIndex++ ) - { - if( inputVolumeTypes[i] == inputVolumeTypes[prevImageIndex] ) - // If it matches a previous found image type, - // then histogram equalize - { - muLogMacro( << "Equalizing image (" << i << ") to image (" << prevImageIndex << ")" << std::endl ); - typedef itk::HistogramMatchingImageFilter HistogramMatchingFilterType; - HistogramMatchingFilterType::Pointer histogramfilter - = HistogramMatchingFilterType::New(); - - histogramfilter->SetInput( imgreader->GetOutput() ); - histogramfilter->SetReferenceImage( intraSubjectNoiseRemovedImageList[prevImageIndex] ); - - histogramfilter->SetNumberOfHistogramLevels( 128 ); - histogramfilter->SetNumberOfMatchPoints( 16 ); - // histogramfilter->ThresholdAtMeanIntensityOn(); - histogramfilter->Update(); - // Overwrite if necessary. - typewiseEqualizedToFirstImage = histogramfilter->GetOutput(); - break; - } - } -#endif // Normalize Image Intensities: muLogMacro( << "Standardizing Intensities: ...\n" ); @@ -842,47 +814,43 @@ int main(int argc, char * *argv) 1, 0.95 * MAX_IMAGE_OUTPUT_VALUE, 0, MAX_IMAGE_OUTPUT_VALUE); muLogMacro( << "done.\n" ); -#if 1 + { + std::vector unused_gridSize; + double localFilterTimeStep = filterTimeStep; + if( localFilterTimeStep <= 0 ) { - std::vector unused_gridSize; - double localFilterTimeStep = filterTimeStep; - if( localFilterTimeStep <= 0 ) + FloatImageType::SpacingType::ValueType minPixelSize = + vcl_numeric_limits::max(); + const FloatImageType::SpacingType & imageSpacing = intraSubjectRawImageList[i]->GetSpacing(); + for( int is = 0; is < FloatImageType::ImageDimension; ++is ) { - FloatImageType::SpacingType::ValueType minPixelSize = - vcl_numeric_limits::max(); - const FloatImageType::SpacingType & imageSpacing = intraSubjectRawImageList[i]->GetSpacing(); - for( int is = 0; is < FloatImageType::ImageDimension; ++is ) - { - minPixelSize = vcl_min( minPixelSize, imageSpacing[is]); - } - localFilterTimeStep = - ( (minPixelSize - vcl_numeric_limits::epsilon() ) - / ( vcl_pow(2.0, FloatImageType::ImageDimension + 1 ) ) - ); - } - intraSubjectNoiseRemovedImageList[i] = - DenoiseFiltering(intraSubjectRawImageList[i], filterMethod, filterIteration, - localFilterTimeStep, - unused_gridSize); - if( debuglevel > 1 ) - { - // DEBUG: This code is for debugging purposes only; - typedef itk::ImageFileWriter WriterType; - WriterType::Pointer writer = WriterType::New(); - writer->UseCompressionOn(); - - std::stringstream template_index_stream(""); - template_index_stream << i; - const std::string fn = outputDir + "/DENOISED_INDEX_" + template_index_stream.str() + ".nii.gz"; - writer->SetInput(intraSubjectNoiseRemovedImageList[i]); - writer->SetFileName(fn.c_str() ); - writer->Update(); - muLogMacro( << "DEBUG: Wrote image " << fn << std::endl); + minPixelSize = vcl_min( minPixelSize, imageSpacing[is]); } + localFilterTimeStep = + ( (minPixelSize - vcl_numeric_limits::epsilon() ) + / ( vcl_pow(2.0, FloatImageType::ImageDimension + 1 ) ) + ); } -#else - intraSubjectNoiseRemovedImageList[i] = intraSubjectRawImageList[i]; -#endif + intraSubjectNoiseRemovedImageList[i] = + DenoiseFiltering(intraSubjectRawImageList[i], filterMethod, filterIteration, + localFilterTimeStep, + unused_gridSize); + if( debuglevel > 1 ) + { + // DEBUG: This code is for debugging purposes only; + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->UseCompressionOn(); + + std::stringstream template_index_stream(""); + template_index_stream << i; + const std::string fn = outputDir + "/DENOISED_INDEX_" + template_index_stream.str() + ".nii.gz"; + writer->SetInput(intraSubjectNoiseRemovedImageList[i]); + writer->SetFileName(fn.c_str() ); + writer->Update(); + muLogMacro( << "DEBUG: Wrote image " << fn << std::endl); + } + } intraSubjectTransformFileNames[i] = outputDir + GetStripedImageFileNameExtension(inputVolumes[i]) + std::string( "_to_") @@ -1278,56 +1246,6 @@ int main(int argc, char * *argv) // Start the segmentation process. for( unsigned int segmentationLevel = 0; segmentationLevel < 1; segmentationLevel++ ) { -#if 0 - // Prior Names - // Prior Names : [ AIR GM BGM CSF NOTCSF NOTGM NOTVB NOTWM VB - // WM ] - // Prior weight scales : [ 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 - // 1.00 0.50 ] - // Prior Label Codes : [ 0 2 3 4 6 7 9 8 5 - // 1 ] - std::map > ParentLabels; - // Air - ParentLabels[100].push_back(0); - // GM-BGM-NOTGM - ParentLabels[101].push_back(2); - ParentLabels[101].push_back(3); - ParentLabels[101].push_back(7); - // CSF-NOTCSF - ParentLabels[102].push_back(4); - ParentLabels[102].push_back(6); - // WM-NOTWM - ParentLabels[103].push_back(1); - ParentLabels[103].push_back(8); - // VB-NOTVB - ParentLabels[104].push_back(5); - ParentLabels[104].push_back(9); - - std::map ParentLabelNames; - for( - std::map >::const_iterator plIter = ParentLabels.begin(); - plIter != ParentLabels.end(); plIter++ ) - { - std::string combinedName = ""; - bool firstpass = true; - for( - std::vector::const_iterator vIter = plIter->second.begin(); - vIter != plIter->second.end(); vIter++ ) - { - if( firstpass ) - { - firstpass = false; - } - else - { - combinedName += "-"; - } - combinedName += PriorNames[*vIter]; - } - ParentLabelNames[plIter->first] = combinedName; - } - -#endif SegFilterType::Pointer segfilter = SegFilterType::New(); // __MAX__PROBS @@ -1396,19 +1314,6 @@ int main(int argc, char * *argv) // TODO: Expose the transform type to the BRAINSABC command line // segfilter->SetAtlasTransformType("SyN"); // atlasTransformType); -#if 0 - // THIS ACTAULLY NEEDS TO BE USED FOR ONLY DOING THE FINAL ITERATION OF THE - // SEGMENTATION PROCESS AND SKIPPING MOST OF EMLoop - // If this "Post" transform is given, then jump over the looping, and - // warp priors,do estimates, and bias correct in one step. - { - ReadGenericTransform; - - segfilter->SetTemplateGenericTransform(); - // Turn off warping, and just use the input given from disk. - } -#endif - if( !atlasWarpingOff ) { segfilter->UpdateTransformationOn(); diff --git a/BRAINSABC/brainseg/BRAINSABCUtilities.hxx b/BRAINSABC/brainseg/BRAINSABCUtilities.hxx index bf0e06d2fc..c0e40e4b49 100644 --- a/BRAINSABC/brainseg/BRAINSABCUtilities.hxx +++ b/BRAINSABC/brainseg/BRAINSABCUtilities.hxx @@ -109,38 +109,6 @@ DuplicateImageList(const std::vector & inputList) return outputList; } -#if 0 -template -void WinnerTakesAll(std::vector & listOfImages) -{ - itk::ImageRegionIteratorWithIndex it(listOfImages[0], listOfImages[0]->GetLargestPossibleRegion() ); - - while( !it.IsAtEnd() ) - { - const typename ImageType::IndexType currIndex = it.GetIndex(); - unsigned int currMaxIndex = 0; - typename ImageType::PixelType currMax = listOfImages[0]->GetPixel(currIndex); - typename ImageType::PixelType sum = currMax; - listOfImages[0]->SetPixel(currIndex, 0); - for( unsigned int im = 1; im < listOfImages.size(); im++ ) - { - const typename ImageType::PixelType currValue = listOfImages[im]->GetPixel(currIndex); - sum += currValue; - if( currValue > currMax ) - { - currMax = currValue; - currMaxIndex = im; - } - listOfImages[im]->SetPixel(currIndex, 0); - } - listOfImages[currMaxIndex]->SetPixel(currIndex, sum); - ++it; - } - - return; -} - -#endif template typename ByteImageType::Pointer ComputeForegroundProbMask( @@ -188,30 +156,6 @@ typename ByteImageType::Pointer ComputeForegroundProbMask( } } } -#if 0 - { - // Pre-Dilate mask - typedef itk::BinaryBallStructuringElement StructElementType; - typedef - itk::BinaryDilateImageFilter DilateType; - - StructElementType structel; - structel.SetRadius(1); - structel.CreateStructuringElement(); - - typename DilateType::Pointer dil = DilateType::New(); - dil->SetDilateValue(50); - dil->SetKernel(structel); - dil->SetInput(currForegroundMask); - - dil->Update(); - - ByteImagePointer dilmask = dil->GetOutput(); - // a simple assignment is probably sufficient, test both ways? - currForegroundMask = CopyImage(dilmask); - } -#endif return currForegroundMask; } diff --git a/BRAINSABC/brainseg/DenoiseFiltering.h b/BRAINSABC/brainseg/DenoiseFiltering.h index 5439cd2377..152e5c5f7c 100644 --- a/BRAINSABC/brainseg/DenoiseFiltering.h +++ b/BRAINSABC/brainseg/DenoiseFiltering.h @@ -75,50 +75,5 @@ typename TInputImageType::Pointer DenoiseFiltering( return img; } -#if 0 -< !--Standard options for the denoising phase of a program-- > -< integer - vector > -< name > medianFilterSize -< longflag > medianFilterSize -< label > Median Filter Size -< description > The radius for the optional MedianImageFilter preprocessing in all 3 directions.< / description > -< default > 0, 0, 0 < / default > -< / integer - vector > -< integer > -< name > filterIteration -< description > Filter iterations -< label > Filter Iterations -< longflag > filterIteration -< default > 10 < / default > -< constraints > -< minimum > 0 < / minimum > -< maximum > 50 < / maximum > -< step > 1 < / step > -< / constraints > -< / integer > -< float > -< name > filterTimeStep -< description > Filter time step -< label > Filter Time Step -< longflag > filterTimeStep -< default > 0.01 < / default > -< constraints > -< minimum > 0 < / minimum > -< maximum > 0.5 < / maximum > -< step > 0.01 < / step > -< / constraints > -< / float > -< string - enumeration > -< name > filterMethod -< label > Filter Method -< longflag > filterMethod -< description > Filter method for preprocessing of registration -< element > None -< element > CurvatureFlow -< element > GradientAnisotropicDiffusion -< element > Median -< default > None -< / string - enumeration > -#endif #endif diff --git a/BRAINSABC/brainseg/EMSegmentationFilter.hxx b/BRAINSABC/brainseg/EMSegmentationFilter.hxx index bc296b9b06..69113c0d5b 100644 --- a/BRAINSABC/brainseg/EMSegmentationFilter.hxx +++ b/BRAINSABC/brainseg/EMSegmentationFilter.hxx @@ -1486,19 +1486,6 @@ EMSegmentationFilter atlasToSubjectRegistrationHelper->Update(); unsigned int actualIterations = atlasToSubjectRegistrationHelper->GetActualNumberOfIterations(); muLogMacro( << "Registration tool " << actualIterations << " iterations." << std::endl ); -#if 0 // ERROR: This is not working correctly, because the proper number of - // iterations is not reportd by the optimizer. - while( actualIterations == 0 ) - { - double newGradientTolerance = atlasToSubjectRegistrationHelper->GetProjectedGradientTolerance() * 0.1; - atlasToSubjectRegistrationHelper->SetProjectedGradientTolerance(newGradientTolerance); - muLogMacro( << "Reducing Gradient Tolerance to " << newGradientTolerance << std::endl ); - atlasToSubjectRegistrationHelper->Update(); - actualIterations = atlasToSubjectRegistrationHelper->GetActualNumberOfIterations(); - muLogMacro( << "Registration tool " << actualIterations << " iterations." << std::endl ); - } - -#endif m_TemplateGenericTransform = atlasToSubjectRegistrationHelper->GetCurrentGenericTransform(); } } @@ -1652,15 +1639,6 @@ EMSegmentationFilter this->m_WarpedPriors, this->m_PriorUseForBiasVector, this->m_SampleSpacing, this->m_DebugLevel, this->m_OutputDebugDir); WriteDebugCorrectedImages(this->m_CorrectedImages, 0); -#if 0 // This is probably overkill since the update of NonAirRegion is likely - // not making much in changes at this level of bias correction! Otsu - // works well even in hightly biased images - this->m_NonAirRegion = ComputeTissueRegion(this->m_CorrectedImages[0], 0); - if( this->m_DebugLevel > 5 ) - { - this->WriteDebugHeadRegion(0); - } -#endif this->m_ListOfClassStatistics = ComputeDistributions(SubjectCandidateRegions, this->m_WarpedPriors, this->m_CorrectedImages, @@ -1842,15 +1820,6 @@ EMSegmentationFilter ComputeDistributions(SubjectCandidateRegions, this->m_Posteriors, this->m_CorrectedImages, this->m_DebugLevel); -#if 0 // This is probably overkill since the update of NonAirRegion is likely - // not making much in changes at this level of bias correction! Otsu - // works well even in hightly biased images - this->m_NonAirRegion = ComputeTissueRegion(this->m_CorrectedImages[0], 0); - if( this->m_DebugLevel > 9 ) - { - this->WriteDebugHeadRegion(0); - } -#endif this->m_RawCorrectedImages = CorrectBias(biasdegree, CurrentEMIteration + 100, SubjectCandidateRegions, this->m_RawInputImages, this->m_CleanedLabels, this->m_NonAirRegion, this->m_Posteriors, this->m_PriorUseForBiasVector, diff --git a/BRAINSABC/brainseg/EMSplitPriorMST.h b/BRAINSABC/brainseg/EMSplitPriorMST.h deleted file mode 100644 index 2dcd5dcf9e..0000000000 --- a/BRAINSABC/brainseg/EMSplitPriorMST.h +++ /dev/null @@ -1,392 +0,0 @@ -#if 0 -#ifndef __EMSplitPriorMST_h -#define __EMSplitPriorMST_h -template -void -EMSegmentationFilter -::SplitPriorMST(unsigned int /*iprior*/) -{ - throw; // don't exit -#if 0 - const unsigned int numChannels = m_InputImages.size(); - const unsigned int numPriors = m_WarpedPriors.size(); - - if( iprior >= numPriors ) - { - itkExceptionMacro(<< "Invalid prior index"); - } - - if( m_PriorGaussianClusterCountVector[iprior] == 1 ) - { - return; // No split necessary - } - - unsigned int numClasses = 0; - for( unsigned int i = 0; i < numPriors; i++ ) - { - numClasses += m_PriorGaussianClusterCountVector[i]; - } - - muLogMacro(<< "Splitting distributions for prior " << iprior + 1 << " (MST)\n"); - - InputImageSizeType size - = m_InputImages[0]->GetLargestPossibleRegion().GetSize(); - - ProbabilityImagePointer probImg = m_WarpedPriors[iprior]; - - // Find max prob - double maxP = 0.0; - { - // #pragma omp parallel for - for( long kk = 0; kk < (long)size[2]; kk++ ) - { - for( long jj = 0; jj < (long)size[1]; jj++ ) - { - for( long ii = 0; ii < (long)size[0]; ii++ ) - { - const ProbabilityImageIndexType currIndex = {{ii, jj, kk}}; - if( m_CurrentIterationForegroundMask->GetPixel(currIndex) == 0 ) - { - continue; - } - if( probImg->GetPixel(currIndex) > maxP ) - { - maxP = probImg->GetPixel(currIndex); - } - } - } - } - } - - // Select uniqueClassAverageSamples by thresholding prior with value above tau - const double tau = 0.85 * maxP; - muLogMacro(<< "Sampling with threshold tau = " << tau << "\n"); - - unsigned int numPossibleSamples = 0; - { - // #pragma omp parallel for - for( long kk = 0; kk < (long)size[2]; kk++ ) - { - for( long jj = 0; jj < (long)size[1]; jj++ ) - { - for( long ii = 0; ii < (long)size[0]; ii++ ) - { - const ProbabilityImageIndexType currIndex = {{ii, jj, kk}}; - if( m_CurrentIterationForegroundMask->GetPixel(currIndex) == 0 ) - { - continue; - } - if( probImg->GetPixel(currIndex) >= tau ) - { - numPossibleSamples++; - } - } - } - } - } - - // Make a mapping of unique type names to a list of corresponding corrected - // images indicies - // The uniqueClassAverageSamples should only have dimension of the number of - // unique class types, - // and not the total number of images being corrected. - std::map > uniqueClassMappings; - { - for( unsigned int u = 0; u < m_InputVolumeTypes.size(); u++ ) - { - uniqueClassMappings[m_InputVolumeTypes[u]].push_back(u); - } - } - // Sample selection mask - std::vector uniqueClassAverageSamples; - std::vector uniqueClassAverageIndicies; - size_t numSamples = numPossibleSamples; - { - std::vector selectMask(numPossibleSamples); - for( unsigned int i = 0; i < numPossibleSamples; i++ ) - { - selectMask[i] = 0; - } - - if( numSamples > 50000 ) - { - numSamples = 50000; - } - - muLogMacro(<< " Selecting " << numSamples << " / " << numPossibleSamples << "\n"); - - itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer - rng( itk::Statistics::MersenneTwisterRandomVariateGenerator::New() ); - static unsigned int rngInitializer = 0; - rng->Initialize(rngInitializer); - rngInitializer++; - - if( numSamples < numPossibleSamples ) - { - unsigned int c = 0; - while( c < numSamples ) - { - const unsigned int which = (unsigned int) - rng->GetIntegerVariate(numPossibleSamples - 1); - if( selectMask[which] != 0 ) - { - continue; - } - selectMask[which] = 1; - c++; - } - } - else - { - for( unsigned int i = 0; i < numSamples; i++ ) - { - selectMask[i] = 1; - } - } - - uniqueClassAverageSamples.clear(); - uniqueClassAverageSamples.reserve(numSamples); - uniqueClassAverageIndicies.clear(); - uniqueClassAverageIndicies.reserve(numSamples); - - muLogMacro(<< " Finding uniqueClassAverageSamples...\n"); - - // Compute the muAccumulator - { - const size_t numUniqueClassMappings = uniqueClassMappings.size(); - unsigned int selectMaskIndex = 0; - { - // #pragma omp parallel for - for( long kk = 0; kk < (long)size[2]; kk++ ) - { - for( long jj = 0; jj < (long)size[1]; jj++ ) - { - for( long ii = 0; ii < (long)size[0]; ii++ ) - { - const ProbabilityImageIndexType currIndex = {{ii, jj, kk}}; - if( ( m_CurrentIterationForegroundMask->GetPixel(currIndex) == 0 ) || - ( probImg->GetPixel(currIndex) < tau ) ) - { - continue; - } - - if( selectMask[selectMaskIndex] != 0 ) - { - QHullMSTClusteringProcess::VertexType x(numUniqueClassMappings); - unsigned int uniqueIndex = 0; - for( std::map >::const_iterator mi = uniqueClassMappings.begin(); - mi != uniqueClassMappings.end(); - mi++ ) - { - const std::list & currentListOfImages = mi->second; - const double invCurrentListSize = 1.0 - / ( static_cast( currentListOfImages.size() ) ); - x[uniqueIndex] = 0.0; - for( std::list::const_iterator li = currentListOfImages.begin(); - li != currentListOfImages.end(); - li++ - ) - { - const typename TInputImage::PixelType currPixel = this->m_CorrectedImages[*li]->GetPixel(currIndex); - x[uniqueIndex] += currPixel * invCurrentListSize; - } - uniqueIndex++; - } - uniqueClassAverageSamples.push_back(x); - uniqueClassAverageIndicies.push_back(currIndex); - } - ++selectMaskIndex; - } - } - } - } - } - } - - muLogMacro(<< " Create MST...\n"); - - QHullMSTClusteringProcess mstProc; - mstProc.SetInputVertices(uniqueClassAverageSamples); - mstProc.SortOn(); - - muLogMacro(<< " Allocate maps...\n"); - - // TODO: Make this an std::vector - unsigned int *mapOfSamplesIntoCluster = new unsigned int[numSamples]; - { - for( double T = 2.0; T > 0; T -= 0.01 ) - { - if( T < 1.0 ) - { - itkExceptionMacro(<< "Failed clustering prior " << iprior + 1); - } - muLogMacro(<< " MST clustering, T = " << T << "\n"); - const unsigned int numClusters = mstProc.GetClusters(mapOfSamplesIntoCluster, T); - - if( numClusters < m_PriorGaussianClusterCountVector[iprior] ) - { - continue; - } - - // Check cluster sizes - bool sizeOK = true; - for( unsigned int currentPriorSubClassCluster = 0; - currentPriorSubClassCluster < m_PriorGaussianClusterCountVector[iprior]; - currentPriorSubClassCluster++ ) - { - unsigned int numSamplesInSubClassCluster = 0; - for( unsigned int i = 0; i < numSamples; i++ ) - { - if( mapOfSamplesIntoCluster[i] == currentPriorSubClassCluster ) - { - numSamplesInSubClassCluster++; - } - } - if( numSamplesInSubClassCluster < 10 ) // TODO: Determine what this - // does? - { - sizeOK = false; - } - } - if( sizeOK ) - { - break; - } - } - } - // TODO: Should be able to clear uniqueClassAverageSamples here. - // =============Clustering Done - - // Find beginning class index for this prior - unsigned int istart = 0; - for( unsigned int j = 0; j < iprior; j++ ) - { - istart += m_PriorGaussianClusterCountVector[j]; - } - // Use the largest clusters to estimate the Gaussian parameters - for( unsigned int currentPriorSubClassCluster = 0; - currentPriorSubClassCluster < m_PriorGaussianClusterCountVector[iprior]; - currentPriorSubClassCluster++ ) - { - const unsigned int iclass = istart + currentPriorSubClassCluster; - - unsigned int numSamplesInSubClassCluster = 0; - for( unsigned int i = 0; i < numSamples; i++ ) - { - if( mapOfSamplesIntoCluster[i] == currentPriorSubClassCluster ) - { - numSamplesInSubClassCluster++; - } - } - - muLogMacro(<< " Estimating Gaussian parameters for prior " << iprior + 1 << ":class " << iclass + 1); - muLogMacro(<< " with " << numSamplesInSubClassCluster << " uniqueClassAverageSamples\n"); - - QHullMSTClusteringProcess::VertexType muAccumulator( m_InputVolumeTypes.size() ); - std::fill( muAccumulator.begin(), muAccumulator.end(), 0.0F ); - for( unsigned int ichan = 0; ichan < numChannels; ichan++ ) - { - muAccumulator[ichan] = 0.0; - for( unsigned int i = 0; i < numSamples; i++ ) - { - if( mapOfSamplesIntoCluster[i] == currentPriorSubClassCluster ) - { - const typename TInputImage::PixelType currPixel - = this->m_CorrectedImages[ichan]->GetPixel(uniqueClassAverageIndicies[i]); - muAccumulator[ichan] += currPixel; - } - } - muAccumulator[ichan] /= static_cast( numSamplesInSubClassCluster ); - this->m_ListOfClassStatistics[iclass].m_Means[ichan] = muAccumulator[ichan]; - } - // Compute symmetric covariance matrix. - MatrixType tmpSubClassSampleCovarianceBetweenChannels(numChannels, numChannels); - for( unsigned int channel1Index = 0; channel1Index < numChannels; channel1Index++ ) - { - const double channel1mu = this->m_ListOfClassStatistics[iclass].m_Means[channel1Index]; - for( unsigned int channel2Index = channel1Index; channel2Index < numChannels; channel2Index++ ) - { - const double channel2mu = this->m_ListOfClassStatistics[iclass].m_Means[channel2Index]; - double tmpSubClassClusterCovarianceBetweenChannels12 = 0.0; - for( unsigned int i = 0; i < numSamples; i++ ) - { - if( mapOfSamplesIntoCluster[i] == currentPriorSubClassCluster ) - { - const typename TInputImage::PixelType channel1Pixel - = this->m_CorrectedImages[channel1Index]->GetPixel(uniqueClassAverageIndicies[i]); - const double diff1 = static_cast( channel1Pixel ) - channel1mu; - - const typename TInputImage::PixelType channel2Pixel - = this->m_CorrectedImages[channel2Index]->GetPixel(uniqueClassAverageIndicies[i]); - const double diff2 = static_cast( channel2Pixel ) - channel2mu; - tmpSubClassClusterCovarianceBetweenChannels12 += ( diff1 * diff2 ); - } - } - tmpSubClassClusterCovarianceBetweenChannels12 - /= ( static_cast( numSamplesInSubClassCluster ) - 1.0 + vnl_math::eps ); - if( channel1Index == channel2Index ) - { - tmpSubClassClusterCovarianceBetweenChannels12 += 1e-20; - } - - tmpSubClassSampleCovarianceBetweenChannels(channel1Index, - channel2Index) = tmpSubClassClusterCovarianceBetweenChannels12; - tmpSubClassSampleCovarianceBetweenChannels(channel2Index, - channel1Index) = tmpSubClassClusterCovarianceBetweenChannels12; - } - } - // Scale the covariance up a bit for softer initialization - tmpSubClassSampleCovarianceBetweenChannels *= 1.2; - this->m_ListOfClassStatistics[iclass].m_Covariance = tmpSubClassSampleCovarianceBetweenChannels; - } - for( unsigned int iclass = 0; iclass < numClasses; iclass++ ) - { - const double detcov = vnl_determinant(this->m_ListOfClassStatistics[iclass].m_Covariance); - if( detcov <= 0.0 ) - { - itkExceptionMacro(<< "Determinant of covariance for class " << iclass - << " is <= 0.0 (" << detcov << "), covariance matrix:\n" - << this->m_ListOfClassStatistics[iclass].m_Covariance); - } - } - - uniqueClassAverageSamples.clear(); - uniqueClassAverageIndicies.clear(); // Free the memory. - - delete[] mapOfSamplesIntoCluster; - - // Special case for background, always set the "darkest" mean - // to be the last and set it to zero - if( iprior == ( numPriors - 1 ) ) - { - unsigned int imin = istart; - double minm = this->m_ListOfClassStatistics[istart].m_Means.squared_magnitude(); - for( unsigned int currentPriorSubClassCluster = 1; - currentPriorSubClassCluster < m_PriorGaussianClusterCountVector[iprior]; - currentPriorSubClassCluster++ ) - { - unsigned int iclass = istart + currentPriorSubClassCluster; - double mag = this->m_ListOfClassStatistics[iclass].m_Means.squared_magnitude(); - if( mag < minm ) - { - imin = iclass; - minm = mag; - } - } - - if( imin != ( numClasses - 1 ) ) - { - muLogMacro(<< " Replacing " << this->m_ListOfClassStatistics[imin].m_Means << " with zero\n"); - VectorType v = this->m_ListOfClassStatistics[numClasses - 1].m_Means; - this->m_ListOfClassStatistics[imin].m_Means = v; - } - for( unsigned int ichan = 0; ichan < numChannels; ichan++ ) - { - this->m_ListOfClassStatistics[numClasses - 1].m_Means[ichan] = 0; - } - } -#endif -} - -#endif -#endif diff --git a/BRAINSABC/brainseg/RegistrationParameters.h b/BRAINSABC/brainseg/RegistrationParameters.h deleted file mode 100644 index d76845b836..0000000000 --- a/BRAINSABC/brainseg/RegistrationParameters.h +++ /dev/null @@ -1,64 +0,0 @@ -#if 0 -/******************************************************************************* - Common point for the registration parameters - - Sequence of parameters for affine: - Translation x, y, z - Rotation x, y, z - Scaling x, y, z - Skew x, y, z - ***************************************************************************** */ - -#ifndef _RegistrationParameters_h -#define _RegistrationParameters_h - -#define MU_REGISTER_MAJORVER "2" -#define MU_REGISTER_MINORVER "0a" - -// -// Line search steps -// - -#define MU_AFFINE_STEP_TRANSLATE 10.0 -#define MU_AFFINE_STEP_ROTATE 0.1 -#define MU_AFFINE_STEP_SCALE 0.1 -#define MU_AFFINE_STEP_SKEW 0.01 - -// -// Optimization order -// - -// Translation, rotation, scale, then skew - -#define MU_AFFINE_ORDER0 0 -#define MU_AFFINE_ORDER1 1 -#define MU_AFFINE_ORDER2 2 -#define MU_AFFINE_ORDER3 3 -#define MU_AFFINE_ORDER4 4 -#define MU_AFFINE_ORDER5 5 -#define MU_AFFINE_ORDER6 6 -#define MU_AFFINE_ORDER7 7 -#define MU_AFFINE_ORDER8 8 -#define MU_AFFINE_ORDER9 9 -#define MU_AFFINE_ORDER10 10 -#define MU_AFFINE_ORDER11 11 - -/* -// Translation, scale, rotation, then skew - -#define MU_AFFINE_ORDER0 0 -#define MU_AFFINE_ORDER1 1 -#define MU_AFFINE_ORDER2 2 -#define MU_AFFINE_ORDER3 6 -#define MU_AFFINE_ORDER4 7 -#define MU_AFFINE_ORDER5 8 -#define MU_AFFINE_ORDER6 3 -#define MU_AFFINE_ORDER7 4 -#define MU_AFFINE_ORDER8 5 -#define MU_AFFINE_ORDER9 9 -#define MU_AFFINE_ORDER10 10 -#define MU_AFFINE_ORDER11 11 -*/ - -#endif -#endif diff --git a/BRAINSCommonLib/genericRegistrationHelper.h b/BRAINSCommonLib/genericRegistrationHelper.h index e119ae31f9..0d950148ae 100644 --- a/BRAINSCommonLib/genericRegistrationHelper.h +++ b/BRAINSCommonLib/genericRegistrationHelper.h @@ -78,30 +78,6 @@ void MakeDebugJointHistogram(const std::string & debugOutputDirectory, const typ << "_" << std::setw( 4 ) << std::setfill( '0' ) << currentIteration << ".png"; - // typename JOINTPDFType::ConstPointer = - // typedef itk::CastImageFilter - // > CasterType; - // typedef itk::RescaleIntensityImageFilter > CasterType; - // typedef itk::ShiftScaleImageFilter > CasterType; - // typename CasterType::Pointer myCaster = CasterType::New(); - // myCaster->SetInput(myHistogram); - // myCaster->SetShift(0); - // myCaster->SetScale(255); - // myCaster->Update(); - // -#if 0 - typedef itk::StatisticsImageFilter StatisticsImageFilterType; - typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New(); - statisticsImageFilter->SetInput(myHistogram); - statisticsImageFilter->Update(); - - const float histMean = statisticsImageFilter->GetMean(); - std::cout << "Mean: " << histMean << std::endl; - std::cout << "Std.: " << statisticsImageFilter->GetSigma() << std::endl; - std::cout << "Min: " << statisticsImageFilter->GetMinimum() << std::endl; - std::cout << "Max: " << statisticsImageFilter->GetMaximum() << std::endl; -#endif itk::ImageRegionConstIterator origIter(myHistogram, myHistogram->GetLargestPossibleRegion() ); origIter.GoToBegin(); diff --git a/BRAINSCommonLib/genericRegistrationHelper.hxx b/BRAINSCommonLib/genericRegistrationHelper.hxx index ce0195687c..66d0a4d949 100644 --- a/BRAINSCommonLib/genericRegistrationHelper.hxx +++ b/BRAINSCommonLib/genericRegistrationHelper.hxx @@ -309,27 +309,6 @@ MultiModal3DMutualRegistrationHelperm_Registration->GetNumberOfThreads() << " of " << itk::MultiThreader::GetGlobalDefaultNumberOfThreads() << std::endl; -#if 0 - // if (globalVerbose) - if( 0 ) - { - std::cout << "Initializer, RelaxationFactor: " << m_RelaxationFactor - << "." << std::endl; - std::cout << "Initializer, MaximumStepLength: " << m_MaximumStepLength - << "." << std::endl; - std::cout << "Initializer, MinimumStepLength: " << m_MinimumStepLength - << "." << std::endl; - std::cout << "Initializer, NumberOfIterations: " << m_NumberOfIterations - << "." << std::endl; - std::cout << "Registration, Transform : " << m_Transform << "." - << std::endl; - std::cout << "Registration, FixedImage : " << m_FixedImage << "." - << std::endl; - std::cout << "Registration, MovingImage : " << m_MovingImage << "." - << std::endl; - } -#endif - // Create the Command observer and register it with the optimizer. // TODO: make this output optional. // diff --git a/BRAINSCommonLib/itkIO.h b/BRAINSCommonLib/itkIO.h index 362a79695f..3f37590bc0 100644 --- a/BRAINSCommonLib/itkIO.h +++ b/BRAINSCommonLib/itkIO.h @@ -295,7 +295,7 @@ ScaleAndCast(const typename InputImageType::Pointer & image, } catch( itk::ExceptionObject & e ) { - throw; + throw e; } typename OutputImageType::Pointer returnScaled = RealToProbMapCast->GetOutput(); return returnScaled; diff --git a/BRAINSConstellationDetector/src/BRAINSHoughEyeDetector.hxx b/BRAINSConstellationDetector/src/BRAINSHoughEyeDetector.hxx index f305bde03d..4be1749777 100644 --- a/BRAINSConstellationDetector/src/BRAINSHoughEyeDetector.hxx +++ b/BRAINSConstellationDetector/src/BRAINSHoughEyeDetector.hxx @@ -296,28 +296,6 @@ BRAINSHoughEyeDetector image->TransformIndexToPhysicalPoint( stopIndex, physicalStopLocation ); } -#if 0 - // Compute the bound of image - double bound[6]; - bound[0] = - physicalStartLocation[0] < physicalStopLocation[0] ? - physicalStartLocation[0] : physicalStopLocation[0]; - bound[1] = - physicalStartLocation[0] >= physicalStopLocation[0] ? - physicalStartLocation[0] : physicalStopLocation[0]; - bound[2] = - physicalStartLocation[1] < physicalStopLocation[1] ? - physicalStartLocation[1] : physicalStopLocation[1]; - bound[3] = - physicalStartLocation[1] >= physicalStopLocation[1] ? - physicalStartLocation[1] : physicalStopLocation[1]; - bound[4] = - physicalStartLocation[2] < physicalStopLocation[2] ? - physicalStartLocation[2] : physicalStopLocation[2]; - bound[5] = - physicalStartLocation[2] >= physicalStopLocation[2] ? - physicalStartLocation[2] : physicalStopLocation[2]; -#endif // Space coordinate origin -> center of eye centers of input image VersorVectorType translation1; diff --git a/BRAINSConstellationDetector/src/BRAINSTransformFromFiducials.cxx b/BRAINSConstellationDetector/src/BRAINSTransformFromFiducials.cxx index 3f2ab68426..46cba63ad6 100644 --- a/BRAINSConstellationDetector/src/BRAINSTransformFromFiducials.cxx +++ b/BRAINSConstellationDetector/src/BRAINSTransformFromFiducials.cxx @@ -262,21 +262,7 @@ int main(int argc, char* argv[]) genericTransform = NULL; return EXIT_FAILURE; } -#if 0 - // Convert into an affine transform for saving to slicer - itk::AffineTransform::Pointer atransform = - itk::AffineTransform::New(); - atransform->SetCenter(genericTransform->GetCenter() ); - atransform->SetMatrix(genericTransform->GetMatrix() ); - atransform->SetTranslation(genericTransform->GetTranslation() ); - - itk::TransformFileWriter::Pointer twriter = itk::TransformFileWriter::New(); - twriter->SetInput(atransform); - twriter->SetFileName(saveTransform); - - twriter->Update(); -#endif WriteTransformToDisk(genericTransform.GetPointer(), saveTransform); return EXIT_SUCCESS; diff --git a/BRAINSConstellationDetector/src/landmarksConstellationCommon.cxx b/BRAINSConstellationDetector/src/landmarksConstellationCommon.cxx index 670b391f83..c7f0b018c0 100644 --- a/BRAINSConstellationDetector/src/landmarksConstellationCommon.cxx +++ b/BRAINSConstellationDetector/src/landmarksConstellationCommon.cxx @@ -312,39 +312,9 @@ void ComputeEulerAnglesFromRotationMatrix(const itk::Matrix & m, initialHeadingAngle = vcl_asin(m[1][0]); } -#if 0 -itk::Matrix CreateRotationMatrixFromAngles(const double alpha, const double beta, const double gamma) -{ - // alpha is rotate the X axis -- Attitude - // beta is rotate the Y axis -- Bank - // gamma is rotate the Z axis -- Heading - const double ca = vcl_cos(alpha); - const double sa = vcl_sin(alpha); - const double cb = vcl_cos(beta); - const double sb = vcl_sin(beta); - const double cg = vcl_cos(gamma); - const double sg = vcl_sin(gamma); - - itk::Matrix R; - - R(0, 0) = cb * cg; R(0, 1) = -ca * sg + sa * sb * cg; R(0, 2) = sa * sg + ca * sb * cg; - R(1, 0) = cb * sg; R(1, 1) = ca * cg + sa * sb * sg; R(1, 2) = -sa * cg + ca * sb * sg; - R(2, 0) = -sb; R(2, 1) = sa * cb; R(2, 2) = ca * cb; - itk::Matrix::InternalMatrixType test = - R.GetVnlMatrix() * R.GetTranspose(); - if( !test.is_identity(1.0e-10) ) - { - std::cout << "Computed matrix is not orthogonal!!!" << std::endl; - std::cout << R << std::endl; - } - return R; -} - -#endif itk::Versor CreateRotationVersorFromAngles(const double alpha, const double beta, const double gamma) { -#if 1 // http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles // psi = alpha is rotate the X axis -- Attitude // theta= beta is rotate the Y axis -- Bank @@ -365,12 +335,6 @@ itk::Versor CreateRotationVersorFromAngles(const double alpha, const dou itk::Versor v; v.Set(q[0], q[1], q[2], q[3]); return v; -#else - itk::Matrix R = CreateRotationMatrixFromAngles(alpha, beta, gamma); - itk::Versor v; - v.Set(R); - return v; -#endif } VersorTransformType::Pointer ConvertToVersorRigid3D(RigidTransformType::Pointer RT) diff --git a/BRAINSConstellationDetector/src/landmarksConstellationCommon.h b/BRAINSConstellationDetector/src/landmarksConstellationCommon.h index c38bccd5aa..a4b4b70a90 100644 --- a/BRAINSConstellationDetector/src/landmarksConstellationCommon.h +++ b/BRAINSConstellationDetector/src/landmarksConstellationCommon.h @@ -168,21 +168,6 @@ extern SImageType::PointType GetCenterOfHeadMass(SImageType::Pointer volume); extern SImageType::Pointer MakeIsoTropicReferenceImage(); -#if 0 -// #include -#include -// TODO: Replace double precision with single precision math. -typedef itk::ConstantBoundaryCondition BoundaryConditionType; -static const unsigned int WindowedSincHammingWindowRadius = 5; -typedef itk::Function::HammingWindowFunction WindowFunctionType; -typedef itk::WindowedSincInterpolateImageFunction< - SImageType, - WindowedSincHammingWindowRadius, - WindowFunctionType, - BoundaryConditionType, - double> WindowedSincInterpolatorType; -#endif - /*****************************************************************************/ template ValuesType vectorNorm(const std::vector & x) diff --git a/BRAINSConstellationDetector/src/tempSphereFromAreaEstimates.cxx b/BRAINSConstellationDetector/src/tempSphereFromAreaEstimates.cxx index c946321816..56f6761595 100644 --- a/BRAINSConstellationDetector/src/tempSphereFromAreaEstimates.cxx +++ b/BRAINSConstellationDetector/src/tempSphereFromAreaEstimates.cxx @@ -138,9 +138,6 @@ double FindCenterOfBrainBasedOnTopOfHead(SImageType::Pointer & foreground, SIma // Now foreground holds both the silouette information and the // distance-from-extermum information. } -#if 0 - return true; -#endif // //////////////////////////////////////////////////////////////////////// // This will populate a histogram to make an increasing volume distribution. @@ -149,13 +146,6 @@ double FindCenterOfBrainBasedOnTopOfHead(SImageType::Pointer & foreground, SIma typedef itk::ImageRegionIteratorWithIndex Iterator; -#if 0 - typedef itk::MinimumMaximumImageCalculator minMaxCalcType; - minMaxCalcType::Pointer minMaxCalc = minMaxCalcType::New(); - minMaxCalc->SetImage(foreground); - minMaxCalc->Compute(); -#endif - double maxval = 0; double minval = vcl_numeric_limits::max(); typedef itk::ImageRegionIteratorWithIndex SImageIteratorType; @@ -189,13 +179,6 @@ double FindCenterOfBrainBasedOnTopOfHead(SImageType::Pointer & foreground, SIma size[0] = numBins; std::cout << "Histo values " << minval << " ... " << maxval << " -> " << numBins << std::endl; -#if 0 - itk::ImageFileWriter::Pointer dbgWriter = itk::ImageFileWriter::New(); - dbgWriter->UseCompressionOn(); - dbgWriter->SetFileName("distmap.nii.gz"); - dbgWriter->SetInput(foreground); - dbgWriter->Update(); -#endif HistogramType::MeasurementVectorType minValVector, maxValVector; minValVector[0] = minval; @@ -376,17 +359,6 @@ double FindCenterOfBrainBasedOnTopOfHead(SImageType::Pointer & foreground, SIma // what // is // radius. -#if 0 - std::cout << "HISTOGRAM: " << index << "," << histoIter.GetFrequency() << "," << minValVector[0] << "," - << maxValVector[0] << std::endl; - // std::cout << " SupInf_thickness " << SupInf_thickness; - // std::cout << " RLbyAP_area " << RLbyAP_area; - std::cout << "ForegroundLevel = " << ForegroundLevel; - std::cout << " CUMULATIVEVolume: " << CummulativeVolume; - std::cout << " MAXCrossSectonalArea: " << MaxCrossSectionalArea << " CurrentCrossSectionalArea: " - << CurrentCrossSectionalArea << " DesiredVolumeToIncludeBeforeClipping: " - << DesiredVolumeToIncludeBeforeClipping << std::endl; -#endif if( ( CurrentDistanceFromTopOfHead > 100.0 ) // // Can // not @@ -417,31 +389,6 @@ double FindCenterOfBrainBasedOnTopOfHead(SImageType::Pointer & foreground, SIma << DesiredVolumeToIncludeBeforeClipping << std::endl; exitLoop = true; } -#if 0 // This just does not work in many cases. - if( ( ( CurrentCrossSectionalArea > MaxCrossSectionalArea * .25 ) - && ( CurrentCrossSectionalArea < 0.65 * MaxCrossSectionalArea ) ) // - // If - // a - // constriction - // of - // 75% - // max - // occurs, - // then - // assume - // that - // the - // neck - // area - // has - // been - // entered. - ) - { - std::cout << "NECK CONSTRICTION CRITERIA MET, so exiting." << std::endl; - exitLoop = true; - } -#endif // Now ForegroundLevel represents where to threshold the head from the // neck. } diff --git a/BRAINSContinuousClass/liblinear-1.8/linear.cpp b/BRAINSContinuousClass/liblinear-1.8/linear.cpp index 9951585b57..039c45de43 100644 --- a/BRAINSContinuousClass/liblinear-1.8/linear.cpp +++ b/BRAINSContinuousClass/liblinear-1.8/linear.cpp @@ -131,25 +131,10 @@ static void print_string_stdout(const char *s) static void (*liblinear_print_string) (const char *) = &print_string_stdout; -#if 0 -static void info(const char *fmt, ...) -{ - char buf[BUFSIZ]; - va_list ap; - - va_start(ap, fmt); - vsprintf(buf, fmt, ap); - va_end(ap); - (*liblinear_print_string)(buf); -} - -#else static void info(const char *, ...) { } -#endif - class l2r_lr_fun : public function { public: diff --git a/BRAINSCut/BRAINSCutConfiguration/LandmarkType.h b/BRAINSCut/BRAINSCutConfiguration/LandmarkType.h deleted file mode 100644 index fe83ad5e5d..0000000000 --- a/BRAINSCut/BRAINSCutConfiguration/LandmarkType.h +++ /dev/null @@ -1,33 +0,0 @@ -#if 0 -#ifndef LandmarkType_h -#define LandmarkType_h -#include "ElementParser.h" -#include "StringValue.h" -class LandmarkType : public ElementParser -{ -public: - typedef ElementParser SuperClass; - virtual int PrintSelf(std::ostream & os, int indent) const - { - indent += SuperClass::PrintSelf(os, indent); - os << this->PrintSpaces(indent) << "=== LandmarkType ===" << std::endl; - return indent + 2; - } - - LandmarkType() : ElementParser("Landmark") - { - this->Add(new StringValue("Type", ""), "Type"); - this->Add(new StringValue("Filename", ""), "Filename"); - } -}; - -class LandmarkList : public ElementParser -{ -public: - LandmarkList() : ElementParser("LandmarkList") - { - } -}; - -#endif // xoLandmarkType_h -#endif diff --git a/BRAINSCut/BRAINSCutConfiguration/ProcessObjectBase.h b/BRAINSCut/BRAINSCutConfiguration/ProcessObjectBase.h index 00380edd5d..31fff3d27d 100644 --- a/BRAINSCut/BRAINSCutConfiguration/ProcessObjectBase.h +++ b/BRAINSCut/BRAINSCutConfiguration/ProcessObjectBase.h @@ -56,11 +56,6 @@ class ProcessObjectBase virtual int PrintSelf(std::ostream & os, int indent) const = 0; -#if 0 - { - os << "===ProcessObjectBase===!" << std::endl; - } -#endif const std::string & GetName() const { return m_Name; diff --git a/BRAINSCut/BRAINSFeatureCreators/ImageRegionPlotter/ImageRegionPlotter.cxx b/BRAINSCut/BRAINSFeatureCreators/ImageRegionPlotter/ImageRegionPlotter.cxx index cdd4d8aa6a..7ac522cb41 100644 --- a/BRAINSCut/BRAINSFeatureCreators/ImageRegionPlotter/ImageRegionPlotter.cxx +++ b/BRAINSCut/BRAINSFeatureCreators/ImageRegionPlotter/ImageRegionPlotter.cxx @@ -169,18 +169,6 @@ MapHistogramToImage( typename TImageType::Pointer inputImage, /* * main function */ -#if 0 -void dummy(void) -{ - int * x = malloc( 10 * sizeof(int) ); - - x[10] = 0; - - // not a freed??? - return; -} - -#endif int main(int argc, char *argv[]) { diff --git a/BRAINSDemonWarp/VValidationInputParser.hxx b/BRAINSDemonWarp/VValidationInputParser.hxx index 1ede0038f1..9e05f412a1 100644 --- a/BRAINSDemonWarp/VValidationInputParser.hxx +++ b/BRAINSDemonWarp/VValidationInputParser.hxx @@ -108,7 +108,7 @@ VValidationInputParser catch( itk::ExceptionObject & err ) { std::cerr << "Caught an ITK exception: " << std::endl; - throw; + throw err; } if( this->GetOutDebug() ) { diff --git a/BRAINSDemonWarp/commandIterationupdate.h b/BRAINSDemonWarp/commandIterationupdate.h index 5cbab03a64..2108f651af 100644 --- a/BRAINSDemonWarp/commandIterationupdate.h +++ b/BRAINSDemonWarp/commandIterationupdate.h @@ -233,19 +233,6 @@ class CommandIterationUpdate : public itk::Command DebugImageDisplaySender.SendImage(deffield, 0, 0); DebugImageDisplaySender.SendImage(deffield, 1, 1); DebugImageDisplaySender.SendImage(deffield, 2, 2); -#if 0 - typedef typename itk::WarpImageFilter WarpFilterType; - typename WarpFilterType::Pointer warper = WarpFilterType::New(); - warper->SetInput(m_MovingImage); - warper->SetOutputSpacing( deffield->GetSpacing() ); - warper->SetOutputOrigin( deffield->GetOrigin() ); - warper->SetOutputDirection( deffield->GetDirection() ); - warper->SetDisplacementField(deffield); - warper->Update(); - typename InternalImageType::Pointer - DeformedMovingImagePtr = warper->GetOutput(); -#else typename InternalImageType::Pointer DeformedMovingImagePtr = TransformWarp( @@ -258,7 +245,6 @@ class CommandIterationUpdate : public itk::Command // std::cerr << std::endl << "************IMAGES // SENT*************" << std::endl; } -#endif #endif // defined(USE_DebugImageViewer) m_HarmonicEnergyCalculator->SetImage(deffield); diff --git a/BRAINSFit/BRAINSFit.cxx b/BRAINSFit/BRAINSFit.cxx index b68ec9f679..f20275a033 100644 --- a/BRAINSFit/BRAINSFit.cxx +++ b/BRAINSFit/BRAINSFit.cxx @@ -480,37 +480,6 @@ int main(int argc, char *argv[]) myHelper->Update(); currentGenericTransform = myHelper->GetCurrentGenericTransform(); MovingVolumeType::ConstPointer preprocessedMovingVolume = myHelper->GetPreprocessedMovingVolume(); -#if 0 - if( interpolationMode == "ResampleInPlace" ) - { - % { - VersorRigid3DTransformType::ConstPointer versor3D = - dynamic_cast(currentGenericTransform.GetPointer() ); - if( versor3D.IsNotNull() ) - { - FixedVolumeType::Pointer tempInPlaceResample = itk::SetRigidTransformInPlace( - versor3D.GetPointer(), extractMovingVolume.GetPointer() ); - resampledImage = itkUtil::TypeCast(tempInPlaceResample); - } - else - { - // This should be an exception thow instead of exit. - std::cout << "could not convert to rigid versor type" << std::endl; - return EXIT_FAILURE; - } - } - else - { - // Remember: the Data is Moving's, the shape is Fixed's. - resampledImage = TransformResample( - preprocessedMovingVolume, - extractFixedVolume, - backgroundFillValue, - GetInterpolatorFromString(interpolationMode), - currentGenericTransform); - } - } -#else { typedef float VectorComponentType; typedef itk::Vector VectorPixelType; @@ -524,7 +493,6 @@ int main(int argc, char *argv[]) interpolationMode, false); } -#endif // actualIterations = myHelper->GetActualNumberOfIterations(); // permittedIterations = myHelper->GetPermittedNumberOfIterations(); /* @@ -616,23 +584,6 @@ int main(int argc, char *argv[]) WriteOutImageType>(resampledImage) ); itkUtil::WriteImage(CastImage, outputVolume); } -#if 0 - else if( outputVolumePixelType == "char" ) - { - // itkUtil::WriteCastImage, - // FixedVolumeType>(resampledImage,outputVolume); - typedef itk::Image WriteOutImageType; - WriteOutImageType::Pointer CastImage = - ( scaleOutputValues == true ) ? - ( itkUtil::PreserveCast(resampledImage) ) : - ( itkUtil::TypeCast(resampledImage) ); - itkUtil::WriteImage(CastImage, outputVolume); - } -#endif else if( outputVolumePixelType == "uchar" ) { // itkUtil::WriteCastImageUpdate(); currentGenericTransform = myHelper->GetCurrentGenericTransform(); MovingVolumeType::ConstPointer preprocessedMovingVolume = myHelper->GetPreprocessedMovingVolume(); -#if 0 - if( interpolationMode == "ResampleInPlace" ) - { - % { - VersorRigid3DTransformType::ConstPointer versor3D = - dynamic_cast(currentGenericTransform.GetPointer() ); - if( versor3D.IsNotNull() ) - { - FixedVolumeType::Pointer tempInPlaceResample = itk::SetRigidTransformInPlace( - versor3D.GetPointer(), extractMovingVolume.GetPointer() ); - resampledImage = itkUtil::TypeCast(tempInPlaceResample); - } - else - { - // This should be an exception thow instead of exit. - std::cout << "could not convert to rigid versor type" << std::endl; - return EXIT_FAILURE; - } - } - else - { - // Remember: the Data is Moving's, the shape is Fixed's. - resampledImage = TransformResample( - preprocessedMovingVolume, - extractFixedVolume, - backgroundFillValue, - GetInterpolatorFromString(interpolationMode), - currentGenericTransform); - } - } -#else { typedef float VectorComponentType; typedef itk::Vector VectorPixelType; @@ -526,7 +495,6 @@ int main(int argc, char *argv[]) interpolationMode, false); } -#endif // actualIterations = myHelper->GetActualNumberOfIterations(); // permittedIterations = myHelper->GetPermittedNumberOfIterations(); /* @@ -618,23 +586,6 @@ int main(int argc, char *argv[]) WriteOutImageType>(resampledImage) ); itkUtil::WriteImage(CastImage, outputVolume); } -#if 0 - else if( outputVolumePixelType == "char" ) - { - // itkUtil::WriteCastImage, - // FixedVolumeType>(resampledImage,outputVolume); - typedef itk::Image WriteOutImageType; - WriteOutImageType::Pointer CastImage = - ( scaleOutputValues == true ) ? - ( itkUtil::PreserveCast(resampledImage) ) : - ( itkUtil::TypeCast(resampledImage) ); - itkUtil::WriteImage(CastImage, outputVolume); - } -#endif else if( outputVolumePixelType == "uchar" ) { // itkUtil::WriteCastImage for( LineIteratorType lineIter(outputPtr, refIndex, targetIndex); !lineIter.IsAtEnd(); ++lineIter ) { -#if 0 - typename TOutputImage::IndexType testIndex = lineIter.GetIndex(); - bool lineIsInside = true; - for( unsigned int j = 0; j < ImageDimension; j++ ) - { - if( ( testIndex[j] < FirstIndex[j] ) || ( testIndex[j] >= OnePastValidIndex[j] ) ) - { - lineIsInside = false; - break; - } - } - if( lineIsInside ) -#endif { lineIter.Set(m_ForegroundValue); } diff --git a/BRAINSTalairach/vtkTalairachGrid.cxx b/BRAINSTalairach/vtkTalairachGrid.cxx index cb94dd8cba..1bd96d4355 100644 --- a/BRAINSTalairach/vtkTalairachGrid.cxx +++ b/BRAINSTalairach/vtkTalairachGrid.cxx @@ -130,344 +130,6 @@ vtkPoints * vtkTalairachGrid::GetBoundingBoxGridPoints() return this->boundingBoxGridPoints; } -#if 0 -void vtkTalairachGrid::EstablishBoundingBoxGrid() -{ - // vtkPoints *points = vtkPoints::New(); - // int* dim = GetBoundingBoxDimensions(); - // points->Allocate(dim[0]*dim[1]*dim[2]); - - double pnt[3]; - - std::vector talPnt; - - pnt[0] = IRPPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(0, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(1, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(2, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(3, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(4, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(5, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(6, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(7, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(8, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(9, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(10, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(IRPPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(11, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(12, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(13, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(14, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(15, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(16, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(17, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(18, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(19, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(20, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(21, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(22, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = ACPoint[2]; - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(23, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(24, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(25, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = IRPPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(26, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(27, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(28, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = PCPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(29, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(30, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(31, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = ACPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(32, pnt); - - pnt[0] = IRPPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(33, pnt); - - pnt[0] = PCPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(34, pnt); - - pnt[0] = SLAPoint[0]; - pnt[1] = SLAPoint[1]; - pnt[2] = floor(SLAPoint[2] + 1.0); - talPnt = ConvertPixelPointToTalairachPoint(pnt); - pnt[0] = talPnt[0]; - pnt[1] = talPnt[1]; - pnt[2] = talPnt[2]; - boundingBoxGridPoints->InsertPoint(35, pnt); -} - -#endif - void vtkTalairachGrid::EstablishBoundingBoxGrid() { vtkPoints *points = vtkPoints::New(); diff --git a/BRAINSTransformConvert/BRAINSTransformConvert.cxx b/BRAINSTransformConvert/BRAINSTransformConvert.cxx index 2ca6ba4b01..200e260156 100644 --- a/BRAINSTransformConvert/BRAINSTransformConvert.cxx +++ b/BRAINSTransformConvert/BRAINSTransformConvert.cxx @@ -351,19 +351,6 @@ int main(int argc, char *argv[]) } outputXfrm = scaleSkewVersorXfrm.GetPointer(); } -#if 0 - else if( outputTransformType == "BSplineDeformable" ) - { - BSplineTransformType::Pointer bsplineXfrm = - BSplineTransformType::New(); - if( ExtractTransform(bsplineXfrm, inputXfrm.GetPointer() ) == false ) - { - TransformConvertError(inputXfrm, "BSplineDeformable Transform"); - return EXIT_FAILURE; - } - outputXfrm = bsplineXfrm.GetPointer(); - } -#endif if( outputTransformType == "Same" ) { typedef itk::TransformFileWriter TransformWriterType; diff --git a/DebugImageViewer/QDebugImageViewerWindow.cxx b/DebugImageViewer/QDebugImageViewerWindow.cxx index dee4abe3cb..13d0198388 100644 --- a/DebugImageViewer/QDebugImageViewerWindow.cxx +++ b/DebugImageViewer/QDebugImageViewerWindow.cxx @@ -130,18 +130,6 @@ QDebugImageViewerWindow::stateChanged(QAbstractSocket::SocketState state) case QAbstractSocket::ClosingState: // 0 The socket is not connected. case QAbstractSocket::UnconnectedState: -#if 0 - if( this->m_Socket != 0 ) - { - disconnect(this->m_Socket, SIGNAL(readyRead() ), - this, SLOT(readImage() ) ); - - disconnect(this->m_Socket, SIGNAL(stateChanged(QAbstractSocket::SocketState) ), - this, SLOT(stateChanged(QAbstractSocket::SocketState) ) ); - delete this->m_Socket; - this->m_Socket = 0; - } -#endif break; // 1 The socket is performing a host name lookup. case QAbstractSocket::HostLookupState: @@ -235,14 +223,6 @@ QDebugImageViewerWindow::readImage() std::cerr << "Error reading socket" << std::endl; exit(1); } -#if 0 - std::cerr << "DebugImageViewer: size = " << imageSize - << " spacing = " << imageSpacing << std::endl - << "orientation = " << orientation - << " viewIndex = " << " bufferSize " << bufferSize - << std::endl; - std::cerr.flush(); -#endif if( viewIndex < this->m_ViewCount ) { this->m_ImageDisplayList[viewIndex]->SetImage(xferImage); diff --git a/GTRACT/Cmdline/gtractCoregBvalues.cxx b/GTRACT/Cmdline/gtractCoregBvalues.cxx index 16cc8d53d6..d8206cfce1 100644 --- a/GTRACT/Cmdline/gtractCoregBvalues.cxx +++ b/GTRACT/Cmdline/gtractCoregBvalues.cxx @@ -222,20 +222,6 @@ int main(int argc, char *argv[]) registerImageFilter->SetDebugLevel(debugLevel); registerImageFilter->SetInitializeTransformMode("useMomentsAlign" ); -#if 0 - if( debug ) - { - registerImageFilter->DebugOn(); - } - if( outputTransform.size() > 0 ) - { - registerImageFilter->SetOutputParameterFile( outputTransform ); - } - if( registerB0Only ) - { - registerImageFilter->SetRegisterB0Only( true ); - } -#endif try { diff --git a/GTRACT/Cmdline/gtractTransformToDisplacementField.cxx b/GTRACT/Cmdline/gtractTransformToDisplacementField.cxx index 0920a6be98..de41d9f2a6 100644 --- a/GTRACT/Cmdline/gtractTransformToDisplacementField.cxx +++ b/GTRACT/Cmdline/gtractTransformToDisplacementField.cxx @@ -93,15 +93,7 @@ int main(int argc, char *argv[]) GenericTransformType::Pointer baseTransform = itk::ReadTransformFromDisk(inputTransform); DisplacementFieldGeneratorType::Pointer defGenerator = DisplacementFieldGeneratorType::New(); -#if 0 - defGenerator->SetOutputSize( image->GetLargestPossibleRegion().GetSize() ); - defGenerator->SetOutputSpacing( image->GetSpacing() ); - defGenerator->SetOutputOrigin( image->GetOrigin() ); - defGenerator->SetOutputIndex( image->GetLargestPossibleRegion().GetIndex() ); - defGenerator->SetOutputDirection( image->GetDirection() ); -#else defGenerator->SetOutputParametersFromImage( image ); -#endif defGenerator->SetTransform(baseTransform); try { diff --git a/GTRACT/Common/itkAnatomicalBSplineFilter.cxx b/GTRACT/Common/itkAnatomicalBSplineFilter.cxx index 4eaa15c69f..8e1015f2f9 100644 --- a/GTRACT/Common/itkAnatomicalBSplineFilter.cxx +++ b/GTRACT/Common/itkAnatomicalBSplineFilter.cxx @@ -94,23 +94,6 @@ void AnatomicalBSplineFilter::Update() // bsplineRegion.SetSize( totalGridSize ); -#if 0 - - TransformSpacingType spacing = m_FixedImage->GetSpacing(); - TransformOriginType origin = m_FixedImage->GetOrigin(); - - RegisterImageSizeType fixedImageSize = fixedImageRegion.GetSize(); - for( unsigned int r = 0; r < 3; r++ ) - { - spacing[r] *= vcl_floor( static_cast( fixedImageSize[r] - 1 ) - / static_cast( gridSizeOnImage[r] - 1 ) ); - origin[r] -= spacing[r]; - } - - m_Output->SetGridSpacing( spacing ); - m_Output->SetGridOrigin( origin ); - m_Output->SetGridRegion( bsplineRegion ); -#endif typedef itk::BSplineDeformableTransformInitializer InitializerType; InitializerType::Pointer transformInitializer = InitializerType::New(); diff --git a/GTRACT/Common/itkAnatomicalBSplineFilter.h b/GTRACT/Common/itkAnatomicalBSplineFilter.h index 9f84a5d427..865b428e10 100644 --- a/GTRACT/Common/itkAnatomicalBSplineFilter.h +++ b/GTRACT/Common/itkAnatomicalBSplineFilter.h @@ -121,30 +121,6 @@ class GTRACT_COMMON_EXPORT AnatomicalBSplineFilter : public itk::Object typedef Transform BulkTransformType; typedef BulkTransformType::ConstPointer BulkTransformPointer; -#if 0 - /* - void SetBulkTransform( BulkTransformPointer &BulkTransform ) - { - m_Output->SetBulkTransform( BulkTransform ); - } - */ - - /** ImageDimension constants * / - itkStaticConstMacro(InputImageDimension, unsigned int, - TInputImage::ImageDimension); - itkStaticConstMacro(OutputImageDimension, unsigned int, - TOutputImage::ImageDimension); - - / ** The dimensions of the input image must equal those of the - output image. * / - itkConceptMacro(SameDimension, - (Concept::SameDimension)); - - / ** The dimension of the input image must be 4. * / - itkConceptMacro(DimensionShouldBe4, - (Concept::SameDimension)); - */ -#endif /** Standard New method. */ itkNewMacro(Self); diff --git a/GTRACT/Common/itkAnatomicalVersorRigidFilter.h b/GTRACT/Common/itkAnatomicalVersorRigidFilter.h index 2d208ecbd5..b1b32f4fa4 100644 --- a/GTRACT/Common/itkAnatomicalVersorRigidFilter.h +++ b/GTRACT/Common/itkAnatomicalVersorRigidFilter.h @@ -102,24 +102,6 @@ class GTRACT_COMMON_EXPORT AnatomicalVersorRigidFilter : public itk::Object typedef RegistrationType::Pointer RegistrationTypePointer; typedef TransformInitializerType::Pointer TransformInitializerTypePointer; -#if 0 - / **ImageDimension constants - * / itkStaticConstMacro(InputImageDimension, unsigned int, - TInputImage::ImageDimension); - itkStaticConstMacro(OutputImageDimension, unsigned int, - TOutputImage::ImageDimension); - - / **The dimensions of the input image must equal those of the - output image. * - / itkConceptMacro( SameDimension, - ( Concept::SameDimension ) ); - - / **The dimension of the input image must be 4. - * / itkConceptMacro( DimensionShouldBe4, - ( Concept::SameDimension ) ); - * / -#endif /** Standard New method. */ itkNewMacro(Self); diff --git a/GTRACT/Common/itkEigenVectorToColorImageFilter.h b/GTRACT/Common/itkEigenVectorToColorImageFilter.h index 269578ddcb..770a08245a 100644 --- a/GTRACT/Common/itkEigenVectorToColorImageFilter.h +++ b/GTRACT/Common/itkEigenVectorToColorImageFilter.h @@ -80,23 +80,6 @@ class EigenVectorToColorImageFilter : public itk::Object typedef OutputImageType::Pointer OutputImagePointer; typedef OutputImageType::RegionType OutputImageRegionType; -#if 0 - /** ImageDimension constants */ - itkStaticConstMacro(InputImageDimension, unsigned int, - TInputImage::ImageDimension); - itkStaticConstMacro(OutputImageDimension, unsigned int, - TOutputImage::ImageDimension); - - /** The dimensions of the input image must equal those of the - output image. */ - itkConceptMacro( SameDimension, - ( Concept::SameDimension ) ); - -/** The dimension of the input image must be 4. */ - itkConceptMacro( DimensionShouldBe4, - ( Concept::SameDimension ) ); -#endif /** Standard New method. */ itkNewMacro(Self); diff --git a/GTRACT/Common/itkTimeSeriesVersorScaleSkewFilter.h b/GTRACT/Common/itkTimeSeriesVersorScaleSkewFilter.h index 74316736c7..f1326d0e7c 100644 --- a/GTRACT/Common/itkTimeSeriesVersorScaleSkewFilter.h +++ b/GTRACT/Common/itkTimeSeriesVersorScaleSkewFilter.h @@ -124,23 +124,6 @@ class GTRACT_COMMON_EXPORT TimeSeriesVersorScaleSkewFilter : public itk::Object typedef TransformInitializerType::Pointer TransformInitializerTypePointer; typedef ResampleFilterType::Pointer ResampleFilterTypePointer; -#if 0 - /** ImageDimension constants */ - itkStaticConstMacro(InputImageDimension, unsigned int, - TInputImage::ImageDimension); - itkStaticConstMacro(OutputImageDimension, unsigned int, - TOutputImage::ImageDimension); - - /** The dimensions of the input image must equal those of the - output image. */ - itkConceptMacro( SameDimension, - ( Concept::SameDimension ) ); - -/** The dimension of the input image must be 4. */ - itkConceptMacro( DimensionShouldBe4, - ( Concept::SameDimension ) ); -#endif /** Standard New method. */ itkNewMacro(Self); diff --git a/GTRACT/Common/itkVectorImageRegisterAffineFilter.h b/GTRACT/Common/itkVectorImageRegisterAffineFilter.h index 3ac9add74b..c4e0525024 100644 --- a/GTRACT/Common/itkVectorImageRegisterAffineFilter.h +++ b/GTRACT/Common/itkVectorImageRegisterAffineFilter.h @@ -90,16 +90,8 @@ class ITK_EXPORT VectorImageRegisterAffineFilter : typedef typename InputImageType::InternalPixelType InputImageValueType; typedef typename InputImagePixelType::ComponentType InputImageVectorComponentType; -#if 0 - /* The pixel type of FixedImageType should be based on the - * InputImageValueType. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded - */ - typedef itk::Image FixedImageType; -#else // typedef typename itk::Image FixedImageType; typedef typename itk::Image FixedImageType; -#endif typedef typename FixedImageType::Pointer FixedImagePointer; typedef typename FixedImageType::PixelType FixedImagePixelType; @@ -138,19 +130,8 @@ class ITK_EXPORT VectorImageRegisterAffineFilter : FixedImageType, FixedImageType> ResampleFilterType; -#if 0 - /* The output image type should be CastImageType for the - * CastFilter. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded to be - * the FixedImageType to make it a scalar - */ - typedef itk::CastImageFilter< - FixedImageType, - FixedImageType> CastFilterType; -#else typedef itk::CastImageFilter CastFilterType; -#endif typedef typename VectorIndexFilterType::Pointer VectorIndexFilterPointer; typedef typename TransformType::Pointer TransformTypePointer; diff --git a/GTRACT/Common/itkVectorImageRegisterAffineFilter.hxx b/GTRACT/Common/itkVectorImageRegisterAffineFilter.hxx index d5c343ea4b..b63b355b9b 100644 --- a/GTRACT/Common/itkVectorImageRegisterAffineFilter.hxx +++ b/GTRACT/Common/itkVectorImageRegisterAffineFilter.hxx @@ -184,39 +184,6 @@ void VectorImageRegisterAffineFilter OptimizerParameterType finalParameters = registration->GetLastTransformParameters(); -#if 0 - const double rotatorXX = finalParameters[0]; - const double rotatorXY = finalParameters[1]; - const double rotatorXZ = finalParameters[2]; - const double rotatorYX = finalParameters[3]; - const double rotatorYY = finalParameters[4]; - const double rotatorYZ = finalParameters[5]; - const double rotatorZX = finalParameters[6]; - const double rotatorZY = finalParameters[7]; - const double rotatorZZ = finalParameters[8]; - const double finalTranslationX = finalParameters[9]; - const double finalTranslationY = finalParameters[10]; - const double finalTranslationZ = finalParameters[11]; - const unsigned int numberOfIterations = optimizer->GetCurrentIteration(); - const double bestValue = optimizer->GetValue(); - - itkDebugMacro("\tResult = "); - itkDebugMacro("\t\trotator XX = " << rotatorXX); - itkDebugMacro("\t\trotator XY = " << rotatorXY); - itkDebugMacro("\t\trotator XZ = " << rotatorXZ); - itkDebugMacro("\t\trotator YX = " << rotatorYX); - itkDebugMacro("\t\trotator YY = " << rotatorYY); - itkDebugMacro("\t\trotator YZ = " << rotatorYZ); - itkDebugMacro("\t\trotator ZX = " << rotatorZX); - itkDebugMacro("\t\trotator ZY = " << rotatorZY); - itkDebugMacro("\t\trotator ZZ = " << rotatorZZ); - itkDebugMacro("\t\tTranslation X = " << finalTranslationX); - itkDebugMacro("\t\tTranslation Y = " << finalTranslationY); - itkDebugMacro("\t\tTranslation Z = " << finalTranslationZ); - itkDebugMacro("\t\tIterations = " << numberOfIterations); - itkDebugMacro("\t\tMetric value = " << bestValue); -#endif - transform->SetParameters( finalParameters ); /* This step can probably be removed */ diff --git a/GTRACT/Common/itkVectorImageRegisterVersorRigidFilter.h b/GTRACT/Common/itkVectorImageRegisterVersorRigidFilter.h index 3ac88ec231..4faac341ae 100644 --- a/GTRACT/Common/itkVectorImageRegisterVersorRigidFilter.h +++ b/GTRACT/Common/itkVectorImageRegisterVersorRigidFilter.h @@ -90,16 +90,9 @@ class ITK_EXPORT VectorImageRegisterVersorRigidFilter : typedef typename InputImageType::InternalPixelType InputImageValueType; typedef typename InputImagePixelType::ComponentType InputImageVectorComponentType; -#if 0 - /* The pixel type of FixedImageType should be based on the - * InputImageValueType. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded - */ - typedef itk::Image FixedImageType; -#else // typedef typename itk::Image FixedImageType; typedef typename itk::Image FixedImageType; -#endif + typedef typename FixedImageType::Pointer FixedImagePointer; typedef typename FixedImageType::PixelType FixedImagePixelType; @@ -138,19 +131,8 @@ class ITK_EXPORT VectorImageRegisterVersorRigidFilter : FixedImageType, FixedImageType> ResampleFilterType; -#if 0 - /* The output image type should be CastImageType for the - * CastFilter. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded to be - * the FixedImageType to make it a scalar - */ - typedef itk::CastImageFilter< - FixedImageType, - FixedImageType> CastFilterType; -#else typedef itk::CastImageFilter CastFilterType; -#endif typedef typename VectorIndexFilterType::Pointer VectorIndexFilterPointer; typedef typename TransformType::Pointer TransformTypePointer; diff --git a/GTRACT/Common/itkVectorImageRegisterVersorScaleSkewFilter.h b/GTRACT/Common/itkVectorImageRegisterVersorScaleSkewFilter.h index d60d324e5d..bf598d4a7e 100644 --- a/GTRACT/Common/itkVectorImageRegisterVersorScaleSkewFilter.h +++ b/GTRACT/Common/itkVectorImageRegisterVersorScaleSkewFilter.h @@ -91,16 +91,8 @@ class ITK_EXPORT VectorImageRegisterVersorScaleSkewFilter : typedef typename InputImageType::InternalPixelType InputImageValueType; typedef typename InputImagePixelType::ComponentType InputImageVectorComponentType; -#if 0 - /* The pixel type of FixedImageType should be based on the - * InputImageValueType. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded - */ - typedef itk::Image FixedImageType; -#else // typedef typename itk::Image FixedImageType; typedef typename itk::Image FixedImageType; -#endif typedef typename FixedImageType::Pointer FixedImagePointer; typedef typename FixedImageType::PixelType FixedImagePixelType; @@ -139,19 +131,8 @@ class ITK_EXPORT VectorImageRegisterVersorScaleSkewFilter : FixedImageType, FixedImageType> ResampleFilterType; -#if 0 - /* The output image type should be CastImageType for the - * CastFilter. However, the Mac compiler 4.0.1 - * does not like this. Currently this is hardcoded to be - * the FixedImageType to make it a scalar - */ - typedef itk::CastImageFilter< - FixedImageType, - FixedImageType> CastFilterType; -#else typedef itk::CastImageFilter CastFilterType; -#endif typedef typename VectorIndexFilterType::Pointer VectorIndexFilterPointer; typedef typename TransformType::Pointer TransformTypePointer; diff --git a/ICCDEF/ApplyWarp.cxx b/ICCDEF/ApplyWarp.cxx new file mode 100644 index 0000000000..221d262af7 --- /dev/null +++ b/ICCDEF/ApplyWarp.cxx @@ -0,0 +1,322 @@ +/*================================================================== + +TODO: NEED TO COMMENT WHAT THIS PROGRAM IS TO BE USED FOR + +==================================================================*/ + +#include +#include "itkVector.h" +#include "itkImage.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkWarpImageFilter.h" +#include "itkLinearInterpolateImageFunction.h" +#include "itkNearestNeighborInterpolateImageFunction.h" +#include "itkBSplineInterpolateImageFunction.h" +#include "itkWindowedSincInterpolateImageFunction.h" +#include "ApplyWarpCLP.h" + +#include "itkBinaryThresholdImageFilter.h" +#include "itkSignedMaurerDistanceMapImageFilter.h" +#include "itkStatisticsImageFilter.h" +#include "itkIO.h" + +#include "GenericTransformImage.h" + +// A filter to debug the min/max values +template +void PrintImageMinAndMax(TImage * inputImage) +{ +// typename TImage::PixelType resultMaximum: +// typename TImage::PixelType resultMinimum; + typedef typename itk::StatisticsImageFilter StatisticsFilterType; + typename StatisticsFilterType::Pointer statsFilter = StatisticsFilterType::New(); + statsFilter->SetInput( inputImage ); + statsFilter->Update(); +// resultMaximum = statsFilter->GetMaximum(); +// resultMinimum = statsFilter->GetMinimum(); + std::cerr << "StatisticsFilter gave Minimum of " << statsFilter->GetMinimum() + << " and Maximum of " << statsFilter->GetMaximum() << std::endl; +} + +int ApplyWarp(int argc, char *argv[]) +{ + PARSE_ARGS; + + const bool useTransform = (warpTransform.size() > 0); + { + const bool useDisplacementField = (deformationVolume.size() > 0); + const bool debug = true; + + if( debug ) + { + std::cout << "=====================================================" << std::endl; + std::cout << "Input Volume: " << inputVolume << std::endl; + std::cout << "Reference Volume: " << referenceVolume << std::endl; + std::cout << "Output Volume: " << outputVolume << std::endl; + std::cout << "Pixel Type: " << pixelType << std::endl; + std::cout << "Orientation to RAI:" << orientationRAI << std::endl; + std::cout << "Interpolation: " << interpolationMode << std::endl; + std::cout << "Background Value: " << defaultValue << std::endl; + if( useDisplacementField ) + { + std::cout << "Warp by Deformation Volume: " << deformationVolume << std::endl; + } + if( useTransform ) + { + std::cout << "Warp By Transform: " << warpTransform << std::endl; + } + std::cout << "=====================================================" << std::endl; + } + + if( useTransformMode.size() > 0 ) + { + std::cout + << + "Scripting 'code rot' note: The useTransformMode parameter will be ignored. Now ApplyWarp infers the warpTransform type from the contents of the .mat file." + << std::endl; + } + + if( useTransform == useDisplacementField ) + { + std::cout + << "Choose one of the two possibilities, a BRAINSFit transform --or-- a high-dimensional deformation field." + << std::endl; + exit(1); + } + } + + typedef itk::Image ImageType; + typedef itk::Image RefImageType; + ImageType::Pointer PrincipalOperandImage; // One name for the image to be warped. + { + + if( orientationRAI ) + { + PrincipalOperandImage = itkUtil::ReadImage(inputVolume); + PrincipalOperandImage = itkUtil::OrientImage(PrincipalOperandImage, + itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI); + } + else + { + PrincipalOperandImage = itkUtil::ReadImage(inputVolume); + } + } + + // Read ReferenceVolume and DeformationVolume + + typedef float VectorComponentType; + typedef itk::Vector VectorPixelType; + typedef itk::Image DisplacementFieldType; + + // An empty SmartPointer constructor sets up someImage.IsNull() to represent a not-supplied state: + DisplacementFieldType::Pointer DisplacementField; + RefImageType::Pointer ReferenceImage; + + if( useTransform ) + { + + if( referenceVolume.size() > 0 ) + { + if( orientationRAI ) + { + ReferenceImage = itkUtil::ReadImage(referenceVolume); + ReferenceImage = itkUtil::OrientImage(ReferenceImage, + itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI); + } + else + { + ReferenceImage = itkUtil::ReadImage(referenceVolume); + } + } + else + { + std::cout << "Alert: missing Reference Volume defaulted to: " << inputVolume << std::endl; + // ReferenceImage = itkUtil::ReadImage( inputVolume ); + if( orientationRAI ) + { + ReferenceImage = itkUtil::ReadImage(inputVolume); + ReferenceImage = itkUtil::OrientImage(ReferenceImage, + itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI); + } + else + { + ReferenceImage = itkUtil::ReadImage(inputVolume); + } + } + } + else if( !useTransform ) // that is, it's a warp by deformation field: + { + if( orientationRAI ) + { + DisplacementField = itkUtil::ReadImage(deformationVolume); + DisplacementField = itkUtil::OrientImage( + DisplacementField, itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI); + } + else + { + DisplacementField = itkUtil::ReadImage(deformationVolume); + } + if( referenceVolume.size() > 0 ) + { + if( orientationRAI ) + { + ReferenceImage = itkUtil::ReadImage(referenceVolume); + ReferenceImage = itkUtil::OrientImage(ReferenceImage, + itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI); + } + else + { + ReferenceImage = itkUtil::ReadImage(referenceVolume); + } + } + } + + // Read optional transform: + + // An empty SmartPointer constructor sets up someTransform.IsNull() to represent a not-supplied state: + BSplineTransformType::Pointer itkBSplineTransform; + AffineTransformType::Pointer ITKAffineTransform; + + if( useTransform ) + { + std::cerr << "Invalid option, not implemented in ITKv4 version yet." << std::endl; + return -1; + } + + ImageType::Pointer TransformedImage + = GenericTransformImage( + PrincipalOperandImage, + ReferenceImage, + DisplacementField, + NULL, + defaultValue, + interpolationMode, + pixelType == "binary"); + + // Write out the output image; threshold it if necessary. + if( pixelType == "binary" ) + { + // A special case for dealing with binary images + // where signed distance maps are warped and thresholds created + typedef short int MaskPixelType; + typedef itk::Image MaskImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + + MaskImageType::Pointer outputImage = castFilter->GetOutput(); + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "uchar" ) + { + typedef unsigned char NewPixelType; + typedef itk::Image NewImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "short" ) + { + typedef signed short NewPixelType; + typedef itk::Image NewImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "ushort" ) + { + typedef unsigned short NewPixelType; + typedef itk::Image NewImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "int" ) + { + typedef int NewPixelType; + typedef itk::Image NewImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "uint" ) + { + typedef unsigned int NewPixelType; + typedef itk::Image NewImageType; + typedef itk::CastImageFilter CastImageFilter; + CastImageFilter::Pointer castFilter = CastImageFilter::New(); + castFilter->SetInput( TransformedImage ); + castFilter->Update(); + ; + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( castFilter->GetOutput() ); + imageWriter->Update(); + } + else if( pixelType == "float" ) + { + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer imageWriter = WriterType::New(); + imageWriter->SetFileName( outputVolume ); + imageWriter->SetInput( TransformedImage ); + imageWriter->Update(); + } + else + { + std::cout << "ERROR: Invalid pixelType" << std::endl; + exit(-1); + } + + return EXIT_SUCCESS; +} + +int main( int argc, char *argv[] ) +{ + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + std::cout << "This program has been replaced by BRAINSResample. PLEASE TRY TO AVOID USING THIS!" << std::endl; + + return ApplyWarp(argc, argv); +} diff --git a/ICCDEF/ICCDEFWarp.h b/ICCDEF/ICCDEFWarp.h index 7e2b468d9a..637bc451cd 100644 --- a/ICCDEF/ICCDEFWarp.h +++ b/ICCDEF/ICCDEFWarp.h @@ -137,26 +137,6 @@ class ICCDEFWarp : public ICCApplicationBase< itkSetStringMacro(MovingBinaryVolume); itkGetStringMacro(MovingBinaryVolume); -#if 0 - /** Set/Get the lower threshold. The default is 0. */ - itkSetMacro(Lower, PixelType); - itkGetMacro(Lower, PixelType); - - /** Set/Get the upper threshold. The default is 70 */ - itkSetMacro(Upper, PixelType); - itkGetMacro(Upper, PixelType); - - /** Set the Seed of the neighborhood used for a mask. */ - itkSetMacro(Seed, IndexType); - /** Get the radius of the neighborhood used to compute the median */ - itkGetConstReferenceMacro(Seed, IndexType); - - /** Set the radius of the neighborhood used for a mask. */ - itkSetMacro(Radius, SizeType); - /** Get the radius of the neighborhood used to compute the median */ - itkGetConstReferenceMacro(Radius, SizeType); -#endif - /** Set/Get value to replace thresholded pixels. Pixels that lie * * within Lower and Upper (inclusive) will be replaced with this * value. The default is 1. */ diff --git a/ICCDEF/IccdefPreprocessor.h b/ICCDEF/IccdefPreprocessor.h index 7252a418d5..debfa6b327 100644 --- a/ICCDEF/IccdefPreprocessor.h +++ b/ICCDEF/IccdefPreprocessor.h @@ -138,25 +138,6 @@ class ITK_EXPORT IccdefPreprocessor : public Object itkSetStringMacro(MovingBinaryVolume ); itkGetStringMacro( MovingBinaryVolume ); -#if 0 - /** Set/Get the lower threshold. The default is 0. */ - itkSetMacro(Lower, PixelType); - itkGetMacro(Lower, PixelType); - - /** Set/Get the upper threshold. The default is 70 */ - itkSetMacro(Upper, PixelType); - itkGetMacro(Upper, PixelType); - - /** Set the radius of the neighborhood used for a mask. */ - itkSetMacro(Radius, SizeType); - /** Get the radius of the neighborhood used to compute the median */ - itkGetConstReferenceMacro(Radius, SizeType); - - /** Set the Seed of the neighborhood used for a mask. */ - itkSetMacro(Seed, IndexType); - /** Get the radius of the neighborhood used to compute the median */ - itkGetConstReferenceMacro(Seed, IndexType); -#endif itkSetMacro(DefaultPixelValue, PixelType); itkGetMacro(DefaultPixelValue, PixelType); diff --git a/ICCDEF/IccdefPreprocessor.hxx b/ICCDEF/IccdefPreprocessor.hxx index 3200448968..7e1deed59c 100644 --- a/ICCDEF/IccdefPreprocessor.hxx +++ b/ICCDEF/IccdefPreprocessor.hxx @@ -272,7 +272,6 @@ IccdefPreprocessor m_InputMovingImage = medianFilter->GetOutput(); } -#if 1 // Create UnNormalized...Images { this->m_UnNormalizedFixedImage @@ -284,7 +283,6 @@ IccdefPreprocessor = itkUtil::PreserveCast( this->m_InputMovingImage); } -#endif m_OutputMovingImage = itkUtil::CopyImage( m_UnNormalizedMovingImage); @@ -299,19 +297,6 @@ IccdefPreprocessor { std::cout << "Performing Histogram Matching \n"; } -#if 0 - typedef MinimumMaximumImageCalculator MinMaxFilterType; - typename MinMaxFilterType::Pointer minmaxfilter = MinMaxFilterType::New(); - minmaxfilter->SetImage(m_UnNormalizedFixedImage); - minmaxfilter->ComputeMaximum(); - minmaxfilter->ComputeMinimum(); - if( ( minmaxfilter->GetMaximum() - minmaxfilter->GetMinimum() ) < - m_NumberOfHistogramLevels ) - { - std::cout << "The intensity of range is less than Histogram levels!!" - << std::endl; - } -#else if( ( vcl_numeric_limits::max() - vcl_numeric_limits::min() ) < m_NumberOfHistogramLevels ) @@ -319,7 +304,6 @@ IccdefPreprocessor std::cout << "The intensity of range is less than Histogram levels!!" << std::endl; } -#endif histogramfilter->SetInput( m_UnNormalizedMovingImage ); histogramfilter->SetReferenceImage( m_UnNormalizedFixedImage); @@ -360,34 +344,6 @@ typename IccdefPreprocessor(MaskName); -#if 0 - if( ( m_UnNormalizedFixedImage->GetLargestPossibleRegion().GetSize() != - Mask->GetLargestPossibleRegion().GetSize() ) - || ( m_UnNormalizedFixedImage->GetSpacing() != Mask->GetSpacing() ) ) - { -#if 0 // DEBUG: No need to resample if everyting is done in physical - // coordinates. - Mask = itkUtil::ResampleImage( - Mask, - m_UnNormalizedFixedImage - ->GetLargestPossibleRegion().GetSize(), - m_UnNormalizedFixedImage - ->GetSpacing(), - m_UnNormalizedFixedImage - ->GetOrigin(), - m_UnNormalizedFixedImage - ->GetDirection(), - this-> - m_DefaultPixelValue - ); -#endif - if( this->GetOutDebug() ) - { - std::cout << "Writing Resampled Output image" << std::endl; - itkUtil::WriteImage(Mask, "Resampled.mask"); - } - } -#endif typedef BOBFFilter BOBFFilterType; typename BOBFFilterType::Pointer BOBFfilter = BOBFFilterType::New(); if( this->GetOutDebug() ) @@ -397,19 +353,7 @@ typename IccdefPreprocessorSetLower(m_Lower ); - BOBFfilter->SetUpper( m_Upper ); - BOBFfilter->SetRadius( m_Radius ); - BOBFfilter->SetSeed( m_Seed ); -#endif BOBFfilter->SetReplaceValue( (InputPixelType)m_DefaultPixelValue ); BOBFfilter->SetInputImage( input ); diff --git a/ICCDEF/IccdefRegistrator.hxx b/ICCDEF/IccdefRegistrator.hxx index 1c95c49cba..1391a72413 100644 --- a/ICCDEF/IccdefRegistrator.hxx +++ b/ICCDEF/IccdefRegistrator.hxx @@ -71,13 +71,6 @@ void IccdefRegistrator(DisplacementComponentImagePtr, CurrentComponentFilename); -#if 0 - typedef itk::ImageFileWriter FileWriterType; - FileWriterType::Pointer DisplacementImageWriter = FileWriterType::New(); - DisplacementImageWriter->SetInput(DisplacementComponentImagePtr); - DisplacementImageWriter->SetFileName( CurrentComponentFilename.c_str() ); - DisplacementImageWriter->Update(); -#endif } } catch( itk::ExceptionObject & e ) @@ -175,16 +168,6 @@ template < class TFieldValue> void IccdefRegistrator::Execute() { -#if 0 - // Setup the image pyramids - this->m_FixedImagePyramid->SetNumberOfLevels(this->m_NumberOfLevels); - this->m_FixedImagePyramid->SetStartingShrinkFactors( - this->m_FixedImageShrinkFactors.GetDataPointer() ); - - this->m_MovingImagePyramid->SetNumberOfLevels(this->m_NumberOfLevels); - this->m_MovingImagePyramid-> - SetStartingShrinkFactors( this->m_MovingImageShrinkFactors.GetDataPointer() ); -#endif this->m_Registration->SetFixedImage(this->m_FixedImage); this->m_Registration->SetMovingImage(this->m_MovingImage); @@ -240,14 +223,6 @@ void IccdefRegistrator::Execute() this->m_Registration->SetInitialMovingDisplacementField(fieldReader->GetOutput() ); } -#if 0 - if( this->m_InitialDisplacementField.IsNotNull() ) - { -// this->m_Registration->SetInitialDisplacementField(this->m_InitialDisplacementField); - this->m_Registration->SetInitialMovingDisplacementField(this->m_InitialDisplacementField); - } -#endif - // Perform the registration. try { diff --git a/ICCDEF/commandIterationupdate.h b/ICCDEF/commandIterationupdate.h index c71c86e567..5b57657ce5 100644 --- a/ICCDEF/commandIterationupdate.h +++ b/ICCDEF/commandIterationupdate.h @@ -392,29 +392,6 @@ class CommandIterationUpdate : public itk::Command std::cout << "d(.,true) " << fieldDist << " - "; std::cout << "d(.,Jac(true)) " << fieldGradDist << " - "; } -#if 0 -#if defined( USE_DEBUG_IMAGE_VIEWER ) - if( DebugImageDisplaySender.Enabled() ) - { - DebugImageDisplaySender.SendImage( backdeffield, 0, 0); - DebugImageDisplaySender.SendImage( backdeffield, 1, 1); - DebugImageDisplaySender.SendImage( backdeffield, 2, 2); - typedef typename itk::WarpImageFilter WarpFilterType; - typename WarpFilterType::Pointer back_warper = WarpFilterType::New(); - back_warper->SetInput(m_FixedImage); - back_warper->SetOutputSpacing( backdeffield->GetSpacing() ); - back_warper->SetOutputOrigin( backdeffield->GetOrigin() ); - back_warper->SetOutputDirection( backdeffield->GetDirection() ); - back_warper->SetDisplacementField( backdeffield); - back_warper->Update(); - typename InternalImageType::Pointer - DeformedFixedImagePtr = back_warper->GetOutput(); - DebugImageDisplaySender.SendImage(DeformedFixedImagePtr, 3); - // std::cerr << std::endl << "************IMAGES SENT*************" << std::endl; - } -#endif // defined(USE_DEBUG_IMAGE_VIEWER) -#endif m_HarmonicEnergyCalculator->SetImage( backdeffield ); m_HarmonicEnergyCalculator->Compute(); const double backharmonicEnergy diff --git a/ICCDEF/iccdefRegistration_New.cxx b/ICCDEF/iccdefRegistration_New.cxx index 6da209704d..b0901d9e7d 100644 --- a/ICCDEF/iccdefRegistration_New.cxx +++ b/ICCDEF/iccdefRegistration_New.cxx @@ -236,20 +236,6 @@ void ThirionFunction(const struct ICCDEFWarpAppParameters & command) { // Read Landmark spatial objects // what type is the landmark file? -#if 0 - typedef itk::SpatialObjectReader SpatialObjectReaderType; - typedef itk::LandmarkSpatialObject LandmarkSpatialObjectType; - typename SpatialObjectReaderType::Pointer fixedLandmarkReader = SpatialObjectReaderType::New(); - fixedLandmarkReader->SetFileName(command.fixedLandmark.c_str() ); - fixedLandmarkReader->Update(); - filter->SetFixedLandmark(dynamic_cast(fixedLandmarkReader->GetGroup().GetPointer() ) ); - - typename SpatialObjectReaderType::Pointer movingLandmarkReader = SpatialObjectReaderType::New(); - movingLandmarkReader->SetFileName(command.movingLandmark.c_str() ); - movingLandmarkReader->Update(); - filter->SetMovingLandmark( - dynamic_cast(movingLandmarkReader->GetGroup().GetPointer() ) ); -#endif // Transform index point to physical point typename TRealImage::Pointer fixedVolume @@ -282,37 +268,6 @@ void ThirionFunction(const struct ICCDEFWarpAppParameters & command) exit(-1); } -#if 0 - if( command.maskProcessingMode == "ROI" ) - { - if( ( command.fixedBinaryVolume == "" ) - || ( command.movingBinaryVolume == "" ) ) - { - std::cout - << "ERROR: Must specify mask file names when ROI is used for the maskProcessingMode" << std::endl; - exit(-1); - } - std::cout << "Registration with Mask!!!!!!!" << std::endl; - typename TRealImage::Pointer fixedVolume - = itkUtil::ReadImageRAI( command.fixedVolume.c_str() ); - typename TRealImage::Pointer movingVolume - = itkUtil::ReadImageRAI( command.movingVolume.c_str() ); - - fixedMask = ReadImageMask( - command.fixedBinaryVolume, - fixedVolume->GetDirection(), - fixedVolume->GetOrigin() ); - movingMask = ReadImageMask( - command.movingBinaryVolume, - movingVolume->GetDirection(), - movingVolume->GetOrigin() ); - - filter->SetFixedImageMask(fixedMask); - filter->SetMovingImageMask(movingMask); - filter->SetBackgroundFilledValue(command.backgroundFillValue); - } - -#endif typename CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New(); @@ -324,14 +279,6 @@ void ThirionFunction(const struct ICCDEFWarpAppParameters & command) app->SetRegistrationFilter(filter); -#if 0 - if( command.initialDisplacementFieldVolume != "" ) - { - app->SetInitialDisplacementFieldFilename( - command.initialDisplacementFieldVolume.c_str() ); - } -#endif - if( command.initialFixedDisplacementFieldVolume != "" ) { app->SetInitialFixedDisplacementFieldFilename( diff --git a/ICCDEF/itkBrains2LandmarkReader.hxx b/ICCDEF/itkBrains2LandmarkReader.hxx index 4acc4da670..bf035c97ac 100644 --- a/ICCDEF/itkBrains2LandmarkReader.hxx +++ b/ICCDEF/itkBrains2LandmarkReader.hxx @@ -23,13 +23,6 @@ Brains2LandmarkReader m_ReferenceImage = ImageType::New(); } -#if 0 -Brains2LandmarkReader -::~Brains2LandmarkReader() -{ -} - -#endif template void Brains2LandmarkReader ::Update() diff --git a/ICCDEF/itkICCDeformableFunction.h b/ICCDEF/itkICCDeformableFunction.h index c06ddef0c0..7885b98660 100644 --- a/ICCDEF/itkICCDeformableFunction.h +++ b/ICCDEF/itkICCDeformableFunction.h @@ -311,18 +311,6 @@ class ITK_EXPORT ICCDeformableFunction : return this->m_LandmarkWeight; } -#if 0 - virtual void SetChi(double sm) - { - this->m_Chi = sm; - } - - virtual double GetChi() const - { - return this->m_Chi; - } - -#endif virtual void SetSmoothFilter(ComplexImagePointer& filter) { m_SmoothFilter = filter; @@ -363,28 +351,6 @@ class ITK_EXPORT ICCDeformableFunction : m_MovingLandmark = lk; } -#if 0 - virtual void SetFixedImageBackground(float bk) - { - m_FixedImageBackground = bk; - } - - virtual float GetFixedImageBackground() const - { - return m_FixedImageBackground; - } - - virtual void SetMovingImageBackground(float bk) - { - m_MovingImageBackground = bk; - } - - virtual float GetMovingImageBackground() const - { - return m_MovingImageBackground; - } - -#endif protected: ICCDeformableFunction(); ~ICCDeformableFunction() diff --git a/ICCDEF/itkICCDeformableFunction.hxx b/ICCDEF/itkICCDeformableFunction.hxx index cbd38530e3..13be8e9cea 100644 --- a/ICCDEF/itkICCDeformableFunction.hxx +++ b/ICCDEF/itkICCDeformableFunction.hxx @@ -165,14 +165,12 @@ void ICCDeformableFunction ::InitializeIteration() { -#if 1 if( !this->GetMovingImage() || !this->GetFixedImage() || !m_MovingImageInterpolator ) { itkExceptionMacro( << "MovingImage, FixedImage and/or Interpolator not set" ); } -#endif m_SumOfSquaredDifference = 0.0; m_NumberOfPixelsProcessed = 0L; @@ -224,41 +222,7 @@ ICCDeformableFunction { typename DerivativeType::Pointer derivative = DerivativeType::New(); -#if 0 - if( this->GetMovingImageMask() && this->GetFixedImageMask() ) - { - typename MovingImageType::Pointer maskedMovingImage = MovingImageType::New(); - maskedMovingImage->SetRegions(m_MovingImageWarper->GetOutput()->GetLargestPossibleRegion() ); - maskedMovingImage->Allocate(); - maskedMovingImage->SetSpacing(m_MovingImageWarper->GetOutput()->GetSpacing() ); - maskedMovingImage->SetOrigin(m_MovingImageWarper->GetOutput()->GetOrigin() ); - maskedMovingImage->SetDirection(m_MovingImageWarper->GetOutput()->GetDirection() ); - - typedef ImageRegionIterator IterationMaskImageType; - IterationImageType maskIter(maskedMovingImage, maskedMovingImage->GetRequestedRegion() ); - IterationImageType mIt(m_MovingImageWarper->GetOutput(), - m_MovingImageWarper->GetOutput()->GetRequestedRegion() ); - IterationMaskImageType dit(m_MovingMaskImageWarper->GetOutput(), - m_MovingMaskImageWarper->GetOutput()->GetRequestedRegion() ); - for( maskIter.GoToBegin(), mIt.GoToBegin(), dit.GoToBegin(); !maskIter.IsAtEnd(); ++maskIter, ++mIt, ++dit ) - { - if( dit.Get() > 0 ) - { - maskIter.Set(mIt.Get() ); - } - else - { - maskIter.Set(m_BackgroundFilledValue); - } - } - // itkUtil::WriteImage(maskedMovingImage,"maskedImage.nii.gz"); - derivative->SetInput(maskedMovingImage); - } - else -#endif - { - derivative->SetInput(m_MovingImageWarper->GetOutput() ); - } + derivative->SetInput(m_MovingImageWarper->GetOutput() ); // compute the derivative of Warped image derivative->SetDirection(0); derivative->Update(); @@ -302,92 +266,17 @@ ICCDeformableFunction IterationImageType zIter(tempZ, tempZ->GetRequestedRegion() ); typename TDisplacementField::PixelType pixel; -#if 0 - if( this->GetMovingImageMask() && this->GetFixedImageMask() ) - { - for( def12Iter.GoToBegin(), + for( def12Iter.GoToBegin(), fit.GoToBegin(), mit.GoToBegin(), xIter.GoToBegin(), yIter.GoToBegin(), zIter.GoToBegin(); - !def12Iter.IsAtEnd(); ++def12Iter, ++fit, ++mit, ++xIter, ++yIter, ++zIter ) - { - typename DisplacementFieldType::IndexType index = def12Iter.GetIndex(); - typename FixedImageType::PointType fixedPoint; - this->GetFixedImage()->TransformIndexToPhysicalPoint(index, fixedPoint); - if( this->GetFixedImageMask()->IsInside(fixedPoint) && - (m_MovingMaskImageWarper->GetOutput()->GetPixel(index) > 0) ) - { - pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * xIter.Get(); - pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * yIter.Get(); - pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * zIter.Get(); - } - else if( !this->GetFixedImageMask()->IsInside(fixedPoint) && - (m_MovingMaskImageWarper->GetOutput()->GetPixel(index) > 0) ) - { - pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (mit.Get() - m_BackgroundFilledValue) * xIter.Get(); - pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (mit.Get() - m_BackgroundFilledValue) * yIter.Get(); - pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (mit.Get() - m_BackgroundFilledValue) * zIter.Get(); - } - else if( this->GetFixedImageMask()->IsInside(fixedPoint) && - (m_MovingMaskImageWarper->GetOutput()->GetPixel(index) <= 0) ) - { - pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_BackgroundFilledValue - fit.Get() ) * xIter.Get(); - pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_BackgroundFilledValue - fit.Get() ) * yIter.Get(); - pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_BackgroundFilledValue - fit.Get() ) * zIter.Get(); - } - else - { - pixel.Fill(0.0); -#if 0 - pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_MovingImageBackground - m_FixedImageBackground) * xIter.Get(); - pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_MovingImageBackground - m_FixedImageBackground) * yIter.Get(); - pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight - * (m_MovingImageBackground - m_FixedImageBackground) * zIter.Get(); -#endif - } - def12Iter.Set(pixel); - } - } - else -#endif + !def12Iter.IsAtEnd(); ++def12Iter, ++fit, ++mit, ++xIter, ++yIter, ++zIter ) { - for( def12Iter.GoToBegin(), - fit.GoToBegin(), mit.GoToBegin(), xIter.GoToBegin(), yIter.GoToBegin(), zIter.GoToBegin(); - !def12Iter.IsAtEnd(); ++def12Iter, ++fit, ++mit, ++xIter, ++yIter, ++zIter ) - { - pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * xIter.Get(); - pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * yIter.Get(); - pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * zIter.Get(); - def12Iter.Set(pixel); - } + pixel[0] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * xIter.Get(); + pixel[1] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * yIter.Get(); + pixel[2] = 4.0 * m_MaximumUpdateStepLength * m_SimilarityWeight * (mit.Get() - fit.Get() ) * zIter.Get(); + def12Iter.Set(pixel); } similarity->Modified(); // Write the similarity measure image -#if 0 - typedef Image RealImageType; - typedef VectorIndexSelectionCastImageFilter SelectImageType; - typename SelectImageType::Pointer caster = SelectImageType::New(); - caster->SetInput(similarity); - for( unsigned int i = 0; i < 3; i++ ) - { - caster->SetIndex(i); - caster->Update(); - std::string name = "Similarity"; - std::stringstream out; - out << i; - std::string name1 = out.str(); - name = name + name1 + ".nii.gz"; - typename FixedImageType::Pointer simi = FixedImageType::New(); - simi = caster->GetOutput(); - itkUtil::WriteImage(simi, name); - } -#endif /* typename DisplacementFieldFFTType::Pointer coeff = DisplacementFieldFFTType::New(); typename FFTWRealToComplexType::Pointer fft12 = FFTWRealToComplexType::New(); @@ -479,11 +368,6 @@ ICCDeformableFunction this->GetFixedImage()->TransformPhysicalPointToIndex(warpedMovingPoint, pixelIndex1); this->GetFixedImage()->TransformIndexToPhysicalPoint(pixelIndex1, warpedMovingPoint); warpedMovingLandmark->SetPoint(pointID, warpedMovingPoint); -#if 0 - std::cout << "Deformation field: " << this->GetDisplacementField()->GetPixel(pixelIndex) << std::endl; - std::cout << "fixed point: " << fixedPoint << std::endl; - std::cout << "warped moving point: " << warpedMovingPoint << std::endl; -#endif // std::cout<<"difference: "<SetPoint(pointID, movingPoint); ++mlit; @@ -491,19 +375,6 @@ ICCDeformableFunction pointID++; } -#if 0 - // typename KernelPointSetType::PointsContainer::ConstIterator pit; - pit = m_FixedLandmark->GetPoints()->Begin(); - pointID = 0; - while( pit != m_FixedLandmark->GetPoints()->End() ) - { - // std::cout<<"fixed Landmark Position: "<SetPoint(pointID, pit.Value() ); - ++pit; - pointID++; - } - -#endif // warpedMovingLandmark->SetPoints(list); // Wondering to use physical space point or index point?????? // Using kernel transform @@ -533,22 +404,6 @@ ICCDeformableFunction defIter.Set(def); } this->GetDisplacementField()->Modified(); -#if 0 - mlit = m_MovingLandmark->GetPoints()->Begin(); - flit = m_FixedLandmark->GetPoints()->Begin(); - while( mlit != m_MovingLandmark->GetPoints()->End() ) - { - this->GetMovingImage()->TransformPhysicalPointToIndex(mlit.Value(), pixelIndex); - this->GetFixedImage()->TransformIndexToPhysicalPoint(pixelIndex, movingPoint); - point = transform->TransformPoint(movingPoint); - std::cout << "moving point: " << movingPoint << std::endl; - std::cout << "warped moving point: " << point << std::endl; - std::cout << "displacement: " << point - movingPoint << std::endl; - ++mlit; - ++flit; - } - -#endif typename FFTWRealToComplexType::Pointer fft_filter = FFTWRealToComplexType::New(); fft_filter->SetInput(this->GetDisplacementField() ); diff --git a/ICCDEF/itkICCDeformableRegistrationFilter.h b/ICCDEF/itkICCDeformableRegistrationFilter.h index e5e899ff5c..8a742004a9 100644 --- a/ICCDEF/itkICCDeformableRegistrationFilter.h +++ b/ICCDEF/itkICCDeformableRegistrationFilter.h @@ -352,17 +352,6 @@ class ITK_EXPORT ICCDeformableRegistrationFilter : ComplexImagePointer m_sqr11, m_sqr12, m_sqr13, m_sqr22, m_sqr23, m_sqr33; ComplexImagePointer m_SmoothFilter; -#if 0 - MaskImagePointer m_WarpedFixedMask; - MaskImagePointer m_WarpedMovingMask; - - MaskPointer m_FixedMask; - MaskPointer m_MovingMask; - - LandmarkPointer m_FixedLandmark; - LandmarkPointer m_MovingLandmark; -#endif - float m_Alpha, m_Gamma, m_Beta; float m_RegularizationWeight; float m_InverseWeight; diff --git a/ICCDEF/itkICCDeformableRegistrationFilter.hxx b/ICCDEF/itkICCDeformableRegistrationFilter.hxx index 01e7def03e..0404165d2e 100644 --- a/ICCDEF/itkICCDeformableRegistrationFilter.hxx +++ b/ICCDEF/itkICCDeformableRegistrationFilter.hxx @@ -67,14 +67,6 @@ ICCDeformableRegistrationFilter m_sqr33 = ComplexImageType::New(); m_SmoothFilter = ComplexImageType::New(); -#if 0 - m_WarpedFixedMask = MaskImageType::New(); - m_WarpedMovingMask = MaskImageType::New(); - - m_FixedMask = MaskType::New(); - m_MovingMask = MaskType::New(); -#endif - // m_FixedLandmark = LandmarkType::New(); // m_MovingLandmark = LandmarkType::New(); m_Alpha = 0.1; @@ -318,23 +310,6 @@ ICCDeformableRegistrationFilter // Set different landmark matching weight at different resolusion f->SetLandmarkWeight(this->GetLandmarkWeight() ); b->SetLandmarkWeight(this->GetLandmarkWeight() ); -#if 0 - if( this->GetOutput(0)->GetLargestPossibleRegion().GetSize()[0] > 32 ) - { - f->SetSigma(this->GetSimilarityWeight() / 2.0); - b->SetSigma(this->GetSimilarityWeight() / 2.0); - } - else if( this->GetOutput(0)->GetLargestPossibleRegion().GetSize()[0] > 64 ) - { - f->SetSigma(this->GetSimilarityWeight() / 4.0); - b->SetSigma(this->GetSimilarityWeight() / 4.0); - } - else - { - f->SetSigma(this->GetSimilarityWeight() ); - b->SetSigma(this->GetSimilarityWeight() ); - } -#endif } // Compute D[k], Initialize harmonic @@ -710,35 +685,6 @@ ICCDeformableRegistrationFilter return drfp->GetRMSChange(); } -#if 0 -template -typename ICCDeformableRegistrationFilter -::GradientType -ICCDeformableRegistrationFilter -::GetUseGradientType() const -{ - const ICCDeformableFunctionType *drfp = this->GetForwardRegistrationFunctionType(); - - return drfp->GetUseGradientType(); -} - -/** - * - */ -template -void -ICCDeformableRegistrationFilter -::SetUseGradientType(GradientType gtype) -{ - ICCDeformableFunctionType *drfpf = this->GetForwardRegistrationFunctionType(); - - drfpf->SetUseGradientType(gtype); - ICCDeformableFunctionType *drfpb = this->GetBackwardRegistrationFunctionType(); - drfpb->SetUseGradientType(gtype); -} - -#endif - /** * */ @@ -804,12 +750,6 @@ ICCDeformableRegistrationFilter { // If we smooth the update buffer before applying it, then the are // approximating a viscuous problem as opposed to an elastic problem -#if 0 - if( this->GetSmoothUpdateField() ) - { - this->SmoothUpdateField(); - } -#endif // Compute inverse consistency typedef SubtractImageFilter this->SetRMSChange( drfp->GetRMSChange() ); - /** - * Smooth the deformation field - */ -#if 0 - if( this->GetSmoothDisplacementField() ) - { - this->SmoothDisplacementField(); - } -#endif } template diff --git a/ICCDEF/itkMultiResolutionICCDeformableRegistration.hxx b/ICCDEF/itkMultiResolutionICCDeformableRegistration.hxx index 313eb1a646..e4c1b2286f 100644 --- a/ICCDEF/itkMultiResolutionICCDeformableRegistration.hxx +++ b/ICCDEF/itkMultiResolutionICCDeformableRegistration.hxx @@ -337,17 +337,6 @@ MultiResolutionICCDeformableRegistrationDisconnectPipeline(); // tempField12->Print(std::cout,6); -#if 0 - typedef WarpImageFilter WarpImageType; - typename WarpImageType::Pointer warper = WarpImageType::New(); - warper->SetOutputOrigin( m_FixedImagePyramid->GetOutput(fixedLevel)->GetOrigin() ); - warper->SetOutputSpacing( m_FixedImagePyramid->GetOutput(fixedLevel)->GetSpacing() ); - warper->SetOutputDirection(m_FixedImagePyramid->GetOutput(fixedLevel)->GetDirection() ); - warper->SetInput( m_MovingImagePyramid->GetOutput(fixedLevel) ); - warper->SetDisplacementField(tempField12 ); - warper->GetOutput()->SetRequestedRegion( tempField12->GetRequestedRegion() ); - warper->Update(); -#endif if( m_DisplacementFieldOutputNamePrefix != "" && m_DisplacementFieldOutputNamePrefix != "none" ) { std::string name = m_DisplacementFieldOutputNamePrefix + "_iteration"; // resolution";