diff --git a/.clang-format b/.clang-format index 7566ce7..7252d3e 100644 --- a/.clang-format +++ b/.clang-format @@ -1,28 +1,111 @@ +## This config file is only relevant for clang-format version 19.1.7 +## +## Examples of each format style can be found on the in the clang-format documentation +## See: https://clang.llvm.org/docs/ClangFormatStyleOptions.html for details of each option +## +## The clang-format binaries can be downloaded as part of the clang binary distributions +## from https://releases.llvm.org/download.html +## +## Use the script Utilities/Maintenance/clang-format.bash to faciliate +## maintaining a consistent code style. +## +## EXAMPLE apply code style enforcement before commit: +# Utilities/Maintenance/clang-format.bash --clang ${PATH_TO_CLANG_FORMAT_19.1.7} --modified +## EXAMPLE apply code style enforcement after commit: +# Utilities/Maintenance/clang-format.bash --clang ${PATH_TO_CLANG_FORMAT_19.1.7} --last --- +# This configuration requires clang-format version 19.1.7 exactly. Language: Cpp -BasedOnStyle: Mozilla AccessModifierOffset: -2 AlignAfterOpenBracket: Align -AlignConsecutiveAssignments: false -AlignConsecutiveDeclarations: true -AlignEscapedNewlinesLeft: true -AlignOperands: true -AlignTrailingComments: true +AlignArrayOfStructures: None +AlignConsecutiveAssignments: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: true +AlignConsecutiveBitFields: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveDeclarations: + Enabled: true + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: true +AlignConsecutiveMacros: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveShortCaseStatements: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCaseArrows: false + AlignCaseColons: false +AlignConsecutiveTableGenBreakingDAGArgColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveTableGenCondOperatorColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveTableGenDefinitionColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignEscapedNewlines: Left +AlignOperands: Align +AlignTrailingComments: + Kind: Always + OverEmptyLines: 0 +AllowAllArgumentsOnNextLine: true AllowAllParametersOfDeclarationOnNextLine: false -AllowShortBlocksOnASingleLine: false +AllowBreakBeforeNoexceptSpecifier: Never +AllowShortBlocksOnASingleLine: Never +AllowShortCaseExpressionOnASingleLine: true AllowShortCaseLabelsOnASingleLine: false -AllowShortFunctionsOnASingleLine: Inline -AllowShortIfStatementsOnASingleLine: false +AllowShortCompoundRequirementOnASingleLine: true +AllowShortEnumsOnASingleLine: true +#AllowShortFunctionsOnASingleLine: Inline Only merge functions defined inside a class. Implies empty. +#AllowShortFunctionsOnASingleLine: None (in configuration: None) Never merge functions into a single line. +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Never +AllowShortLambdasOnASingleLine: All AllowShortLoopsOnASingleLine: false -AlwaysBreakAfterReturnType: All -AlwaysBreakBeforeMultilineStrings: true -AlwaysBreakTemplateDeclarations: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AttributeMacros: + - __capability BinPackArguments: false -BinPackParameters: true +BinPackParameters: false +BitFieldColonSpacing: Both BraceWrapping: + AfterCaseLabel: true AfterClass: true - AfterControlStatement: true + AfterControlStatement: Always AfterEnum: true + AfterExternBlock: true AfterFunction: true AfterNamespace: true AfterObjCDeclaration: true @@ -30,65 +113,206 @@ BraceWrapping: AfterUnion: true BeforeCatch: true BeforeElse: true - IndentBraces: true + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: false + SplitEmptyRecord: false + SplitEmptyNamespace: false +BreakAdjacentStringLiterals: true +BreakAfterAttributes: Leave +BreakAfterJavaFieldAnnotations: false +BreakAfterReturnType: All +BreakArrays: true BreakBeforeBinaryOperators: None -BreakBeforeBraces: GNU +BreakBeforeConceptDeclarations: Always +BreakBeforeBraces: Custom +BreakBeforeInlineASMColon: OnlyMultiline BreakBeforeTernaryOperators: true -BreakConstructorInitializersBeforeComma: true -BreakAfterJavaFieldAnnotations: false +BreakConstructorInitializers: BeforeComma +BreakFunctionDefinitionParameters: false +BreakInheritanceList: BeforeComma BreakStringLiterals: true -ColumnLimit: 200 -# ColumnLimit: 80 +BreakTemplateDeclarations: Yes +## The following line allows larger lines in non-documentation code +ColumnLimit: 120 CommentPragmas: '^ IWYU pragma:' -ConstructorInitializerAllOnOneLineOrOnePerLine: false +CompactNamespaces: false ConstructorInitializerIndentWidth: 2 ContinuationIndentWidth: 2 Cpp11BracedListStyle: false DerivePointerAlignment: false DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock ExperimentalAutoDetectBinPacking: false -ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ] +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve IncludeCategories: - Regex: '^"(llvm|llvm-c|clang|clang-c)/' Priority: 2 - - Regex: '^(<|"(gtest|isl|json)/)' + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' Priority: 3 + SortPriority: 0 + CaseSensitive: false - Regex: '.*' Priority: 1 -IncludeIsMainRegex: '$' + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseBlocks: false IndentCaseLabels: true +IndentExternBlock: AfterExternBlock +IndentGotoLabels: true +IndentPPDirectives: AfterHash +IndentRequiresClause: true IndentWidth: 2 IndentWrappedFunctionNames: false +InsertBraces: false +InsertNewlineAtEOF: false +InsertTrailingCommas: None +IntegerLiteralSeparator: + Binary: 0 + BinaryMinDigits: 0 + Decimal: 0 + DecimalMinDigits: 0 + Hex: 0 + HexMinDigits: 0 JavaScriptQuotes: Leave JavaScriptWrapImports: true -KeepEmptyLinesAtTheStartOfBlocks: false +KeepEmptyLines: + AtEndOfFile: false + AtStartOfBlock: true + AtStartOfFile: true +LambdaBodyIndentation: Signature +LineEnding: DeriveLF MacroBlockBegin: '' MacroBlockEnd: '' +MainIncludeChar: Quote MaxEmptyLinesToKeep: 2 NamespaceIndentation: None +ObjCBinPackProtocolList: Auto ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true ObjCSpaceAfterProperty: true ObjCSpaceBeforeProtocolList: false +PackConstructorInitializers: BinPack +PenaltyBreakAssignment: 2 PenaltyBreakBeforeFirstCallParameter: 19 PenaltyBreakComment: 300 +## The following line allows larger lines in non-documentation code PenaltyBreakFirstLessLess: 120 +PenaltyBreakOpenParenthesis: 0 +PenaltyBreakScopeResolution: 500 PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 PenaltyExcessCharacter: 1000000 -PenaltyReturnTypeOnItsOwnLine: 0 -PointerAlignment: Left +PenaltyIndentedWhitespace: 0 +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Middle +PPIndentWidth: -1 +QualifierAlignment: Custom +QualifierOrder: + - friend + - static + - inline + - constexpr + - const + - type +ReferenceAlignment: Pointer ReflowComments: true -SortIncludes: true +RemoveBracesLLVM: false +RemoveParentheses: Leave +RemoveSemicolon: false +RequiresClausePosition: OwnLine +RequiresExpressionIndentation: OuterScope +SeparateDefinitionBlocks: Leave +ShortNamespaceLines: 1 +SkipMacroDefinitionBody: false +# We may want to sort the includes as a separate pass +SortIncludes: Never +SortJavaStaticImport: Before +# We may want to revisit this later +SortUsingDeclarations: Never SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceAroundPointerQualifiers: Default SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeJsonColon: false SpaceBeforeParens: ControlStatements -SpaceInEmptyParentheses: false +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterForeachMacros: true + AfterFunctionDefinitionName: false + AfterFunctionDeclarationName: false + AfterIfMacros: true + AfterOverloadedOperator: false + AfterPlacementOperator: true + AfterRequiresInClause: false + AfterRequiresInExpression: false + BeforeNonEmptyParentheses: false +SpaceBeforeRangeBasedForLoopColon: true +SpaceBeforeSquareBrackets: false +SpaceInEmptyBlock: false SpacesBeforeTrailingComments: 1 -SpacesInAngles: true -SpacesInCStyleCastParentheses: false -SpacesInContainerLiterals: true -SpacesInParentheses: true +SpacesInAngles: Never +SpacesInContainerLiterals: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParens: Never +SpacesInParensOptions: + ExceptDoubleParentheses: false + InCStyleCasts: false + InConditionalStatements: false + InEmptyParentheses: false + Other: false SpacesInSquareBrackets: false -Standard: Cpp03 +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION + - ITK_GCC_PRAGMA_PUSH + - ITK_GCC_PRAGMA_POP + - ITK_GCC_SUPPRESS_Wfloat_equal + - ITK_GCC_SUPPRESS_Wformat_nonliteral + - ITK_GCC_SUPPRESS_Warray_bounds + - ITK_CLANG_PRAGMA_PUSH + - ITK_CLANG_PRAGMA_POP + - ITK_CLANG_SUPPRESS_Wzero_as_null_pointer_constant + - CLANG_PRAGMA_PUSH + - CLANG_PRAGMA_POP + - CLANG_SUPPRESS_Wfloat_equal + - INTEL_PRAGMA_WARN_PUSH + - INTEL_PRAGMA_WARN_POP + - INTEL_SUPPRESS_warning_1292 + - itkTemplateFloatingToIntegerMacro + - itkLegacyMacro +TableGenBreakInsideDAGArg: DontBreak TabWidth: 2 UseTab: Never +VerilogBreakBetweenInstancePorts: true +WhitespaceSensitiveMacros: + - BOOST_PP_STRINGIZE + - CF_SWIFT_NAME + - NS_SWIFT_NAME + - PP_STRINGIZE + - STRINGIZE ... diff --git a/.github/workflows/clang-format-linter.yml b/.github/workflows/clang-format-linter.yml index 69166d9..c8e1681 100644 --- a/.github/workflows/clang-format-linter.yml +++ b/.github/workflows/clang-format-linter.yml @@ -7,7 +7,9 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v4 with: fetch-depth: 1 - uses: InsightSoftwareConsortium/ITKClangFormatLinterAction@master + with: + itk-branch: master diff --git a/examples/mciExample.cxx b/examples/mciExample.cxx index 334c00e..7312b82 100644 --- a/examples/mciExample.cxx +++ b/examples/mciExample.cxx @@ -6,49 +6,49 @@ int -main( int argc, char* argv[] ) +main(int argc, char * argv[]) { - if ( argc < 3 ) - { - std::cout << "Syntax: " << argv[0] << " sliceImage interpolatedImage [smoothingRadius=2]" << std::endl; - return EXIT_FAILURE; - } + if (argc < 3) + { + std::cout << "Syntax: " << argv[0] << " sliceImage interpolatedImage [smoothingRadius=2]" << std::endl; + return EXIT_FAILURE; + } try - { - using ImageType = itk::Image< unsigned char, 3 >; - using ReaderType = itk::ImageFileReader< ImageType >; - ReaderType::Pointer imageReader = ReaderType::New(); - imageReader->SetFileName( argv[1] ); - - using mciType = itk::MorphologicalContourInterpolator< ImageType >; - mciType::Pointer mci = mciType::New(); - mci->SetInput( imageReader->GetOutput() ); - mci->Update(); + { + using ImageType = itk::Image; + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer imageReader = ReaderType::New(); + imageReader->SetFileName(argv[1]); - int smoothingRadius = 2; - if ( argc >= 4 ) - { - smoothingRadius = std::stoi( argv[3] ); - } + using mciType = itk::MorphologicalContourInterpolator; + mciType::Pointer mci = mciType::New(); + mci->SetInput(imageReader->GetOutput()); + mci->Update(); - using MedianType = itk::MedianImageFilter< ImageType, ImageType >; - MedianType::Pointer medF = MedianType::New(); - medF->SetInput( mci->GetOutput() ); - medF->SetRadius( smoothingRadius ); - medF->Update(); - - using WriterType = itk::ImageFileWriter< ImageType >; - WriterType::Pointer writer = WriterType::New(); - writer->SetFileName( argv[2] ); - writer->SetInput( medF->GetOutput() ); - writer->SetUseCompression( true ); - writer->Update(); - } - catch ( itk::ExceptionObject& exc ) + int smoothingRadius = 2; + if (argc >= 4) { - std::cerr << exc << std::endl; - return EXIT_FAILURE; + smoothingRadius = std::stoi(argv[3]); } + + using MedianType = itk::MedianImageFilter; + MedianType::Pointer medF = MedianType::New(); + medF->SetInput(mci->GetOutput()); + medF->SetRadius(smoothingRadius); + medF->Update(); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(argv[2]); + writer->SetInput(medF->GetOutput()); + writer->SetUseCompression(true); + writer->Update(); + } + catch (itk::ExceptionObject & exc) + { + std::cerr << exc << std::endl; + return EXIT_FAILURE; + } return EXIT_SUCCESS; } diff --git a/include/itkMorphologicalContourInterpolator.h b/include/itkMorphologicalContourInterpolator.h index 8696e3c..9ba6882 100644 --- a/include/itkMorphologicalContourInterpolator.h +++ b/include/itkMorphologicalContourInterpolator.h @@ -61,108 +61,108 @@ namespace itk * * \ingroup MorphologicalContourInterpolation */ -template< typename TImage > -class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TImage > +template +class MorphologicalContourInterpolator : public ImageToImageFilter { - template< typename T > + template friend class MorphologicalContourInterpolatorParallelInvoker; public: - ITK_DISALLOW_COPY_AND_MOVE( MorphologicalContourInterpolator ); + ITK_DISALLOW_COPY_AND_MOVE(MorphologicalContourInterpolator); /** Standard class type alias. */ using Self = MorphologicalContourInterpolator; - using Superclass = ImageToImageFilter< TImage, TImage >; - using Pointer = SmartPointer< Self >; - using SliceType = Image< typename TImage::PixelType, TImage::ImageDimension - 1 >; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using SliceType = Image; /** Method for creation through the object factory. */ - itkNewMacro( Self ); + itkNewMacro(Self); /** Run-time type information (and related methods). */ - itkTypeMacro( MorphologicalContourInterpolator, ImageToImageFilter ); + itkOverrideGetNameOfClassMacro(MorphologicalContourInterpolator); /** Interpolate only this label. Interpolates all labels if set to 0 (default). */ - itkSetMacro( Label, typename TImage::PixelType ); + itkSetMacro(Label, typename TImage::PixelType); /** Which label is interpolated. 0 means all labels (default). */ - itkGetMacro( Label, typename TImage::PixelType ); + itkGetMacro(Label, typename TImage::PixelType); /** Which label is interpolated. 0 means all labels (default). */ - itkGetConstMacro( Label, typename TImage::PixelType ); + itkGetConstMacro(Label, typename TImage::PixelType); /** Interpolate only along this axis. Interpolates along all axes if set to -1 (default). */ - itkSetMacro( Axis, int ); + itkSetMacro(Axis, int); /** Axis of interpolation. -1 means interpolation along all axes (default). */ - itkGetMacro( Axis, int ); + itkGetMacro(Axis, int); /** Axis of interpolation. -1 means interpolation along all axes (default). */ - itkGetConstMacro( Axis, int ); + itkGetConstMacro(Axis, int); /** Heuristic alignment of regions for interpolation is faster than optimal alignment. * Default is heuristic. */ - itkSetMacro( HeuristicAlignment, bool ); + itkSetMacro(HeuristicAlignment, bool); /** Heuristic alignment of regions for interpolation is faster than optimal alignment. * Default is heuristic. */ - itkGetMacro( HeuristicAlignment, bool ); + itkGetMacro(HeuristicAlignment, bool); /** Heuristic alignment of regions for interpolation is faster than optimal alignment. * Default is heuristic. */ - itkGetConstMacro( HeuristicAlignment, bool ); + itkGetConstMacro(HeuristicAlignment, bool); /** Using distance transform instead of repeated dilations to calculate * the median contour is slightly faster, but produces lower quality interpolations. * Default is OFF (that is, use repeated dilations). */ - itkSetMacro( UseDistanceTransform, bool ); + itkSetMacro(UseDistanceTransform, bool); /** Using distance transform instead of repeated dilations to calculate * the median contour is slightly faster, but produces lower quality interpolations. * Default is OFF (that is, use repeated dilations). */ - itkGetMacro( UseDistanceTransform, bool ); + itkGetMacro(UseDistanceTransform, bool); /** Using distance transform instead of repeated dilations to calculate * the median contour is slightly faster, but produces lower quality interpolations. * Default is OFF (that is, use repeated dilations). */ - itkGetConstMacro( UseDistanceTransform, bool ); + itkGetConstMacro(UseDistanceTransform, bool); /** Use custom slice positions (not slice auto-detection). * SetLabeledSliceIndices has to be called prior to Update(). */ - itkSetMacro( UseCustomSlicePositions, bool ); + itkSetMacro(UseCustomSlicePositions, bool); /** Use custom slice positions (not slice auto-detection). * SetLabeledSliceIndices has to be called prior to Update(). */ - itkGetMacro( UseCustomSlicePositions, bool ); + itkGetMacro(UseCustomSlicePositions, bool); /** Use custom slice positions (not slice auto-detection). * SetLabeledSliceIndices has to be called prior to Update(). */ - itkGetConstMacro( UseCustomSlicePositions, bool ); + itkGetConstMacro(UseCustomSlicePositions, bool); /** Perform extrapolation for branch extremities. * Branch extremities are defined as regions having no overlap with any region in the next slice. * Extrapolation helps generate smooth surface closings. * Default is ON * */ - itkSetMacro( UseExtrapolation, bool ); + itkSetMacro(UseExtrapolation, bool); /** Use ball instead of default cross structuring element for repeated dilations. */ void - SetUseBallStructuringElement( bool useBall ) + SetUseBallStructuringElement(bool useBall) { - if ( useBall != m_UseBallStructuringElement ) - { - m_UseBallStructuringElement = useBall; - m_ConnectedComponents->SetFullyConnected( useBall ); - this->Modified(); - } + if (useBall != m_UseBallStructuringElement) + { + m_UseBallStructuringElement = useBall; + m_ConnectedComponents->SetFullyConnected(useBall); + this->Modified(); + } } /** Use ball instead of default cross structuring element for repeated dilations. */ - itkGetMacro( UseBallStructuringElement, bool ); + itkGetMacro(UseBallStructuringElement, bool); /** Use ball instead of default cross structuring element for repeated dilations. */ - itkGetConstMacro( UseBallStructuringElement, bool ); + itkGetConstMacro(UseBallStructuringElement, bool); /** If there is a pixel whose all 4-way neighbors belong the the same label except along one axis, and along that axis its neighbors are 0 (background), @@ -172,24 +172,26 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma DetermineSliceOrientations(); /** An std::set of slice indices which need to be interpolated. */ - using SliceSetType = std::set< typename TImage::IndexValueType >; + using SliceSetType = std::set; /** Clears all custom slice positions. */ void ClearLabeledSliceIndices() { m_LabeledSlices.clear(); - m_LabeledSlices.resize( TImage::ImageDimension ); + m_LabeledSlices.resize(TImage::ImageDimension); this->Modified(); } /** If default slice detection is not wanted, slice indices * between which interpolation is done can be set using this method. */ void - SetLabeledSliceIndices( unsigned int axis, typename TImage::PixelType label, const std::vector< typename TImage::IndexValueType >& indices ) + SetLabeledSliceIndices(unsigned int axis, + typename TImage::PixelType label, + const std::vector & indices) { SliceSetType sliceSet; - sliceSet.insert( indices.begin(), indices.end() ); + sliceSet.insert(indices.begin(), indices.end()); m_LabeledSlices[axis][label] = sliceSet; this->Modified(); } @@ -197,7 +199,7 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma /** If default slice detection is not wanted, slice indices * between which interpolation is done can be set using this method. */ void - SetLabeledSliceIndices( unsigned int axis, typename TImage::PixelType label, const SliceSetType& indices ) + SetLabeledSliceIndices(unsigned int axis, typename TImage::PixelType label, const SliceSetType & indices) { m_LabeledSlices[axis][label] = indices; this->Modified(); @@ -205,14 +207,14 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma /** Slice indices between which interpolation is done. */ SliceSetType - GetLabeledSliceIndices( unsigned int axis, typename TImage::PixelType label ) + GetLabeledSliceIndices(unsigned int axis, typename TImage::PixelType label) { return m_LabeledSlices[axis][label]; } // each label gets a set of slices in which it is present - using LabeledSlicesType = std::unordered_map< typename TImage::PixelType, SliceSetType >; - using SliceIndicesType = std::vector< LabeledSlicesType >; + using LabeledSlicesType = std::unordered_map; + using SliceIndicesType = std::vector; /** Slice indices between which interpolation is done. */ SliceIndicesType @@ -237,13 +239,13 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma SliceIndicesType m_LabeledSlices; // one for each axis /** Derived image type alias. */ - using BoolImageType = Image< bool, TImage::ImageDimension >; - using FloatSliceType = Image< float, TImage::ImageDimension - 1 >; - using BoolSliceType = Image< bool, TImage::ImageDimension - 1 >; + using BoolImageType = Image; + using FloatSliceType = Image; + using BoolSliceType = Image; /** Are these two slices equal? */ bool - ImagesEqual( typename BoolSliceType::Pointer& a, typename BoolSliceType::Pointer& b ); + ImagesEqual(typename BoolSliceType::Pointer & a, typename BoolSliceType::Pointer & b); /** Does the real work. */ void @@ -251,8 +253,13 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma /** Determines correspondances between two slices and calls apropriate methods. */ void - InterpolateBetweenTwo( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, typename SliceType::Pointer& iconn, - typename SliceType::Pointer& jconn ); + InterpolateBetweenTwo(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iconn, + typename SliceType::Pointer & jconn); /** If interpolation is done along more than one axis, the interpolations are merged using a modified "or" rule: @@ -260,50 +267,78 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma -if just one image has a non-zero label, then that label is chosen -if more than one image has a non-zero label, last written label is chosen */ void - InterpolateAlong( int axis, TImage* out, float startProgress, float endProgress ); + InterpolateAlong(int axis, TImage * out, float startProgress, float endProgress); /** Slice i has a region, slice j does not */ void - Extrapolate( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, typename SliceType::Pointer& iConn, - typename TImage::PixelType iRegionId ); + Extrapolate(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId); /** Creates a signed distance field image. */ typename FloatSliceType::Pointer - MaurerDM( typename BoolSliceType::Pointer& inImage ); + MaurerDM(typename BoolSliceType::Pointer & inImage); /** A sequence of conditional dilations starting with begin and reaching end. begin and end must cover the same region (size and index the same) */ - std::vector< typename BoolSliceType::Pointer > - GenerateDilationSequence( typename BoolSliceType::Pointer& begin, typename BoolSliceType::Pointer& end ); + std::vector + GenerateDilationSequence(typename BoolSliceType::Pointer & begin, typename BoolSliceType::Pointer & end); /** Finds an interpolating mask for these two aligned masks */ typename BoolSliceType::Pointer - FindMedianImageDilations( typename BoolSliceType::Pointer& intersection, typename BoolSliceType::Pointer& iMask, typename BoolSliceType::Pointer& jMask ); + FindMedianImageDilations(typename BoolSliceType::Pointer & intersection, + typename BoolSliceType::Pointer & iMask, + typename BoolSliceType::Pointer & jMask); /** Finds an interpolating mask for these two aligned masks */ typename BoolSliceType::Pointer - FindMedianImageDistances( typename BoolSliceType::Pointer& intersection, typename BoolSliceType::Pointer& iMask, typename BoolSliceType::Pointer& jMask ); + FindMedianImageDistances(typename BoolSliceType::Pointer & intersection, + typename BoolSliceType::Pointer & iMask, + typename BoolSliceType::Pointer & jMask); /** Build transition sequence and pick the median */ void - Interpolate1to1( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, typename SliceType::Pointer& iConn, - typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, typename TImage::PixelType jRegionId, const typename SliceType::IndexType& translation, bool recursive ); - - using PixelList = std::vector< typename TImage::PixelType >; + Interpolate1to1(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + typename TImage::PixelType jRegionId, + const typename SliceType::IndexType & translation, + bool recursive); + + using PixelList = std::vector; /** Splits the bigger region and does N 1-to-1 interpolations */ void - Interpolate1toN( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, typename SliceType::Pointer& iConn, - typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds, const typename SliceType::IndexType& translation ); + Interpolate1toN(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds, + const typename SliceType::IndexType & translation); /** Crates a translated copy of part of the image which fits in the newRegion. */ typename SliceType::Pointer - TranslateImage( typename SliceType::Pointer& image, const typename SliceType::IndexType& translation, typename SliceType::RegionType newRegion ); + TranslateImage(typename SliceType::Pointer & image, + const typename SliceType::IndexType & translation, + typename SliceType::RegionType newRegion); /** The returns cardingal of the symmetric distance between images. The images must cover the same region. */ IdentifierType - CardinalSymmetricDifference( typename BoolSliceType::Pointer& shape1, typename BoolSliceType::Pointer& shape2 ); + CardinalSymmetricDifference(typename BoolSliceType::Pointer & shape1, typename BoolSliceType::Pointer & shape2); /** Copied from ImageSource and changed to allocate a cleared buffer. */ void @@ -311,57 +346,67 @@ class MorphologicalContourInterpolator : public ImageToImageFilter< TImage, TIma /** Returns the centroid of given regions */ typename SliceType::IndexType - Centroid( typename SliceType::Pointer& conn, const PixelList& regionIds ); + Centroid(typename SliceType::Pointer & conn, const PixelList & regionIds); /** Calculates maximum intersection region for both slices given a translation. Both inputs are modified so that jRegion is a translated version of iRegion. */ void - IntersectionRegions( const typename SliceType::IndexType& translation, typename SliceType::RegionType& iRegion, typename SliceType::RegionType& jRegion ); + IntersectionRegions(const typename SliceType::IndexType & translation, + typename SliceType::RegionType & iRegion, + typename SliceType::RegionType & jRegion); /** Returns number of intersecting pixels */ IdentifierType - Intersection( typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds, - const typename SliceType::IndexType& translation ); + Intersection(typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds, + const typename SliceType::IndexType & translation); /** How much j needs to be translated to best align with i */ typename SliceType::IndexType - Align( typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds ); + Align(typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds); - using BoundingBoxesType = std::unordered_map< typename TImage::PixelType, typename TImage::RegionType >; + using BoundingBoxesType = std::unordered_map; BoundingBoxesType m_BoundingBoxes; // bounding box for each label /** Calculates a bounding box of non-zero pixels. */ typename SliceType::RegionType - BoundingBox( itk::SmartPointer< SliceType > image ); + BoundingBox(itk::SmartPointer image); /** Expands a region to incorporate the provided index. * Assumes both a valid region and a valid index. * It can be invoked with 2D or 3D region, hence the additional template parameter. */ - template< typename T2 > + template void - ExpandRegion( typename T2::RegionType& region, const typename T2::IndexType& index ); + ExpandRegion(typename T2::RegionType & region, const typename T2::IndexType & index); /** Connected components of a specified region. */ typename SliceType::Pointer - RegionedConnectedComponents( const typename TImage::RegionType& region, typename TImage::PixelType label, IdentifierType& objectCount ); + RegionedConnectedComponents(const typename TImage::RegionType & region, + typename TImage::PixelType label, + IdentifierType & objectCount); /** Seed and mask must cover the same region (size and index the same). */ typename BoolSliceType::Pointer - Dilate1( typename BoolSliceType::Pointer& seed, typename BoolSliceType::Pointer& mask ); + Dilate1(typename BoolSliceType::Pointer & seed, typename BoolSliceType::Pointer & mask); - using RoiType = ExtractImageFilter< TImage, SliceType >; + using RoiType = ExtractImageFilter; typename RoiType::Pointer m_RoI; - using BinarizerType = BinaryThresholdImageFilter< SliceType, BoolSliceType >; + using BinarizerType = BinaryThresholdImageFilter; typename BinarizerType::Pointer m_Binarizer; - using ConnectedComponentsType = ConnectedComponentImageFilter< BoolSliceType, SliceType >; + using ConnectedComponentsType = ConnectedComponentImageFilter; typename ConnectedComponentsType::Pointer m_ConnectedComponents; }; } // namespace itk #ifndef ITK_MANUAL_INSTANTIATION -#include "itkMorphologicalContourInterpolator.hxx" +# include "itkMorphologicalContourInterpolator.hxx" #endif #endif // itkMorphologicalContourInterpolator_h diff --git a/include/itkMorphologicalContourInterpolator.hxx b/include/itkMorphologicalContourInterpolator.hxx index 737d967..216d075 100644 --- a/include/itkMorphologicalContourInterpolator.hxx +++ b/include/itkMorphologicalContourInterpolator.hxx @@ -41,505 +41,531 @@ namespace itk { -template< typename TImage > +template struct SegmentBetweenTwo { - int axis; - TImage* out; - int label, i, j; - typename MorphologicalContourInterpolator< TImage >::SliceType::Pointer iconn, jconn; + int axis; + TImage * out; + int label, i, j; + typename MorphologicalContourInterpolator::SliceType::Pointer iconn, jconn; }; -template< typename TImage > +template bool -MorphologicalContourInterpolator< TImage >::ImagesEqual( typename BoolSliceType::Pointer& a, typename BoolSliceType::Pointer& b ) +MorphologicalContourInterpolator::ImagesEqual(typename BoolSliceType::Pointer & a, + typename BoolSliceType::Pointer & b) { - ImageRegionConstIterator< BoolSliceType > ita( a, a->GetLargestPossibleRegion() ); - ImageRegionConstIterator< BoolSliceType > itb( b, b->GetLargestPossibleRegion() ); + ImageRegionConstIterator ita(a, a->GetLargestPossibleRegion()); + ImageRegionConstIterator itb(b, b->GetLargestPossibleRegion()); - while ( !ita.IsAtEnd() ) + while (!ita.IsAtEnd()) + { + if (ita.Get() != itb.Get()) { - if ( ita.Get() != itb.Get() ) - { - break; - } - ++ita; - ++itb; + break; } + ++ita; + ++itb; + } - if ( ita.IsAtEnd() ) - { - return true; - } + if (ita.IsAtEnd()) + { + return true; + } else - { - return false; - } + { + return false; + } } // >::ImagesEqual -template< typename TImage > -MorphologicalContourInterpolator< TImage >::MorphologicalContourInterpolator() - : m_Label( 0 ) +template +MorphologicalContourInterpolator::MorphologicalContourInterpolator() + : m_Label(0) , - m_MinAlignIters( std::pow( 2., static_cast< int >( TImage::ImageDimension ) ) ) + m_MinAlignIters(std::pow(2., static_cast(TImage::ImageDimension))) , // smaller of this and pixel count of the search image - m_MaxAlignIters( std::pow( 6., static_cast< int >( TImage::ImageDimension ) ) ) + m_MaxAlignIters(std::pow(6., static_cast(TImage::ImageDimension))) , // bigger of this and root of pixel count of the search image - m_ThreadCount( MultiThreaderBase::GetGlobalDefaultNumberOfThreads() ) - , m_LabeledSlices( TImage::ImageDimension ) // initialize with empty sets + m_ThreadCount(MultiThreaderBase::GetGlobalDefaultNumberOfThreads()) + , m_LabeledSlices(TImage::ImageDimension) // initialize with empty sets { // set up pipeline for regioned connected components m_RoI = RoiType::New(); m_RoI->SetDirectionCollapseToIdentity(); m_Binarizer = BinarizerType::New(); - m_Binarizer->SetInput( m_RoI->GetOutput() ); + m_Binarizer->SetInput(m_RoI->GetOutput()); m_ConnectedComponents = ConnectedComponentsType::New(); - m_ConnectedComponents->SetInput( m_Binarizer->GetOutput() ); + m_ConnectedComponents->SetInput(m_Binarizer->GetOutput()); // FullyConnected is related to structuring element used // true for ball, false for cross - m_ConnectedComponents->SetFullyConnected( m_UseBallStructuringElement ); - m_ConnectedComponents->SetBackgroundValue( NumericTraits< typename TImage::PixelType >::ZeroValue() ); + m_ConnectedComponents->SetFullyConnected(m_UseBallStructuringElement); + m_ConnectedComponents->SetBackgroundValue(NumericTraits::ZeroValue()); } -template< typename TImage > -template< typename T2 > +template +template void -MorphologicalContourInterpolator< TImage >::ExpandRegion( typename T2::RegionType& region, const typename T2::IndexType& index ) +MorphologicalContourInterpolator::ExpandRegion(typename T2::RegionType & region, + const typename T2::IndexType & index) { - for ( unsigned int a = 0; a < T2::ImageDimension; ++a ) + for (unsigned int a = 0; a < T2::ImageDimension; ++a) + { + if (region.GetIndex(a) > index[a]) { - if ( region.GetIndex( a ) > index[a] ) - { - region.SetSize( a, region.GetSize( a ) + region.GetIndex( a ) - index[a] ); - region.SetIndex( a, index[a] ); - } - else if ( region.GetIndex( a ) + (typename T2::IndexValueType)region.GetSize( a ) <= index[a] ) - { - region.SetSize( a, index[a] - region.GetIndex( a ) + 1 ); - } - // else it is already within + region.SetSize(a, region.GetSize(a) + region.GetIndex(a) - index[a]); + region.SetIndex(a, index[a]); } + else if (region.GetIndex(a) + (typename T2::IndexValueType)region.GetSize(a) <= index[a]) + { + region.SetSize(a, index[a] - region.GetIndex(a) + 1); + } + // else it is already within + } } -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::DetermineSliceOrientations() +MorphologicalContourInterpolator::DetermineSliceOrientations() { m_LabeledSlices.clear(); - m_LabeledSlices.resize( TImage::ImageDimension ); // initialize with empty sets + m_LabeledSlices.resize(TImage::ImageDimension); // initialize with empty sets m_BoundingBoxes.clear(); typename TImage::ConstPointer m_Input = this->GetInput(); typename TImage::Pointer m_Output = this->GetOutput(); - typename TImage::RegionType region = m_Output->GetRequestedRegion(); - ImageRegionConstIteratorWithIndex< TImage > it( m_Input, region ); - while ( !it.IsAtEnd() ) + typename TImage::RegionType region = m_Output->GetRequestedRegion(); + ImageRegionConstIteratorWithIndex it(m_Input, region); + while (!it.IsAtEnd()) + { + typename TImage::IndexType indPrev, indNext; + const typename TImage::IndexType ind = it.GetIndex(); + const typename TImage::PixelType val = m_Input->GetPixel(ind); + if (val != 0) { - typename TImage::IndexType indPrev, indNext; - const typename TImage::IndexType ind = it.GetIndex(); - const typename TImage::PixelType val = m_Input->GetPixel( ind ); - if ( val != 0 ) + typename TImage::RegionType boundingBox1; + boundingBox1.SetIndex(ind); + for (unsigned int a = 0; a < TImage::ImageDimension; ++a) + { + boundingBox1.SetSize(a, 1); + } + std::pair resBB = + m_BoundingBoxes.insert(std::make_pair(val, boundingBox1)); + if (!resBB.second) // include this index in existing BB + { + ExpandRegion(resBB.first->second, ind); + } + + unsigned int cTrue = 0; + unsigned int cAdjacent = 0; + unsigned int axis = 0; + for (unsigned int a = 0; a < TImage::ImageDimension; ++a) + { + indPrev = ind; + indPrev[a]--; + indNext = ind; + indNext[a]++; + typename TImage::PixelType prev = 0; + if (region.IsInside(indPrev)) { - typename TImage::RegionType boundingBox1; - boundingBox1.SetIndex( ind ); - for ( unsigned int a = 0; a < TImage::ImageDimension; ++a ) - { - boundingBox1.SetSize( a, 1 ); - } - std::pair< typename BoundingBoxesType::iterator, bool > resBB = m_BoundingBoxes.insert( std::make_pair( val, boundingBox1 ) ); - if ( !resBB.second ) // include this index in existing BB - { - ExpandRegion< TImage >( resBB.first->second, ind ); - } - - unsigned int cTrue = 0; - unsigned int cAdjacent = 0; - unsigned int axis = 0; - for ( unsigned int a = 0; a < TImage::ImageDimension; ++a ) - { - indPrev = ind; - indPrev[a]--; - indNext = ind; - indNext[a]++; - typename TImage::PixelType prev = 0; - if ( region.IsInside( indPrev ) ) - { - prev = m_Input->GetPixel( indPrev ); - } - typename TImage::PixelType next = 0; - if ( region.IsInside( indNext ) ) - { - next = m_Input->GetPixel( indNext ); - } - if ( prev == 0 && next == 0 ) // && - isolated slices only, || - flat edges too - { - axis = a; - ++cTrue; - } - else if ( prev == val && next == val ) - { - ++cAdjacent; - } - } - if ( cTrue == 1 && cAdjacent == TImage::ImageDimension - 1 ) - // slice has empty adjacent space only along one axis - { - if ( m_Axis == -1 || m_Axis == int( axis ) ) - { - m_LabeledSlices[axis][val].insert( ind[axis] ); - } - } + prev = m_Input->GetPixel(indPrev); } - ++it; + typename TImage::PixelType next = 0; + if (region.IsInside(indNext)) + { + next = m_Input->GetPixel(indNext); + } + if (prev == 0 && next == 0) // && - isolated slices only, || - flat edges too + { + axis = a; + ++cTrue; + } + else if (prev == val && next == val) + { + ++cAdjacent; + } + } + if (cTrue == 1 && cAdjacent == TImage::ImageDimension - 1) + // slice has empty adjacent space only along one axis + { + if (m_Axis == -1 || m_Axis == int(axis)) + { + m_LabeledSlices[axis][val].insert(ind[axis]); + } + } } + ++it; + } } // >::DetermineSliceOrientations -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::Extrapolate( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, - typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId ) +MorphologicalContourInterpolator::Extrapolate(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId) { PixelList jRegionIds; - jRegionIds.push_back( iRegionId ); - typename SliceType::IndexType centroid = Centroid( iConn, jRegionIds ); + jRegionIds.push_back(iRegionId); + typename SliceType::IndexType centroid = Centroid(iConn, jRegionIds); typename SliceType::RegionType reg3; typename SliceType::SizeType size; - size.Fill( 3 ); - reg3.SetSize( size ); + size.Fill(3); + reg3.SetSize(size); typename SliceType::IndexType phIndex; - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - phIndex[d] = centroid.GetIndex()[d] - 1; - } - reg3.SetIndex( phIndex ); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + phIndex[d] = centroid.GetIndex()[d] - 1; + } + reg3.SetIndex(phIndex); // create a phantom small slice centered around centroid typename SliceType::Pointer phSlice = SliceType::New(); - phSlice->CopyInformation( iConn ); - phSlice->SetRegions( reg3 ); - phSlice->Allocate( true ); + phSlice->CopyInformation(iConn); + phSlice->SetRegions(reg3); + phSlice->Allocate(true); // add a phantom point to the center of a newly constructed slice - phSlice->SetPixel( centroid, iRegionId ); + phSlice->SetPixel(centroid, iRegionId); // do alignment in case the iShape is concave and centroid is not within the iShape - typename SliceType::IndexType translation = Align( iConn, iRegionId, phSlice, jRegionIds ); + typename SliceType::IndexType translation = Align(iConn, iRegionId, phSlice, jRegionIds); // now translate the phantom slice for best alignment - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - phIndex[d] -= translation[d]; - } - reg3.SetIndex( phIndex ); - phSlice->SetRegions( reg3 ); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + phIndex[d] -= translation[d]; + } + reg3.SetIndex(phIndex); + phSlice->SetRegions(reg3); typename SliceType::IndexType t0; - t0.Fill( 0 ); - Interpolate1to1( axis, out, label, i, j, iConn, iRegionId, phSlice, iRegionId, t0, false ); + t0.Fill(0); + Interpolate1to1(axis, out, label, i, j, iConn, iRegionId, phSlice, iRegionId, t0, false); } // >::Extrapolate -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::Dilate1( typename BoolSliceType::Pointer& seed, typename BoolSliceType::Pointer& mask ) -> typename BoolSliceType::Pointer +MorphologicalContourInterpolator::Dilate1(typename BoolSliceType::Pointer & seed, + typename BoolSliceType::Pointer & mask) -> + typename BoolSliceType::Pointer { // set up structuring element for dilation - using CrossStructuringElementType = BinaryCrossStructuringElement< typename BoolSliceType::PixelType, BoolSliceType::ImageDimension >; - using BallStructuringElementType = BinaryBallStructuringElement< typename BoolSliceType::PixelType, BoolSliceType::ImageDimension >; - using CrossDilateType = BinaryDilateImageFilter< BoolSliceType, BoolSliceType, CrossStructuringElementType >; - using BallDilateType = BinaryDilateImageFilter< BoolSliceType, BoolSliceType, BallStructuringElementType >; + using CrossStructuringElementType = + BinaryCrossStructuringElement; + using BallStructuringElementType = + BinaryBallStructuringElement; + using CrossDilateType = BinaryDilateImageFilter; + using BallDilateType = BinaryDilateImageFilter; thread_local bool initialized = false; thread_local typename CrossDilateType::Pointer crossDilator = CrossDilateType::New(); thread_local typename BallDilateType::Pointer ballDilator = BallDilateType::New(); thread_local CrossStructuringElementType crossStructuringElement; thread_local BallStructuringElementType ballStructuringElement; - using AndFilterType = AndImageFilter< BoolSliceType, BoolSliceType, BoolSliceType >; + using AndFilterType = AndImageFilter; thread_local typename AndFilterType::Pointer andFilter = AndFilterType::New(); - if ( !initialized ) // make sure these non-trivial operations are executed only once per thread - { - andFilter->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - using SizeType = Size< BoolSliceType::ImageDimension >; - SizeType size; - size.Fill( 1 ); - - crossDilator->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - crossStructuringElement.SetRadius( size ); - crossStructuringElement.CreateStructuringElement(); - crossDilator->SetKernel( crossStructuringElement ); - - ballDilator->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - ballStructuringElement.SetRadius( size ); - ballStructuringElement.CreateStructuringElement(); - ballDilator->SetKernel( ballStructuringElement ); - - initialized = true; - } + if (!initialized) // make sure these non-trivial operations are executed only once per thread + { + andFilter->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + using SizeType = Size; + SizeType size; + size.Fill(1); + + crossDilator->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + crossStructuringElement.SetRadius(size); + crossStructuringElement.CreateStructuringElement(); + crossDilator->SetKernel(crossStructuringElement); + + ballDilator->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + ballStructuringElement.SetRadius(size); + ballStructuringElement.CreateStructuringElement(); + ballDilator->SetKernel(ballStructuringElement); + + initialized = true; + } typename BoolSliceType::Pointer temp; - if ( m_UseBallStructuringElement ) - { - ballDilator->SetInput( seed ); - ballDilator->GetOutput()->SetRegions( seed->GetRequestedRegion() ); - ballDilator->Update(); - temp = ballDilator->GetOutput(); - } + if (m_UseBallStructuringElement) + { + ballDilator->SetInput(seed); + ballDilator->GetOutput()->SetRegions(seed->GetRequestedRegion()); + ballDilator->Update(); + temp = ballDilator->GetOutput(); + } else - { - crossDilator->SetInput( seed ); - crossDilator->GetOutput()->SetRegions( seed->GetRequestedRegion() ); - crossDilator->Update(); - temp = crossDilator->GetOutput(); - } + { + crossDilator->SetInput(seed); + crossDilator->GetOutput()->SetRegions(seed->GetRequestedRegion()); + crossDilator->Update(); + temp = crossDilator->GetOutput(); + } temp->DisconnectPipeline(); // temp->SetRegions(mask->GetLargestPossibleRegion()); //not needed when seed and mask have same regions - andFilter->SetInput( 0, mask ); - andFilter->SetInput( 1, temp ); - andFilter->GetOutput()->SetRegions( seed->GetRequestedRegion() ); + andFilter->SetInput(0, mask); + andFilter->SetInput(1, temp); + andFilter->GetOutput()->SetRegions(seed->GetRequestedRegion()); andFilter->Update(); typename BoolSliceType::Pointer result = andFilter->GetOutput(); result->DisconnectPipeline(); return result; } // >::Dilate1 -template< typename TImage > -std::vector< typename MorphologicalContourInterpolator< TImage >::BoolSliceType::Pointer > -MorphologicalContourInterpolator< TImage >::GenerateDilationSequence( typename BoolSliceType::Pointer& begin, typename BoolSliceType::Pointer& end ) +template +std::vector::BoolSliceType::Pointer> +MorphologicalContourInterpolator::GenerateDilationSequence(typename BoolSliceType::Pointer & begin, + typename BoolSliceType::Pointer & end) { - std::vector< typename BoolSliceType::Pointer > seq; - seq.push_back( Dilate1( begin, end ) ); + std::vector seq; + seq.push_back(Dilate1(begin, end)); do - { - seq.back()->DisconnectPipeline(); - seq.push_back( Dilate1( seq.back(), end ) ); - } - while ( !ImagesEqual( seq.back(), seq[seq.size() - 2] ) ); + { + seq.back()->DisconnectPipeline(); + seq.push_back(Dilate1(seq.back(), end)); + } while (!ImagesEqual(seq.back(), seq[seq.size() - 2])); seq.pop_back(); // remove duplicate image return seq; } -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::FindMedianImageDilations( typename BoolSliceType::Pointer& intersection, typename BoolSliceType::Pointer& iMask, typename BoolSliceType::Pointer& jMask ) -> +MorphologicalContourInterpolator::FindMedianImageDilations(typename BoolSliceType::Pointer & intersection, + typename BoolSliceType::Pointer & iMask, + typename BoolSliceType::Pointer & jMask) -> typename BoolSliceType::Pointer { - std::vector< typename BoolSliceType::Pointer > iSeq = GenerateDilationSequence( intersection, iMask ); - std::vector< typename BoolSliceType::Pointer > jSeq = GenerateDilationSequence( intersection, jMask ); - std::reverse( iSeq.begin(), iSeq.end() ); // we want to start from i and end at intersection - if ( iSeq.size() < jSeq.size() ) - { - iSeq.swap( jSeq ); // swap so iSeq.size() >= jSeq.size() - } - float ratio = float( jSeq.size() ) / iSeq.size(); + std::vector iSeq = GenerateDilationSequence(intersection, iMask); + std::vector jSeq = GenerateDilationSequence(intersection, jMask); + std::reverse(iSeq.begin(), iSeq.end()); // we want to start from i and end at intersection + if (iSeq.size() < jSeq.size()) + { + iSeq.swap(jSeq); // swap so iSeq.size() >= jSeq.size() + } + float ratio = float(jSeq.size()) / iSeq.size(); // generate union of transition sequences - using OrType = OrImageFilter< BoolSliceType >; + using OrType = OrImageFilter; thread_local typename OrType::Pointer orFilter = OrType::New(); - orFilter->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive + orFilter->SetNumberOfWorkUnits(1); // excessive threading is counterproductive - std::vector< typename BoolSliceType::Pointer > seq; - for ( unsigned x = 0; x < iSeq.size(); x++ ) - { - orFilter->SetInput( 0, iSeq[x] ); - unsigned xj = ratio * x; - orFilter->SetInput( 1, jSeq[xj] ); - orFilter->GetOutput()->SetRegions( iMask->GetRequestedRegion() ); - orFilter->Update(); - seq.push_back( orFilter->GetOutput() ); - seq.back()->DisconnectPipeline(); - } + std::vector seq; + for (unsigned x = 0; x < iSeq.size(); x++) + { + orFilter->SetInput(0, iSeq[x]); + unsigned xj = ratio * x; + orFilter->SetInput(1, jSeq[xj]); + orFilter->GetOutput()->SetRegions(iMask->GetRequestedRegion()); + orFilter->Update(); + seq.push_back(orFilter->GetOutput()); + seq.back()->DisconnectPipeline(); + } // find median unsigned minIndex = 0; IdentifierType min = iMask->GetRequestedRegion().GetNumberOfPixels(); - for ( unsigned x = 0; x < iSeq.size(); x++ ) + for (unsigned x = 0; x < iSeq.size(); x++) + { + IdentifierType iS = CardinalSymmetricDifference(seq[x], iMask); + IdentifierType jS = CardinalSymmetricDifference(seq[x], jMask); + IdentifierType xScore = iS >= jS ? iS - jS : jS - iS; // itk::Math::abs(iS-jS) + if (xScore < min) { - IdentifierType iS = CardinalSymmetricDifference( seq[x], iMask ); - IdentifierType jS = CardinalSymmetricDifference( seq[x], jMask ); - IdentifierType xScore = iS >= jS ? iS - jS : jS - iS; // itk::Math::abs(iS-jS) - if ( xScore < min ) - { - min = xScore; - minIndex = x; - } + min = xScore; + minIndex = x; } + } return seq[minIndex]; } // >::FindMedianImageDilations -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::MaurerDM( typename BoolSliceType::Pointer& mask ) -> typename FloatSliceType::Pointer +MorphologicalContourInterpolator::MaurerDM(typename BoolSliceType::Pointer & mask) -> + typename FloatSliceType::Pointer { - using FilterType = itk::SignedMaurerDistanceMapImageFilter< BoolSliceType, FloatSliceType >; + using FilterType = itk::SignedMaurerDistanceMapImageFilter; thread_local typename FilterType::Pointer filter = FilterType::New(); - filter->SetUseImageSpacing( false ); // interpolation algorithm calls for working in index space - filter->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - filter->SetInput( mask ); - filter->GetOutput()->SetRequestedRegion( mask->GetRequestedRegion() ); + filter->SetUseImageSpacing(false); // interpolation algorithm calls for working in index space + filter->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + filter->SetInput(mask); + filter->GetOutput()->SetRequestedRegion(mask->GetRequestedRegion()); filter->Update(); return filter->GetOutput(); } -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::FindMedianImageDistances( typename BoolSliceType::Pointer& intersection, typename BoolSliceType::Pointer& iMask, typename BoolSliceType::Pointer& jMask ) -> +MorphologicalContourInterpolator::FindMedianImageDistances(typename BoolSliceType::Pointer & intersection, + typename BoolSliceType::Pointer & iMask, + typename BoolSliceType::Pointer & jMask) -> typename BoolSliceType::Pointer { // calculate distance field - typename FloatSliceType::Pointer sdf = MaurerDM( intersection ); + typename FloatSliceType::Pointer sdf = MaurerDM(intersection); // create histograms of distances and union typename BoolSliceType::Pointer orImage = BoolSliceType::New(); - orImage->CopyInformation( intersection ); - orImage->SetRegions( iMask->GetRequestedRegion() ); - orImage->Allocate( true ); - std::vector< long long > iHist; - std::vector< long long > jHist; - ImageRegionConstIterator< BoolSliceType > iti( iMask, iMask->GetRequestedRegion() ); - ImageRegionConstIterator< BoolSliceType > itj( jMask, iMask->GetRequestedRegion() ); - ImageRegionIterator< BoolSliceType > ito( orImage, iMask->GetRequestedRegion() ); - ImageRegionConstIterator< FloatSliceType > itsdf( sdf, iMask->GetRequestedRegion() ); - constexpr short fractioning = 10; // how many times more precise distance than rounding to int - while ( !itsdf.IsAtEnd() ) + orImage->CopyInformation(intersection); + orImage->SetRegions(iMask->GetRequestedRegion()); + orImage->Allocate(true); + std::vector iHist; + std::vector jHist; + ImageRegionConstIterator iti(iMask, iMask->GetRequestedRegion()); + ImageRegionConstIterator itj(jMask, iMask->GetRequestedRegion()); + ImageRegionIterator ito(orImage, iMask->GetRequestedRegion()); + ImageRegionConstIterator itsdf(sdf, iMask->GetRequestedRegion()); + constexpr short fractioning = 10; // how many times more precise distance than rounding to int + while (!itsdf.IsAtEnd()) + { + bool iM = iti.Get(); + bool jM = itj.Get(); + typename TImage::PixelType dist = fractioning * itsdf.Get(); + if (iM && !jM) { - bool iM = iti.Get(); - bool jM = itj.Get(); - typename TImage::PixelType dist = fractioning * itsdf.Get(); - if ( iM && !jM ) - { - if ( size_t( dist ) >= iHist.size() ) - { - iHist.resize( dist + 1, 0 ); - } - iHist[dist]++; - ito.Set( true ); - } - else if ( jM && !iM ) - { - if ( size_t( dist ) >= jHist.size() ) - { - jHist.resize( dist + 1, 0 ); - } - jHist[dist]++; - ito.Set( true ); - } - else if ( iM && jM ) - { - ito.Set( true ); - } - - ++iti; - ++itj; - ++ito; - ++itsdf; + if (size_t(dist) >= iHist.size()) + { + iHist.resize(dist + 1, 0); + } + iHist[dist]++; + ito.Set(true); } - - // sum of histogram bins for i and j and - std::vector< long long >::size_type maxSize = std::max( iHist.size(), jHist.size() ); - if ( maxSize == 0 ) + else if (jM && !iM) { - return intersection; + if (size_t(dist) >= jHist.size()) + { + jHist.resize(dist + 1, 0); + } + jHist[dist]++; + ito.Set(true); } - iHist.resize( maxSize, 0 ); - jHist.resize( maxSize, 0 ); - std::vector< long long > iSum( maxSize, 0 ); - std::vector< long long > jSum( maxSize, 0 ); - iSum[0] = iHist[0]; - jSum[0] = jHist[0]; - for ( unsigned b = 1; b < maxSize; b++ ) + else if (iM && jM) { - iSum[b] = iSum[b - 1] + iHist[b]; - jSum[b] = jSum[b - 1] + jHist[b]; + ito.Set(true); } + + ++iti; + ++itj; + ++ito; + ++itsdf; + } + + // sum of histogram bins for i and j and + std::vector::size_type maxSize = std::max(iHist.size(), jHist.size()); + if (maxSize == 0) + { + return intersection; + } + iHist.resize(maxSize, 0); + jHist.resize(maxSize, 0); + std::vector iSum(maxSize, 0); + std::vector jSum(maxSize, 0); + iSum[0] = iHist[0]; + jSum[0] = jHist[0]; + for (unsigned b = 1; b < maxSize; b++) + { + iSum[b] = iSum[b - 1] + iHist[b]; + jSum[b] = jSum[b - 1] + jHist[b]; + } long long iTotal = iSum[maxSize - 1]; long long jTotal = jSum[maxSize - 1]; // find minimum of differences of sums int bestBin = 0; long long bestDiff = LLONG_MAX; - for ( unsigned b = 0; b < maxSize; b++ ) + for (unsigned b = 0; b < maxSize; b++) + { + long long iS = itk::Math::abs(iTotal - iSum[b] + jSum[b]); + long long jS = itk::Math::abs(jTotal - jSum[b] + iSum[b]); + long long diff = itk::Math::abs(iS - jS); + if (diff < bestDiff) { - long long iS = itk::Math::abs( iTotal - iSum[b] + jSum[b] ); - long long jS = itk::Math::abs( jTotal - jSum[b] + iSum[b] ); - long long diff = itk::Math::abs( iS - jS ); - if ( diff < bestDiff ) - { - bestDiff = diff; - bestBin = b; - } + bestDiff = diff; + bestBin = b; } + } // threshold at distance bestBin is the median intersection - using FloatBinarizerType = BinaryThresholdImageFilter< FloatSliceType, BoolSliceType >; - using AndFilterType = AndImageFilter< BoolSliceType, BoolSliceType, BoolSliceType >; + using FloatBinarizerType = BinaryThresholdImageFilter; + using AndFilterType = AndImageFilter; thread_local typename FloatBinarizerType::Pointer threshold = FloatBinarizerType::New(); thread_local typename AndFilterType::Pointer andFilter = AndFilterType::New(); // excessive threading is counterproductive - threshold->SetNumberOfWorkUnits( 1 ); - andFilter->SetNumberOfWorkUnits( 1 ); + threshold->SetNumberOfWorkUnits(1); + andFilter->SetNumberOfWorkUnits(1); - threshold->SetInput( sdf ); - threshold->SetUpperThreshold( float( bestBin ) / fractioning ); - threshold->GetOutput()->SetRequestedRegion( sdf->GetRequestedRegion() ); + threshold->SetInput(sdf); + threshold->SetUpperThreshold(float(bestBin) / fractioning); + threshold->GetOutput()->SetRequestedRegion(sdf->GetRequestedRegion()); threshold->Update(); - andFilter->SetInput( threshold->GetOutput() ); - andFilter->SetInput( 1, orImage ); - andFilter->GetOutput()->SetRequestedRegion( orImage->GetRequestedRegion() ); + andFilter->SetInput(threshold->GetOutput()); + andFilter->SetInput(1, orImage); + andFilter->GetOutput()->SetRequestedRegion(orImage->GetRequestedRegion()); andFilter->Update(); typename BoolSliceType::Pointer median = andFilter->GetOutput(); return median; } // >::FindMedianImageDistances -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::BoundingBox( itk::SmartPointer< SliceType > image ) -> typename SliceType::RegionType +MorphologicalContourInterpolator::BoundingBox(itk::SmartPointer image) -> + typename SliceType::RegionType { - typename SliceType::RegionType newRegion = image->GetLargestPossibleRegion(); - typename SliceType::IndexType minInd = newRegion.GetIndex() + newRegion.GetSize(); - typename SliceType::IndexType maxInd = newRegion.GetIndex(); - ImageRegionConstIteratorWithIndex< SliceType > iIt( image, newRegion ); + typename SliceType::RegionType newRegion = image->GetLargestPossibleRegion(); + typename SliceType::IndexType minInd = newRegion.GetIndex() + newRegion.GetSize(); + typename SliceType::IndexType maxInd = newRegion.GetIndex(); + ImageRegionConstIteratorWithIndex iIt(image, newRegion); - while ( !iIt.IsAtEnd() ) + while (!iIt.IsAtEnd()) + { + if (iIt.Get()) { - if ( iIt.Get() ) + typename SliceType::IndexType ind = iIt.GetIndex(); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + if (ind[d] < minInd[d]) { - typename SliceType::IndexType ind = iIt.GetIndex(); - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - if ( ind[d] < minInd[d] ) - { - minInd[d] = ind[d]; - } - if ( ind[d] > maxInd[d] ) - { - maxInd[d] = ind[d]; - } - } + minInd[d] = ind[d]; } - ++iIt; + if (ind[d] > maxInd[d]) + { + maxInd[d] = ind[d]; + } + } } + ++iIt; + } - newRegion.SetIndex( minInd ); - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - newRegion.SetSize( d, maxInd[d] - minInd[d] + 1 ); - } + newRegion.SetIndex(minInd); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + newRegion.SetSize(d, maxInd[d] - minInd[d] + 1); + } return newRegion; } // >::BoundingBox -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::Interpolate1to1( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, - typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, - typename TImage::PixelType jRegionId, const typename SliceType::IndexType& translation, bool recursive ) +MorphologicalContourInterpolator::Interpolate1to1(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + typename TImage::PixelType jRegionId, + const typename SliceType::IndexType & translation, + bool recursive) { // translate iConn by t/2 and jConn by -t/2 typename SliceType::IndexType iTrans; @@ -548,261 +574,274 @@ MorphologicalContourInterpolator< TImage >::Interpolate1to1( int axis, TImage* o typename SliceType::RegionType jRegion = jConn->GetLargestPossibleRegion(); typename SliceType::IndexType jBottom = jRegion.GetIndex(); bool carry = false; - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + if (!carry) { - if ( !carry ) - { - iTrans[d] = translation[d] / 2; - carry = translation[d] % 2; - } - else if ( translation[d] % 2 == 0 ) - { - iTrans[d] = translation[d] / 2; - } - else // use carry - { - if ( translation[d] > 0 ) - { - iTrans[d] = translation[d] / 2 + 1; - } - else - { - iTrans[d] = translation[d] / 2 - 1; - } - carry = false; - } - jTrans[d] = iTrans[d] - translation[d]; - iRegion.SetIndex( d, iRegion.GetIndex()[d] + iTrans[d] ); - jRegion.SetIndex( d, jRegion.GetIndex()[d] + jTrans[d] ); - jBottom[d] = jRegion.GetIndex()[d] + jRegion.GetSize( d ) - 1; + iTrans[d] = translation[d] / 2; + carry = translation[d] % 2; + } + else if (translation[d] % 2 == 0) + { + iTrans[d] = translation[d] / 2; + } + else // use carry + { + if (translation[d] > 0) + { + iTrans[d] = translation[d] / 2 + 1; + } + else + { + iTrans[d] = translation[d] / 2 - 1; + } + carry = false; } + jTrans[d] = iTrans[d] - translation[d]; + iRegion.SetIndex(d, iRegion.GetIndex()[d] + iTrans[d]); + jRegion.SetIndex(d, jRegion.GetIndex()[d] + jTrans[d]); + jBottom[d] = jRegion.GetIndex()[d] + jRegion.GetSize(d) - 1; + } typename SliceType::RegionType newRegion = iRegion; - ExpandRegion< SliceType >( newRegion, jRegion.GetIndex() ); - ExpandRegion< SliceType >( newRegion, jBottom ); - typename SliceType::IndexValueType mid = ( i + j + ( carry ? 1 : 0 ) ) / 2; // index of middle slice + ExpandRegion(newRegion, jRegion.GetIndex()); + ExpandRegion(newRegion, jBottom); + typename SliceType::IndexValueType mid = (i + j + (carry ? 1 : 0)) / 2; // index of middle slice - typename SliceType::Pointer iConnT = TranslateImage( iConn, iTrans, newRegion ); - typename SliceType::Pointer jConnT = TranslateImage( jConn, jTrans, newRegion ); + typename SliceType::Pointer iConnT = TranslateImage(iConn, iTrans, newRegion); + typename SliceType::Pointer jConnT = TranslateImage(jConn, jTrans, newRegion); - if ( !recursive ) // reduce newRegion to bounding box so we deal with less pixels + if (!recursive) // reduce newRegion to bounding box so we deal with less pixels + { + newRegion = BoundingBox(iConnT); + typename SliceType::RegionType jBB = BoundingBox(jConnT); + typename SliceType::IndexType i2 = jBB.GetIndex(); + ExpandRegion(newRegion, i2); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) { - newRegion = BoundingBox( iConnT ); - typename SliceType::RegionType jBB = BoundingBox( jConnT ); - typename SliceType::IndexType i2 = jBB.GetIndex(); - ExpandRegion< SliceType >( newRegion, i2 ); - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - i2[d] += jBB.GetSize( d ) - 1; - } - ExpandRegion< SliceType >( newRegion, i2 ); + i2[d] += jBB.GetSize(d) - 1; } + ExpandRegion(newRegion, i2); + } // create and allocate slice binary masks typename BoolSliceType::Pointer iSlice = BoolSliceType::New(); - iSlice->CopyInformation( iConnT ); - iSlice->SetRegions( newRegion ); - iSlice->Allocate( true ); + iSlice->CopyInformation(iConnT); + iSlice->SetRegions(newRegion); + iSlice->Allocate(true); typename BoolSliceType::Pointer jSlice = BoolSliceType::New(); - jSlice->CopyInformation( jConnT ); - jSlice->SetRegions( newRegion ); - jSlice->Allocate( true ); + jSlice->CopyInformation(jConnT); + jSlice->SetRegions(newRegion); + jSlice->Allocate(true); // convert to binary by iteration - ImageRegionConstIterator< SliceType > iIt( iConnT, newRegion ); - ImageRegionConstIterator< SliceType > jIt( jConnT, newRegion ); - ImageRegionIterator< BoolSliceType > ibIt( iSlice, newRegion ); - ImageRegionIterator< BoolSliceType > jbIt( jSlice, newRegion ); - while ( !iIt.IsAtEnd() ) + ImageRegionConstIterator iIt(iConnT, newRegion); + ImageRegionConstIterator jIt(jConnT, newRegion); + ImageRegionIterator ibIt(iSlice, newRegion); + ImageRegionIterator jbIt(jSlice, newRegion); + while (!iIt.IsAtEnd()) + { + if (iIt.Get() == iRegionId) { - if ( iIt.Get() == iRegionId ) - { - ibIt.Set( true ); - } - if ( jIt.Get() == jRegionId ) - { - jbIt.Set( true ); - } - ++iIt; - ++jIt; - ++ibIt; - ++jbIt; + ibIt.Set(true); + } + if (jIt.Get() == jRegionId) + { + jbIt.Set(true); } + ++iIt; + ++jIt; + ++ibIt; + ++jbIt; + } // create intersection - using AndSliceType = AndImageFilter< BoolSliceType >; + using AndSliceType = AndImageFilter; thread_local typename AndSliceType::Pointer sAnd = AndSliceType::New(); - sAnd->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - sAnd->SetInput( 0, iSlice ); - sAnd->SetInput( 1, jSlice ); - sAnd->GetOutput()->SetRegions( iSlice->GetRequestedRegion() ); + sAnd->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + sAnd->SetInput(0, iSlice); + sAnd->SetInput(1, jSlice); + sAnd->GetOutput()->SetRegions(iSlice->GetRequestedRegion()); sAnd->Update(); typename BoolSliceType::Pointer intersection = sAnd->GetOutput(); intersection->DisconnectPipeline(); typename BoolSliceType::Pointer median; - if ( m_UseDistanceTransform ) - { - median = FindMedianImageDistances( intersection, iSlice, jSlice ); - } + if (m_UseDistanceTransform) + { + median = FindMedianImageDistances(intersection, iSlice, jSlice); + } else - { - median = FindMedianImageDilations( intersection, iSlice, jSlice ); - } + { + median = FindMedianImageDilations(intersection, iSlice, jSlice); + } // finally write it out into the output image pointer typename TImage::RegionType outRegion = this->GetOutput()->GetRequestedRegion(); typename SliceType::RegionType sliceRegion; - for ( int d = 0; d < int( TImage::ImageDimension ) - 1; d++ ) + for (int d = 0; d < int(TImage::ImageDimension) - 1; d++) + { + if (d < axis) { - if ( d < axis ) - { - sliceRegion.SetIndex( d, outRegion.GetIndex( d ) ); - sliceRegion.SetSize( d, outRegion.GetSize( d ) ); - } - else - { - sliceRegion.SetIndex( d, outRegion.GetIndex( d + 1 ) ); - sliceRegion.SetSize( d, outRegion.GetSize( d + 1 ) ); - } + sliceRegion.SetIndex(d, outRegion.GetIndex(d)); + sliceRegion.SetSize(d, outRegion.GetSize(d)); + } + else + { + sliceRegion.SetIndex(d, outRegion.GetIndex(d + 1)); + sliceRegion.SetSize(d, outRegion.GetSize(d + 1)); } + } typename SliceType::IndexType t0; - t0.Fill( 0 ); - IntersectionRegions( t0, sliceRegion, newRegion ); // clips new region to output region + t0.Fill(0); + IntersectionRegions(t0, sliceRegion, newRegion); // clips new region to output region // sliceRegion possibly shrunk, copy it into outRegion - for ( int d = 0; d < int( TImage::ImageDimension ) - 1; d++ ) + for (int d = 0; d < int(TImage::ImageDimension) - 1; d++) + { + if (d < axis) { - if ( d < axis ) - { - outRegion.SetIndex( d, sliceRegion.GetIndex( d ) ); - outRegion.SetSize( d, sliceRegion.GetSize( d ) ); - } - else - { - outRegion.SetIndex( d + 1, sliceRegion.GetIndex( d ) ); - outRegion.SetSize( d + 1, sliceRegion.GetSize( d ) ); - } + outRegion.SetIndex(d, sliceRegion.GetIndex(d)); + outRegion.SetSize(d, sliceRegion.GetSize(d)); + } + else + { + outRegion.SetIndex(d + 1, sliceRegion.GetIndex(d)); + outRegion.SetSize(d + 1, sliceRegion.GetSize(d)); } - outRegion.SetIndex( axis, mid ); - outRegion.SetSize( axis, 1 ); + } + outRegion.SetIndex(axis, mid); + outRegion.SetSize(axis, 1); typename SliceType::Pointer midConn = SliceType::New(); - midConn->CopyInformation( iConnT ); - midConn->SetRegions( newRegion ); - midConn->Allocate( true ); + midConn->CopyInformation(iConnT); + midConn->SetRegions(newRegion); + midConn->Allocate(true); - ImageRegionConstIterator< BoolSliceType > seqIt( median, newRegion ); - ImageRegionIterator< SliceType > midIt( midConn, newRegion ); - while ( !seqIt.IsAtEnd() ) + ImageRegionConstIterator seqIt(median, newRegion); + ImageRegionIterator midIt(midConn, newRegion); + while (!seqIt.IsAtEnd()) + { + if (seqIt.Get()) { - if ( seqIt.Get() ) - { - midIt.Set( 1 ); - } - ++seqIt; - ++midIt; + midIt.Set(1); } + ++seqIt; + ++midIt; + } bool withinReq = true; typename TImage::RegionType reqRegion = this->GetOutput()->GetRequestedRegion(); - for ( unsigned d = 0; d < TImage::ImageDimension; d++ ) + for (unsigned d = 0; d < TImage::ImageDimension; d++) + { + if (outRegion.GetIndex(d) < reqRegion.GetIndex(d) || + outRegion.GetIndex(d) + outRegion.GetSize(d) > reqRegion.GetIndex(d) + reqRegion.GetSize(d)) { - if ( outRegion.GetIndex( d ) < reqRegion.GetIndex( d ) || outRegion.GetIndex( d ) + outRegion.GetSize( d ) > reqRegion.GetIndex( d ) + reqRegion.GetSize( d ) ) - { - withinReq = false; - break; - } + withinReq = false; + break; } + } -#if !defined( __wasi__ ) && !defined( __EMSCRIPTEN__ ) +#if !defined(__wasi__) && !defined(__EMSCRIPTEN__) static std::mutex mutexLock; #endif - if ( withinReq ) // else we should not write it - { - seqIt.GoToBegin(); - // writing through one RLEImage iterator invalidates all the others - // so this whole writing loop needs to be serialized -#if !defined( __wasi__ ) && !defined( __EMSCRIPTEN__ ) - std::lock_guard< std::mutex > mutexHolder( mutexLock ); + if (withinReq) // else we should not write it + { + seqIt.GoToBegin(); + // writing through one RLEImage iterator invalidates all the others + // so this whole writing loop needs to be serialized +#if !defined(__wasi__) && !defined(__EMSCRIPTEN__) + std::lock_guard mutexHolder(mutexLock); #endif - ImageRegionIterator< TImage > outIt( out, outRegion ); - while ( !outIt.IsAtEnd() ) - { - if ( seqIt.Get() && outIt.Get() < label ) - { - outIt.Set( label ); - } - ++seqIt; - ++outIt; - } - } // iterator destroyed here + ImageRegionIterator outIt(out, outRegion); + while (!outIt.IsAtEnd()) + { + if (seqIt.Get() && outIt.Get() < label) + { + outIt.Set(label); + } + ++seqIt; + ++outIt; + } + } // iterator destroyed here // recurse if needed - if ( itk::Math::abs( i - j ) > 2 ) + if (itk::Math::abs(i - j) > 2) + { + PixelList regionIDs; + regionIDs.push_back(1); + + int iReq = i < reqRegion.GetIndex(axis) + ? -1 + : (i > reqRegion.GetIndex(axis) + IndexValueType(reqRegion.GetSize(axis)) ? +1 : 0); + int jReq = j < reqRegion.GetIndex(axis) + ? -1 + : (j > reqRegion.GetIndex(axis) + IndexValueType(reqRegion.GetSize(axis)) ? +1 : 0); + int mReq = mid < reqRegion.GetIndex(axis) + ? -1 + : (mid > reqRegion.GetIndex(axis) + IndexValueType(reqRegion.GetSize(axis)) ? +1 : 0); + bool first = itk::Math::abs(i - mid) > 1 && itk::Math::abs(iReq + mReq) <= 1; // i-mid? + bool second = itk::Math::abs(j - mid) > 1 && itk::Math::abs(jReq + mReq) <= 1; // j-mid? + + if (first) { - PixelList regionIDs; - regionIDs.push_back( 1 ); - - int iReq = i < reqRegion.GetIndex( axis ) ? -1 : ( i > reqRegion.GetIndex( axis ) + IndexValueType( reqRegion.GetSize( axis ) ) ? +1 : 0 ); - int jReq = j < reqRegion.GetIndex( axis ) ? -1 : ( j > reqRegion.GetIndex( axis ) + IndexValueType( reqRegion.GetSize( axis ) ) ? +1 : 0 ); - int mReq = mid < reqRegion.GetIndex( axis ) ? -1 : ( mid > reqRegion.GetIndex( axis ) + IndexValueType( reqRegion.GetSize( axis ) ) ? +1 : 0 ); - bool first = itk::Math::abs( i - mid ) > 1 && itk::Math::abs( iReq + mReq ) <= 1; // i-mid? - bool second = itk::Math::abs( j - mid ) > 1 && itk::Math::abs( jReq + mReq ) <= 1; // j-mid? - - if ( first ) - { - Interpolate1to1( axis, out, label, i, mid, iConn, iRegionId, midConn, 1, iTrans, true ); - } - if ( second ) - { - Interpolate1to1( axis, out, label, j, mid, jConn, jRegionId, midConn, 1, jTrans, true ); - } + Interpolate1to1(axis, out, label, i, mid, iConn, iRegionId, midConn, 1, iTrans, true); + } + if (second) + { + Interpolate1to1(axis, out, label, j, mid, jConn, jRegionId, midConn, 1, jTrans, true); } + } } // >::Interpolate1to1 -template< typename TImage > +template class MatchesID { typename TImage::PixelType m_ID; public: MatchesID() = default; - MatchesID( typename TImage::PixelType id ) - : m_ID( id ) - { - } + MatchesID(typename TImage::PixelType id) + : m_ID(id) + {} bool - operator!=( const MatchesID& other ) + operator!=(const MatchesID & other) { return m_ID != other.m_ID; } bool - operator==( const MatchesID& other ) + operator==(const MatchesID & other) { return m_ID == other.m_ID; } inline bool - operator()( const typename TImage::PixelType& val ) const + operator()(const typename TImage::PixelType & val) const { return val == m_ID; } }; -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::Interpolate1toN( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, - typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds, - const typename SliceType::IndexType& translation ) +MorphologicalContourInterpolator::Interpolate1toN(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds, + const typename SliceType::IndexType & translation) { // first convert iConn into binary mask - MatchesID< TImage > matchesID( iRegionId ); + MatchesID matchesID(iRegionId); - using CastType = UnaryFunctorImageFilter< SliceType, BoolSliceType, MatchesID< TImage > >; + using CastType = UnaryFunctorImageFilter>; typename CastType::Pointer caster = CastType::New(); - caster->SetNumberOfWorkUnits( 1 ); // excessive threading is counterproductive - caster->SetFunctor( matchesID ); - caster->SetInput( iConn ); + caster->SetNumberOfWorkUnits(1); // excessive threading is counterproductive + caster->SetFunctor(matchesID); + caster->SetInput(iConn); caster->Update(); typename BoolSliceType::Pointer mask = caster->GetOutput(); @@ -810,753 +849,791 @@ MorphologicalContourInterpolator< TImage >::Interpolate1toN( int axis, TImage* o iRegion = iConn->GetLargestPossibleRegion(); jRegion = jConn->GetLargestPossibleRegion(); newjRegion = jRegion; - newjRegion.SetSize( iRegion.GetSize() ); + newjRegion.SetSize(iRegion.GetSize()); // construct n empty images - std::vector< typename BoolSliceType::Pointer > blobs; - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - typename BoolSliceType::Pointer temp = BoolSliceType::New(); - temp->CopyInformation( jConn ); - temp->SetRegions( iRegion ); - temp->Allocate( true ); - blobs.push_back( temp ); - } + std::vector blobs; + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + typename BoolSliceType::Pointer temp = BoolSliceType::New(); + temp->CopyInformation(jConn); + temp->SetRegions(iRegion); + temp->Allocate(true); + blobs.push_back(temp); + } // fill the n images with intersections - these are seeds typename SliceType::Pointer belongs = SliceType::New(); - belongs->CopyInformation( mask ); - belongs->SetRegions( iRegion ); - belongs->Allocate( true ); // initialize to zero (false) - ImageRegionIterator< SliceType > belongIt( belongs, iRegion ); - IntersectionRegions( translation, iRegion, jRegion ); - ImageRegionConstIterator< BoolSliceType > maskIt( mask, iRegion ); - ImageRegionConstIteratorWithIndex< SliceType > jIt( jConn, jRegion ); - ImageRegionIterator< SliceType > belongInit( belongs, iRegion ); + belongs->CopyInformation(mask); + belongs->SetRegions(iRegion); + belongs->Allocate(true); // initialize to zero (false) + ImageRegionIterator belongIt(belongs, iRegion); + IntersectionRegions(translation, iRegion, jRegion); + ImageRegionConstIterator maskIt(mask, iRegion); + ImageRegionConstIteratorWithIndex jIt(jConn, jRegion); + ImageRegionIterator belongInit(belongs, iRegion); // convert jConn into n blobs, translating them into the index space of iConn - while ( !maskIt.IsAtEnd() ) + while (!maskIt.IsAtEnd()) + { + if (maskIt.Get()) { - if ( maskIt.Get() ) - { - typename TImage::PixelType jVal = jIt.Get(); - auto res = std::find( jRegionIds.begin(), jRegionIds.end(), jVal ); - if ( res != jRegionIds.end() ) - { - blobs[res - jRegionIds.begin()]->SetPixel( maskIt.GetIndex(), true ); - belongInit.Set( res - jRegionIds.begin() + 1 ); - } - } - ++maskIt; - ++jIt; - ++belongInit; + typename TImage::PixelType jVal = jIt.Get(); + auto res = std::find(jRegionIds.begin(), jRegionIds.end(), jVal); + if (res != jRegionIds.end()) + { + blobs[res - jRegionIds.begin()]->SetPixel(maskIt.GetIndex(), true); + belongInit.Set(res - jRegionIds.begin() + 1); + } } + ++maskIt; + ++jIt; + ++belongInit; + } // prepare dilation filter iRegion = iConn->GetLargestPossibleRegion(); // expand to full i image - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - blobs[x]->SetRegions( iRegion ); - } - ImageRegionConstIterator< BoolSliceType > maskIt2( mask, iRegion ); - ImageRegionConstIteratorWithIndex< BoolSliceType > jIt2( blobs[0], iRegion ); + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + blobs[x]->SetRegions(iRegion); + } + ImageRegionConstIterator maskIt2(mask, iRegion); + ImageRegionConstIteratorWithIndex jIt2(blobs[0], iRegion); bool hollowedMaskEmpty; do // while hollowed mask is not empty + { + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + blobs[x] = Dilate1(blobs[x], mask); + blobs[x]->DisconnectPipeline(); + } + + hollowedMaskEmpty = true; + maskIt2.GoToBegin(); + jIt2.GoToBegin(); + belongIt.GoToBegin(); + while (!maskIt2.IsAtEnd()) // hollow out the big mask with dilated seeds while avoiding conflicts { - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) + if (maskIt2.Get()) + { + if (!belongIt.Get()) { - blobs[x] = Dilate1( blobs[x], mask ); - blobs[x]->DisconnectPipeline(); - } + unsigned x = 0; + while (x < jRegionIds.size()) + { + if (blobs[x]->GetPixel(jIt2.GetIndex())) + { + break; + } + ++x; + } - hollowedMaskEmpty = true; - maskIt2.GoToBegin(); - jIt2.GoToBegin(); - belongIt.GoToBegin(); - while ( !maskIt2.IsAtEnd() ) // hollow out the big mask with dilated seeds while avoiding conflicts + if (x < jRegionIds.size()) // covered by a blob, hollow it out + { + belongIt.Set(x + 1); + for (x++; x < jRegionIds.size(); x++) + { + // pixel does not belong to this blob + blobs[x]->SetPixel(jIt2.GetIndex(), false); + } + } + else // keep it + { + hollowedMaskEmpty = false; + } + } + else // the pixel already belongs to some blob { - if ( maskIt2.Get() ) + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + if (unsigned(belongIt.Get()) != x + 1) { - if ( !belongIt.Get() ) - { - unsigned x = 0; - while ( x < jRegionIds.size() ) - { - if ( blobs[x]->GetPixel( jIt2.GetIndex() ) ) - { - break; - } - ++x; - } - - if ( x < jRegionIds.size() ) // covered by a blob, hollow it out - { - belongIt.Set( x + 1 ); - for ( x++; x < jRegionIds.size(); x++ ) - { - // pixel does not belong to this blob - blobs[x]->SetPixel( jIt2.GetIndex(), false ); - } - } - else // keep it - { - hollowedMaskEmpty = false; - } - } - else // the pixel already belongs to some blob - { - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - if ( unsigned( belongIt.Get() ) != x + 1 ) - { - // pixel does not belong to this blob - blobs[x]->SetPixel( jIt2.GetIndex(), false ); - } - } - } + // pixel does not belong to this blob + blobs[x]->SetPixel(jIt2.GetIndex(), false); } - ++maskIt2; - ++jIt2; - ++belongIt; + } } + } + ++maskIt2; + ++jIt2; + ++belongIt; } - while ( !hollowedMaskEmpty ); + } while (!hollowedMaskEmpty); blobs.clear(); // deallocates the images // convert the belongs into n Conn-style images - std::vector< typename SliceType::Pointer > conns; - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - typename SliceType::Pointer temp2 = SliceType::New(); - temp2->CopyInformation( iConn ); - temp2->SetRegions( iConn->GetLargestPossibleRegion() ); - temp2->Allocate( true ); - conns.push_back( temp2 ); - } - ImageRegionConstIteratorWithIndex< SliceType > belongIt2( belongs, iRegion ); - while ( !belongIt2.IsAtEnd() ) + std::vector conns; + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + typename SliceType::Pointer temp2 = SliceType::New(); + temp2->CopyInformation(iConn); + temp2->SetRegions(iConn->GetLargestPossibleRegion()); + temp2->Allocate(true); + conns.push_back(temp2); + } + ImageRegionConstIteratorWithIndex belongIt2(belongs, iRegion); + while (!belongIt2.IsAtEnd()) + { + const typename SliceType::PixelType belong = belongIt2.Get(); + if (belong > 0) { - const typename SliceType::PixelType belong = belongIt2.Get(); - if ( belong > 0 ) - { - conns[belong - 1]->SetPixel( belongIt2.GetIndex(), iRegionId ); - } - ++belongIt2; + conns[belong - 1]->SetPixel(belongIt2.GetIndex(), iRegionId); } + ++belongIt2; + } // make n 1-to-1 interpolations - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - Interpolate1to1( axis, out, label, i, j, conns[x], iRegionId, jConn, jRegionIds[x], translation, false ); - } + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + Interpolate1to1(axis, out, label, i, j, conns[x], iRegionId, jConn, jRegionIds[x], translation, false); + } } // >::Interpolate1toN -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::TranslateImage( typename SliceType::Pointer& image, const typename SliceType::IndexType& translation, typename SliceType::RegionType newRegion ) -> +MorphologicalContourInterpolator::TranslateImage(typename SliceType::Pointer & image, + const typename SliceType::IndexType & translation, + typename SliceType::RegionType newRegion) -> typename SliceType::Pointer { typename SliceType::Pointer result = SliceType::New(); - result->CopyInformation( image ); - result->SetRegions( newRegion ); - result->Allocate( true ); // initialize to zero (false) + result->CopyInformation(image); + result->SetRegions(newRegion); + result->Allocate(true); // initialize to zero (false) typename SliceType::RegionType inRegion = image->GetLargestPossibleRegion(); - IntersectionRegions( translation, inRegion, newRegion ); - ImageAlgorithm::Copy< SliceType, SliceType >( image.GetPointer(), result.GetPointer(), inRegion, newRegion ); + IntersectionRegions(translation, inRegion, newRegion); + ImageAlgorithm::Copy(image.GetPointer(), result.GetPointer(), inRegion, newRegion); return result; } -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::IntersectionRegions( const typename SliceType::IndexType& translation, typename SliceType::RegionType& iRegion, typename SliceType::RegionType& jRegion ) +MorphologicalContourInterpolator::IntersectionRegions(const typename SliceType::IndexType & translation, + typename SliceType::RegionType & iRegion, + typename SliceType::RegionType & jRegion) { typename SliceType::IndexType iBegin = iRegion.GetIndex(); typename SliceType::IndexType jBegin = jRegion.GetIndex(); - for ( IdentifierType d = 0; d < SliceType::ImageDimension; d++ ) - { - IndexValueType iSize = iRegion.GetSize( d ); - IndexValueType jSize = jRegion.GetSize( d ); - iBegin[d] += translation[d]; - IndexValueType t = std::max( iBegin[d], jBegin[d] ); - iRegion.SetSize( d, std::min( IndexValueType( iSize ) - ( t - iBegin[d] ), IndexValueType( jSize ) - ( t - jBegin[d] ) ) ); - iRegion.SetIndex( d, t - translation[d] ); - jRegion.SetIndex( d, t ); - } - jRegion.SetSize( iRegion.GetSize() ); // size is the same + for (IdentifierType d = 0; d < SliceType::ImageDimension; d++) + { + IndexValueType iSize = iRegion.GetSize(d); + IndexValueType jSize = jRegion.GetSize(d); + iBegin[d] += translation[d]; + IndexValueType t = std::max(iBegin[d], jBegin[d]); + iRegion.SetSize(d, std::min(IndexValueType(iSize) - (t - iBegin[d]), IndexValueType(jSize) - (t - jBegin[d]))); + iRegion.SetIndex(d, t - translation[d]); + jRegion.SetIndex(d, t); + } + jRegion.SetSize(iRegion.GetSize()); // size is the same } -template< typename TImage > +template IdentifierType -MorphologicalContourInterpolator< TImage >::Intersection( typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds, - const typename SliceType::IndexType& translation ) +MorphologicalContourInterpolator::Intersection(typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds, + const typename SliceType::IndexType & translation) { typename SliceType::RegionType iRegion, jRegion; iRegion = iConn->GetLargestPossibleRegion(); jRegion = jConn->GetLargestPossibleRegion(); - IntersectionRegions( translation, iRegion, jRegion ); + IntersectionRegions(translation, iRegion, jRegion); - std::vector< IdentifierType > counts( jRegionIds.size() ); - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) - { - counts[x] = 0; - } - ImageRegionConstIterator< SliceType > iIt( iConn, iRegion ); - ImageRegionConstIterator< SliceType > jIt( jConn, jRegion ); - while ( !iIt.IsAtEnd() ) + std::vector counts(jRegionIds.size()); + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + counts[x] = 0; + } + ImageRegionConstIterator iIt(iConn, iRegion); + ImageRegionConstIterator jIt(jConn, jRegion); + while (!iIt.IsAtEnd()) + { + if (iIt.Get() == iRegionId) { - if ( iIt.Get() == iRegionId ) - { - typename TImage::PixelType jVal = jIt.Get(); - auto res = std::find( jRegionIds.begin(), jRegionIds.end(), jVal ); - if ( res != jRegionIds.end() ) - { - ++counts[res - jRegionIds.begin()]; - } - } - ++iIt; - ++jIt; + typename TImage::PixelType jVal = jIt.Get(); + auto res = std::find(jRegionIds.begin(), jRegionIds.end(), jVal); + if (res != jRegionIds.end()) + { + ++counts[res - jRegionIds.begin()]; + } } + ++iIt; + ++jIt; + } IdentifierType sum = 0; - for ( unsigned x = 0; x < jRegionIds.size(); x++ ) + for (unsigned x = 0; x < jRegionIds.size(); x++) + { + if (counts[x] == 0) { - if ( counts[x] == 0 ) - { - return 0; // iConn must intersect all subregions of jConn - } - sum += counts[x]; + return 0; // iConn must intersect all subregions of jConn } + sum += counts[x]; + } return sum; } // >::Intersection -template< typename TImage > +template IdentifierType -MorphologicalContourInterpolator< TImage >::CardinalSymmetricDifference( typename BoolSliceType::Pointer& iShape, typename BoolSliceType::Pointer& jShape ) +MorphologicalContourInterpolator::CardinalSymmetricDifference(typename BoolSliceType::Pointer & iShape, + typename BoolSliceType::Pointer & jShape) { - typename BoolSliceType::RegionType region = iShape->GetLargestPossibleRegion(); - IdentifierType count = 0; - ImageRegionConstIterator< BoolSliceType > iIt( iShape, region ); - ImageRegionConstIterator< BoolSliceType > jIt( jShape, region ); - while ( !iIt.IsAtEnd() ) + typename BoolSliceType::RegionType region = iShape->GetLargestPossibleRegion(); + IdentifierType count = 0; + ImageRegionConstIterator iIt(iShape, region); + ImageRegionConstIterator jIt(jShape, region); + while (!iIt.IsAtEnd()) + { + if (iIt.Get() != jIt.Get()) { - if ( iIt.Get() != jIt.Get() ) - { - count++; - } - ++iIt; - ++jIt; + count++; } + ++iIt; + ++jIt; + } return count; } -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::Centroid( typename SliceType::Pointer& conn, const PixelList& regionIds ) -> typename SliceType::IndexType +MorphologicalContourInterpolator::Centroid(typename SliceType::Pointer & conn, const PixelList & regionIds) -> + typename SliceType::IndexType { - ImageRegionConstIteratorWithIndex< SliceType > it( conn, conn->GetLargestPossibleRegion() ); - IndexValueType ind[SliceType::ImageDimension] = { 0 }; // all components are initialized to zero - IdentifierType pixelCount = 0; - while ( !it.IsAtEnd() ) + ImageRegionConstIteratorWithIndex it(conn, conn->GetLargestPossibleRegion()); + IndexValueType ind[SliceType::ImageDimension] = { 0 }; // all components are initialized to zero + IdentifierType pixelCount = 0; + while (!it.IsAtEnd()) + { + typename TImage::PixelType val = it.Get(); + if (val) { - typename TImage::PixelType val = it.Get(); - if ( val ) + auto res = std::find(regionIds.begin(), regionIds.end(), val); + if (res != regionIds.end()) + { + ++pixelCount; + typename SliceType::IndexType pInd = it.GetIndex(); + for (unsigned d = 0; d < SliceType::ImageDimension; d++) { - auto res = std::find( regionIds.begin(), regionIds.end(), val ); - if ( res != regionIds.end() ) - { - ++pixelCount; - typename SliceType::IndexType pInd = it.GetIndex(); - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - ind[d] += pInd[d]; - } - } + ind[d] += pInd[d]; } - ++it; + } } + ++it; + } typename SliceType::IndexType retVal; - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - retVal[d] = ind[d] / pixelCount; - } + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + retVal[d] = ind[d] / pixelCount; + } return retVal; } // >::Centroid -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::Align( typename SliceType::Pointer& iConn, typename TImage::PixelType iRegionId, typename SliceType::Pointer& jConn, const PixelList& jRegionIds ) -> - typename SliceType::IndexType +MorphologicalContourInterpolator::Align(typename SliceType::Pointer & iConn, + typename TImage::PixelType iRegionId, + typename SliceType::Pointer & jConn, + const PixelList & jRegionIds) -> typename SliceType::IndexType { // calculate centroids PixelList iRegionIds; - iRegionIds.push_back( iRegionId ); - typename SliceType::IndexType iCentroid = Centroid( iConn, iRegionIds ); - typename SliceType::IndexType jCentroid = Centroid( jConn, jRegionIds ); + iRegionIds.push_back(iRegionId); + typename SliceType::IndexType iCentroid = Centroid(iConn, iRegionIds); + typename SliceType::IndexType jCentroid = Centroid(jConn, jRegionIds); typename SliceType::IndexType ind; - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - ind[d] = jCentroid[d] - iCentroid[d]; - } + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + ind[d] = jCentroid[d] - iCentroid[d]; + } // construct an image with all possible translations typename SliceType::RegionType searchRegion; typename SliceType::RegionType iLPR = iConn->GetLargestPossibleRegion(); typename SliceType::RegionType jLPR = jConn->GetLargestPossibleRegion(); - for ( IdentifierType d = 0; d < SliceType::ImageDimension; d++ ) - { - searchRegion.SetIndex( d, jLPR.GetIndex()[d] - iLPR.GetIndex()[d] - iLPR.GetSize( d ) + 1 ); - searchRegion.SetSize( d, iLPR.GetSize( d ) + jLPR.GetSize( d ) - 1 ); - } + for (IdentifierType d = 0; d < SliceType::ImageDimension; d++) + { + searchRegion.SetIndex(d, jLPR.GetIndex()[d] - iLPR.GetIndex()[d] - iLPR.GetSize(d) + 1); + searchRegion.SetSize(d, iLPR.GetSize(d) + jLPR.GetSize(d) - 1); + } typename BoolSliceType::Pointer searched = BoolSliceType::New(); - searched->SetRegions( searchRegion ); - searched->Allocate( true ); // initialize to zero (false) + searched->SetRegions(searchRegion); + searched->Allocate(true); // initialize to zero (false) // breadth first search starting from centroid, implicitly: // when intersection scores are equal, chooses the one closer to centroid - std::queue< typename SliceType::IndexType > uncomputed; - typename SliceType::IndexType t0; - t0.Fill( 0 ); - uncomputed.push( t0 ); // no translation - guaranteed to find a non-zero intersection - uncomputed.push( ind ); // this introduces movement, and possibly has the same score - searched->SetPixel( t0, true ); - searched->SetPixel( ind, true ); + std::queue uncomputed; + typename SliceType::IndexType t0; + t0.Fill(0); + uncomputed.push(t0); // no translation - guaranteed to find a non-zero intersection + uncomputed.push(ind); // this introduces movement, and possibly has the same score + searched->SetPixel(t0, true); + searched->SetPixel(ind, true); IdentifierType score, maxScore = 0; typename SliceType::IndexType bestIndex; IdentifierType iter = 0; - IdentifierType minIter = std::min( m_MinAlignIters, searchRegion.GetNumberOfPixels() ); - IdentifierType maxIter = std::max( m_MaxAlignIters, static_cast< IdentifierType >( std::sqrt( static_cast< double >( searchRegion.GetNumberOfPixels() ) ) ) ); + IdentifierType minIter = std::min(m_MinAlignIters, searchRegion.GetNumberOfPixels()); + IdentifierType maxIter = std::max( + m_MaxAlignIters, static_cast(std::sqrt(static_cast(searchRegion.GetNumberOfPixels())))); - while ( !uncomputed.empty() ) + while (!uncomputed.empty()) + { + ind = uncomputed.front(); + uncomputed.pop(); + score = Intersection(iConn, iRegionId, jConn, jRegionIds, ind); + if (score > maxScore) { - ind = uncomputed.front(); - uncomputed.pop(); - score = Intersection( iConn, iRegionId, jConn, jRegionIds, ind ); - if ( score > maxScore ) + maxScore = score; + bestIndex = ind; + } + + // we breadth this search + if (!m_HeuristicAlignment || maxScore == 0 || iter <= minIter || (score > maxScore * 0.9 && iter <= maxIter)) + { + for (unsigned d = 0; d < SliceType::ImageDimension; d++) + { + ind[d] -= 1; // "left" + if (searchRegion.IsInside(ind) && !searched->GetPixel(ind)) { - maxScore = score; - bestIndex = ind; + uncomputed.push(ind); + searched->SetPixel(ind, true); + ++iter; } - - // we breadth this search - if ( !m_HeuristicAlignment || maxScore == 0 || iter <= minIter || ( score > maxScore * 0.9 && iter <= maxIter ) ) + ind[d] += 2; // "right" + if (searchRegion.IsInside(ind) && !searched->GetPixel(ind)) { - for ( unsigned d = 0; d < SliceType::ImageDimension; d++ ) - { - ind[d] -= 1; // "left" - if ( searchRegion.IsInside( ind ) && !searched->GetPixel( ind ) ) - { - uncomputed.push( ind ); - searched->SetPixel( ind, true ); - ++iter; - } - ind[d] += 2; // "right" - if ( searchRegion.IsInside( ind ) && !searched->GetPixel( ind ) ) - { - uncomputed.push( ind ); - searched->SetPixel( ind, true ); - ++iter; - } - ind[d] -= 1; // return to initial - } + uncomputed.push(ind); + searched->SetPixel(ind, true); + ++iter; } + ind[d] -= 1; // return to initial + } } + } return bestIndex; } // >::Align -template< typename TImage > +template auto -MorphologicalContourInterpolator< TImage >::RegionedConnectedComponents( const typename TImage::RegionType& region, typename TImage::PixelType label, IdentifierType& objectCount ) -> +MorphologicalContourInterpolator::RegionedConnectedComponents(const typename TImage::RegionType & region, + typename TImage::PixelType label, + IdentifierType & objectCount) -> typename SliceType::Pointer { - m_RoI->SetExtractionRegion( region ); - m_RoI->SetInput( this->GetInput() ); - m_Binarizer->SetLowerThreshold( label ); - m_Binarizer->SetUpperThreshold( label ); + m_RoI->SetExtractionRegion(region); + m_RoI->SetInput(this->GetInput()); + m_Binarizer->SetLowerThreshold(label); + m_Binarizer->SetUpperThreshold(label); m_ConnectedComponents->Update(); objectCount = m_ConnectedComponents->GetObjectCount(); return m_ConnectedComponents->GetOutput(); } -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::InterpolateBetweenTwo( int axis, TImage* out, typename TImage::PixelType label, typename TImage::IndexValueType i, typename TImage::IndexValueType j, - typename SliceType::Pointer& iconn, typename SliceType::Pointer& jconn ) +MorphologicalContourInterpolator::InterpolateBetweenTwo(int axis, + TImage * out, + typename TImage::PixelType label, + typename TImage::IndexValueType i, + typename TImage::IndexValueType j, + typename SliceType::Pointer & iconn, + typename SliceType::Pointer & jconn) { // go through comparison image and create correspondence pairs - using PairSet = std::set< std::pair< typename TImage::PixelType, typename TImage::PixelType > >; - PairSet pairs, unwantedPairs, uncleanPairs; - typename SliceType::RegionType ri = iconn->GetRequestedRegion(); - typename SliceType::RegionType rj = jconn->GetRequestedRegion(); - ImageRegionConstIterator< SliceType > iti( iconn, ri ); - ImageRegionConstIterator< SliceType > itj( jconn, rj ); - while ( !iti.IsAtEnd() ) + using PairSet = std::set>; + PairSet pairs, unwantedPairs, uncleanPairs; + typename SliceType::RegionType ri = iconn->GetRequestedRegion(); + typename SliceType::RegionType rj = jconn->GetRequestedRegion(); + ImageRegionConstIterator iti(iconn, ri); + ImageRegionConstIterator itj(jconn, rj); + while (!iti.IsAtEnd()) + { + if (iti.Get() != 0 || itj.Get() != 0) { - if ( iti.Get() != 0 || itj.Get() != 0 ) - { - uncleanPairs.insert( std::make_pair( iti.Get(), itj.Get() ) ); - // std::cout << " iti:" << iti.GetIndex() << iti.Get() << - // " itj:" << itj.GetIndex() << itj.Get() << std::endl; - if ( iti.Get() != 0 && itj.Get() != 0 ) - { - unwantedPairs.insert( std::make_pair( 0, itj.Get() ) ); - unwantedPairs.insert( std::make_pair( iti.Get(), 0 ) ); - } - } - ++iti; - ++itj; + uncleanPairs.insert(std::make_pair(iti.Get(), itj.Get())); + // std::cout << " iti:" << iti.GetIndex() << iti.Get() << + // " itj:" << itj.GetIndex() << itj.Get() << std::endl; + if (iti.Get() != 0 && itj.Get() != 0) + { + unwantedPairs.insert(std::make_pair(0, itj.Get())); + unwantedPairs.insert(std::make_pair(iti.Get(), 0)); + } } + ++iti; + ++itj; + } - std::set_difference( uncleanPairs.begin(), uncleanPairs.end(), unwantedPairs.begin(), unwantedPairs.end(), std::inserter( pairs, pairs.end() ) ); + std::set_difference(uncleanPairs.begin(), + uncleanPairs.end(), + unwantedPairs.begin(), + unwantedPairs.end(), + std::inserter(pairs, pairs.end())); // first do extrapolation for components without overlaps auto p = pairs.begin(); - while ( p != pairs.end() ) + while (p != pairs.end()) + { + if (p->second == 0) { - if ( p->second == 0 ) - { - if ( m_UseExtrapolation ) - { - Extrapolate( axis, out, label, i, j, iconn, p->first ); - } - pairs.erase( p++ ); - } - else if ( p->first == 0 ) - { - if ( m_UseExtrapolation ) - { - Extrapolate( axis, out, label, j, i, jconn, p->second ); - } - pairs.erase( p++ ); - } - else - { - ++p; - } + if (m_UseExtrapolation) + { + Extrapolate(axis, out, label, i, j, iconn, p->first); + } + pairs.erase(p++); } + else if (p->first == 0) + { + if (m_UseExtrapolation) + { + Extrapolate(axis, out, label, j, i, jconn, p->second); + } + pairs.erase(p++); + } + else + { + ++p; + } + } // count ocurrances of each component - using CountMap = std::map< typename TImage::PixelType, IdentifierType >; + using CountMap = std::map; CountMap iCounts, jCounts; - for ( p = pairs.begin(); p != pairs.end(); ++p ) - { - iCounts[p->first]++; - jCounts[p->second]++; - } + for (p = pairs.begin(); p != pairs.end(); ++p) + { + iCounts[p->first]++; + jCounts[p->second]++; + } // now handle 1 to 1 correspondences p = pairs.begin(); - while ( p != pairs.end() ) + while (p != pairs.end()) + { + if (iCounts[p->first] == 1 && jCounts[p->second] == 1) { - if ( iCounts[p->first] == 1 && jCounts[p->second] == 1 ) - { - PixelList regionIDs; - regionIDs.push_back( p->second ); - typename SliceType::IndexType translation = Align( iconn, p->first, jconn, regionIDs ); - Interpolate1to1( axis, out, label, i, j, iconn, p->first, jconn, p->second, translation, false ); - iCounts.erase( p->first ); - jCounts.erase( p->second ); - pairs.erase( p++ ); - } - else - { - ++p; - } + PixelList regionIDs; + regionIDs.push_back(p->second); + typename SliceType::IndexType translation = Align(iconn, p->first, jconn, regionIDs); + Interpolate1to1(axis, out, label, i, j, iconn, p->first, jconn, p->second, translation, false); + iCounts.erase(p->first); + jCounts.erase(p->second); + pairs.erase(p++); } + else + { + ++p; + } + } - PixelList regionIDs( pairs.size() ); // preallocate + PixelList regionIDs(pairs.size()); // preallocate // now do 1-to-N and M-to-1 cases p = pairs.begin(); - while ( p != pairs.end() ) - { - regionIDs.clear(); + while (p != pairs.end()) + { + regionIDs.clear(); - if ( iCounts[p->first] == 1 ) // M-to-1 + if (iCounts[p->first] == 1) // M-to-1 + { + for (auto rest = pairs.begin(); rest != pairs.end(); ++rest) + { + if (rest->second == p->second) { - for ( auto rest = pairs.begin(); rest != pairs.end(); ++rest ) - { - if ( rest->second == p->second ) - { - regionIDs.push_back( rest->first ); - } - } + regionIDs.push_back(rest->first); + } + } - typename SliceType::IndexType translation = Align( jconn, p->second, iconn, regionIDs ); - Interpolate1toN( axis, out, label, j, i, jconn, p->second, iconn, regionIDs, translation ); + typename SliceType::IndexType translation = Align(jconn, p->second, iconn, regionIDs); + Interpolate1toN(axis, out, label, j, i, jconn, p->second, iconn, regionIDs, translation); - auto rest = pairs.begin(); - while ( rest != pairs.end() ) - { - if ( rest != p && rest->second == p->second ) - { - --iCounts[rest->first]; - --jCounts[rest->second]; - pairs.erase( rest++ ); - } - else - { - ++rest; - } - } - - --iCounts[p->first]; - --jCounts[p->second]; - pairs.erase( p++ ); - } // M-to-1 - else if ( jCounts[p->second] == 1 ) // 1-to-N + auto rest = pairs.begin(); + while (rest != pairs.end()) + { + if (rest != p && rest->second == p->second) { - for ( auto rest = pairs.begin(); rest != pairs.end(); ++rest ) - { - if ( rest->first == p->first ) - { - regionIDs.push_back( rest->second ); - } - } + --iCounts[rest->first]; + --jCounts[rest->second]; + pairs.erase(rest++); + } + else + { + ++rest; + } + } - typename SliceType::IndexType translation = Align( iconn, p->first, jconn, regionIDs ); - Interpolate1toN( axis, out, label, i, j, iconn, p->first, jconn, regionIDs, translation ); + --iCounts[p->first]; + --jCounts[p->second]; + pairs.erase(p++); + } // M-to-1 + else if (jCounts[p->second] == 1) // 1-to-N + { + for (auto rest = pairs.begin(); rest != pairs.end(); ++rest) + { + if (rest->first == p->first) + { + regionIDs.push_back(rest->second); + } + } - auto rest = pairs.begin(); - ++rest; - while ( rest != pairs.end() ) - { - if ( rest != p && rest->first == p->first ) - { - --iCounts[rest->first]; - --jCounts[rest->second]; - pairs.erase( rest++ ); - } - else - { - ++rest; - } - } + typename SliceType::IndexType translation = Align(iconn, p->first, jconn, regionIDs); + Interpolate1toN(axis, out, label, i, j, iconn, p->first, jconn, regionIDs, translation); - --iCounts[p->first]; - --jCounts[p->second]; - pairs.erase( p++ ); - } // 1-to-N - else + auto rest = pairs.begin(); + ++rest; + while (rest != pairs.end()) + { + if (rest != p && rest->first == p->first) { - ++p; + --iCounts[rest->first]; + --jCounts[rest->second]; + pairs.erase(rest++); } - } // 1-to-N and M-to-1 + else + { + ++rest; + } + } + + --iCounts[p->first]; + --jCounts[p->second]; + pairs.erase(p++); + } // 1-to-N + else + { + ++p; + } + } // 1-to-N and M-to-1 // only M-to-N correspondences remain // we turn each M-to-N case into m 1-to-N cases p = pairs.begin(); - while ( p != pairs.end() ) + while (p != pairs.end()) + { + regionIDs.clear(); + for (auto rest = p; rest != pairs.end(); ++rest) { - regionIDs.clear(); - for ( auto rest = p; rest != pairs.end(); ++rest ) - { - if ( rest->first == p->first ) - { - regionIDs.push_back( rest->second ); - } - } + if (rest->first == p->first) + { + regionIDs.push_back(rest->second); + } + } - typename SliceType::IndexType translation = Align( iconn, p->first, jconn, regionIDs ); - Interpolate1toN( axis, out, label, i, j, iconn, p->first, jconn, regionIDs, translation ); + typename SliceType::IndexType translation = Align(iconn, p->first, jconn, regionIDs); + Interpolate1toN(axis, out, label, i, j, iconn, p->first, jconn, regionIDs, translation); - auto rest = p; - ++rest; - while ( rest != pairs.end() ) - { - if ( rest->first == p->first ) - { - pairs.erase( rest++ ); - } - else - { - ++rest; - } - } + auto rest = p; + ++rest; + while (rest != pairs.end()) + { + if (rest->first == p->first) + { + pairs.erase(rest++); + } + else + { + ++rest; + } + } - // counts no longer matter, do not waste time deleting them - pairs.erase( p++ ); - } // M-to-N + // counts no longer matter, do not waste time deleting them + pairs.erase(p++); + } // M-to-N } // void MorphologicalContourInterpolator::InterpolateBetweenTwo() -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::InterpolateAlong( int axis, TImage* out, float startProgress, float endProgress ) +MorphologicalContourInterpolator::InterpolateAlong(int axis, + TImage * out, + float startProgress, + float endProgress) { // a list of segments which need to be interpolated - std::vector< SegmentBetweenTwo< TImage > > segments; - typename TImage::RegionType reqRegion = this->GetOutput()->GetRequestedRegion(); - for ( typename LabeledSlicesType::iterator it = m_LabeledSlices[axis].begin(); it != m_LabeledSlices[axis].end(); ++it ) + std::vector> segments; + typename TImage::RegionType reqRegion = this->GetOutput()->GetRequestedRegion(); + for (typename LabeledSlicesType::iterator it = m_LabeledSlices[axis].begin(); it != m_LabeledSlices[axis].end(); ++it) + { + if (m_Label == 0 || m_Label == it->first) // label needs to be interpolated { - if ( m_Label == 0 || m_Label == it->first ) // label needs to be interpolated + auto prev = it->second.begin(); + if (prev == it->second.end()) + { + continue; // nothing to do for this label + } + + typename TImage::RegionType ri = reqRegion; + if (!m_UseCustomSlicePositions) + { + typename BoundingBoxesType::iterator iBB = m_BoundingBoxes.find(it->first); + if (iBB == m_BoundingBoxes.end()) { - auto prev = it->second.begin(); - if ( prev == it->second.end() ) - { - continue; // nothing to do for this label - } - - typename TImage::RegionType ri = reqRegion; - if ( !m_UseCustomSlicePositions ) - { - typename BoundingBoxesType::iterator iBB = m_BoundingBoxes.find( it->first ); - if ( iBB == m_BoundingBoxes.end() ) - { - continue; // this label not present in requested region - } - else - { - ri = iBB->second; - } - } - ri.SetSize( axis, 0 ); - ri.SetIndex( axis, *prev ); - IdentifierType xCount; - typename SliceType::Pointer iconn = this->RegionedConnectedComponents( ri, it->first, xCount ); - iconn->DisconnectPipeline(); - int iReq = *prev < reqRegion.GetIndex( axis ) ? -1 : ( *prev > reqRegion.GetIndex( axis ) + IndexValueType( reqRegion.GetSize( axis ) ) ? +1 : 0 ); - - auto next = it->second.begin(); - for ( ++next; next != it->second.end(); ++next ) - { - typename TImage::RegionType rj = ri; - rj.SetIndex( axis, *next ); - typename SliceType::Pointer jconn = this->RegionedConnectedComponents( rj, it->first, xCount ); - jconn->DisconnectPipeline(); - int jReq = *next < reqRegion.GetIndex( axis ) ? -1 : ( *next > reqRegion.GetIndex( axis ) + IndexValueType( reqRegion.GetSize( axis ) ) ? +1 : 0 ); - - if ( *prev + 1 < *next // only if they are not adjacent slices - && itk::Math::abs( iReq + jReq ) <= 1 ) // and not out of the requested region - // unless they are on opposite ends - { - SegmentBetweenTwo< TImage > s; - s.axis = axis; - s.out = out; - s.label = it->first; - s.i = *prev; - s.j = *next; - s.iconn = iconn; - s.jconn = jconn; - segments.push_back( s ); - } - iconn = jconn; - iReq = jReq; - prev = next; - } + continue; // this label not present in requested region + } + else + { + ri = iBB->second; } + } + ri.SetSize(axis, 0); + ri.SetIndex(axis, *prev); + IdentifierType xCount; + typename SliceType::Pointer iconn = this->RegionedConnectedComponents(ri, it->first, xCount); + iconn->DisconnectPipeline(); + int iReq = *prev < reqRegion.GetIndex(axis) + ? -1 + : (*prev > reqRegion.GetIndex(axis) + IndexValueType(reqRegion.GetSize(axis)) ? +1 : 0); + + auto next = it->second.begin(); + for (++next; next != it->second.end(); ++next) + { + typename TImage::RegionType rj = ri; + rj.SetIndex(axis, *next); + typename SliceType::Pointer jconn = this->RegionedConnectedComponents(rj, it->first, xCount); + jconn->DisconnectPipeline(); + int jReq = *next < reqRegion.GetIndex(axis) + ? -1 + : (*next > reqRegion.GetIndex(axis) + IndexValueType(reqRegion.GetSize(axis)) ? +1 : 0); + + if (*prev + 1 < *next // only if they are not adjacent slices + && itk::Math::abs(iReq + jReq) <= 1) // and not out of the requested region + // unless they are on opposite ends + { + SegmentBetweenTwo s; + s.axis = axis; + s.out = out; + s.label = it->first; + s.i = *prev; + s.j = *next; + s.iconn = iconn; + s.jconn = jconn; + segments.push_back(s); + } + iconn = jconn; + iReq = jReq; + prev = next; + } } + } - ProgressTransformer pt( startProgress, endProgress, this ); - MultiThreaderBase* mt = this->GetMultiThreader(); + ProgressTransformer pt(startProgress, endProgress, this); + MultiThreaderBase * mt = this->GetMultiThreader(); mt->ParallelizeArray( 0, segments.size(), - [&]( SizeValueType ii ) { this->InterpolateBetweenTwo( segments[ii].axis, segments[ii].out, segments[ii].label, segments[ii].i, segments[ii].j, segments[ii].iconn, segments[ii].jconn ); }, - pt.GetProcessObject() ); + [&](SizeValueType ii) { + this->InterpolateBetweenTwo(segments[ii].axis, + segments[ii].out, + segments[ii].label, + segments[ii].i, + segments[ii].j, + segments[ii].iconn, + segments[ii].jconn); + }, + pt.GetProcessObject()); } // >::InterpolateAlong -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::AllocateOutputs() +MorphologicalContourInterpolator::AllocateOutputs() { - using ImageBaseType = ImageBase< TImage::ImageDimension >; + using ImageBaseType = ImageBase; typename ImageBaseType::Pointer outputPtr; - for ( OutputDataObjectIterator it( this ); !it.IsAtEnd(); it++ ) + for (OutputDataObjectIterator it(this); !it.IsAtEnd(); it++) + { + outputPtr = dynamic_cast(it.GetOutput()); + if (outputPtr) { - outputPtr = dynamic_cast< ImageBaseType* >( it.GetOutput() ); - if ( outputPtr ) - { - outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() ); - outputPtr->Allocate( true ); - } + outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion()); + outputPtr->Allocate(true); } + } } -template< typename TImage > +template void -MorphologicalContourInterpolator< TImage >::GenerateData() +MorphologicalContourInterpolator::GenerateData() { typename TImage::ConstPointer m_Input = this->GetInput(); typename TImage::Pointer m_Output = this->GetOutput(); this->AllocateOutputs(); - if ( m_UseCustomSlicePositions ) - { - SliceIndicesType t = m_LabeledSlices; - this->DetermineSliceOrientations(); // calculates bounding boxes - m_LabeledSlices = t; // restore custom positions - } + if (m_UseCustomSlicePositions) + { + SliceIndicesType t = m_LabeledSlices; + this->DetermineSliceOrientations(); // calculates bounding boxes + m_LabeledSlices = t; // restore custom positions + } else - { - this->DetermineSliceOrientations(); - } + { + this->DetermineSliceOrientations(); + } - if ( m_BoundingBoxes.empty() && !m_UseCustomSlicePositions ) - { - ImageAlgorithm::Copy< TImage, TImage >( m_Input.GetPointer(), m_Output.GetPointer(), m_Output->GetRequestedRegion(), m_Output->GetRequestedRegion() ); - return; // no contours detected - } + if (m_BoundingBoxes.empty() && !m_UseCustomSlicePositions) + { + ImageAlgorithm::Copy( + m_Input.GetPointer(), m_Output.GetPointer(), m_Output->GetRequestedRegion(), m_Output->GetRequestedRegion()); + return; // no contours detected + } - if ( m_Axis == -1 ) // interpolate along all axes + if (m_Axis == -1) // interpolate along all axes + { + FixedArray aggregate; + aggregate.Fill(false); + for (unsigned i = 0; i < TImage::ImageDimension; i++) { - FixedArray< bool, TImage::ImageDimension > aggregate; - aggregate.Fill( false ); - for ( unsigned i = 0; i < TImage::ImageDimension; i++ ) + if (this->m_Label == 0) // examine all labels + { + for (unsigned l = 0; l < m_LabeledSlices[i].size(); l++) { - if ( this->m_Label == 0 ) // examine all labels - { - for ( unsigned l = 0; l < m_LabeledSlices[i].size(); l++ ) - { - if ( m_LabeledSlices[i][l].size() > 1 ) - { - aggregate[i] = true; - } - } - } - else // we only care about this label - { - if ( m_LabeledSlices[i][m_Label].size() > 1 ) - { - aggregate[i] = true; - } - } + if (m_LabeledSlices[i][l].size() > 1) + { + aggregate[i] = true; + } } - - for ( unsigned int a = 0; a < TImage::ImageDimension; ++a ) + } + else // we only care about this label + { + if (m_LabeledSlices[i][m_Label].size() > 1) { - if ( aggregate[a] ) - { - this->InterpolateAlong( a, m_Output, a / float( TImage::ImageDimension ), ( a + 1 ) / float( TImage::ImageDimension ) ); - } + aggregate[i] = true; } - } // interpolate along all axes - else // interpolate along the specified axis + } + } + + for (unsigned int a = 0; a < TImage::ImageDimension; ++a) { - this->InterpolateAlong( m_Axis, m_Output, 0.0f, 1.0f ); + if (aggregate[a]) + { + this->InterpolateAlong(a, m_Output, a / float(TImage::ImageDimension), (a + 1) / float(TImage::ImageDimension)); + } } + } // interpolate along all axes + else // interpolate along the specified axis + { + this->InterpolateAlong(m_Axis, m_Output, 0.0f, 1.0f); + } // Overwrites m_Output with non-zeroes from m_Input - ImageRegionIterator< TImage > itO( this->GetOutput(), this->GetOutput()->GetBufferedRegion() ); - ImageRegionConstIterator< TImage > itI( this->GetInput(), this->GetOutput()->GetBufferedRegion() ); - while ( !itI.IsAtEnd() ) + ImageRegionIterator itO(this->GetOutput(), this->GetOutput()->GetBufferedRegion()); + ImageRegionConstIterator itI(this->GetInput(), this->GetOutput()->GetBufferedRegion()); + while (!itI.IsAtEnd()) + { + typename TImage::PixelType val = itI.Get(); + if (val != 0) { - typename TImage::PixelType val = itI.Get(); - if ( val != 0 ) - { - itO.Set( val ); - } - - ++itI; - ++itO; + itO.Set(val); } + + ++itI; + ++itO; + } } // >::GenerateData } // end namespace itk diff --git a/test/dscComparison.cxx b/test/dscComparison.cxx index b0a70ee..e6000bd 100644 --- a/test/dscComparison.cxx +++ b/test/dscComparison.cxx @@ -27,168 +27,172 @@ using TestPixelType = unsigned char; constexpr unsigned int testDim = 3; -using TestImageType = itk::Image< TestPixelType, testDim >; +using TestImageType = itk::Image; TestImageType::Pointer -createSparseCopy( const TestImageType::Pointer& inImage, TestImageType::IndexType nth ) +createSparseCopy(const TestImageType::Pointer & inImage, TestImageType::IndexType nth) { - const TestImageType::RegionType& lpr = inImage->GetLargestPossibleRegion(); + const TestImageType::RegionType & lpr = inImage->GetLargestPossibleRegion(); TestImageType::Pointer outImage = TestImageType::New(); - outImage->CopyInformation( inImage ); - outImage->SetRegions( lpr ); - outImage->Allocate( true ); + outImage->CopyInformation(inImage); + outImage->SetRegions(lpr); + outImage->Allocate(true); - itk::ImageRegionConstIterator< TestImageType > iIt( inImage, lpr ); - itk::ImageRegionIteratorWithIndex< TestImageType > oIt( outImage, lpr ); + itk::ImageRegionConstIterator iIt(inImage, lpr); + itk::ImageRegionIteratorWithIndex oIt(outImage, lpr); - while ( !oIt.IsAtEnd() ) + while (!oIt.IsAtEnd()) + { + const TestPixelType & val = iIt.Get(); + if (val) { - const TestPixelType& val = iIt.Get(); - if ( val ) + const TestImageType::IndexType & ind = oIt.GetIndex(); + bool write = false; + for (unsigned axis = 0; axis < testDim; axis++) + { + if (ind[axis] % nth[axis] == 0) { - const TestImageType::IndexType& ind = oIt.GetIndex(); - bool write = false; - for ( unsigned axis = 0; axis < testDim; axis++ ) - { - if ( ind[axis] % nth[axis] == 0 ) - { - write = true; - break; - } - } - - if ( write ) - { - oIt.Set( val ); - } + write = true; + break; } - ++iIt; - ++oIt; + } + + if (write) + { + oIt.Set(val); + } } + ++iIt; + ++oIt; + } return outImage; } void -calcOverlap( const TestImageType::Pointer& autoSeg, const TestImageType::Pointer& groundTruth, unsigned long long& tpCount, unsigned long long& fpCount, unsigned long long& fnCount ) +calcOverlap(const TestImageType::Pointer & autoSeg, + const TestImageType::Pointer & groundTruth, + unsigned long long & tpCount, + unsigned long long & fpCount, + unsigned long long & fnCount) { - const TestImageType::RegionType& lpr = groundTruth->GetLargestPossibleRegion(); + const TestImageType::RegionType & lpr = groundTruth->GetLargestPossibleRegion(); - itk::ImageRegionConstIterator< TestImageType > itAS( autoSeg, lpr ); - itk::ImageRegionConstIterator< TestImageType > itGT( groundTruth, lpr ); + itk::ImageRegionConstIterator itAS(autoSeg, lpr); + itk::ImageRegionConstIterator itGT(groundTruth, lpr); tpCount = 0; fpCount = 0; fnCount = 0; - while ( !itAS.IsAtEnd() ) + while (!itAS.IsAtEnd()) + { + if (itAS.Get() != 0 && itAS.Get() == itGT.Get()) { - if ( itAS.Get() != 0 && itAS.Get() == itGT.Get() ) - { - tpCount++; // true positive - } - else if ( itAS.Get() != 0 ) - { - fpCount++; // false positive - } - if ( itAS.Get() == 0 && itGT.Get() != 0 ) - { - fnCount++; // false negative - } - - ++itAS; - ++itGT; + tpCount++; // true positive + } + else if (itAS.Get() != 0) + { + fpCount++; // false positive } + if (itAS.Get() == 0 && itGT.Get() != 0) + { + fnCount++; // false negative + } + + ++itAS; + ++itGT; + } } int -main( int argc, char* argv[] ) +main(int argc, char * argv[]) { - if ( argc < 2 ) - { - std::cerr << "Usage: " << argv[0]; - std::cerr << " inputImage [outFilenameBase] [saveIntermediateImages]\n"; - return 1; - } + if (argc < 2) + { + std::cerr << "Usage: " << argv[0]; + std::cerr << " inputImage [outFilenameBase] [saveIntermediateImages]\n"; + return 1; + } bool saveImages = argc > 3; std::string outFilenameBase; std::fstream fout; - if ( argc > 2 ) - { - outFilenameBase = argv[2]; - fout.open( ( outFilenameBase + ".csv" ).c_str(), std::ios::out ); - } + if (argc > 2) + { + outFilenameBase = argv[2]; + fout.open((outFilenameBase + ".csv").c_str(), std::ios::out); + } RegisterRequiredFactories(); - using ReaderType = itk::ImageFileReader< TestImageType >; + using ReaderType = itk::ImageFileReader; ReaderType::Pointer reader = ReaderType::New(); - reader->SetFileName( argv[1] ); + reader->SetFileName(argv[1]); reader->Update(); TestImageType::Pointer inImage = reader->GetOutput(); inImage->DisconnectPipeline(); - using mciType = itk::MorphologicalContourInterpolator< TestImageType >; + using mciType = itk::MorphologicalContourInterpolator; mciType::Pointer mci = mciType::New(); - mci->SetUseBallStructuringElement( true ); // test cross? + mci->SetUseBallStructuringElement(true); // test cross? - using WriterType = itk::ImageFileWriter< TestImageType >; + using WriterType = itk::ImageFileWriter; WriterType::Pointer writer = WriterType::New(); - writer->SetUseCompression( true ); - - const TestImageType::RegionType& lpr = inImage->GetLargestPossibleRegion(); - TestImageType::IndexType maxInd; - for ( unsigned axis = 0; axis < testDim; axis++ ) - { - maxInd[axis] = itk::IndexValueType( lpr.GetSize( axis ) ); - } + writer->SetUseCompression(true); + + const TestImageType::RegionType & lpr = inImage->GetLargestPossibleRegion(); + TestImageType::IndexType maxInd; + for (unsigned axis = 0; axis < testDim; axis++) + { + maxInd[axis] = itk::IndexValueType(lpr.GetSize(axis)); + } fout << "Image, nth, axis, time, TP, FP, FN, TN\n"; // the big for loop which does the work - for ( int sparsity = 2; sparsity <= 8; sparsity++ ) + for (int sparsity = 2; sparsity <= 8; sparsity++) + { + unsigned long long tpCount, fpCount, fnCount; + for (int axis = -1; axis < int(testDim); axis++) { - unsigned long long tpCount, fpCount, fnCount; - for ( int axis = -1; axis < int( testDim ); axis++ ) + mci->SetAxis(axis); + TestImageType::IndexType axisSparsity = maxInd; + if (axis < 0) // all axes + { + for (unsigned a = 0; a < testDim; a++) { - mci->SetAxis( axis ); - TestImageType::IndexType axisSparsity = maxInd; - if ( axis < 0 ) // all axes - { - for ( unsigned a = 0; a < testDim; a++ ) - { - axisSparsity[a] = sparsity; - } - } - else // just one axis - { - axisSparsity[axis] = sparsity; - } - - TestImageType::Pointer sparseImage = createSparseCopy( inImage, axisSparsity ); - mci->SetInput( sparseImage ); - itk::TimeProbe timeProbe; - timeProbe.Start(); - mci->Update(); - timeProbe.Stop(); - TestImageType::Pointer result = mci->GetOutput(); - result->DisconnectPipeline(); - calcOverlap( result, inImage, tpCount, fpCount, fnCount ); - - fout << argv[1] << ", " << sparsity << ", " << axis << ", " << timeProbe.GetMean(); - fout << ", " << tpCount << ", " << fpCount << ", " << fnCount << ", "; - fout << ( lpr.GetNumberOfPixels() - tpCount - fpCount - fnCount ) << std::endl; - - if ( saveImages ) - { - std::cout << outFilenameBase + '_' + char( sparsity + '0' ) + char( axis + 'X' ) + "\nWriting sparse... "; - writer->SetInput( sparseImage ); - writer->SetFileName( outFilenameBase + "_in" + char( sparsity + '0' ) + char( axis + 'X' ) + ".mha" ); - writer->Update(); - std::cout << "interpolated... "; - writer->SetInput( result ); - writer->SetFileName( outFilenameBase + "_out" + char( sparsity + '0' ) + char( axis + 'X' ) + ".mha" ); - writer->Update(); - std::cout << "finished." << std::endl; - } + axisSparsity[a] = sparsity; } + } + else // just one axis + { + axisSparsity[axis] = sparsity; + } + + TestImageType::Pointer sparseImage = createSparseCopy(inImage, axisSparsity); + mci->SetInput(sparseImage); + itk::TimeProbe timeProbe; + timeProbe.Start(); + mci->Update(); + timeProbe.Stop(); + TestImageType::Pointer result = mci->GetOutput(); + result->DisconnectPipeline(); + calcOverlap(result, inImage, tpCount, fpCount, fnCount); + + fout << argv[1] << ", " << sparsity << ", " << axis << ", " << timeProbe.GetMean(); + fout << ", " << tpCount << ", " << fpCount << ", " << fnCount << ", "; + fout << (lpr.GetNumberOfPixels() - tpCount - fpCount - fnCount) << std::endl; + + if (saveImages) + { + std::cout << outFilenameBase + '_' + char(sparsity + '0') + char(axis + 'X') + "\nWriting sparse... "; + writer->SetInput(sparseImage); + writer->SetFileName(outFilenameBase + "_in" + char(sparsity + '0') + char(axis + 'X') + ".mha"); + writer->Update(); + std::cout << "interpolated... "; + writer->SetInput(result); + writer->SetFileName(outFilenameBase + "_out" + char(sparsity + '0') + char(axis + 'X') + ".mha"); + writer->Update(); + std::cout << "finished." << std::endl; + } } + } return 0; } diff --git a/test/itkMorphologicalContourInterpolationTest.cxx b/test/itkMorphologicalContourInterpolationTest.cxx index f21681f..960ee56 100644 --- a/test/itkMorphologicalContourInterpolationTest.cxx +++ b/test/itkMorphologicalContourInterpolationTest.cxx @@ -22,109 +22,111 @@ #include #include -template< typename ImageType > +template void -doTest( std::string inFilename, std::string outFilename, bool UseDistanceTransform, bool ball, int axis, int label ) +doTest(std::string inFilename, std::string outFilename, bool UseDistanceTransform, bool ball, int axis, int label) { - using ReaderType = itk::ImageFileReader< ImageType >; + using ReaderType = itk::ImageFileReader; typename ReaderType::Pointer reader = ReaderType::New(); - reader->SetFileName( inFilename ); + reader->SetFileName(inFilename); reader->Update(); typename ImageType::Pointer test = reader->GetOutput(); - using mciType = itk::MorphologicalContourInterpolator< ImageType >; + using mciType = itk::MorphologicalContourInterpolator; typename mciType::Pointer mci = mciType::New(); - mci->SetInput( test ); - mci->SetUseDistanceTransform( UseDistanceTransform ); - mci->SetUseBallStructuringElement( ball ); - mci->SetAxis( axis ); - mci->SetLabel( label ); + mci->SetInput(test); + mci->SetUseDistanceTransform(UseDistanceTransform); + mci->SetUseBallStructuringElement(ball); + mci->SetAxis(axis); + mci->SetLabel(label); - using WriterType = itk::ImageFileWriter< ImageType >; + using WriterType = itk::ImageFileWriter; typename WriterType::Pointer writer = WriterType::New(); - writer->SetFileName( outFilename ); - writer->SetInput( mci->GetOutput() ); - writer->SetUseCompression( true ); + writer->SetFileName(outFilename); + writer->SetInput(mci->GetOutput()); + writer->SetUseCompression(true); writer->Update(); } int -itkMorphologicalContourInterpolationTest( int argc, char* argv[] ) +itkMorphologicalContourInterpolationTest(int argc, char * argv[]) { - if ( argc < 3 ) + if (argc < 3) + { + std::cerr << "Usage: " << argv[0]; + std::cerr << " inputImage outputImage [algorithm] [axis] [label]\n"; + std::cerr << " algorithms:\n"; + std::cerr << " B = repeated dilations with ball structuring element"; + std::cerr << " C = repeated dilations with cross structuring element"; + std::cerr << " T = distance transform (not geodesic!)"; + std::cerr << " defaults: algo B, axis -1 (all axes), label 0 (all labels)"; + std::cerr << std::endl; + return EXIT_FAILURE; + } + const char * inputImageFileName = argv[1]; + const char * outputImageFileName = argv[2]; + bool dt = false; // DistanceTransform + bool ball = true; + int axis = -1, label = 0; + if (argc >= 4) + { + char algo = toupper(argv[3][0]); + if (algo == 'T') { - std::cerr << "Usage: " << argv[0]; - std::cerr << " inputImage outputImage [algorithm] [axis] [label]\n"; - std::cerr << " algorithms:\n"; - std::cerr << " B = repeated dilations with ball structuring element"; - std::cerr << " C = repeated dilations with cross structuring element"; - std::cerr << " T = distance transform (not geodesic!)"; - std::cerr << " defaults: algo B, axis -1 (all axes), label 0 (all labels)"; - std::cerr << std::endl; - return EXIT_FAILURE; + dt = true; } - const char* inputImageFileName = argv[1]; - const char* outputImageFileName = argv[2]; - bool dt = false; // DistanceTransform - bool ball = true; - int axis = -1, label = 0; - if ( argc >= 4 ) + else if (algo == 'C') { - char algo = toupper( argv[3][0] ); - if ( algo == 'T' ) - { - dt = true; - } - else if ( algo == 'C' ) - { - ball = false; - } - // else B - } - if ( argc >= 5 ) - { - axis = strtol( argv[4], nullptr, 10 ); - } - if ( argc >= 6 ) - { - label = strtol( argv[5], nullptr, 10 ); + ball = false; } + // else B + } + if (argc >= 5) + { + axis = strtol(argv[4], nullptr, 10); + } + if (argc >= 6) + { + label = strtol(argv[5], nullptr, 10); + } using ScalarPixelType = itk::CommonEnums::IOComponent; - itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO( inputImageFileName, itk::CommonEnums::IOFileMode::ReadMode ); - if ( !imageIO ) - { - std::cerr << "Could not CreateImageIO for: " << inputImageFileName << std::endl; - return EXIT_FAILURE; - } - imageIO->SetFileName( inputImageFileName ); + itk::ImageIOBase::Pointer imageIO = + itk::ImageIOFactory::CreateImageIO(inputImageFileName, itk::CommonEnums::IOFileMode::ReadMode); + if (!imageIO) + { + std::cerr << "Could not CreateImageIO for: " << inputImageFileName << std::endl; + return EXIT_FAILURE; + } + imageIO->SetFileName(inputImageFileName); imageIO->ReadImageInformation(); const ScalarPixelType pixelType = imageIO->GetComponentType(); const size_t numDimensions = imageIO->GetNumberOfDimensions(); try + { + // unused cases are not instantiated because they greatly increase compile time + if (numDimensions == 3 && + (pixelType == itk::CommonEnums::IOComponent::SHORT || pixelType == itk::CommonEnums::IOComponent::USHORT)) { - // unused cases are not instantiated because they greatly increase compile time - if ( numDimensions == 3 && ( pixelType == itk::CommonEnums::IOComponent::SHORT || pixelType == itk::CommonEnums::IOComponent::USHORT ) ) - { - doTest< itk::Image< short, 3 > >( inputImageFileName, outputImageFileName, dt, ball, axis, label ); - return EXIT_SUCCESS; - } - if ( numDimensions == 4 && pixelType == itk::CommonEnums::IOComponent::UCHAR ) - { - doTest< itk::Image< unsigned char, 4 > >( inputImageFileName, outputImageFileName, dt, ball, axis, label ); - return EXIT_SUCCESS; - } - - std::cerr << "Unsupported image type:\n Dimensions: " << numDimensions; - std::cerr << "\n Pixel type:" << imageIO->GetComponentTypeAsString( pixelType ) << std::endl; - return EXIT_FAILURE; + doTest>(inputImageFileName, outputImageFileName, dt, ball, axis, label); + return EXIT_SUCCESS; } - catch ( itk::ExceptionObject& error ) + if (numDimensions == 4 && pixelType == itk::CommonEnums::IOComponent::UCHAR) { - std::cerr << "Error: " << error << std::endl; - return EXIT_FAILURE; + doTest>(inputImageFileName, outputImageFileName, dt, ball, axis, label); + return EXIT_SUCCESS; } + + std::cerr << "Unsupported image type:\n Dimensions: " << numDimensions; + std::cerr << "\n Pixel type:" << imageIO->GetComponentTypeAsString(pixelType) << std::endl; + return EXIT_FAILURE; + } + catch (itk::ExceptionObject & error) + { + std::cerr << "Error: " << error << std::endl; + return EXIT_FAILURE; + } return EXIT_SUCCESS; } diff --git a/test/itkMorphologicalContourInterpolationTestWithRLEImage.cxx b/test/itkMorphologicalContourInterpolationTestWithRLEImage.cxx index 7277e7e..aecea4e 100644 --- a/test/itkMorphologicalContourInterpolationTestWithRLEImage.cxx +++ b/test/itkMorphologicalContourInterpolationTestWithRLEImage.cxx @@ -23,20 +23,20 @@ #include #include -template< typename ImageType > +template void -doTest( std::string inFilename, std::string outFilename, bool UseDistanceTransform, bool ball, int axis, int label ) +doTest(std::string inFilename, std::string outFilename, bool UseDistanceTransform, bool ball, int axis, int label) { - using ReaderType = itk::ImageFileReader< ImageType >; + using ReaderType = itk::ImageFileReader; typename ReaderType::Pointer reader = ReaderType::New(); - reader->SetFileName( inFilename ); + reader->SetFileName(inFilename); reader->Update(); - using myRLEImage = itk::RLEImage< typename ImageType::PixelType, ImageType::ImageDimension >; - using inConverterType = itk::RegionOfInterestImageFilter< ImageType, myRLEImage >; + using myRLEImage = itk::RLEImage; + using inConverterType = itk::RegionOfInterestImageFilter; typename inConverterType::Pointer inConv = inConverterType::New(); - inConv->SetInput( reader->GetOutput() ); - inConv->SetRegionOfInterest( reader->GetOutput()->GetLargestPossibleRegion() ); + inConv->SetInput(reader->GetOutput()); + inConv->SetRegionOfInterest(reader->GetOutput()->GetLargestPossibleRegion()); inConv->Update(); typename myRLEImage::Pointer test = inConv->GetOutput(); @@ -49,104 +49,106 @@ doTest( std::string inFilename, std::string outFilename, bool UseDistanceTransfo // reg.SetSize(i, (reg.GetSize(i) + 1) / 2); // } - using mciType = itk::MorphologicalContourInterpolator< myRLEImage >; + using mciType = itk::MorphologicalContourInterpolator; typename mciType::Pointer mci = mciType::New(); - mci->SetInput( test ); - mci->SetUseDistanceTransform( UseDistanceTransform ); - mci->SetUseBallStructuringElement( ball ); - mci->SetAxis( axis ); - mci->SetLabel( label ); + mci->SetInput(test); + mci->SetUseDistanceTransform(UseDistanceTransform); + mci->SetUseBallStructuringElement(ball); + mci->SetAxis(axis); + mci->SetLabel(label); - using outConverterType = itk::RegionOfInterestImageFilter< myRLEImage, ImageType >; + using outConverterType = itk::RegionOfInterestImageFilter; typename outConverterType::Pointer outConv = outConverterType::New(); - outConv->SetInput( mci->GetOutput() ); - outConv->SetRegionOfInterest( reg ); + outConv->SetInput(mci->GetOutput()); + outConv->SetRegionOfInterest(reg); outConv->Update(); - using WriterType = itk::ImageFileWriter< ImageType >; + using WriterType = itk::ImageFileWriter; typename WriterType::Pointer writer = WriterType::New(); - writer->SetFileName( outFilename ); - writer->SetInput( outConv->GetOutput() ); - writer->SetUseCompression( true ); + writer->SetFileName(outFilename); + writer->SetInput(outConv->GetOutput()); + writer->SetUseCompression(true); writer->Update(); } int -itkMorphologicalContourInterpolationTestWithRLEImage( int argc, char* argv[] ) +itkMorphologicalContourInterpolationTestWithRLEImage(int argc, char * argv[]) { - if ( argc < 3 ) + if (argc < 3) + { + std::cerr << "Usage: " << argv[0]; + std::cerr << " inputImage outputImage [algorithm] [axis] [label]\n"; + std::cerr << " algorithms:\n"; + std::cerr << " B = repeated dilations with ball structuring element"; + std::cerr << " C = repeated dilations with cross structuring element"; + std::cerr << " T = distance transform (not geodesic!)"; + std::cerr << " defaults: algo B, axis -1 (all axes), label 0 (all labels)"; + std::cerr << std::endl; + return EXIT_FAILURE; + } + const char * inputImageFileName = argv[1]; + const char * outputImageFileName = argv[2]; + bool dt = false; // DistanceTransform + bool ball = true; + int axis = -1, label = 0; + if (argc >= 4) + { + char algo = toupper(argv[3][0]); + if (algo == 'T') { - std::cerr << "Usage: " << argv[0]; - std::cerr << " inputImage outputImage [algorithm] [axis] [label]\n"; - std::cerr << " algorithms:\n"; - std::cerr << " B = repeated dilations with ball structuring element"; - std::cerr << " C = repeated dilations with cross structuring element"; - std::cerr << " T = distance transform (not geodesic!)"; - std::cerr << " defaults: algo B, axis -1 (all axes), label 0 (all labels)"; - std::cerr << std::endl; - return EXIT_FAILURE; + dt = true; } - const char* inputImageFileName = argv[1]; - const char* outputImageFileName = argv[2]; - bool dt = false; // DistanceTransform - bool ball = true; - int axis = -1, label = 0; - if ( argc >= 4 ) + else if (algo == 'C') { - char algo = toupper( argv[3][0] ); - if ( algo == 'T' ) - { - dt = true; - } - else if ( algo == 'C' ) - { - ball = false; - } - // else B - } - if ( argc >= 5 ) - { - axis = strtol( argv[4], nullptr, 10 ); - } - if ( argc >= 6 ) - { - label = strtol( argv[5], nullptr, 10 ); + ball = false; } + // else B + } + if (argc >= 5) + { + axis = strtol(argv[4], nullptr, 10); + } + if (argc >= 6) + { + label = strtol(argv[5], nullptr, 10); + } using ScalarPixelType = itk::CommonEnums::IOComponent; - itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO( inputImageFileName, itk::CommonEnums::IOFileMode::ReadMode ); - if ( !imageIO ) - { - std::cerr << "Could not CreateImageIO for: " << inputImageFileName << std::endl; - return EXIT_FAILURE; - } - imageIO->SetFileName( inputImageFileName ); + itk::ImageIOBase::Pointer imageIO = + itk::ImageIOFactory::CreateImageIO(inputImageFileName, itk::CommonEnums::IOFileMode::ReadMode); + if (!imageIO) + { + std::cerr << "Could not CreateImageIO for: " << inputImageFileName << std::endl; + return EXIT_FAILURE; + } + imageIO->SetFileName(inputImageFileName); imageIO->ReadImageInformation(); const ScalarPixelType pixelType = imageIO->GetComponentType(); const size_t numDimensions = imageIO->GetNumberOfDimensions(); try + { + // unused cases are not instantiated because they greatly increase compile time + if (numDimensions == 3 && + (pixelType == itk::CommonEnums::IOComponent::SHORT || pixelType == itk::CommonEnums::IOComponent::USHORT)) { - // unused cases are not instantiated because they greatly increase compile time - if ( numDimensions == 3 && ( pixelType == itk::CommonEnums::IOComponent::SHORT || pixelType == itk::CommonEnums::IOComponent::USHORT ) ) - { - doTest< itk::Image< short, 3 > >( inputImageFileName, outputImageFileName, dt, ball, axis, label ); - return EXIT_SUCCESS; - } - if ( numDimensions == 4 && pixelType == itk::CommonEnums::IOComponent::UCHAR ) - { - doTest< itk::Image< unsigned char, 4 > >( inputImageFileName, outputImageFileName, dt, ball, axis, label ); - return EXIT_SUCCESS; - } - - std::cerr << "Unsupported image type:\n Dimensions: " << numDimensions; - std::cerr << "\n Pixel type:" << imageIO->GetComponentTypeAsString( pixelType ) << std::endl; - return EXIT_FAILURE; + doTest>(inputImageFileName, outputImageFileName, dt, ball, axis, label); + return EXIT_SUCCESS; } - catch ( itk::ExceptionObject& error ) + if (numDimensions == 4 && pixelType == itk::CommonEnums::IOComponent::UCHAR) { - std::cerr << "Error: " << error << std::endl; - return EXIT_FAILURE; + doTest>(inputImageFileName, outputImageFileName, dt, ball, axis, label); + return EXIT_SUCCESS; } + + std::cerr << "Unsupported image type:\n Dimensions: " << numDimensions; + std::cerr << "\n Pixel type:" << imageIO->GetComponentTypeAsString(pixelType) << std::endl; + return EXIT_FAILURE; + } + catch (itk::ExceptionObject & error) + { + std::cerr << "Error: " << error << std::endl; + return EXIT_FAILURE; + } return EXIT_SUCCESS; } diff --git a/test/manualTest.cxx b/test/manualTest.cxx index 713a1b0..5b6d80f 100644 --- a/test/manualTest.cxx +++ b/test/manualTest.cxx @@ -17,17 +17,17 @@ *=========================================================================*/ #ifdef _MSC_VER -#pragma warning( disable : 4996 ) /* deprecation */ +# pragma warning(disable : 4996) /* deprecation */ #endif #include "itkTestDriverIncludeRequiredIOFactories.h" extern int -itkMorphologicalContourInterpolationTest( int argc, char* argv[] ); +itkMorphologicalContourInterpolationTest(int argc, char * argv[]); int -main( int argc, char* argv[] ) +main(int argc, char * argv[]) { RegisterRequiredFactories(); - itkMorphologicalContourInterpolationTest( argc, argv ); + itkMorphologicalContourInterpolationTest(argc, argv); } diff --git a/wasm/morphological-contour-interpolation.cxx b/wasm/morphological-contour-interpolation.cxx index 09eb68d..d09b2ed 100644 --- a/wasm/morphological-contour-interpolation.cxx +++ b/wasm/morphological-contour-interpolation.cxx @@ -23,103 +23,119 @@ #include "itkMorphologicalContourInterpolator.h" -template< typename TImage > +template int -MorphologicalContourInterpolation( itk::wasm::Pipeline& pipeline, const TImage* inputImage ) +MorphologicalContourInterpolation(itk::wasm::Pipeline & pipeline, const TImage * inputImage) { using ImageType = TImage; - pipeline.get_option( "input-image" )->required()->type_name( "INPUT_IMAGE" ); + pipeline.get_option("input-image")->required()->type_name("INPUT_IMAGE"); - using OutputImageType = itk::wasm::OutputImage< ImageType >; + using OutputImageType = itk::wasm::OutputImage; OutputImageType outputImage; - pipeline.add_option( "output-image", outputImage, "The output image" )->required()->type_name( "OUTPUT_IMAGE" ); + pipeline.add_option("output-image", outputImage, "The output image")->required()->type_name("OUTPUT_IMAGE"); typename ImageType::IndexValueType label = 0; - pipeline.add_option( "-l,--label", label, "The label to interpolate. Interpolates all labels if set to 0 (default)." ); + pipeline.add_option("-l,--label", label, "The label to interpolate. Interpolates all labels if set to 0 (default)."); int axis = -1; - pipeline.add_option( "-a,--axis", axis, "Interpolate only along this axis. Interpolates along all axes if set to -1 (default)." ); + pipeline.add_option( + "-a,--axis", axis, "Interpolate only along this axis. Interpolates along all axes if set to -1 (default)."); bool noHeuristicAlignment{ false }; - pipeline.add_flag( "--no-heuristic-alignment", noHeuristicAlignment, "Heuristic alignment of regions for interpolation is faster than optimal alignment." ); + pipeline.add_flag("--no-heuristic-alignment", + noHeuristicAlignment, + "Heuristic alignment of regions for interpolation is faster than optimal alignment."); bool noUseDistanceTransform{ false }; - pipeline.add_flag( "--no-use-distance-transform", - noUseDistanceTransform, - "Using distance transform instead of repeated dilations to calculate the median contour is slightly faster, but produces lower quality interpolations." ); + pipeline.add_flag("--no-use-distance-transform", + noUseDistanceTransform, + "Using distance transform instead of repeated dilations to calculate the median contour is " + "slightly faster, but produces lower quality interpolations."); bool useCustomSlicePositions{ false }; - pipeline.add_flag( "--use-custom-slice-positions", useCustomSlicePositions, "Use custom slice positions (not slice auto-detection)." ); + pipeline.add_flag( + "--use-custom-slice-positions", useCustomSlicePositions, "Use custom slice positions (not slice auto-detection)."); bool noUseExtrapolation{ false }; - pipeline.add_flag( "--no-use-extrapolation", - noUseExtrapolation, - "Perform extrapolation for branch extremities. Branch extremities are defined as regions having no overlap with any region in the next slice. Extrapolation helps generate smooth " - "surface closings." ); + pipeline.add_flag("--no-use-extrapolation", + noUseExtrapolation, + "Perform extrapolation for branch extremities. Branch extremities are defined as regions having no " + "overlap with any region in the next slice. Extrapolation helps generate smooth " + "surface closings."); bool useBallStructuringElement; - pipeline.add_flag( "--use-ball-structuring-element", useBallStructuringElement, "Use ball instead of default cross structuring element for repeated dilations." ); + pipeline.add_flag("--use-ball-structuring-element", + useBallStructuringElement, + "Use ball instead of default cross structuring element for repeated dilations."); int labeledSliceIndicesAxis = -1; - pipeline.add_option( "--labeled-slice-indices-axis", labeledSliceIndicesAxis, "Axis along which the labeled slice indices are defined. Default is -1 (that is, auto-detection)." ); + pipeline.add_option( + "--labeled-slice-indices-axis", + labeledSliceIndicesAxis, + "Axis along which the labeled slice indices are defined. Default is -1 (that is, auto-detection)."); int labeledSliceIndicesLabel = 1; - pipeline.add_option( "--labeled-slice-indices-label", labeledSliceIndicesLabel, "Label of the slice indices. Default is 1." ); + pipeline.add_option( + "--labeled-slice-indices-label", labeledSliceIndicesLabel, "Label of the slice indices. Default is 1."); - std::vector< typename ImageType::IndexValueType > labeledSliceIndices; - pipeline.add_option( "--labeled-slice-indices", labeledSliceIndices, "List of labeled slice indices. Default is empty." ); + std::vector labeledSliceIndices; + pipeline.add_option( + "--labeled-slice-indices", labeledSliceIndices, "List of labeled slice indices. Default is empty."); - ITK_WASM_PARSE( pipeline ); + ITK_WASM_PARSE(pipeline); - using InterpolatorType = itk::MorphologicalContourInterpolator< ImageType >; + using InterpolatorType = itk::MorphologicalContourInterpolator; auto interpolator = InterpolatorType::New(); - interpolator->SetInput( inputImage ); - - interpolator->SetLabel( label ); - interpolator->SetAxis( axis ); - interpolator->SetHeuristicAlignment( !noHeuristicAlignment ); - interpolator->SetUseDistanceTransform( !noUseDistanceTransform ); - interpolator->SetUseCustomSlicePositions( useCustomSlicePositions ); - interpolator->SetUseExtrapolation( !noUseExtrapolation ); - interpolator->SetUseBallStructuringElement( !useBallStructuringElement ); - if ( labeledSliceIndicesAxis != -1 && !labeledSliceIndices.empty() ) - { - interpolator->SetLabeledSliceIndices( labeledSliceIndicesAxis, labeledSliceIndicesLabel, labeledSliceIndices ); - } - - ITK_WASM_CATCH_EXCEPTION( pipeline, interpolator->UpdateLargestPossibleRegion() ); + interpolator->SetInput(inputImage); + + interpolator->SetLabel(label); + interpolator->SetAxis(axis); + interpolator->SetHeuristicAlignment(!noHeuristicAlignment); + interpolator->SetUseDistanceTransform(!noUseDistanceTransform); + interpolator->SetUseCustomSlicePositions(useCustomSlicePositions); + interpolator->SetUseExtrapolation(!noUseExtrapolation); + interpolator->SetUseBallStructuringElement(!useBallStructuringElement); + if (labeledSliceIndicesAxis != -1 && !labeledSliceIndices.empty()) + { + interpolator->SetLabeledSliceIndices(labeledSliceIndicesAxis, labeledSliceIndicesLabel, labeledSliceIndices); + } + + ITK_WASM_CATCH_EXCEPTION(pipeline, interpolator->UpdateLargestPossibleRegion()); typename ImageType::ConstPointer output = interpolator->GetOutput(); - outputImage.Set( output ); + outputImage.Set(output); return EXIT_SUCCESS; } -template< typename TImage > +template class PipelineFunctor { public: int - operator()( itk::wasm::Pipeline& pipeline ) + operator()(itk::wasm::Pipeline & pipeline) { using ImageType = TImage; - using InputImageType = itk::wasm::InputImage< ImageType >; + using InputImageType = itk::wasm::InputImage; InputImageType inputImage; - pipeline.add_option( "input-image", inputImage, "The input image" ); + pipeline.add_option("input-image", inputImage, "The input image"); - ITK_WASM_PRE_PARSE( pipeline ); + ITK_WASM_PRE_PARSE(pipeline); typename ImageType::ConstPointer image = inputImage.Get(); - return MorphologicalContourInterpolation< ImageType >( pipeline, image ); + return MorphologicalContourInterpolation(pipeline, image); } }; int -main( int argc, char* argv[] ) +main(int argc, char * argv[]) { - itk::wasm::Pipeline pipeline( "morphological-contour-interpolation", "Interpolates contours between slices.", argc, argv ); + itk::wasm::Pipeline pipeline( + "morphological-contour-interpolation", "Interpolates contours between slices.", argc, argv); - return itk::wasm::SupportInputImageTypes< PipelineFunctor, uint8_t, int8_t, uint16_t, int16_t, uint32_t, int32_t, uint64_t, int64_t >::Dimensions< 3U, 4U >( "input-image", pipeline ); + return itk::wasm:: + SupportInputImageTypes:: + Dimensions<3U, 4U>("input-image", pipeline); }