diff --git a/.circleci/config.yml b/.circleci/config.yml new file mode 100644 index 0000000..c078b1b --- /dev/null +++ b/.circleci/config.yml @@ -0,0 +1,63 @@ +version: 2 +jobs: + build-and-test: + working_directory: /ITKPrincipalComponentsAnalysis-build + docker: + - image: insighttoolkit/module-ci:latest + steps: + - checkout: + path: /ITKPrincipalComponentsAnalysis + - run: + name: Fetch CTest driver script + command: | + curl -L https://raw.githubusercontent.com/InsightSoftwareConsortium/ITK/dashboard/itk_common.cmake -O + - run: + name: Configure CTest script + command: | + SHASNIP=$(echo $CIRCLE_SHA1 | cut -c1-7) + cat > dashboard.cmake << EOF + set(CTEST_SITE "CircleCI") + set(CTEST_BUILD_NAME "External-ITKPrincipalComponentsAnalysis-${CIRCLE_BRANCH}-${CIRCLE_BUILD_NUM}-${SHASNIP}") + set(CTEST_BUILD_CONFIGURATION "MinSizeRel") + set(CTEST_CMAKE_GENERATOR "Unix Makefiles") + set(CTEST_BUILD_FLAGS: "-j5") + set(CTEST_SOURCE_DIRECTORY /ITKPrincipalComponentsAnalysis) + set(CTEST_BINARY_DIRECTORY /ITKPrincipalComponentsAnalysis-build) + set(dashboard_model Experimental) + set(dashboard_no_clean 1) + set(dashboard_cache " + ITK_DIR:PATH=/ITK-build + BUILD_TESTING:BOOL=ON + ") + include(\${CTEST_SCRIPT_DIRECTORY}/itk_common.cmake) + EOF + - run: + name: Build and Test + no_output_timeout: 1.0h + command: | + ctest -j 2 -VV -S dashboard.cmake + package: + working_directory: ~/ITKPrincipalComponentsAnalysis + machine: true + steps: + - checkout + - run: + name: Fetch build script + command: | + curl -L https://rawgit.com/InsightSoftwareConsortium/ITKPythonPackage/master/scripts/dockcross-manylinux-download-cache-and-build-module-wheels.sh -O + chmod u+x dockcross-manylinux-download-cache-and-build-module-wheels.sh + - run: + name: Build Python packages + no_output_timeout: 1.0h + command: | + ./dockcross-manylinux-download-cache-and-build-module-wheels.sh + - store_artifacts: + path: dist + destination: dist + +workflows: + version: 2 + build-test-package: + jobs: + - build-and-test + - package diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..c0a5d09 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,15 @@ +sudo: true +language: cpp +os: +- osx +compiler: +- gcc +cache: + directories: + - "$HOME/Library/Caches/Homebrew" +script: +- curl -L https://raw.githubusercontent.com/InsightSoftwareConsortium/ITKPythonPackage/master/scripts/macpython-download-cache-and-build-module-wheels.sh -O +- chmod u+x macpython-download-cache-and-build-module-wheels.sh +- ./macpython-download-cache-and-build-module-wheels.sh +- tar -zcvf dist.tar.gz dist/ +- curl --upload-file dist.tar.gz https://transfer.sh/dist.tar.gz diff --git a/CMakeLists.txt b/CMakeLists.txt index 2231f55..0edc2e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,12 +1,12 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.9) +CMAKE_MINIMUM_REQUIRED(VERSION 3.10.2) IF(COMMAND CMAKE_POLICY) CMAKE_POLICY(SET CMP0003 NEW) ENDIF(COMMAND CMAKE_POLICY) -SET(kit PrincipalComponentsAnalysis) -PROJECT(${kit}) +SET(PCA PrincipalComponentsAnalysis) +PROJECT(${PCA}) -IF (NOT ITK_SOURCE_DIR) +IF(NOT ITK_SOURCE_DIR) FIND_PACKAGE(ITK REQUIRED) LIST(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR}) INCLUDE(ITKModuleExternal) @@ -14,13 +14,14 @@ ELSE() itk_module_impl() ENDIF() -OPTION(${kit}_BUILD_EXAMPLES "Build the examples" OFF) -IF(${kit}_BUILD_EXAMPLES) - ADD_SUBDIRECTORY( Examples ) +CMAKE_DEPENDENT_OPTION(Module_${PCA}_BUILD_EXAMPLES "Build the examples" OFF "BUILD_EXAMPLES" OFF) +if(Module_${PCA}_BUILD_EXAMPLES) + ADD_SUBDIRECTORY( examples ) ENDIF() -OPTION(${kit}_GENERATE_REPORT "Generate the PDF report from latex files, source code examples and result screenshots" OFF) -IF(${kit}_GENERATE_REPORT) - ADD_SUBDIRECTORY( Documents ) +CMAKE_DEPENDENT_OPTION(Module_${PCA}_BUILD_DOCUMENTATION "Generate documentation from LaTeX files, source code examples and result screenshots" OFF + "BUILD_DOCUMENTATION" OFF) +IF(Module_${PCA}_BUILD_DOCUMENTATION) + ADD_SUBDIRECTORY( doc ) ENDIF() diff --git a/CTestConfig.cmake b/CTestConfig.cmake new file mode 100644 index 0000000..62f68ac --- /dev/null +++ b/CTestConfig.cmake @@ -0,0 +1,7 @@ +set(CTEST_PROJECT_NAME "ITK") +set(CTEST_NIGHTLY_START_TIME "1:00:00 UTC") + +set(CTEST_DROP_METHOD "http") +set(CTEST_DROP_SITE "open.cdash.org") +set(CTEST_DROP_LOCATION "/submit.php?project=Insight") +set(CTEST_DROP_SITE_CDASH TRUE) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d645695 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/README.rst b/README.rst index 9f8e0ca..b7b48a7 100644 --- a/README.rst +++ b/README.rst @@ -1,19 +1,46 @@ PrincipalComponentsAnalysis -============================ +=========================== -.. ADD circle ci badge once this module is merge to ITK's repository +.. |CircleCI| image:: https://circleci.com/gh/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis.svg?style=shield + :target: https://circleci.com/gh/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis + +.. |TravisCI| image:: https://travis-ci.org/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis.svg?branch=master + :target: https://travis-ci.org/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis + +.. |AppVeyor| image:: https://img.shields.io/appveyor/ci/itkrobot/itkprincipalcomponentsanalysis.svg + :target: https://ci.appveyor.com/project/itkrobot/itkprincipalcomponentsanalysis + +=========== =========== =========== + Linux macOS Windows +=========== =========== =========== +|CircleCI| |TravisCI| |AppVeyor| +=========== =========== =========== + +Overview +-------- An `ITK `_-based implementation of principal components analysis. -A more detailed description can be found in the Insight Journal article:: - Bowers M., Younes L. "Principal Components Analysis of Scalar, Vector, and Mesh Vertex Data." - http://hdl.handle.net/10380/3386 - http://www.insight-journal.org/browse/publication/878 - August, 2013. +For more information, see the `Insight Journal article `_:: + + Bowers M., Younes L. + Principal Components Analysis of Scalar, Vector, and Mesh Vertex Data + The Insight Journal. August. 2013. + http://hdl.handle.net/10380/3386 + http://www.insight-journal.org/browse/publication/878 + +Installation +------------ Since ITK 4.10.0, this module is available in the ITK source tree as a Remote -module. To enable it, set:: +module. To enable it, set:: Module_PrincipalComponentsAnalysis:BOOL=ON in ITK's CMake build configuration. + +License +------- + +This software is distributed under the Apache 2.0 license. Please see +the *LICENSE* file for details. diff --git a/appveyor.yml b/appveyor.yml new file mode 100644 index 0000000..e18ea36 --- /dev/null +++ b/appveyor.yml @@ -0,0 +1,21 @@ +branches: + only: + - master + +version: "0.0.1.{build}" + +install: + + - curl -L https://raw.githubusercontent.com/InsightSoftwareConsortium/ITKPythonPackage/master/scripts/windows-download-cache-and-build-module-wheels.ps1 -O + - ps: .\windows-download-cache-and-build-module-wheels.ps1 + +build: off + +test: off + +artifacts: + + # pushing entire folder as a zip archive + - path: dist\* + +deploy: off \ No newline at end of file diff --git a/Documents/CMakeLists.txt b/doc/CMakeLists.txt similarity index 100% rename from Documents/CMakeLists.txt rename to doc/CMakeLists.txt diff --git a/Documents/Report001/CMakeLists.txt b/doc/Report001/CMakeLists.txt similarity index 100% rename from Documents/Report001/CMakeLists.txt rename to doc/Report001/CMakeLists.txt diff --git a/Documents/Report001/CaudateFirstPrincipalComponentOnMomentum.jpg b/doc/Report001/CaudateFirstPrincipalComponentOnMomentum.jpg similarity index 100% rename from Documents/Report001/CaudateFirstPrincipalComponentOnMomentum.jpg rename to doc/Report001/CaudateFirstPrincipalComponentOnMomentum.jpg diff --git a/Documents/Report001/InsightArticle.cls b/doc/Report001/InsightArticle.cls similarity index 100% rename from Documents/Report001/InsightArticle.cls rename to doc/Report001/InsightArticle.cls diff --git a/Documents/Report001/InsightJournal.bib b/doc/Report001/InsightJournal.bib similarity index 100% rename from Documents/Report001/InsightJournal.bib rename to doc/Report001/InsightJournal.bib diff --git a/Documents/Report001/InsightJournal.ist b/doc/Report001/InsightJournal.ist similarity index 100% rename from Documents/Report001/InsightJournal.ist rename to doc/Report001/InsightJournal.ist diff --git a/Documents/Report001/InsightJournal.sty b/doc/Report001/InsightJournal.sty similarity index 100% rename from Documents/Report001/InsightJournal.sty rename to doc/Report001/InsightJournal.sty diff --git a/Documents/Report001/PrincipalComponentsAnalysis.tex b/doc/Report001/PrincipalComponentsAnalysis.tex similarity index 100% rename from Documents/Report001/PrincipalComponentsAnalysis.tex rename to doc/Report001/PrincipalComponentsAnalysis.tex diff --git a/Documents/Report001/algorithm.sty b/doc/Report001/algorithm.sty similarity index 100% rename from Documents/Report001/algorithm.sty rename to doc/Report001/algorithm.sty diff --git a/Documents/Report001/algorithmic.sty b/doc/Report001/algorithmic.sty similarity index 100% rename from Documents/Report001/algorithmic.sty rename to doc/Report001/algorithmic.sty diff --git a/Documents/Report001/amssymb.sty b/doc/Report001/amssymb.sty similarity index 100% rename from Documents/Report001/amssymb.sty rename to doc/Report001/amssymb.sty diff --git a/Documents/Report001/fancyhdr.sty b/doc/Report001/fancyhdr.sty similarity index 100% rename from Documents/Report001/fancyhdr.sty rename to doc/Report001/fancyhdr.sty diff --git a/Documents/Report001/floatflt.sty b/doc/Report001/floatflt.sty similarity index 100% rename from Documents/Report001/floatflt.sty rename to doc/Report001/floatflt.sty diff --git a/Documents/Report001/fncychap.sty b/doc/Report001/fncychap.sty similarity index 100% rename from Documents/Report001/fncychap.sty rename to doc/Report001/fncychap.sty diff --git a/Documents/Report001/picins.sty b/doc/Report001/picins.sty similarity index 100% rename from Documents/Report001/picins.sty rename to doc/Report001/picins.sty diff --git a/Documents/Report001/times.sty b/doc/Report001/times.sty similarity index 100% rename from Documents/Report001/times.sty rename to doc/Report001/times.sty diff --git a/Examples/CMakeLists.txt b/examples/CMakeLists.txt similarity index 100% rename from Examples/CMakeLists.txt rename to examples/CMakeLists.txt diff --git a/Examples/VectorKernelPCA.cxx b/examples/VectorKernelPCA.cxx similarity index 75% rename from Examples/VectorKernelPCA.cxx rename to examples/VectorKernelPCA.cxx index 991c16e..cfb6ab8 100644 --- a/Examples/VectorKernelPCA.cxx +++ b/examples/VectorKernelPCA.cxx @@ -25,10 +25,10 @@ int showUsage(const char* programName) { - vcl_cerr << "USAGE: " << programName << " " << vcl_endl; - vcl_cerr << " " << vcl_endl; - vcl_cerr << "\t\tpcaCount : number of principal components to calculate" << vcl_endl; - vcl_cerr << "\t\tkernelSigma : KernelSigma" << vcl_endl; + std::cerr << "USAGE: " << programName << " " << std::endl; + std::cerr << " " << std::endl; + std::cerr << "\t\tpcaCount : number of principal components to calculate" << std::endl; + std::cerr << "\t\tkernelSigma : KernelSigma" << std::endl; return EXIT_FAILURE; } @@ -40,41 +40,39 @@ int main( int argc, char *argv[] ) if(argc < MIN_ARG_COUNT) return (showUsage(argv[0])); - int pcaCount = atoi(argv[1]); - double kernelSigma = atof(argv[2]); + int pcaCount = std::stoi(argv[1]); + double kernelSigma = std::stod(argv[2]); - typedef double PointDataType; - typedef itk::Array PointDataVectorType; - typedef PointDataVectorType PixelType; - typedef itk::Vector OutPointDataVectorType; - typedef OutPointDataVectorType OutPixelType; + using PointDataType = double; + using PointDataVectorType = itk::Array; + using PixelType = PointDataVectorType; + using OutPointDataVectorType = itk::Vector; + using OutPixelType = OutPointDataVectorType; - typedef double DDataType; - typedef double CoordRep; - typedef double InterpRep; + using CoordRep = double; const unsigned int Dimension = 3; -// typedef float PCAResultsType; - typedef double PCAResultsType; +// using PCAResultsType = float; + using PCAResultsType = double; // Declare the type of the input mesh - typedef itk::Mesh InMeshType; - typedef itk::Mesh OutMeshType; + using InMeshType = itk::Mesh; + using OutMeshType = itk::Mesh; // Declare the type of the kernel function class - typedef itk::GaussianDistanceKernel KernelType; + using KernelType = itk::GaussianDistanceKernel; // Declare the type of the PCA calculator - typedef itk::VectorFieldPCA< PointDataType, PCAResultsType, - PixelType, CoordRep, KernelType, InMeshType > PCACalculatorType; + using PCACalculatorType = itk::VectorFieldPCA< PointDataType, PCAResultsType, + PixelType, CoordRep, KernelType, InMeshType >; // Here we recover the file names from the command line arguments const char* inMeshFile = argv[3]; const char* outFileNameBase = argv[4]; // We can now instantiate the types of the reader/writer. - typedef itk::MeshFileReader< InMeshType > ReaderType; - typedef itk::MeshFileWriter< OutMeshType > WriterType; + using ReaderType = itk::MeshFileReader< InMeshType >; + using WriterType = itk::MeshFileWriter< OutMeshType >; // create readers/writers ReaderType::Pointer meshReader = ReaderType::New(); @@ -90,17 +88,17 @@ int main( int argc, char *argv[] ) } catch( itk::ExceptionObject & excp ) { - vcl_cerr << "Error reading mesh file " << inMeshFile << vcl_endl; - vcl_cerr << excp << vcl_endl; + std::cerr << "Error reading mesh file " << inMeshFile << std::endl; + std::cerr << excp << std::endl; } // get the objects InMeshType::Pointer mesh = meshReader->GetOutput(); - vcl_cout << "Vertex Count: " << - mesh->GetNumberOfPoints() << vcl_endl; - vcl_cout << "Cell Count: " << - mesh->GetNumberOfCells() << vcl_endl; + std::cout << "Vertex Count: " << + mesh->GetNumberOfPoints() << std::endl; + std::cout << "Cell Count: " << + mesh->GetNumberOfCells() << std::endl; const char* vectorFieldName; PCACalculatorType::VectorFieldType vectorField; @@ -129,8 +127,8 @@ int main( int argc, char *argv[] ) } catch( itk::ExceptionObject & excp ) { - vcl_cerr << "Error reading mesh field file " << vectorFieldName << vcl_endl; - vcl_cerr << excp << vcl_endl; + std::cerr << "Error reading mesh field file " << vectorFieldName << std::endl; + std::cerr << excp << std::endl; } // get the objects @@ -148,9 +146,9 @@ int main( int argc, char *argv[] ) } if (vectorFieldCount != meshWithField->GetNumberOfPoints()) { - vcl_cerr << "Vector field count (" << vectorFieldCount << + std::cerr << "Vector field count (" << vectorFieldCount << ") doesn't match mesh vertext count (" << - meshWithField->GetNumberOfPoints() << ")." << vcl_endl; + meshWithField->GetNumberOfPoints() << ")." << std::endl; exit (EXIT_FAILURE); } vectorField.set_size(vectorFieldCount, vectorFieldDimension); @@ -160,9 +158,9 @@ int main( int argc, char *argv[] ) if (vectorFieldDimension != vectorField.cols() || vectorFieldCount != vectorField.rows()) { - vcl_cerr << "Unexpected dimensions in vector field file " << vectorFieldName << vcl_endl; - vcl_cerr << "\tExpected " << vectorFieldCount << " x " << vectorFieldDimension; - vcl_cerr << "\t, got " << vectorField.rows() << " x " << vectorField.cols() << vcl_endl; + std::cerr << "Unexpected dimensions in vector field file " << vectorFieldName << std::endl; + std::cerr << "\tExpected " << vectorFieldCount << " x " << vectorFieldDimension; + std::cerr << "\t, got " << vectorField.rows() << " x " << vectorField.cols() << std::endl; exit(1); } } @@ -198,7 +196,7 @@ int main( int argc, char *argv[] ) } catch( itk::ExceptionObject & excp ) { - vcl_cerr << excp << vcl_endl; + std::cerr << excp << std::endl; return EXIT_FAILURE; } @@ -222,8 +220,8 @@ int main( int argc, char *argv[] ) InMeshType::CellsContainerConstIterator itCellsEnd = cells->End(); outMesh->GetCells()->Reserve( cells->Size() ); - typedef OutMeshType::CellType::CellAutoPointer CellAutoPointer; - typedef itk::TriangleCell< OutMeshType::CellType > TriangleType; + using CellAutoPointer = OutMeshType::CellType::CellAutoPointer; + using TriangleType = itk::TriangleCell< OutMeshType::CellType >; while ( itCells != itCellsEnd ) { @@ -268,8 +266,8 @@ int main( int argc, char *argv[] ) // The name of the file to be read or written is passed with the // SetFileName() method. char outFileName[2048]; - vcl_strcpy(outFileName, outFileNameBase); - vcl_strcat(outFileName, "Ave.vtk"); + std::strcpy(outFileName, outFileNameBase); + std::strcat(outFileName, "Ave.vtk"); meshWriter->SetFileName(outFileName); meshWriter->SetInput(outMesh); @@ -280,7 +278,7 @@ int main( int argc, char *argv[] ) } catch( itk::ExceptionObject & excp ) { - vcl_cerr << excp << vcl_endl; + std::cerr << excp << std::endl; return EXIT_FAILURE; } @@ -300,12 +298,12 @@ int main( int argc, char *argv[] ) } // The name of the file to be read or written is passed with the // SetFileName() method. - vcl_strcpy(outFileName, outFileNameBase); - vcl_strcat(outFileName, "Basis"); + std::strcpy(outFileName, outFileNameBase); + std::strcat(outFileName, "Basis"); char fileCount[8]; sprintf(fileCount, "%d", j + 1); - vcl_strcat(outFileName, fileCount); - vcl_strcat(outFileName, ".vtk"); + std::strcat(outFileName, fileCount); + std::strcat(outFileName, ".vtk"); meshWriter->SetFileName(outFileName); meshWriter->SetInput(outMesh); @@ -316,7 +314,7 @@ int main( int argc, char *argv[] ) } catch( itk::ExceptionObject & excp ) { - vcl_cerr << excp << vcl_endl; + std::cerr << excp << std::endl; return EXIT_FAILURE; } } diff --git a/include/itkVectorFieldPCA.h b/include/itkVectorFieldPCA.h index 8b32218..75f5af4 100644 --- a/include/itkVectorFieldPCA.h +++ b/include/itkVectorFieldPCA.h @@ -55,11 +55,11 @@ template< typename TRealValueType = double > class ITK_EXPORT GaussianDistanceKernel : public KernelFunctionBase { public: - /** Standard class typedefs. */ - typedef GaussianDistanceKernel Self; - typedef KernelFunctionBase Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; + /** Standard class type alias. */ + using Self = GaussianDistanceKernel; + using Superclass = KernelFunctionBase; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; /** Run-time type information (and related methods). */ itkTypeMacro(GaussianDistanceKernel, KernelFunction); @@ -78,13 +78,13 @@ class ITK_EXPORT GaussianDistanceKernel : public KernelFunctionBase Pointer; - typedef SmartPointer ConstPointer; + ITK_DISALLOW_COPY_AND_ASSIGN(VectorFieldPCA); + + /** Standard class type alias. */ + using Self = VectorFieldPCA; + using Superclass = Object; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -117,20 +119,20 @@ class ITK_EXPORT VectorFieldPCA : public Object itkTypeMacro(VectorFieldPCA, Object); /** Type definitions for the PointSet. */ - typedef TPointSetType InputPointSetType; + using InputPointSetType = TPointSetType; /** Definitions for points of the PointSet. */ - typedef typename InputPointSetType::PointType InputPointType; + using InputPointType = typename InputPointSetType::PointType; /** Definitions for the PointsContainer. */ - typedef typename InputPointSetType::PointsContainer PointsContainer; - typedef typename PointsContainer::Iterator PointsContainerIterator; + using PointsContainer = typename InputPointSetType::PointsContainer; + using PointsContainerIterator = typename PointsContainer::Iterator; /** Pointer types for the PointSet. */ - typedef typename InputPointSetType::Pointer InputPointSetPointer; + using InputPointSetPointer = typename InputPointSetType::Pointer; /** Const Pointer type for the PointSet. */ - typedef typename InputPointSetType::ConstPointer InputPointSetConstPointer; + using InputPointSetConstPointer = typename InputPointSetType::ConstPointer; /** * \brief Input PointSet dimension @@ -139,21 +141,21 @@ class ITK_EXPORT VectorFieldPCA : public Object TPointSetType::PointDimension); /** type for the vector fields. */ - typedef vnl_matrix< TVectorFieldElementType > VectorFieldType; - typedef VectorContainer< unsigned int, VectorFieldType > VectorFieldSetType; + using VectorFieldType = vnl_matrix< TVectorFieldElementType >; + using VectorFieldSetType = VectorContainer< unsigned int, VectorFieldType >; - typedef typename VectorFieldSetType::Pointer VectorFieldSetTypePointer; - typedef typename VectorFieldSetType::ConstPointer VectorFieldSetTypeConstPointer; + using VectorFieldSetTypePointer = typename VectorFieldSetType::Pointer; + using VectorFieldSetTypeConstPointer = typename VectorFieldSetType::ConstPointer; /** types for the output. */ - typedef vnl_matrix< TPCType > MatrixType; - typedef vnl_vector< TPCType > VectorType; + using MatrixType = vnl_matrix< TPCType >; + using VectorType = vnl_vector< TPCType >; - typedef VectorContainer< unsigned int, MatrixType > BasisSetType; - typedef VectorContainer< unsigned int, VectorType > ResSetType; + using BasisSetType = VectorContainer< unsigned int, MatrixType >; + using ResSetType = VectorContainer< unsigned int, VectorType >; - typedef typename BasisSetType::Pointer BasisSetTypePointer; - typedef typename KernelFunctionType::Pointer KernelFunctionPointer; + using BasisSetTypePointer = typename BasisSetType::Pointer; + using KernelFunctionPointer = typename KernelFunctionType::Pointer; /** * \brief Set and get the input point set. @@ -194,17 +196,16 @@ class ITK_EXPORT VectorFieldPCA : public Object protected: VectorFieldPCA(); - virtual ~VectorFieldPCA() {}; - void PrintSelf(std::ostream& os, Indent indent) const; - - void KernelPCA(void); - void computeMomentumSCP(void); + ~VectorFieldPCA() override {}; + void PrintSelf(std::ostream& os, Indent indent) const override; -private: + /** Kernel PCA. */ + void KernelPCA(); - VectorFieldPCA(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented + /** Compute Momentum SCP. */ + void ComputeMomentumSCP(); +private: VectorType m_PCAEigenValues; BasisSetTypePointer m_BasisVectors; @@ -212,7 +213,7 @@ class ITK_EXPORT VectorFieldPCA : public Object InputPointSetPointer m_PointSet; KernelFunctionPointer m_KernelFunction; - // problem dimensions + // Problem dimensions unsigned int m_ComponentCount; unsigned int m_SetSize; unsigned int m_VectorDimCount; @@ -229,6 +230,8 @@ class ITK_EXPORT VectorFieldPCA : public Object } // end namespace itk +#ifndef ITK_MANUAL_INSTANTIATION #include "itkVectorFieldPCA.hxx" +#endif #endif diff --git a/include/itkVectorFieldPCA.hxx b/include/itkVectorFieldPCA.hxx index 95510e9..9a6871f 100644 --- a/include/itkVectorFieldPCA.hxx +++ b/include/itkVectorFieldPCA.hxx @@ -22,14 +22,11 @@ #include "itkVectorFieldPCA.h" #include "vnl/algo/vnl_symmetric_eigensystem.h" #include "vnl/vnl_c_vector.h" -#include "vnl/vnl_math.h" +#include "itkMath.h" namespace itk { -/** - * Constructor - */ template< typename TVectorFieldElementType, typename TPCType, @@ -39,22 +36,17 @@ template< class TPointSetType > VectorFieldPCA -::VectorFieldPCA() +::VectorFieldPCA() : + m_BasisVectors( BasisSetType::New() ), + m_ComponentCount( 0 ), + m_SetSize( 0 ), + m_VectorDimCount( 0 ), + m_VertexCount( 0 ), + m_PointDim( 0 ), + m_PCACalculated( false ) { - m_PCACalculated = false; - m_SetSize = 0; - m_VectorDimCount = 0; - m_PointDim = 0; - m_ComponentCount = 0; - m_KernelFunction = NULL; - m_BasisVectors = BasisSetType::New(); - m_VectorFieldSet = NULL; - m_PointSet = NULL; } -/** - * - */ template< typename TVectorFieldElementType, typename TPCType, @@ -65,71 +57,9 @@ template< > void VectorFieldPCA -::PrintSelf(std::ostream& os, Indent indent) const +::Compute() { - Superclass::PrintSelf(os,indent); - - if (this->m_VectorFieldSet.IsNull()) - { - os << indent << "Vector Field Set Empty" << std::endl; - } - else - { - os << indent << "Vector Field Set Count: " << this->m_VectorFieldSet->Size() << std::endl; - } - - os << indent << "Component Count: " << this->m_ComponentCount << std::endl; - os << indent << "Eigenvalues: " << this->m_PCAEigenValues << std::endl; - if (this->m_BasisVectors.IsNull()) - { - os << indent << "Basis Vector Empty" << std::endl; - } - else - { - os << indent << "Basis Vector Count: " << this->m_BasisVectors->Size() << std::endl; - if (this->m_BasisVectors->Size()) - { - MatrixType& thisBasis = this->m_BasisVectors->ElementAt(0); - os << indent << "Basis Vector Dimensions: " << thisBasis.rows() << "x" << thisBasis.cols() << std::endl; - } - } - - if (this->m_PointSet.IsNull()) - { - os << indent << "PointSet Empty" << std::endl; - } - else - { - os << indent << "PointSet is " << m_PointSet->GetNumberOfPoints() << - "x" << this->m_PointSet->PointDimension << std::endl; - } - - if (this->m_KernelFunction.IsNull()) - { - os << indent << "Kernel Function not set" << std::endl; - } - else - { - os << indent << "KernelFunction is set." << std::endl; - } -} - -/** - * Compute the principal components - */ -template< - typename TVectorFieldElementType, - typename TPCType, - typename TPointSetPixelType, - typename TPointSetCoordRepType, - typename KernelFunctionType, - class TPointSetType - > -void -VectorFieldPCA -::Compute(void) -{ - // check parameters + // Check parameters if (!m_VectorFieldSet || !m_VectorFieldSet->Size()) { itkExceptionMacro("Vector Field Set not specified."); @@ -144,13 +74,13 @@ VectorFieldPCAElementAt(0); m_VectorDimCount = firstField.rows(); m_PointDim = firstField.cols(); - // check all vector dimensions in the set + // Check all vector dimensions in the set for (unsigned int i = 1; i < m_VectorFieldSet->Size(); i++) { VectorFieldType& thisField = m_VectorFieldSet->ElementAt(i); @@ -190,13 +120,13 @@ VectorFieldPCAComputeMomentumSCP(); + this->KernelPCA(); - // save only the desired eigenvalues + // Save only the desired eigenvalues m_PCAEigenValues = m_PCAEigenValues.extract(m_ComponentCount); - // save only the desired eigenvectors + // Save only the desired eigenvectors m_V0 = m_V0.extract(m_V0.rows(), m_ComponentCount); m_BasisVectors->Reserve(m_ComponentCount); @@ -225,10 +155,6 @@ VectorFieldPCA void VectorFieldPCA -::computeMomentumSCP(void) +::ComputeMomentumSCP() { - VectorFieldType accum; accum.set_size(m_VectorDimCount, m_PointDim); accum = 0.0; - // determine the average of the vector field over the set + // Determine the average of the vector field over the set for (unsigned k = 0; k < m_SetSize; k++) { accum += m_VectorFieldSet->ElementAt(k); @@ -261,7 +186,7 @@ VectorFieldPCA void VectorFieldPCA -::KernelPCA(void) +::KernelPCA() { VectorType rowMeans(m_SetSize); - unsigned int k, l; - for (k = 0; k < m_SetSize; k++) + for (unsigned int k = 0; k < m_SetSize; k++) { rowMeans(k) = m_K.get_row(k).mean(); } TPCType meanOfMeans = rowMeans.mean(); MatrixType K0(m_K - meanOfMeans); - for (k = 0; k < m_SetSize; k++) + for (unsigned int k = 0; k < m_SetSize; k++) { - for (l = 0; l < m_SetSize; l++) + for (unsigned int l = 0; l < m_SetSize; l++) { K0(k, l) -= rowMeans(k) + rowMeans(l); } @@ -354,20 +275,78 @@ VectorFieldPCA eigs(K0); m_PCAEigenValues = eigs.D.diagonal(); - // eigs come out in ascending order, reorder them + + // Eigenvalues come out in ascending order, reorder them m_PCAEigenValues.flip(); - // reorder eigenvectors + // Reorder eigenvectors m_V0 = eigs.V; m_V0.fliplr(); const double eigenvalue_epsilon = 1.0e-10; - for (k = 0; k < m_SetSize; k++) + for (unsigned int k = 0; k < m_SetSize; k++) { - m_V0.scale_column(k, 1.0 / vcl_sqrt(m_PCAEigenValues(k) + eigenvalue_epsilon)); + m_V0.scale_column(k, 1.0 / std::sqrt(m_PCAEigenValues(k) + eigenvalue_epsilon)); } } +template< + typename TVectorFieldElementType, + typename TPCType, + typename TPointSetPixelType, + typename TPointSetCoordRepType, + typename KernelFunctionType, + class TPointSetType + > +void +VectorFieldPCA +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf( os, indent ); + + os << indent << "PCAEigenvalues: " << this->m_PCAEigenValues << std::endl; + + if( this->m_BasisVectors.IsNotNull() ) + { + os << indent << "Basis Vector count: " << this->m_BasisVectors->Size() + << std::endl; + if( this->m_BasisVectors->Size() ) + { + MatrixType& thisBasis = this->m_BasisVectors->ElementAt(0); + os << indent << "Basis Vector dimensions: " << thisBasis.rows() << "x" + << thisBasis.cols() << std::endl; + } + } + itkPrintSelfObjectMacro( BasisVectors ); + + if( this->m_VectorFieldSet.IsNotNull() ) + { + os << indent << "Vector Field Set count: " << this->m_VectorFieldSet->Size() + << std::endl; + } + itkPrintSelfObjectMacro( VectorFieldSet ); + + if( this->m_PointSet.IsNotNull() ) + { + os << indent << "PointSet dimensions: " << m_PointSet->GetNumberOfPoints() + << "x" << this->m_PointSet->PointDimension << std::endl; + } + itkPrintSelfObjectMacro( PointSet ); + + itkPrintSelfObjectMacro( KernelFunction ); + + os << indent << "ComponentCount: " << this->m_ComponentCount << std::endl; + os << indent << "SetSize: " << this->m_SetSize << std::endl; + os << indent << "VectorDimCount: " << this->m_VectorDimCount << std::endl; + os << indent << "VertexCount: " << this->m_VertexCount << std::endl; + os << indent << "PointDim: " << this->m_PointDim << std::endl; + + os << indent << "V0 : " << this->m_V0 << std::endl; + os << indent << "AveVectorField: " << this->m_AveVectorField << std::endl; + os << indent << "K: " << this->m_K << std::endl; + + os << indent << "PCACalculated: " << this->m_PCACalculated << std::endl; +} } // end namespace itk #endif diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..3e9df8e --- /dev/null +++ b/setup.py @@ -0,0 +1,54 @@ +# -*- coding: utf-8 -*- +from __future__ import print_function +from os import sys + +try: + from skbuild import setup +except ImportError: + print('scikit-build is required to build from source.', file=sys.stderr) + print('Please run:', file=sys.stderr) + print('', file=sys.stderr) + print(' python -m pip install scikit-build') + sys.exit(1) + +setup( + name='itk-principalcomponentsanalysis', + version='1.0.0', + author='Michael Bowers and Laurent Younes', + author_email='mbowers@cis.jhu.edu', + packages=['itk'], + package_dir={'itk': 'itk'}, + download_url=r'https://github.com/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis', + description=r'ITK classes for the analysis of the principal components of data sets, optionally point data associated with the vertices of a mesh.', + long_description='ITKPrincipalComponentsAnalysis provides an' + 'implementation of standard principal components analysis' + 'algorithms for use on scalar or vector data sets.\n' + 'Please refer to:\n' + 'M. Bowers and L. Younes, “Principal Components Analysis of Scalar, Vector, and Mesh Vertex Data.”,' + 'Insight Journal, January-December 2012, http://hdl.handle.net/10380/3386', + classifiers=[ + "License :: OSI Approved :: Apache Software License", + "Programming Language :: Python", + "Programming Language :: C++", + "Development Status :: 4 - Beta", + "Intended Audience :: Developers", + "Intended Audience :: Education", + "Intended Audience :: Healthcare Industry", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering", + "Topic :: Scientific/Engineering :: Medical Science Apps.", + "Topic :: Scientific/Engineering :: Information Analysis", + "Topic :: Software Development :: Libraries", + "Operating System :: Android", + "Operating System :: Microsoft :: Windows", + "Operating System :: POSIX", + "Operating System :: Unix", + "Operating System :: MacOS" + ], + license='Apache', + keywords='ITK PCA Mesh', + url=r'https://github.com/InsightSoftwareConsortium/ITKPrincipalComponentsAnalysis', + install_requires=[ + r'itk' + ] + ) \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8e8ff2d..b446f4f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,16 +1,16 @@ itk_module_test() -SET(kit PrincipalComponentsAnalysis) -SET(${kit}Tests +SET(PCA PrincipalComponentsAnalysis) +SET(${PCA}Tests itkVectorKernelPCATest.cxx ) -SET(TEST_DATA_ROOT ${${kit}_SOURCE_DIR}/Data) +SET(TEST_DATA_ROOT ${${PCA}_SOURCE_DIR}/Data) -CreateTestDriver(${kit} "${${kit}-Test_LIBRARIES}" "${${kit}Tests}") +CreateTestDriver(${PCA} "${${PCA}-Test_LIBRARIES}" "${${PCA}Tests}") itk_add_test(NAME itkVectorKernelPCATest - COMMAND ${kit}TestDriver itkVectorKernelPCATest + COMMAND ${PCA}TestDriver itkVectorKernelPCATest DATA{${TEST_DATA_ROOT}/PCATestSurface.vtk} DATA{${TEST_DATA_ROOT}/PCATestSurface_alpha0_01.vtk} DATA{${TEST_DATA_ROOT}/PCATestSurface_alpha0_02.vtk} diff --git a/test/itkVectorKernelPCATest.cxx b/test/itkVectorKernelPCATest.cxx index 62d608d..9763ffd 100644 --- a/test/itkVectorKernelPCATest.cxx +++ b/test/itkVectorKernelPCATest.cxx @@ -19,274 +19,265 @@ #include "itkMesh.h" #include "itkMeshFileReader.h" #include "itkVectorFieldPCA.h" +#include "itkTestingMacros.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_vector.h" -int showUsage(const char* programName) -{ - std::cerr << "USAGE: " << programName << std::endl; - std::cerr << " ... " << std::endl; - std::cerr << "\t\tN must be greater than 3" << std::endl; - return EXIT_FAILURE; -} -int itkVectorKernelPCATest(int argc, char *argv[]) +template< typename TPixel, typename TMesh, typename TVectorContainer > +int ParseVectorFields( std::vector< std::string > vectorFieldFilenames, typename TVectorContainer::Pointer vectorFieldSet ) { -#define MIN_ARG_COUNT 6 -#define FIRST_VECTOR_FIELD_ARG 2 - if(argc < MIN_ARG_COUNT) - return (showUsage(argv[0])); - - unsigned int pcaCount = 3; - double kernelSigma = 6.25; - - typedef double PointDataType; - typedef itk::Array PointDataVectorType; - typedef PointDataVectorType PixelType; - typedef double CoordRep; - const unsigned int Dimension = 3; - -// typedef float PCAResultsType; - typedef double PCAResultsType; + int testStatus = EXIT_SUCCESS; - // Declare the type of the input mesh - typedef itk::Mesh InMeshType; - - // Declare the type of the kernel function class - // typedef itk::GaussianDistanceKernel KernelType; - typedef itk::GaussianDistanceKernel KernelType; + using ReaderType = itk::MeshFileReader< TMesh >; + typename ReaderType::Pointer meshReader = ReaderType::New(); - // Declare the type of the PCA calculator - typedef itk::VectorFieldPCA< PointDataType, PCAResultsType, - PixelType, CoordRep, KernelType, InMeshType > - PCACalculatorType; + unsigned int fieldSetCount = vectorFieldFilenames.size(); + vectorFieldSet->Reserve( fieldSetCount ); - // Here we recover the file names from the command line arguments - const char* inMeshFile = argv[1]; - - // We can now instantiate the types of the reader/writer. - typedef itk::MeshFileReader< InMeshType > ReaderType; - - // create readers/writers - ReaderType::Pointer meshReader = ReaderType::New(); - // The name of the file to be read or written is passed with the - // SetFileName() method. - meshReader->SetFileName( inMeshFile ); + typename TVectorContainer::Element vectorField; - try - { - meshReader->Update(); - } - catch( itk::ExceptionObject & excp ) - { - std::cerr << "Error reading mesh file " << inMeshFile << std::endl; - std::cerr << excp << std::endl; - } - - // get the objects - InMeshType::Pointer mesh = meshReader->GetOutput(); - - std::cout << "Vertex Count: " << - mesh->GetNumberOfPoints() << std::endl; - std::cout << "Cell Count: " << - mesh->GetNumberOfCells() << std::endl; - - const char* vectorFieldName; - PCACalculatorType::VectorFieldType vectorField; - - // should know vector field dimensions now unsigned int vectorFieldDim = 0; unsigned int vectorFieldCount = 0; - // how many vector field sets? - unsigned int fieldSetCount = argc - FIRST_VECTOR_FIELD_ARG; - - PCACalculatorType::VectorFieldSetTypePointer vectorFieldSet = - PCACalculatorType::VectorFieldSetType::New(); - - vectorFieldSet->Reserve(fieldSetCount); - unsigned int setIx = 0; - for (int i = FIRST_VECTOR_FIELD_ARG; i < argc; i++) + for( unsigned int i = 0; i < fieldSetCount; i++ ) { - vectorFieldName = argv[i]; - // The name of the file to be read or written is passed with the - // SetFileName() method. - meshReader->SetFileName( vectorFieldName ); + std::string vectorFieldName = vectorFieldFilenames[i]; - try - { - meshReader->Update(); - } - catch( itk::ExceptionObject & excp ) - { - std::cerr << "Error reading mesh field file " << vectorFieldName << std::endl; - std::cerr << excp << std::endl; - } + meshReader->SetFileName( vectorFieldName ); + + TRY_EXPECT_NO_EXCEPTION( meshReader->Update() ); - // get the objects - InMeshType::Pointer meshWithField = meshReader->GetOutput(); + // Get the objects + typename TMesh::Pointer meshWithField = meshReader->GetOutput(); - InMeshType::PointDataContainerPointer pointData = meshWithField->GetPointData(); - if (setIx == 0) + typename TMesh::PointDataContainerPointer pointData = meshWithField->GetPointData(); + if( setIx == 0 ) { vectorFieldCount = pointData->Size(); - if (vectorFieldCount) + if( vectorFieldCount ) { - PixelType oneDataSetVal = pointData->GetElement(0); + TPixel oneDataSetVal = pointData->GetElement( 0 ); vectorFieldDim = oneDataSetVal.size(); } - if (vectorFieldCount != meshWithField->GetNumberOfPoints()) + if( vectorFieldCount != meshWithField->GetNumberOfPoints() ) { + std::cerr << "Test failed!" << std::endl; std::cerr << "Vector field count (" << vectorFieldCount << - ") doesn't match mesh vertext count (" << - meshWithField->GetNumberOfPoints() << ")." << std::endl; - exit (EXIT_FAILURE); + ") doesn't match mesh vertext count (" << + meshWithField->GetNumberOfPoints() << ")." << std::endl; + testStatus = EXIT_FAILURE; } - vectorField.set_size(vectorFieldCount, vectorFieldDim); + vectorField.set_size( vectorFieldCount, vectorFieldDim ); } else { - if (vectorFieldDim != vectorField.cols() || - vectorFieldCount != meshWithField->GetNumberOfPoints()) + if( vectorFieldDim != vectorField.cols() || + vectorFieldCount != meshWithField->GetNumberOfPoints() ) { - std::cerr << "Unexpected dimensions in vector field file " << vectorFieldName << std::endl; - std::cerr << "\tExpected " << vectorFieldCount << " x " << vectorFieldDim; - std::cerr << "\t, got " << meshWithField->GetNumberOfPoints() << " x " << vectorField.cols() << std::endl; - exit(1); + std::cerr << "Test failed!" << std::endl; + std::cerr << "Unexpected dimensions in vector field file " + << vectorFieldName << std::endl; + std::cerr << "Expected: " << vectorFieldCount << " x " << vectorFieldDim + << ", but got " << meshWithField->GetNumberOfPoints() << " x " + << vectorField.cols() << std::endl; + testStatus = EXIT_FAILURE; } } - for (unsigned int k = 0; k < pointData->Size(); k++) + for( unsigned int k = 0; k < pointData->Size(); k++ ) { - PixelType oneDataSetVal = pointData->GetElement(k); - vectorField.set_row(k, oneDataSetVal); + TPixel oneDataSetVal = pointData->GetElement(k); + vectorField.set_row( k, oneDataSetVal ); } - vectorFieldSet->SetElement(setIx++, vectorField); + vectorFieldSet->SetElement( setIx++, vectorField ); } - PCACalculatorType::Pointer pcaCalc = PCACalculatorType::New(); - - std::cout << "Name of Class = " << pcaCalc->GetNameOfClass() << std::endl; + return testStatus; +} - std::cout << "Test Print() = " << std::endl; - pcaCalc->Print( std::cout ); - // try to Compute before setting much of anything - expect failure - try +int itkVectorKernelPCATest( int argc, char *argv[] ) +{ + if( argc < 6 ) { - pcaCalc->Compute(); - std::cerr << "Failed to throw expected exception" << std::endl; + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << argv[0] + << " ... " << std::endl; + std::cerr << "where N must be greater than 3" << std::endl; return EXIT_FAILURE; } - catch( itk::ExceptionObject & excp ) - { - std::cout << "SUCCESSFULLY caught expected exception" << std::endl; - std::cout << excp << std::endl; - } - // set user variables + + int testStatus = EXIT_SUCCESS; + + const unsigned int Dimension = 3; + + using PointDataType = double; + using PointDataVectorType = itk::Array< PointDataType >; + using PixelType = PointDataVectorType; + using CoordRep = double; + + using PCAResultsType = double; + + // Declare the type of the input mesh + using MeshType = itk::Mesh< PixelType, Dimension >; + + // Declare the type of the kernel function class + using KernelType = itk::GaussianDistanceKernel< CoordRep >; + + // Declare the type of the PCA calculator + using PCACalculatorType = itk::VectorFieldPCA< PointDataType, PCAResultsType, PixelType, + CoordRep, KernelType, MeshType >; + + // Instantiate the reader + using ReaderType = itk::MeshFileReader< MeshType >; + ReaderType::Pointer meshReader = ReaderType::New(); + + meshReader->SetFileName( argv[1] ); + + TRY_EXPECT_NO_EXCEPTION( meshReader->Update() ); + + + // Get the input mesh + MeshType::Pointer mesh = meshReader->GetOutput(); + + + PCACalculatorType::Pointer pcaCalc = PCACalculatorType::New(); + + EXERCISE_BASIC_OBJECT_METHODS( pcaCalc, VectorFieldPCA, Object ); + + + // Test exception when trying to compute before setting much of anything + TRY_EXPECT_EXCEPTION( pcaCalc->Compute() ); + + + // Set user variables + unsigned int pcaCount = 3; pcaCalc->SetComponentCount( pcaCount ); - // - // Now connect the input. - // + // Connect the input pcaCalc->SetPointSet( mesh ); - // set vector fields - pcaCalc->SetVectorFieldSet( vectorFieldSet ); - // - // Now verify that it runs fine. - // - try - { - pcaCalc->Compute(); - std::cout << "SUCCESSFULLY ran non-Kernel version" << std::endl; - } - catch( itk::ExceptionObject & excp ) + TEST_SET_GET_VALUE( mesh, pcaCalc->GetPointSet() ); + + // Set vector fields + + PCACalculatorType::VectorFieldType vectorField; + + // Should know vector field dimensions now + unsigned int vectorFieldDim = 0; + + // how many vector field sets? + std::vector< std::string > vectorFieldFilenames; + unsigned int firstVectorFieldIdx = 2; + unsigned int fieldSetCount = argc - firstVectorFieldIdx; + for( unsigned int i = firstVectorFieldIdx; i < firstVectorFieldIdx + fieldSetCount; i++ ) { - std::cout << "Failed to run non-Kernel version" << std::endl; - std::cerr << excp << std::endl; - return EXIT_FAILURE; + vectorFieldFilenames.push_back( argv[i] ); } + PCACalculatorType::VectorFieldSetTypePointer vectorFieldSet = + PCACalculatorType::VectorFieldSetType::New(); + + testStatus = ParseVectorFields< PixelType, MeshType, PCACalculatorType::VectorFieldSetType >( + vectorFieldFilenames, vectorFieldSet ); + + pcaCalc->SetVectorFieldSet( vectorFieldSet ); + TEST_SET_GET_VALUE( vectorFieldSet, pcaCalc->GetVectorFieldSet() ); + + + // Execute the PCA calculator + TRY_EXPECT_NO_EXCEPTION( pcaCalc->Compute() ); + + + double kernelSigma = 6.25; KernelType::Pointer distKernel = KernelType::New(); distKernel->SetKernelSigma( kernelSigma ); pcaCalc->SetKernelFunction( distKernel ); - std::cout << "PCA Calculator All Set Up: Print() = " << std::endl; - pcaCalc->Print( std::cout ); std::ofstream debugOut; - debugOut.precision(15); + debugOut.precision( 15 ); - // get the output and perform basic checks - for (unsigned int j = 0; j < pcaCalc->GetComponentCount(); j++) + // Get the output and perform basic checks + unsigned int computedNumberOfAverageVectorFieldCols = + pcaCalc->GetAveVectorField().cols(); + + unsigned int computedNumberOfAverageVectorFieldRows = + pcaCalc->GetAveVectorField().rows(); + + for( unsigned int j = 0; j < pcaCalc->GetComponentCount(); j++ ) { - if ((pcaCalc->GetBasisVectors()->GetElement(j)).cols() != vectorFieldDim) + unsigned int expectedBasisVectorCols = + pcaCalc->GetVectorFieldSet()->GetElement(j).cols(); + unsigned int computedBasisVectorCols = + pcaCalc->GetBasisVectors()->GetElement(j).cols(); + if( computedBasisVectorCols != expectedBasisVectorCols ) { - std::cout << "Basis Vector Results Failed Dimension check:" << std::endl; - std::cout << "Expected: " << vectorFieldDim << std::endl - << " columns. Got: " - << (pcaCalc->GetBasisVectors()->GetElement(j)).cols() - << std::endl; - return EXIT_FAILURE; + std::cout << "Test failed!" << std::endl; + std::cout << "Error in GetBasisVectors() dimension check at index [" << j << "]" + << std::endl; + std::cout << "Expected: " << expectedBasisVectorCols << " columns, but got: " + << computedBasisVectorCols << std::endl; + testStatus = EXIT_FAILURE; } - if ((pcaCalc->GetBasisVectors()->GetElement(j)).rows() != vectorFieldCount) + + unsigned int expectedBasisVectorRows = + pcaCalc->GetVectorFieldSet()->GetElement(j).rows(); + unsigned int computedBasisVectorRows = + pcaCalc->GetBasisVectors()->GetElement(j).rows(); + if( computedBasisVectorRows != expectedBasisVectorRows ) { - std::cout << "Basis Vector Results Failed Dimension check:" << std::endl; - std::cout << "Expected: " << vectorFieldDim << std::endl - << " rows. Got: " - << pcaCalc->GetBasisVectors()->GetElement(j).rows() - << std::endl; - return EXIT_FAILURE; + std::cout << "Test failed!" << std::endl; + std::cout << "Error in GetBasisVectors() dimension check at index [" << j << "]" + << std::endl; + std::cout << "Expected: " << expectedBasisVectorRows << " row, but got: " + << computedBasisVectorRows << std::endl; + testStatus = EXIT_FAILURE; } - } - if (pcaCalc->GetPCAEigenValues().size() != pcaCount) - { - std::cout << "Eigenvalue Vector Results Failed Dimension check:" << std::endl; - std::cout << "Expected: " << pcaCount << std::endl - << ". Got: " - << pcaCalc->GetPCAEigenValues().size() - << std::endl; - return EXIT_FAILURE; - } + if( computedNumberOfAverageVectorFieldCols != expectedBasisVectorCols ) + { + std::cout << "Test failed!" << std::endl; + std::cout << "Error in GetAveVectorField() dimension check at index [" << j << "]" + << std::endl; + std::cout << "Expected: " << vectorFieldDim << " columns, but got: " + << computedNumberOfAverageVectorFieldCols << std::endl; + testStatus = EXIT_FAILURE; + } - if (pcaCalc->GetAveVectorField().cols() != vectorFieldDim) - { - std::cout << "Average Vector Field Results Failed Dimension check:" << std::endl; - std::cout << "Expected: " << vectorFieldDim << std::endl - << " columns. Got: " - << pcaCalc->GetAveVectorField().cols() - << std::endl; - return EXIT_FAILURE; + if( computedNumberOfAverageVectorFieldRows != expectedBasisVectorRows ) + { + std::cout << "Test failed!" << std::endl; + std::cout << "Error in GetAveVectorField() dimension check at index [" << j << "]" + << std::endl; + std::cout << "Expected: " << vectorFieldDim << " rows, but got: " + << computedNumberOfAverageVectorFieldRows << std::endl; + testStatus = EXIT_FAILURE; + } } - if (pcaCalc->GetAveVectorField().rows() != vectorFieldCount) + + unsigned int computedNumberOfEigenValueVectors = + pcaCalc->GetPCAEigenValues().size(); + if( computedNumberOfEigenValueVectors != pcaCount ) { - std::cout << "Average Vector Field Failed Dimension check:" << std::endl; - std::cout << "Expected: " << vectorFieldDim << std::endl - << " rows. Got: " - << pcaCalc->GetAveVectorField().rows() - << std::endl; - return EXIT_FAILURE; + std::cout << "Test failed!" << std::endl; + std::cout << "Error in GetPCAEigenValues() dimension check." << std::endl; + std::cout << "Expected: " << pcaCount << ", but got: " + << computedNumberOfEigenValueVectors << std::endl; + testStatus = EXIT_FAILURE; } - // set the requested input count greater than the number of vector field sets + // Test exception when trying to compute with a requested input count greater + // than the number of vector field sets pcaCalc->SetComponentCount( fieldSetCount + 1 ); - // try to Compute - expect failure - try - { - pcaCalc->Compute(); - std::cerr << "Failed to throw expected exception" << std::endl; - return EXIT_FAILURE; - } - catch( itk::ExceptionObject & excp ) - { - std::cout << "SUCCESSFULLY caught expected exception" << std::endl; - std::cout << excp << std::endl; - } - return EXIT_SUCCESS; + TRY_EXPECT_EXCEPTION( pcaCalc->Compute() ); + + + std::cout << "Test finished." << std::endl; + return testStatus; }