diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..96e9401 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*~ +build +*.DS_Store +.svn diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..5fea7e0 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,25 @@ +cmake_minimum_required(VERSION 2.8.9) + +project(PETCPhantom) + +#----------------------------------------------------------------------------- +# Extension meta-information +set(EXTENSION_HOMEPAGE "http://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Extensions/PETCPhantom") +set(EXTENSION_CATEGORY "Quantification") +set(EXTENSION_CONTRIBUTORS "Christian Bauer (University of Iowa), Ethan Ulrich (University of Iowa), Andrey Fedorov (SPL), Reinhard R. Beichel (University of Iowa), John Buatti (University of Iowa)") +set(EXTENSION_DESCRIPTION "The PETCPhantom Extension allows measurement of calibration and uniformity in a cylinder phantom PET scan.") +set(EXTENSION_ICONURL "https://raw.githubusercontent.com/QIICR/PETPhantomAnalysis/master/PETCPhantom.png") +set(EXTENSION_SCREENSHOTURLS "") + +#----------------------------------------------------------------------------- +# Extension dependencies +find_package(Slicer REQUIRED) +include(${Slicer_USE_FILE}) + +#----------------------------------------------------------------------------- +# Extension modules +add_subdirectory(PETPhantomAnalysisCLI) +add_subdirectory(PETPhantomAnalysis) + +#----------------------------------------------------------------------------- +include(${Slicer_EXTENSION_CPACK}) diff --git a/License.txt b/License.txt new file mode 100644 index 0000000..a640153 --- /dev/null +++ b/License.txt @@ -0,0 +1,199 @@ + +For more information, please see: + + http://www.slicer.org + +The 3D Slicer license below is a BSD style license, with extensions +to cover contributions and other issues specific to 3D Slicer. + + +3D Slicer Contribution and Software License Agreement ("Agreement") +Version 1.0 (December 20, 2005) + +This Agreement covers contributions to and downloads from the 3D +Slicer project ("Slicer") maintained by The Brigham and Women's +Hospital, Inc. ("Brigham"). Part A of this Agreement applies to +contributions of software and/or data to Slicer (including making +revisions of or additions to code and/or data already in Slicer). Part +B of this Agreement applies to downloads of software and/or data from +Slicer. Part C of this Agreement applies to all transactions with +Slicer. If you distribute Software (as defined below) downloaded from +Slicer, all of the paragraphs of Part B of this Agreement must be +included with and apply to such Software. + +Your contribution of software and/or data to Slicer (including prior +to the date of the first publication of this Agreement, each a +"Contribution") and/or downloading, copying, modifying, displaying, +distributing or use of any software and/or data from Slicer +(collectively, the "Software") constitutes acceptance of all of the +terms and conditions of this Agreement. If you do not agree to such +terms and conditions, you have no right to contribute your +Contribution, or to download, copy, modify, display, distribute or use +the Software. + +PART A. CONTRIBUTION AGREEMENT - License to Brigham with Right to +Sublicense ("Contribution Agreement"). + +1. As used in this Contribution Agreement, "you" means the individual + contributing the Contribution to Slicer and the institution or + entity which employs or is otherwise affiliated with such + individual in connection with such Contribution. + +2. This Contribution Agreement applies to all Contributions made to + Slicer, including without limitation Contributions made prior to + the date of first publication of this Agreement. If at any time you + make a Contribution to Slicer, you represent that (i) you are + legally authorized and entitled to make such Contribution and to + grant all licenses granted in this Contribution Agreement with + respect to such Contribution; (ii) if your Contribution includes + any patient data, all such data is de-identified in accordance with + U.S. confidentiality and security laws and requirements, including + but not limited to the Health Insurance Portability and + Accountability Act (HIPAA) and its regulations, and your disclosure + of such data for the purposes contemplated by this Agreement is + properly authorized and in compliance with all applicable laws and + regulations; and (iii) you have preserved in the Contribution all + applicable attributions, copyright notices and licenses for any + third party software or data included in the Contribution. + +3. Except for the licenses granted in this Agreement, you reserve all + right, title and interest in your Contribution. + +4. You hereby grant to Brigham, with the right to sublicense, a + perpetual, worldwide, non-exclusive, no charge, royalty-free, + irrevocable license to use, reproduce, make derivative works of, + display and distribute the Contribution. If your Contribution is + protected by patent, you hereby grant to Brigham, with the right to + sublicense, a perpetual, worldwide, non-exclusive, no-charge, + royalty-free, irrevocable license under your interest in patent + rights covering the Contribution, to make, have made, use, sell and + otherwise transfer your Contribution, alone or in combination with + any other code. + +5. You acknowledge and agree that Brigham may incorporate your + Contribution into Slicer and may make Slicer available to members + of the public on an open source basis under terms substantially in + accordance with the Software License set forth in Part B of this + Agreement. You further acknowledge and agree that Brigham shall + have no liability arising in connection with claims resulting from + your breach of any of the terms of this Agreement. + +6. YOU WARRANT THAT TO THE BEST OF YOUR KNOWLEDGE YOUR CONTRIBUTION + DOES NOT CONTAIN ANY CODE THAT REQURES OR PRESCRIBES AN "OPEN + SOURCE LICENSE" FOR DERIVATIVE WORKS (by way of non-limiting + example, the GNU General Public License or other so-called + "reciprocal" license that requires any derived work to be licensed + under the GNU General Public License or other "open source + license"). + +PART B. DOWNLOADING AGREEMENT - License from Brigham with Right to +Sublicense ("Software License"). + +1. As used in this Software License, "you" means the individual + downloading and/or using, reproducing, modifying, displaying and/or + distributing the Software and the institution or entity which + employs or is otherwise affiliated with such individual in + connection therewith. The Brigham and Women?s Hospital, + Inc. ("Brigham") hereby grants you, with right to sublicense, with + respect to Brigham's rights in the software, and data, if any, + which is the subject of this Software License (collectively, the + "Software"), a royalty-free, non-exclusive license to use, + reproduce, make derivative works of, display and distribute the + Software, provided that: + +(a) you accept and adhere to all of the terms and conditions of this +Software License; + +(b) in connection with any copy of or sublicense of all or any portion +of the Software, all of the terms and conditions in this Software +License shall appear in and shall apply to such copy and such +sublicense, including without limitation all source and executable +forms and on any user documentation, prefaced with the following +words: "All or portions of this licensed product (such portions are +the "Software") have been obtained under license from The Brigham and +Women's Hospital, Inc. and are subject to the following terms and +conditions:" + +(c) you preserve and maintain all applicable attributions, copyright +notices and licenses included in or applicable to the Software; + +(d) modified versions of the Software must be clearly identified and +marked as such, and must not be misrepresented as being the original +Software; and + +(e) you consider making, but are under no obligation to make, the +source code of any of your modifications to the Software freely +available to others on an open source basis. + +2. The license granted in this Software License includes without + limitation the right to (i) incorporate the Software into + proprietary programs (subject to any restrictions applicable to + such programs), (ii) add your own copyright statement to your + modifications of the Software, and (iii) provide additional or + different license terms and conditions in your sublicenses of + modifications of the Software; provided that in each case your use, + reproduction or distribution of such modifications otherwise + complies with the conditions stated in this Software License. + +3. This Software License does not grant any rights with respect to + third party software, except those rights that Brigham has been + authorized by a third party to grant to you, and accordingly you + are solely responsible for (i) obtaining any permissions from third + parties that you need to use, reproduce, make derivative works of, + display and distribute the Software, and (ii) informing your + sublicensees, including without limitation your end-users, of their + obligations to secure any such required permissions. + +4. The Software has been designed for research purposes only and has + not been reviewed or approved by the Food and Drug Administration + or by any other agency. YOU ACKNOWLEDGE AND AGREE THAT CLINICAL + APPLICATIONS ARE NEITHER RECOMMENDED NOR ADVISED. Any + commercialization of the Software is at the sole risk of the party + or parties engaged in such commercialization. You further agree to + use, reproduce, make derivative works of, display and distribute + the Software in compliance with all applicable governmental laws, + regulations and orders, including without limitation those relating + to export and import control. + +5. The Software is provided "AS IS" and neither Brigham nor any + contributor to the software (each a "Contributor") shall have any + obligation to provide maintenance, support, updates, enhancements + or modifications thereto. BRIGHAM AND ALL CONTRIBUTORS SPECIFICALLY + DISCLAIM ALL EXPRESS AND IMPLIED WARRANTIES OF ANY KIND INCLUDING, + BUT NOT LIMITED TO, ANY WARRANTIES OF MERCHANTABILITY, FITNESS FOR + A PARTICULAR PURPOSE AND NON-INFRINGEMENT. IN NO EVENT SHALL + BRIGHAM OR ANY CONTRIBUTOR BE LIABLE TO ANY PARTY FOR DIRECT, + INDIRECT, SPECIAL, INCIDENTAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES + HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY ARISING IN ANY WAY + RELATED TO THE SOFTWARE, EVEN IF BRIGHAM OR ANY CONTRIBUTOR HAS + BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. TO THE MAXIMUM + EXTENT NOT PROHIBITED BY LAW OR REGULATION, YOU FURTHER ASSUME ALL + LIABILITY FOR YOUR USE, REPRODUCTION, MAKING OF DERIVATIVE WORKS, + DISPLAY, LICENSE OR DISTRIBUTION OF THE SOFTWARE AND AGREE TO + INDEMNIFY AND HOLD HARMLESS BRIGHAM AND ALL CONTRIBUTORS FROM AND + AGAINST ANY AND ALL CLAIMS, SUITS, ACTIONS, DEMANDS AND JUDGMENTS + ARISING THEREFROM. + +6. None of the names, logos or trademarks of Brigham or any of + Brigham's affiliates or any of the Contributors, or any funding + agency, may be used to endorse or promote products produced in + whole or in part by operation of the Software or derived from or + based on the Software without specific prior written permission + from the applicable party. + +7. Any use, reproduction or distribution of the Software which is not + in accordance with this Software License shall automatically revoke + all rights granted to you under this Software License and render + Paragraphs 1 and 2 of this Software License null and void. + +8. This Software License does not grant any rights in or to any + intellectual property owned by Brigham or any Contributor except + those rights expressly granted hereunder. + +PART C. MISCELLANEOUS + +This Agreement shall be governed by and construed in accordance with +the laws of The Commonwealth of Massachusetts without regard to +principles of conflicts of law. This Agreement shall supercede and +replace any license terms that you may have agreed to previously with +respect to Slicer. diff --git a/PETCphantom.png b/PETCphantom.png new file mode 100644 index 0000000..f795463 Binary files /dev/null and b/PETCphantom.png differ diff --git a/PETPhantomAnalysis/CMakeLists.txt b/PETPhantomAnalysis/CMakeLists.txt new file mode 100644 index 0000000..22340b8 --- /dev/null +++ b/PETPhantomAnalysis/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required(VERSION 2.8.6) + +#----------------------------------------------------------------------------- +set(MODULE_NAME PETPhantomAnalysis) + +#----------------------------------------------------------------------------- +set(MODULE_PYTHON_SCRIPTS + PETPhantomAnalysis.py + ) + +set(MODULE_PYTHON_RESOURCES + ) + +#----------------------------------------------------------------------------- +SlicerMacroBuildScriptedModule( + NAME PETPhantomAnalysis + SCRIPTS "${MODULE_PYTHON_SCRIPTS}" + RESOURCES "${MODULE_PYTHON_RESOURCES}" + ) + +#----------------------------------------------------------------------------- +if(BUILD_TESTING) + + # Register the unittest subclass in the main script as a ctest. + # Note that the test will also be available at runtime. + slicer_add_python_unittest(SCRIPT ${MODULE_NAME}.py + SLICER_ARGS --additional-module-paths + ${CMAKE_BINARY_DIR}/${Slicer_QTSCRIPTEDMODULES_LIB_DIR} + ${DEPENDENCIES_ADDITIONAL_MODULE_PATHS}) + + #add_subdirectory(Testing) +endif() + diff --git a/PETPhantomAnalysis/PETPhantomAnalysis.py b/PETPhantomAnalysis/PETPhantomAnalysis.py new file mode 100644 index 0000000..829b004 --- /dev/null +++ b/PETPhantomAnalysis/PETPhantomAnalysis.py @@ -0,0 +1,464 @@ +import os +import unittest +import vtk, qt, ctk, slicer +from slicer.ScriptedLoadableModule import * +import logging +import getpass +from datetime import datetime +import json +import shutil +import tempfile +import dicom + +# +# PETPhantomAnalysis +# + +class PETPhantomAnalysis(ScriptedLoadableModule): + """Uses ScriptedLoadableModule base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def __init__(self, parent): + ScriptedLoadableModule.__init__(self, parent) + self.parent.title = "PET Cylinder Phantom Analysis" + self.parent.categories = ["Quantification"] + self.parent.dependencies = [] + self.parent.contributors = ["Christian Bauer (University of Iowa)"] + self.parent.helpText = """ + Measurement of calibration and uniformity in a cylinder phantom PET scan. \ + Documentation. + """ + self.parent.acknowledgementText = """ + This work was partially funded by NIH grants U01-CA140206 and U24-CA180918. +""" # replace with organization, grant and thanks. + +# +# PETPhantomAnalysisWidget +# + +class PETPhantomAnalysisWidget(ScriptedLoadableModuleWidget): + """Uses ScriptedLoadableModuleWidget base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def __init__(self,parent=None): + ScriptedLoadableModuleWidget.__init__(self, parent) + self.slicerTempDir = slicer.util.tempDirectory() + + def setup(self): + self.measurementsLogic = PETPhantomAnalysisLogic() + + ScriptedLoadableModuleWidget.setup(self) + + # Instantiate and connect widgets ... + + # + # Parameters Area + # + parametersCollapsibleButton = ctk.ctkCollapsibleButton() + parametersCollapsibleButton.text = "PET Cylinder Phantom Analysis" + self.layout.addWidget(parametersCollapsibleButton) + parametersFormLayout = qt.QFormLayout(parametersCollapsibleButton) + + # + # input volume selector + # + self.inputSelector = slicer.qMRMLNodeComboBox() + self.inputSelector.nodeTypes = ["vtkMRMLScalarVolumeNode"] + #self.inputSelector.addAttribute("vtkMRMLScalarVolumeNode", "DICOM.instanceUIDs", None) # add this to restrict to dicom datasets + self.inputSelector.selectNodeUponCreation = True + self.inputSelector.addEnabled = False + self.inputSelector.removeEnabled = False + self.inputSelector.noneEnabled = False + self.inputSelector.showHidden = False + self.inputSelector.showChildNodeTypes = False + self.inputSelector.setMRMLScene( slicer.mrmlScene ) + self.inputSelector.setToolTip( "Input SUVbw normalized DICOM PET volume") + self.inputSelector.connect("currentNodeChanged(bool)",self.refreshUIElements) + self.inputSelector.connect("currentNodeChanged(bool)",self.inputVolumeSelected) + parametersFormLayout.addRow("Input Volume", self.inputSelector) + + self.normalizationFactorLineEdit = qt.QLineEdit("1.0") + self.normalizationFactorValidator = qt.QDoubleValidator() + self.normalizationFactorValidator.bottom = 0.0 + self.normalizationFactorLineEdit.setValidator(self.normalizationFactorValidator) + self.normalizationFactorLineEdit.setToolTip( "Normalization factor for input dataset") + self.normalizationFactorLineEdit.connect("textChanged(QString)",self.refreshUIElements) + parametersFormLayout.addRow("Normalization Factor",self.normalizationFactorLineEdit) + + self.segmentationSelector = slicer.qMRMLNodeComboBox() + self.segmentationSelector.nodeTypes = ["vtkMRMLLabelMapVolumeNode"] + self.segmentationSelector.selectNodeUponCreation = True + self.segmentationSelector.addEnabled = True + self.segmentationSelector.removeEnabled = True + self.segmentationSelector.noneEnabled = False + self.segmentationSelector.showHidden = False + self.segmentationSelector.showChildNodeTypes = False + self.segmentationSelector.setMRMLScene( slicer.mrmlScene ) + self.segmentationSelector.setToolTip( "Output liver reference region volume") + self.segmentationSelector.connect("currentNodeChanged(bool)",self.refreshUIElements) + parametersFormLayout.addRow("Output Volume", self.segmentationSelector) + + self.segmentButton = qt.QPushButton("Analyze cylinder phantom") + self.segmentButton.toolTip = "Analyze cylinder phantom and show results" + self.segmentButton.connect('clicked(bool)', self.onApplyButton) + parametersFormLayout.addRow("Analysis",self.segmentButton) + + # Normalization form + normalizationCollapsibleButton = ctk.ctkCollapsibleButton() + normalizationCollapsibleButton.text = "Normalization Metadata" + self.layout.addWidget(normalizationCollapsibleButton) + normalizationFormLayout = qt.QFormLayout(normalizationCollapsibleButton) + + self.inputTypeLabel = qt.QLabel('unspecified') + normalizationFormLayout.addRow("Input Volume Type",self.inputTypeLabel) + + self.halfLifeLineEdit = qt.QLineEdit(float('inf')) + self.halfLifeLineEditValidator = qt.QDoubleValidator() + self.halfLifeLineEditValidator.bottom = 0.0 + self.halfLifeLineEdit.setValidator(self.halfLifeLineEditValidator) + self.halfLifeLineEdit.setToolTip( "Radionuclide half life time") + self.halfLifeLineEdit.connect("textChanged(QString)",self.updateNormalizationFactor) + normalizationFormLayout.addRow("Radionuclide Half Life (s)",self.halfLifeLineEdit) + + self.activityLineEdit = qt.QLineEdit("1.0") + self.activityLineEditValidator = qt.QDoubleValidator() + self.activityLineEditValidator.bottom = 0.0 + self.activityLineEdit.setValidator(self.activityLineEditValidator) + self.activityLineEdit.setToolTip( "Normalization factor for input dataset") + self.activityLineEdit.connect("textChanged(QString)",self.updateNormalizationFactor) + normalizationFormLayout.addRow(u"Syringe Activity (\u00B5Ci)",self.activityLineEdit) + + self.activityTimeEdit = qt.QTimeEdit() + self.activityTimeEdit.connect("timeChanged(QTime)",self.updateNormalizationFactor) + normalizationFormLayout.addRow("Syringe Assay Time",self.activityTimeEdit) + + self.residualActivityLineEdit = qt.QLineEdit("0.0") + self.residualActivityLineEditValidator = qt.QDoubleValidator() + self.residualActivityLineEditValidator.bottom = 0.0 + self.residualActivityLineEdit.setValidator(self.residualActivityLineEditValidator) + self.residualActivityLineEdit.setToolTip( "Normalization factor for input dataset") + self.residualActivityLineEdit.connect("textChanged(QString)",self.updateNormalizationFactor) + normalizationFormLayout.addRow(u"Syringe Residual Activity (\u00B5Ci)",self.residualActivityLineEdit) + + self.residualactivityTimeEdit = qt.QTimeEdit() + self.residualactivityTimeEdit.connect("timeChanged(QTime)",self.updateNormalizationFactor) + normalizationFormLayout.addRow("Syringe Residual Activity Time",self.residualactivityTimeEdit) + + self.scanTimeEdit = qt.QTimeEdit() + self.scanTimeEdit.connect("timeChanged(QTime)",self.updateNormalizationFactor) + normalizationFormLayout.addRow("Scan Time",self.scanTimeEdit) + + self.weightLineEdit = qt.QLineEdit("1000.0") + self.weightLineEditValidator = qt.QDoubleValidator() + self.weightLineEditValidator.bottom = 0.0 + self.weightLineEdit.setValidator(self.weightLineEditValidator) + self.weightLineEdit.setToolTip( "Normalization factor for input dataset") + self.weightLineEdit.connect("textChanged(QString)",self.updateNormalizationFactor) + normalizationFormLayout.addRow("Fill Weight (g)",self.weightLineEdit) + + self.inputVolumeSUVNormalizationLabel = qt.QLabel("1.0") + self.inputVolumeSUVNormalizationLabel.setToolTip( "Normalization factor applied to input volume") + normalizationFormLayout.addRow("Input Volume Normalization Factor",self.inputVolumeSUVNormalizationLabel) + + # Measurement results + measurementsCollapsibleButton = ctk.ctkCollapsibleButton() + measurementsCollapsibleButton.text = "Measurements" + self.layout.addWidget(measurementsCollapsibleButton) + parametersFormLayout = qt.QFormLayout(measurementsCollapsibleButton) + + self.meanValueLineEdit = qt.QLineEdit("-") + self.meanValueLineEdit.setReadOnly(True) + self.meanValueLineEdit.setToolTip( "Mean value in measurement region") + parametersFormLayout.addRow("Mean",self.meanValueLineEdit) + self.stdValueLineEdit = qt.QLineEdit("-") + self.stdValueLineEdit.setReadOnly(True) + self.stdValueLineEdit.setToolTip( "Standard deviation in measurement region") + parametersFormLayout.addRow("Standard Deviation",self.stdValueLineEdit) + self.maxRelDiffValueLineEdit = qt.QLineEdit("-") + self.maxRelDiffValueLineEdit.setReadOnly(True) + self.maxRelDiffValueLineEdit.setToolTip( "Maximum relative difference in axial direction in reference region") + parametersFormLayout.addRow("Maximum Relative Difference",self.maxRelDiffValueLineEdit) + + # Add vertical spacer + self.layout.addStretch(1) + + self.inputVolumeSelected() + self.refreshUIElements() + + def enter(self): + pass + + def exit(self): + pass + + def inputVolumeSelected(self): + self.inputTypeLabel.text = 'unspecified' + self.halfLifeLineEdit.text = float('inf') + self.activityLineEdit.text = 1 + self.activityTimeEdit.time = qt.QTime(0,0) + self.residualActivityLineEdit.text = 0 + self.residualactivityTimeEdit.time = qt.QTime(0,0) + self.scanTimeEdit.time = qt.QTime(0,0) + self.weightLineEdit.text = 1000.0 + self.inputVolumeSUVNormalizationLabel.text = 1.0 + + inputVolume = self.inputSelector.currentNode() + if not inputVolume: return + dicomInstanceUIDs = inputVolume.GetAttribute('DICOM.instanceUIDs') + rwvmInstanceUID = inputVolume.GetAttribute('DICOM.RWV.instanceUID') + if dicomInstanceUIDs: + self.inputTypeLabel.text = "DICOM volume" + sourceFileName = slicer.dicomDatabase.fileForInstance(dicomInstanceUIDs.split(" ")[0]) + d = dicom.read_file(sourceFileName) + if d.Modality=='PT': + self.inputTypeLabel.text = "PET DICOM volume" + self.halfLifeLineEdit.text = d.RadiopharmaceuticalInformationSequence[0].RadionuclideHalfLife + self.activityLineEdit.text = self._doseInMicroCi(d.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose, d.Units) + self.activityTimeEdit.time = self._getTime(d.RadiopharmaceuticalInformationSequence[0].RadiopharmaceuticalStartTime) + self.residualactivityTimeEdit.time = self.activityTimeEdit.time + self.scanTimeEdit.time = self._getTime(d.SeriesTime) + self.weightLineEdit.text = d.PatientWeight*1000.0 + if inputVolume.GetVoxelValueUnits() and inputVolume.GetVoxelValueUnits().GetCodeValue()=='{SUVbw}g/ml': # calculate SUV normalization factor based on PET DICOM + self.inputTypeLabel.text = "SUVbw normalized PET DICOM volume" + halfLife = d.RadiopharmaceuticalInformationSequence[0].RadionuclideHalfLife + injectedDose = self._doseInMicroCi(d.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose, d.Units) + weight = d.PatientWeight + decayTime = 0.001*(\ + self._getTime(d.SeriesTime).msecsSinceStartOfDay()-\ + self._getTime(d.RadiopharmaceuticalInformationSequence[0].RadiopharmaceuticalStartTime).msecsSinceStartOfDay()) + decayedDose = injectedDose * pow(2.0, -decayTime/halfLife) + Ci2BqFactor = 37000000000.0 + decayedDose = decayedDose * Ci2BqFactor * 1e-9 # convert to kBq/mL + SUVNormalizationFactor = weight / decayedDose + self.inputVolumeSUVNormalizationLabel.text = SUVNormalizationFactor + elif rwvmInstanceUID: # load SUV normalization factor from RWVM file + self.inputTypeLabel.text = "SUV normalized PET DICOM volume" + sourceFileName = slicer.dicomDatabase.fileForInstance(rwvmInstanceUID) + d = dicom.read_file(sourceFileName) + self.inputVolumeSUVNormalizationLabel.text = d.ReferencedImageRealWorldValueMappingSequence[0].RealWorldValueMappingSequence[0].RealWorldValueSlope + self.updateNormalizationFactor() + + def updateNormalizationFactor(self): + if self.inputTypeLabel.text == 'unspecified': + return + try: + halfLife = float(self.halfLifeLineEdit.text) + seriesTime = float(self.scanTimeEdit.time.msecsSinceStartOfDay())/1000.0 + activity = float(self.activityLineEdit.text) + assayTime = float(self.activityTimeEdit.time.msecsSinceStartOfDay())/1000.0 + residualActivity = float(self.residualActivityLineEdit.text) + residualActivityTime = float(self.residualactivityTimeEdit.time.msecsSinceStartOfDay())/1000.0 + weight = float(self.weightLineEdit.text)/1000.0 + t1 = residualActivityTime-assayTime + injectedDose = activity * pow(2.0, -t1/halfLife) - residualActivity + t2 = seriesTime - residualActivityTime + decayedDose = injectedDose * pow(2.0, -t2/halfLife) + Ci2BqFactor = 37000000000.0 + decayedDose = decayedDose * Ci2BqFactor * 1e-9 # convert to kBq/mL + SUVNormalizationFactor = weight / decayedDose + normalizationFactor = SUVNormalizationFactor/float(self.inputVolumeSUVNormalizationLabel.text) + self.normalizationFactorLineEdit.text = normalizationFactor + except: + pass + + def _getTime(self, time): + dot = time.find('.') if time.find('.')>=0 else len(time) + hours = int(time[:dot-4]) + minutes = int(time[dot-4:dot-2]) + seconds = int(round(float(time[dot-2:dot]))) + mseconds = int(round(float(time[dot:-1])*1000.0)) if dot2 and units[0] in magnitudeFactors: + dose = dose*magnitudeFactors[units[0]] + return dose*1e6 + + def refreshUIElements(self): + self.meanValueLineEdit.text = "-" + self.stdValueLineEdit.text = "-" + self.maxRelDiffValueLineEdit.text = "-" + self.segmentButton.enabled = self.inputSelector.currentNode()!=None + + def onApplyButton(self): + inputVolume = self.inputSelector.currentNode() + outputVolume = self.segmentationSelector.currentNode() + jsonFile = tempfile.NamedTemporaryFile() + normalizationFactor = 1.0 + try: + normalizationFactor = float(self.normalizationFactorLineEdit.text) + except: + pass + + cliParams = {'inputVolume': inputVolume.GetID(), \ + 'normalizationFactor': normalizationFactor, \ + 'measurementsData': jsonFile.name, + } + if outputVolume: cliParams['outputVolume'] = outputVolume.GetID() + + pd = qt.QProgressDialog('Running PET Cylinder Phantom Analysis...', 'Cancel', 0, 100, slicer.util.mainWindow()) + pd.setModal(True) + pd.setMinimumDuration(0) + pd.show() + pd.setValue(30) + cliNode = None + cliNode = slicer.cli.run(slicer.modules.petphantomanalysiscli, cliNode, cliParams, wait_for_completion=False) + while cliNode.IsBusy(): + slicer.app.processEvents() + if pd.wasCanceled: + cliNode.Cancel() + pd.setValue(100) + if pd.wasCanceled: + return + + if cliNode.GetStatusString() != 'Completed': + qt.QMessageBox().warning(None,"Warning","Analysis of PET cylinder phantom was not successful.") + return + + for i in range(cliNode.GetNumberOfParametersInGroup(1)): + name = cliNode.GetParameterName(1,i) + value = cliNode.GetParameterDefault(1,i) + if name=='Mean_s': self.meanValueLineEdit.setText(value) + if name=='Std_s': self.stdValueLineEdit.setText(value) + if name=='MaxRelDiff_s': self.maxRelDiffValueLineEdit.setText(value) + + if self.meanValueLineEdit.text=='--': + qt.QMessageBox().warning(None,"Warning","Analysis of PET cylinder phantom was not successful.") + return + + try: + measurements = json.load(open(jsonFile.name)) + except: + qt.QMessageBox().warning(None,"Warning","Analysis of PET cylinder phantom was not successful.") + return + + # plot normalized slice measurements + mean = measurements['CylinderMean'] + sliceOffsets = measurements['SliceOffsets'] + sliceMeasurements = measurements['SliceMeasurements'] + normalizedSliceMeasurements = [v/mean for v in sliceMeasurements] + self._createPlot(normalizedSliceMeasurements, sliceOffsets) + + # visualize matched cylinder + cylinderModelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode") + cylinderModelNode.CreateDefaultDisplayNodes() + cylinderModelNode.GetDisplayNode().SetColor(0,1,0) + cylinderModelNode.GetDisplayNode().SliceIntersectionVisibilityOn() + cylinderModelNode.SetAndObserveMesh(self._createCylinderMesh( + measurements['CylinderCenter'], measurements['CylinderDirection'])) + + def _createPlot(self, normalizedSliceMeasurements, sliceOffsets): + tableNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode") + table = tableNode.GetTable() + + arrX = vtk.vtkFloatArray() + arrX.SetName("offset") + table.AddColumn(arrX) + + arrY = vtk.vtkFloatArray() + arrY.SetName("intensity") + table.AddColumn(arrY) + + numPoints = len(normalizedSliceMeasurements) + table.SetNumberOfRows(numPoints) + for i in range(numPoints): + table.SetValue(i, 0, sliceOffsets[i] ) + table.SetValue(i, 1, normalizedSliceMeasurements[i]) + + plotSeriesNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLPlotSeriesNode", "Axial uniformity") + plotSeriesNode.SetAndObserveTableNodeID(tableNode.GetID()) + plotSeriesNode.SetXColumnName("offset") + plotSeriesNode.SetYColumnName("intensity") + plotSeriesNode.SetPlotType(slicer.vtkMRMLPlotSeriesNode.PlotTypeScatter) + plotSeriesNode.SetUniqueColor() + + # Create plot chart node + plotChartNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLPlotChartNode") + plotChartNode.AddAndObservePlotSeriesNodeID(plotSeriesNode.GetID()) + plotChartNode.SetTitle('Mean normalized intensity profile') + plotChartNode.SetXAxisTitle('Offset from center in axial direction (mm)') + plotChartNode.SetYAxisTitle('Mean normalized intensity') + + # Switch to a layout that contains a plot view to create a plot widget + layoutManager = slicer.app.layoutManager() + layoutWithPlot = slicer.modules.plots.logic().GetLayoutWithPlot(layoutManager.layout) + layoutManager.setLayout(layoutWithPlot) + + # Select chart in plot view + plotWidget = layoutManager.plotWidget(0) + plotViewNode = plotWidget.mrmlPlotViewNode() + plotViewNode.SetPlotChartNodeID(plotChartNode.GetID()) + + def _createCylinderMesh(self, center, direction, radius=100.0, height=200.0): + startPoint = [center[i]-direction[i]*height/2.0 for i in range(3)] + cylinderSource = vtk.vtkCylinderSource() + cylinderSource.SetResolution(365) + cylinderSource.SetRadius(radius) + cylinderSource.Update() + normalizedX = direction + normalizedY = [0,0,0] + normalizedZ = [0,0,0] + vtk.vtkMath.Normalize(normalizedX) + arbitrary = [1,1,1] + vtk.vtkMath.Cross(normalizedX, arbitrary, normalizedZ) + vtk.vtkMath.Normalize(normalizedZ) + vtk.vtkMath.Cross(normalizedZ, normalizedX, normalizedY) + matrix = vtk.vtkMatrix4x4() + matrix.Identity() + for i in range(3): + matrix.SetElement(i, 0, normalizedX[i]) + matrix.SetElement(i, 1, normalizedY[i]) + matrix.SetElement(i, 2, normalizedZ[i]) + transform = vtk.vtkTransform() + transform.Translate(startPoint) + transform.Concatenate(matrix) + transform.RotateZ(-90.0) + transform.Scale(1.0, height, 1.0) + transform.Translate(0, .5, 0) + transformPD = vtk.vtkTransformPolyDataFilter() + transformPD.SetTransform(transform) + transformPD.SetInputConnection(cylinderSource.GetOutputPort()) + transformPD.Update() + return transformPD.GetOutputDataObject(0) + + +# +# PETPhantomAnalysisLogic +# + +class PETPhantomAnalysisLogic(ScriptedLoadableModuleLogic): + + def __init__(self): + pass + + def __del__(self): + pass + +# +# PETPhantomAnalysisTest +# + +from DICOMLib import DICOMUtils +class PETPhantomAnalysisTest(ScriptedLoadableModuleTest): + """ + This is the test case for your scripted module. + Uses ScriptedLoadableModuleTest base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def runTest(self): + """Run as few or as many tests as needed here. + """ + self.test_PETPhantomAnalysis() + + def test_PETPhantomAnalysis(self): + """ test standard processing + """ + self.delayDisplay('Test passed!') + diff --git a/PETPhantomAnalysisCLI/CMakeLists.txt b/PETPhantomAnalysisCLI/CMakeLists.txt new file mode 100644 index 0000000..db45538 --- /dev/null +++ b/PETPhantomAnalysisCLI/CMakeLists.txt @@ -0,0 +1,31 @@ + +#----------------------------------------------------------------------------- +set(MODULE_NAME PETPhantomAnalysisCLI) + +#----------------------------------------------------------------------------- +set(MODULE_INCLUDE_DIRECTORIES + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${RapidJSON_INCLUDE_DIR} + ) + +set(MODULE_SRCS + ) + +set(MODULE_TARGET_LIBRARIES + ${ITK_LIBRARIES} + ${RapidJSON_LIBRARY} + ) + +#----------------------------------------------------------------------------- +SEMMacroBuildCLI( + NAME ${MODULE_NAME} + TARGET_LIBRARIES ${MODULE_TARGET_LIBRARIES} + INCLUDE_DIRECTORIES ${MODULE_INCLUDE_DIRECTORIES} + ADDITIONAL_SRCS ${MODULE_SRCS} + ) + +#----------------------------------------------------------------------------- +if(BUILD_TESTING) + #add_subdirectory(Testing) +endif() + diff --git a/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.cxx b/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.cxx new file mode 100644 index 0000000..a29af80 --- /dev/null +++ b/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.cxx @@ -0,0 +1,192 @@ +/*============================================================================== + + Program: PETPhantomAnlysisCLI + + (c) Copyright University of Iowa All Rights Reserved. + + See COPYRIGHT.txt + or http://www.slicer.org/copyright/copyright.txt for details. + + 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. + + ==============================================================================*/ + +#include "itkImageFileWriter.h" + +#include "itkPluginUtilities.h" + +#include "PETPhantomAnalysisCLICLP.h" + +#include "itkCylinderMatchingImageFilter.h" +#include "itkCylinderUniformityMeasurementImageFilter.h" +#include "itkShiftScaleImageFilter.h" +#include + +// rapidjson +#include "rapidjson/rapidjson.h" +#include "rapidjson/document.h" +#include "rapidjson/prettywriter.h" +#include "rapidjson/filewritestream.h" +#include "rapidjson/ostreamwrapper.h" + +// Use an anonymous namespace to keep class types and function names +// from colliding when module is used as shared object module. Every +// thing should be in an anonymous namespace except for the module +// entry point, e.g. main() +// +namespace +{ + +template +itk::FixedArray lps2ras(const itk::FixedArray& lps) +{ + itk::FixedArray ras; + ras[0] = -lps[0]; + ras[1] = -lps[1]; + ras[2] = lps[2]; + return ras; +} + +template +void jsonAddMember(rapidjson::Value& element, rapidjson::Document::AllocatorType& allocator, + const std::string& name, const itk::FixedArray& values) +{ + using namespace rapidjson; + Value r_values(kArrayType); + for (size_t i=0; i +void jsonAddMember(rapidjson::Value& element, rapidjson::Document::AllocatorType& allocator, + const std::string& name, const std::vector& values) +{ + using namespace rapidjson; + Value r_values(kArrayType); + for (std::vector::const_iterator it=values.begin(); it!=values.end(); ++it) + r_values.PushBack(*it, allocator); + Value elementName(name.c_str(), allocator); + element.AddMember(elementName.Move(), r_values.Move(), allocator); +} + +} // end of anonymous namespace + + +int main( int argc, char * argv[] ) +{ + PARSE_ARGS; + + try + { + typedef double PixelType; + typedef short LabelPixelType; + typedef itk::Image ImageType; + typedef itk::Image LabelImageType; + + // read PET scan + typedef itk::ImageFileReader ReaderType; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName( inputVolume.c_str() ); + reader->Update(); + ImageType::Pointer petScan = reader->GetOutput(); + + // normalize input data, if necessary + if (normalizationFactor!=1.0) + { + typedef itk::ShiftScaleImageFilter NormalizerFilterType; + NormalizerFilterType::Pointer normalizer = NormalizerFilterType::New(); + normalizer->SetShift(0.0); + normalizer->SetScale(normalizationFactor); + normalizer->SetInput(petScan); + normalizer->Update(); + petScan = normalizer->GetOutput(); + } + + // find cylinder center and orientation + typedef itk::CylinderMatchingImageFilter CylinderMatchingType; + CylinderMatchingType::Pointer cylinderMatching = CylinderMatchingType::New(); + cylinderMatching->SetInput(petScan); + cylinderMatching->SetSmoothingSigma(12.0); + cylinderMatching->Update(); + CylinderMatchingType::PointType cylinderCenter = cylinderMatching->GetCenter(); + CylinderMatchingType::VectorType cylinderDirection = cylinderMatching->GetDirection(); + if (cylinderDirection[2]<0) cylinderDirection *= -1.0; + + // measure calibration and uniformity + typedef itk::CylinderUniformityMeasurementImageFilter UniformityMeasurementFilterType; + UniformityMeasurementFilterType::Pointer uniformityMeasurements = + UniformityMeasurementFilterType::New(); + uniformityMeasurements->SetInput(petScan); + uniformityMeasurements->SetCenter(cylinderCenter); + uniformityMeasurements->SetDirection(cylinderDirection); + uniformityMeasurements->SetRadius(73.0); // todo: or make configurable? + uniformityMeasurements->SetHeight(160.0); // todo: or make configurable? + uniformityMeasurements->SetLabelInsideRadiusLimit(0); + uniformityMeasurements->SetLabelInsideHeightLimit(0); + uniformityMeasurements->Update(); + LabelImageType::Pointer measurementRegion = uniformityMeasurements->GetOutput(); + + // write measurements + std::ofstream writeFile; + writeFile.open( returnParameterFile.c_str() ); + writeFile << "Mean_s = " << uniformityMeasurements->GetCylinderMean() << std::endl; + writeFile << "Std_s = " << uniformityMeasurements->GetCylinderStd() << std::endl; + writeFile << "MaxRelDiff_s = " << uniformityMeasurements->GetMaxRelativeDifference() << std::endl; + writeFile.close(); + + // write measurement region volume, if requested + if (outputVolume.size()>0) + { + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName( outputVolume.c_str() ); + writer->SetInput( measurementRegion ); + writer->SetUseCompression(1); + writer->Update(); + } + + // write JSON report, if requested + if (measurementsData.size()>0) + { + // prepare JSON document + using namespace rapidjson; + Document document; + document.SetObject(); + Document::AllocatorType& allocator = document.GetAllocator(); + + // add information + document.AddMember("NormalizationFactor", normalizationFactor, allocator); + jsonAddMember(document, allocator, "CylinderCenter", lps2ras(uniformityMeasurements->GetCenter())); + jsonAddMember(document, allocator, "CylinderDirection", lps2ras(uniformityMeasurements->GetDirection())); + document.AddMember("CylinderRadius", uniformityMeasurements->GetRadius(), allocator); + document.AddMember("CylinderHeight", uniformityMeasurements->GetHeight(), allocator); + document.AddMember("CylinderMean", uniformityMeasurements->GetCylinderMean(), allocator); + document.AddMember("CylinderStd", uniformityMeasurements->GetCylinderStd(), allocator); + document.AddMember("MaxRelDiff", uniformityMeasurements->GetMaxRelativeDifference(), allocator); + jsonAddMember(document, allocator, "SliceMeasurements", + uniformityMeasurements->GetSliceMeasurements()->CastToSTLConstContainer()); + jsonAddMember(document, allocator, "SliceOffsets", + uniformityMeasurements->GetSliceOffsets()->CastToSTLConstContainer()); + + // write JSON document to file + std::ofstream outputFile(measurementsData.c_str()); + OStreamWrapper osw(outputFile); + PrettyWriter writer(osw); + document.Accept(writer); + } + + } + catch( itk::ExceptionObject & excep ) + { + std::cerr << argv[0] << ": exception caught!" << std::endl; + std::cerr << excep << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.xml b/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.xml new file mode 100644 index 0000000..c40ec28 --- /dev/null +++ b/PETPhantomAnalysisCLI/PETPhantomAnalysisCLI.xml @@ -0,0 +1,65 @@ + + + Quantification + PET Cylinder Phantom Analysis CLI + + 0.0.1 + https://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modules/PETCPhantom + Slicer + Christian Bauer (University of Iowa), Ethan Ulrich (University of Iowa) + This work was partially funded by NIH grants U01-CA140206 and U24-CA180918. + + + + + inputVolume + + input + 0 + + + + normalizationFactor + --normalizationFactor + + + 1.0 + + + outputVolume + + --outputVolume + output + + + + measurementsData + + output + outputMetadata + File name of the JSON file that will keep the metadata and measurements information. + + + + + + + Mean_s + + output + + + + Std_s + + output + + + + MaxRelDiff_s + + output + + + + diff --git a/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.cxx b/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.cxx new file mode 100644 index 0000000..9ba4d0b --- /dev/null +++ b/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.cxx @@ -0,0 +1,112 @@ +#ifndef __itkCylinderMatchingImageFilter_cxx +#define __itkCylinderMatchingImageFilter_cxx + +#include "itkCylinderMatchingImageFilter.h" + +#include "itkOtsuThresholdImageFilter.h" +#include "itkConnectedComponentImageFilter.h" +#include "itkRelabelComponentImageFilter.h" +#include "itkSmoothingRecursiveGaussianImageFilter.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "vnl/algo/vnl_symmetric_eigensystem.h" + +namespace itk +{ +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +CylinderMatchingImageFilter< TInputImage, TOutputImage > +::CylinderMatchingImageFilter() : + m_SmoothingSigma(1.0), m_Volume(0.0) +{ + m_Center.Fill(0); + m_Direction.Fill(0); +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +void CylinderMatchingImageFilter< TInputImage, TOutputImage > +::GenerateData() +{ + //---------------------------------------------------------------------------- + // processing steps: smoothing, find and apply Otsu threshold, label objects, sort + // them by size, use calculate centroid and eigenvector of largest component + // as cylinder parameters + + typename InputImageType::ConstPointer image = this->GetInput(); + + typedef itk::SmoothingRecursiveGaussianImageFilter SmoothingImageFilterType; + typename SmoothingImageFilterType::Pointer smoothingFilter = SmoothingImageFilterType::New(); + smoothingFilter->SetSigma(m_SmoothingSigma); + smoothingFilter->SetInput(image); + + typedef itk::OtsuThresholdImageFilter ThresholdImageFilterType; + typename ThresholdImageFilterType::Pointer otsuFilter = ThresholdImageFilterType::New(); + otsuFilter->SetInput(smoothingFilter->GetOutput()); + otsuFilter->SetInsideValue(0); + otsuFilter->SetOutsideValue(1); + otsuFilter->Update(); + //double otsuThreshold = otsuFilter->GetThreshold(); // just for debugging + + typedef itk::ConnectedComponentImageFilter LabelingImageFilterType; + typename LabelingImageFilterType::Pointer labeingFilter = LabelingImageFilterType::New(); + labeingFilter->SetInput(otsuFilter->GetOutput()); + + typedef itk::RelabelComponentImageFilter RelabelingImageFilterType; + typename RelabelingImageFilterType::Pointer relablingFilter = RelabelingImageFilterType::New(); + relablingFilter->SetInput(labeingFilter->GetOutput()); + relablingFilter->Update(); + + typename OutputImageType::Pointer segmentation = relablingFilter->GetOutput(); + + PointType tmpPoint; + typedef std::vector< PointType > PointList; + PointList segPoints; + ImageRegionConstIteratorWithIndex it(segmentation, + segmentation->GetLargestPossibleRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + if (it.Value()==1) + { + segmentation->TransformIndexToPhysicalPoint(it.GetIndex(), tmpPoint); + segPoints.push_back(tmpPoint); + } + } + VectorType centerVector; centerVector.Fill(0); + for (typename PointList::const_iterator it=segPoints.begin(); it!=segPoints.end(); ++it) + centerVector += it->GetVectorFromOrigin(); + m_Center = (centerVector/double(segPoints.size())); + + vnl_matrix< double > normalizedSecondOrderCentralMoments(3, 3, 0.0); + for (typename PointList::const_iterator it=segPoints.begin(); it!=segPoints.end(); ++it) + { + vnl_vector v = (*it - m_Center).GetVnlVector(); + normalizedSecondOrderCentralMoments += outer_product(v,v); + } + vnl_symmetric_eigensystem< double > eig(normalizedSecondOrderCentralMoments); + + m_Direction[0] = eig.get_eigenvector(2)[0]; + m_Direction[1] = eig.get_eigenvector(2)[1]; + m_Direction[2] = eig.get_eigenvector(2)[2]; + + // get volume (for plausibility check) + typename InputImageType::SpacingType spacing = image->GetSpacing(); + double voxelVolume = spacing[0]*spacing[1]*spacing[2]; + m_Volume = double(segPoints.size())*voxelVolume; + + this->GetOutput()->Graft(relablingFilter->GetOutput()); +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +void CylinderMatchingImageFilter< TInputImage, TOutputImage > +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); + os << indent << "Center: " << this->m_Center << std::endl; + os << indent << "Direction: " << this->m_Direction << std::endl; + os << indent << "Volume: " << this->m_Volume << std::endl; +} + +}// end namespace + +#endif diff --git a/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.h b/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.h new file mode 100644 index 0000000..f64fc08 --- /dev/null +++ b/PETPhantomAnalysisCLI/include/itkCylinderMatchingImageFilter.h @@ -0,0 +1,74 @@ +#ifndef __itkCylinderMatchingImageFilter_h +#define __itkCylinderMatchingImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkVectorContainer.h" +#include "itkVector.h" + +namespace itk +{ + +/**\class CylinderMatchingImageFilter + * \brief short description + * More detailed desciption + */ +template< class TInputImage, class TOutputImage=TInputImage > +class CylinderMatchingImageFilter : public ImageToImageFilter< TInputImage, TOutputImage > +{ +public: + /** Standard class typedefs. */ + typedef CylinderMatchingImageFilter Self; + typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass; + typedef SmartPointer< Self > Pointer; + + /** Useful class typedefs*/ + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + typedef Image LabelImageType; + typedef typename TInputImage::PixelType PixelType; + typedef typename TInputImage::PointType PointType; + typedef typename TInputImage::IndexType IndexType; + typedef Vector VectorType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(CylinderMatchingImageFilter, ImageToImageFilter); + + /** Dimension of the underlying image. */ + itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension); + + /** Set/Get Macros */ + itkSetMacro(SmoothingSigma, double); + itkGetMacro(SmoothingSigma, double); + itkGetMacro(Center, PointType); + itkGetMacro(Direction, VectorType); + itkGetMacro(Volume, double); + +protected: + CylinderMatchingImageFilter(); + ~CylinderMatchingImageFilter(){} + + virtual void PrintSelf(std::ostream& os, Indent indent) const; + + /** Does the real work. */ + virtual void GenerateData(); + +private: + CylinderMatchingImageFilter(const Self &); //purposely not implemented + void operator=(const Self &); //purposely not implemented + + double m_SmoothingSigma; + PointType m_Center; + VectorType m_Direction; + double m_Volume; +}; +} //namespace ITK + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkCylinderMatchingImageFilter.cxx" +#endif + +#endif // __itkCylinderMatchingImageFilter_h diff --git a/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.cxx b/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.cxx new file mode 100644 index 0000000..42a88bb --- /dev/null +++ b/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.cxx @@ -0,0 +1,169 @@ +#ifndef __itkCylinderUniformityMeasurementImageFilter_cxx +#define __itkCylinderUniformityMeasurementImageFilter_cxx + +#include "itkCylinderUniformityMeasurementImageFilter.h" + +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageRegionConstIteratorWithIndex.h" + +#include + +namespace itk +{ +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::CylinderUniformityMeasurementImageFilter() : + m_Radius(0.0), m_Height(0.0), m_CylinderMean(0.0), m_CylinderStd(0.0), + m_MaxRelativeDifference(0.0), m_LabelInside(1), m_LabelInsideRadiusLimit(2), + m_LabelInsideHeightLimit(3), m_LabelOutside(0) +{ + m_Center.Fill(0); + m_Direction.Fill(0); + m_SliceMeasurements = MeasurementVectorType::New(); + m_SliceOffsets = MeasurementVectorType::New(); +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +void CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::GenerateData() +{ + typename InputImageType::ConstPointer image = this->GetInput(); + size_t nSlices = image->GetLargestPossibleRegion().GetSize()[2]; + + // initialize output measurements + m_SliceMeasurements->Initialize(); + m_SliceMeasurements->Reserve(nSlices); + m_SliceOffsets->Initialize(); + m_SliceOffsets->Reserve(nSlices); + + // initialize output image + typename OutputImageType::Pointer measurementRegion = this->GetOutput(); + measurementRegion->CopyInformation(image); + measurementRegion->SetRegions(image->GetLargestPossibleRegion()); + measurementRegion->Allocate(); + measurementRegion->FillBuffer(0); + + // temporary variables + std::vector cylinderVoxels; + std::vector sliceVoxelsInRegion(nSlices, 0.0); + + // collect measurements and create output measurement region + itk::ImageRegionConstIteratorWithIndex imageIt(image, image->GetLargestPossibleRegion()); + itk::ImageRegionIteratorWithIndex measurementRegionIt(measurementRegion, measurementRegion->GetLargestPossibleRegion()); + imageIt.GoToBegin(); + measurementRegionIt.GoToBegin(); + PointType p; + while (!imageIt.IsAtEnd()) + { + image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), p); + double distanceFromAxis = this->DistanceFromAxis(p); + double distanceAlongAxis = this->DistanceAlongAxis(p); + double isInside = this->IsInside(p); + if (distanceFromAxisSetElement(z, + m_SliceMeasurements->GetElement(z)+imageIt.Value()); + sliceVoxelsInRegion[z]+=1; + } + if (isInside) + { + cylinderVoxels.push_back(imageIt.Value()); + } + + OutputPixelType regionLabel = m_LabelOutside; + if (isInside) + regionLabel = m_LabelInside; + else if (distanceFromAxis0) + m_SliceMeasurements->SetElement(z, + m_SliceMeasurements->GetElement(z)/double(sliceVoxelsInRegion[z])); + + // calculate mean and std + m_CylinderMean = std::accumulate(cylinderVoxels.begin(), + cylinderVoxels.end(), 0.0)/double(cylinderVoxels.size()); + m_CylinderStd = 0.0; + for (std::vector::const_iterator it=cylinderVoxels.begin(); + it!=cylinderVoxels.end(); ++it) + m_CylinderStd += pow(*it-m_CylinderMean,2.0); + m_CylinderStd = sqrt(m_CylinderStd/double(cylinderVoxels.size())); + + // calculate slice offsets + IndexType idx; idx.Fill(0); + for (size_t z=0; zTransformIndexToPhysicalPoint(idx, p); + m_SliceOffsets->SetElement(z, p[2]-m_Center[2]); + } + + // calculate maximum relative difference + m_MaxRelativeDifference = 0.0; + for (size_t z=0; zGetElement(z); + if (fabs(oz)<(m_Direction[2]*m_Height/2.0)) + { + double relativeDifference = + (m_SliceMeasurements->GetElement(z)-m_CylinderMean)/m_CylinderMean; + if (fabs(relativeDifference)>fabs(m_MaxRelativeDifference)) + m_MaxRelativeDifference = relativeDifference; + } + } +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +bool CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::IsInside(PointType p) const +{ + return this->DistanceFromAxis(p)DistanceAlongAxis(p)) +double CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::DistanceFromAxis(PointType p) const +{ + double dCenter = (p-m_Center).GetNorm(); + double dAlongAxis = DistanceAlongAxis(p); + return std::sqrt(dCenter*dCenter-dAlongAxis*dAlongAxis); +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +double CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::DistanceAlongAxis(PointType p) const +{ + return (p-m_Center)*m_Direction; +} + +//---------------------------------------------------------------------------- +template< class TInputImage, class TOutputImage > +void CylinderUniformityMeasurementImageFilter< TInputImage, TOutputImage > +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); + os << indent << "Center: " << this->m_Center << std::endl; + os << indent << "Direction: " << this->m_Direction << std::endl; + os << indent << "Radius: " << this->m_Radius << std::endl; + os << indent << "Height: " << this->m_Height << std::endl; +} + +}// end namespace + +#endif diff --git a/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.h b/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.h new file mode 100644 index 0000000..522c649 --- /dev/null +++ b/PETPhantomAnalysisCLI/include/itkCylinderUniformityMeasurementImageFilter.h @@ -0,0 +1,105 @@ +#ifndef __itkCylinderUniformityMeasurementImageFilter_h +#define __itkCylinderUniformityMeasurementImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkVectorContainer.h" +#include "itkVector.h" + +namespace itk +{ + +/**\class CylinderUniformityMeasurementImageFilter + * \brief short description + * More detailed desciption + */ +template< class TInputImage, class TOutputImage=TInputImage > +class CylinderUniformityMeasurementImageFilter : public ImageToImageFilter< TInputImage, TOutputImage > +{ +public: + /** Standard class typedefs. */ + typedef CylinderUniformityMeasurementImageFilter Self; + typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass; + typedef SmartPointer< Self > Pointer; + + /** Useful class typedefs*/ + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + typedef typename TInputImage::PixelType PixelType; + typedef typename TOutputImage::PixelType OutputPixelType; + typedef typename TInputImage::PointType PointType; + typedef typename TInputImage::IndexType IndexType; + typedef Vector VectorType; + typedef VectorContainer MeasurementVectorType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(CylinderUniformityMeasurementImageFilter, ImageToImageFilter); + + /** Dimension of the underlying image. */ + itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension); + + /** Set/Get Macros */ + itkSetMacro(Center, PointType); + itkGetMacro(Center, PointType); + itkSetMacro(Direction, VectorType); + itkGetMacro(Direction, VectorType); + itkSetMacro(Radius, double); + itkGetMacro(Radius, double); + itkSetMacro(Height, double); + itkGetMacro(Height, double); + itkGetMacro(CylinderMean, double); + itkGetMacro(CylinderStd, double); + itkGetMacro(MaxRelativeDifference, double); + itkGetObjectMacro(SliceMeasurements, MeasurementVectorType); + itkGetObjectMacro(SliceOffsets, MeasurementVectorType); + + itkSetMacro(LabelInside, OutputPixelType); + itkGetMacro(LabelInside, OutputPixelType); + itkSetMacro(LabelInsideRadiusLimit, OutputPixelType); + itkGetMacro(LabelInsideRadiusLimit, OutputPixelType); + itkSetMacro(LabelInsideHeightLimit, OutputPixelType); + itkGetMacro(LabelInsideHeightLimit, OutputPixelType); + itkSetMacro(LabelOutside, OutputPixelType); + itkGetMacro(LabelOutside, OutputPixelType); + + bool IsInside(PointType p) const; + double DistanceFromAxis(PointType p) const; + double DistanceAlongAxis(PointType p) const; + +protected: + CylinderUniformityMeasurementImageFilter(); + ~CylinderUniformityMeasurementImageFilter(){} + + virtual void PrintSelf(std::ostream& os, Indent indent) const; + + /** Does the real work. */ + virtual void GenerateData(); + +private: + CylinderUniformityMeasurementImageFilter(const Self &); //purposely not implemented + void operator=(const Self &); //purposely not implemented + + PointType m_Center; + VectorType m_Direction; + double m_Radius; + double m_Height; + MeasurementVectorType::Pointer m_SliceMeasurements; + MeasurementVectorType::Pointer m_SliceOffsets; + double m_CylinderMean; + double m_CylinderStd; + double m_MaxRelativeDifference; + OutputPixelType m_LabelInside; + OutputPixelType m_LabelInsideRadiusLimit; + OutputPixelType m_LabelInsideHeightLimit; + OutputPixelType m_LabelOutside; +}; +} //namespace ITK + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkCylinderUniformityMeasurementImageFilter.cxx" +#endif + +#endif // __itkCylinderUniformityMeasurementImageFilter_h diff --git a/README.md b/README.md new file mode 100644 index 0000000..b45b7d7 --- /dev/null +++ b/README.md @@ -0,0 +1,11 @@ +SlicerPETPhantomAnalysis Extension +===================================== +3D Slicer Extension for Analysis of Uniform Cylinder Phantoms in PET scans. + +Acknowledgements +-------- + +This work is funded in part by [Quantitative Imaging to Assess Response in Cancer Therapy Trials][] NIH grant U01-CA140206 (PIs John Buatti, Tom Casavant, Michael Graham, Milan Sonka) and [Quantitative Image Informatics for Cancer Research][] (QIICR) NIH grant U24 CA180918. + +[Quantitative Imaging to Assess Response in Cancer Therapy Trials]: http://imaging.cancer.gov/programsandresources/specializedinitiatives/qin/iowa +[Quantitative Image Informatics for Cancer Research]: http://qiicr.org