https://github.com/ElsevierSoftwareX/SOFTX_2019_219
Raw File
Tip revision: 99932db9e71d31a76425d68eef7d5859523ca0e4 authored by Konrad Werys on 04 November 2019, 09:33:32 UTC
doc: changes in the documentation
Tip revision: 99932db
itkReadFileListFilter.h
//
//  ReadFileListFilter
//  TomatoLib
//
//  Created by Konrad Werys on 19/08/17.
//  Copyright © 2017 Konrad Werys. All rights reserved.
//

#ifndef TomatoLIB_ITKREADFILELISTFILTER_H
#define TomatoLIB_ITKREADFILELISTFILTER_H

#include "CmakeConfigForTomato.h"

#ifdef USE_ITK

#include <iostream>
#include "itkImageSource.h"
#include "itkImageFileReader.h"
#include "itkTileImageFilter.h"
#include "itkGDCMImageIO.h"
#include "itkMetaDataObject.h"
#include "itkStatisticsImageFilter.h"
#include "KWUtil.h"

namespace itk {
    template< typename TImage>
    class ReadFileListFilter :public ImageSource< TImage > {
    public:

        /** Standard class typedefs. */
        typedef ReadFileListFilter Self;
        typedef ImageSource <TImage> Superclass;
        typedef SmartPointer <Self> Pointer;

        /** Method for creation through the object factory. */
        itkNewMacro(Self);

        /** Run-time type information (and related methods). */

        typedef TImage ImageType3D;
        typedef typename ImageType3D::PixelType PixelType;
        typedef itk::Image< PixelType, 2 > ImageType2D;
        typedef itk::ImageFileReader< ImageType3D > ReaderType;
        typedef itk::TileImageFilter< ImageType3D, ImageType3D > TileImageType;

        typedef std::vector<std::string> FileListType;
        typedef std::vector<MetaDataDictionary> MetaDataDictionaryArrayType;

        // setters
        void SetFileList(std::vector<std::string> fileList){
            m_FileList = fileList;
        }
        //itkSetMacro(MetaDataDictionaryArray, MetaDataDictionaryArrayType);
        itkSetMacro(Verbose, bool);

        template< typename TYPE >
        vnl_vector<TYPE> GetVnlVectorFromStdVector(std::vector<TYPE> stdVector){
            if (stdVector.size() == 0){
                return vnl_vector<TYPE>();
            } else {
                return vnl_vector<TYPE>(&stdVector.front(), stdVector.size());
            }
        }

        // getters
        vnl_vector<double> GetInvTimes(){
            vnl_vector<double> invTimeVnl = GetVnlVectorFromStdVector(m_InvTimes);
            vnl_vector<double> invTimeFromImageCommentsVnl = GetVnlVectorFromStdVector(m_InvTimesFromImageComments);
            vnl_vector<double> triggerTimesVnl = GetVnlVectorFromStdVector(m_TriggerTimes);
            //vnl_vector<double> invTimeVnl(m_InvTimes.data(),m_InvTimes.size());
            //vnl_vector<double> invTimeFromImageCommentsVnl(m_InvTimesFromImageComments.data(),m_InvTimesFromImageComments.size());
            //vnl_vector<double> triggerTimesVnl(m_TriggerTimes.data(),m_TriggerTimes.size());
            if (invTimeVnl.min_value() != invTimeVnl.max_value()) {
                return invTimeVnl;
            } else if (invTimeFromImageCommentsVnl.min_value() != invTimeFromImageCommentsVnl.max_value()) {
                return invTimeFromImageCommentsVnl;
            } else {
                return triggerTimesVnl;
            }
        }

        vnl_vector<double> GetRepTimes(){
            return GetVnlVectorFromStdVector(m_RepTimes);
            //return vnl_vector<double>(m_RepTimes.data(),m_RepTimes.size());
        }

        vnl_vector<double> GetEchoTimes(){
            return GetVnlVectorFromStdVector(m_EchoTimes);
            //return vnl_vector<double>(m_EchoTimes.data(),m_EchoTimes.size());
        }

        vnl_vector<double> GetTriggerTimes(){
            return GetVnlVectorFromStdVector(m_TriggerTimes);
            //return vnl_vector<double>(m_TriggerTimes.data(),m_TriggerTimes.size());
        }

        vnl_vector<double> GetAcqTimes(){
            return GetVnlVectorFromStdVector(m_AcqTimes);
            //return vnl_vector<double>(m_RelAcqTimes.data(),m_RelAcqTimes.size());
        }

        vnl_vector<double> GetRelAcqTimes(){
            vnl_vector<double> temp = GetAcqTimes();
            double mymin = temp.min_value();
            for (unsigned int i = 0; i < temp.size(); ++i){
                temp[i] = temp[i] - mymin;
            }
            return temp;
        }

        itkGetMacro(MetaDataDictionaryArray, MetaDataDictionaryArrayType);
        itkGetObjectMacro(DicomIO, GDCMImageIO);
        itkGetMacro(Verbose, bool);

    protected:
        ReadFileListFilter(){
            m_Verbose = false;
        };
        ~ReadFileListFilter(){};

        /** Does the real work. */
        virtual void GenerateData() ITK_OVERRIDE;

    private:
        ITK_DISALLOW_COPY_AND_ASSIGN(ReadFileListFilter);

        std::vector<std::string> m_FileList;
        std::vector<double> m_InvTimes;
        std::vector<double> m_InvTimesFromImageComments;
        std::vector<double> m_RepTimes;
        std::vector<double> m_EchoTimes;
        std::vector<double> m_TriggerTimes;
        std::vector<double> m_AcqTimes;
        MetaDataDictionaryArrayType m_MetaDataDictionaryArray;
        GDCMImageIO::Pointer m_DicomIO;
        bool m_Verbose;

        double FindInversionTime(ReaderType* reader);
        double FindInversionTimeFromImageComments(ReaderType* reader);
        double FindRepetitionTime(ReaderType* reader);
        double FindEchoTime(ReaderType* reader);
        double FindTriggerTime(ReaderType* reader);
        double FindAcqTime(ReaderType* reader);
        double AcqTime2Seconds(std::string acqTimes);
    };

    template< class TImage >
    void
    ReadFileListFilter< TImage >
    ::GenerateData() {
        // see https://itk.org/Wiki/ITK/Examples/ImageProcessing/TileImageFilter_CreateVolume
        typename ImageType3D::Pointer outputImage;

        m_DicomIO = GDCMImageIO::New();

        itk::FixedArray< PixelType, 3 > layout;
        layout[0] = 1;
        layout[1] = 1;
        layout[2] = m_FileList.size();

        // tiler
        typename TileImageType::Pointer tiler = TileImageType::New();
        tiler->SetLayout( layout );

        typename ReaderType::Pointer reader;

        // the main loop
        unsigned int inputImageNumber = 0;
        for ( unsigned int i = 0; i < m_FileList.size(); ++i ){

            //reader
            reader = ReaderType::New();
            reader->SetImageIO(m_DicomIO);
            reader->SetFileName( m_FileList[i] );
            //reader->MetaDataDictionaryArrayUpdateOn();
            tiler->SetInput(i,reader->GetOutput());

            try {
                tiler->Update();
                //std::cout << std::setprecision(9) << reader->GetOutput()->GetOrigin() << std::endl;
                m_MetaDataDictionaryArray.push_back(reader->GetMetaDataDictionary());
                m_InvTimes.push_back(FindInversionTime(reader));
                m_InvTimesFromImageComments.push_back(FindInversionTimeFromImageComments(reader));
                m_RepTimes.push_back(FindRepetitionTime(reader));
                m_EchoTimes.push_back(FindEchoTime(reader));
                m_TriggerTimes.push_back(FindTriggerTime(reader));
                m_AcqTimes.push_back(FindAcqTime(reader));
                inputImageNumber++;
            } catch( itk::ExceptionObject & err ) {
                std::cerr << "Unable to read file: " << m_FileList[i] << std::endl;
                std::cerr << "ExceptionObject caught!" << std::endl << err << std::endl;
            }
        }

        if (inputImageNumber == 0) return;

        tiler->UpdateLargestPossibleRegion();
        // the final update
        outputImage = tiler->GetOutput();
        //outputImage->SetOrigin(reader->GetOutput()->GetOrigin());
        //outputImage->SetDirection(reader->GetOutput()->GetDirection());

        try {
            outputImage->Update();
        } catch( itk::ExceptionObject & err ) {
            std::cerr << "ExceptionObject caught!" << std::endl;
            std::cerr << err << std::endl;
        }

        if (GetVerbose()) {
            //
            typename ImageType3D::PixelType rescaleIntercept = m_DicomIO->GetRescaleIntercept();
            typename ImageType3D::PixelType rescaleSlope = m_DicomIO->GetRescaleSlope();

            typedef itk::StatisticsImageFilter<ImageType3D> StatisticsImageFilterType;
            typename StatisticsImageFilterType::Pointer statisticsImageFilter
                    = StatisticsImageFilterType::New();
            statisticsImageFilter->SetInput(outputImage);
            statisticsImageFilter->Update();

            std::cout << "Mean: " << statisticsImageFilter->GetMean();
            std::cout << ", Std.: " << statisticsImageFilter->GetSigma();
            std::cout << ", Min: " << statisticsImageFilter->GetMinimum();
            std::cout << ", Max: " << statisticsImageFilter->GetMaximum();

            std::cout << ", RescaleIntercept: " << rescaleIntercept << ", RescaleSlope: " << rescaleSlope << std::endl;

            //std::cout << outputImage << std::endl;

            vcl_cout << "Inversion Times: " << GetInvTimes() << vcl_endl;
            vcl_cout << "Repetition Times: " << GetRepTimes() << vcl_endl;
            vcl_cout << "Echo Times: " << GetEchoTimes() << vcl_endl;
            vcl_cout << "Trigger Times: " << GetTriggerTimes() << vcl_endl;
            vcl_cout << "Rel Acq Times: " << GetRelAcqTimes() << vcl_endl;
        }

        this->GraftOutput(outputImage);
    }

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindInversionTime(ReaderType* reader) {

        typedef itk::MetaDataDictionary DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        double invTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0018|0082";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);

        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr = dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                invTime = KWUtil::StringToNumber<double>(tagvalue);
                //invTime = (std::stof(tagvalue));
            }
        }
        return invTime;
    };

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindInversionTimeFromImageComments(ReaderType* reader) {

        typedef itk::MetaDataDictionary DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        int invTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0020|4000";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);

        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr = dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                std::sscanf(tagvalue.c_str(), "TIeff %i ms ", &invTime);
            }
        }
        return (double)invTime;
    };

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindRepetitionTime(ReaderType* reader) {

        typedef itk::MetaDataDictionary   DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        double repTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0018|0080";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr =
                    dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                repTime = KWUtil::StringToNumber<double>(tagvalue);
                //repTime = (std::stof(tagvalue));
            }
        }
        return repTime;
    };

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindEchoTime(ReaderType* reader) {

        typedef itk::MetaDataDictionary   DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        double echoTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0018|0081";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr =
                    dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                echoTime = KWUtil::StringToNumber<double>(tagvalue);
                //echoTime = (std::stof(tagvalue));
            }
        }
        return echoTime;
    };

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindTriggerTime(ReaderType* reader) {

        typedef itk::MetaDataDictionary   DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        double triggerTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0018|1060";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr =
                    dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                triggerTime = KWUtil::StringToNumber<double>(tagvalue);
                //triggerTime = (std::stof(tagvalue));
            }
        }
        return triggerTime;
    };

    template< class TImage>
    double
    ReadFileListFilter< TImage>
    ::FindAcqTime(ReaderType* reader) {

        typedef itk::MetaDataDictionary   DictionaryType;
        typedef itk::MetaDataObject< std::string > MetaDataStringType;
        double acqTime = 0;

        const  DictionaryType & dictionary = reader->GetMetaDataDictionary();
        DictionaryType::ConstIterator end = dictionary.End();

        std::string entryId = "0008|0032";
        DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
        if ( tagItr != end ) {
            MetaDataStringType::ConstPointer entryvalueStr =
                    dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());

            if (entryvalueStr) {
                std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
                acqTime = KWUtil::dicomTime2Seconds<double>(tagvalue);
            }
        }
        return acqTime;
    };

}

#endif // USE_ITK

#endif //TomatoLIB_ITKREADFILELISTFILTER_H
back to top