Merge pull request #1409 from sarthakpati/dicom

Some dicom-related fixes

CMakeLists.txt

@@ -323,7 +323,6 @@ IF( NOT BUILD_DOCUMENTATION_ONLY )
     FIND_PACKAGE(Qt5 COMPONENTS X11Extras REQUIRED)
   ENDIF()
 
-  LINK_DIRECTORIES(${QT_LIBRARY_DIR})
   #INCLUDE(${QT_USE_FILE})
   SET(QT_USE_QTGUI TRUE)
   SET(QT_USE_QTMAIN TRUE)

src/applications/FeatureExtraction/src/FeatureExtraction.hxx

@@ -753,6 +753,16 @@ template< class TImage >
 void FeatureExtraction< TImage >::CalculateIntensity(std::vector< typename TImage::PixelType >& nonZeroVoxels, std::map< std::string, double >& featurevec, bool latticePatch)
 {
   cbica::Statistics< typename TImage::PixelType > statisticsCalculatorToUse;
+  if (m_QuantizationExtent.empty())
+  {
+    std::cerr << "'Quantization_Extent' needs to be defined under 'Generic'.\n";
+    exit(EXIT_FAILURE);
+  }
+  if (m_histogramBinningType == -1)
+  {
+    std::cerr << "'Quantization_Type' needs to be defined under 'Generic'.\n";
+    exit(EXIT_FAILURE);
+  }
   if (m_QuantizationExtent == "Image")
   {
     statisticsCalculatorToUse = m_statistics_global[m_currentROIValue];

src/applications/Preprocessing/src/Preprocessing.cxx

@@ -47,7 +47,7 @@ registrationAffineTransformInput, registrationDeformableTransformInput;
 
 int histoMatchQuantiles = 40, histoMatchBins = 100,
   registrationTypeInt, registrationRigidDof = 12;
-bool registrationIntermediate = false, registrationSegmentationMoving = false, registrationJacobian = false;
+bool registrationIntermediate = false, registrationSegmentationMoving = false, registrationJacobian = false, registrationInvert = false;
 float zNormCutLow = 3, zNormCutHigh = 3, zNormQuantLow = 5, zNormQuantHigh = 95,
   bias_fwhm = BiasCorrection::default_fwhm, rescaleLower = 0, rescaleUpper = 1000,
   ssSigma = 0.5, ssIntensityThreshold = 80,
@@ -513,6 +513,10 @@ int algorithmsRunner()
       }
       default: // we shall always assume affine/rigid
       {
+        if (registrationInvert)
+        {
+          interimFiles_affineTransform += ",-1";
+        }
         if (registrationSegmentationMoving)
         {
           commandToCall = greedyPathAndDim +
@@ -706,6 +710,7 @@ int main(int argc, char** argv)
   parser.addOptionalParameter("rIS", "regInterSave", cbica::Parameter::BOOLEAN, "0 or 1", "Whether the intermediate files are to be saved or not", "Defaults to " + std::to_string(registrationIntermediate));
   parser.addOptionalParameter("rSg", "regSegMoving", cbica::Parameter::BOOLEAN, "0 or 1", "Whether the Moving Image(s) is a segmentation file", "If 1, the 'Nearest Label' Interpolation is applied", "Defaults to " + std::to_string(registrationSegmentationMoving));
   parser.addOptionalParameter("rJC", "regJacobian", cbica::Parameter::BOOLEAN, "0 or 1", "Whether Jacobian will be generated or not");
+  parser.addOptionalParameter("rIv", "regInverse", cbica::Parameter::BOOLEAN, "0 or 1", "Whether the affine matrix should be inverted or not, defaults to false");
   parser.addOptionalParameter("rIA", "regInterAffn", cbica::Parameter::FILE, "mat", "The path to the affine transformation to apply to moving image", "If this is present, the Affine registration step will be skipped", "Also used for rigid transformation");
   parser.addOptionalParameter("rID", "regInterDefm", cbica::Parameter::FILE, "NIfTI", "The path to the deformable transformation to apply to moving image", "If this is present, the Deformable registration step will be skipped");
   parser.addOptionalParameter("rsc", "rescaleImage", cbica::Parameter::STRING, "Output Intensity range", "The output intensity range after image rescaling", "Defaults to " + std::to_string(rescaleLower) + ":" + std::to_string(rescaleUpper), "If multiple inputs are passed (comma-separated), the rescaling is done in a cumulative manner,", "i.e., stats from all images are considered for the scaling");
@@ -974,6 +979,10 @@ int main(int argc, char** argv)
     {
       parser.getParameterValue("rJC", registrationJacobian);
     }
+    if (parser.isPresent("rIv"))
+    {
+      parser.getParameterValue("rIv", registrationInvert);
+    }
     if (parser.isPresent("rIA"))
     {
       parser.getParameterValue("rIA", registrationAffineTransformInput);

src/applications/Utilities/src/Utilities.cxx

@@ -1180,7 +1180,7 @@ int algorithmsRunner_imageStack2join(std::vector< std::string >& inputImageFiles
   {
     inputImages[i] = cbica::ReadImage< TImageType >(inputImageFiles[i]);
   }
-  auto output = cbica::GetJoinedImage< TImageType, TOutputImageType >(inputImages, imageStack2JoinSpacing);
+  auto output = cbica::GetJoinedImage< TImageType, TOutputImageType >(inputImages, imageStack2JoinSpacing, nifti2dicomSpacingTolerance);
   cbica::WriteImage< TOutputImageType >(output, outputImageFile);
   return EXIT_SUCCESS;
 }
@@ -1262,6 +1262,7 @@ int main(int argc, char** argv)
   parser.addOptionalParameter("w2i", "world2image", cbica::Parameter::STRING, "i,j,k", "The world coordinates that will be converted to image coordinates for the input image", "Example: '-w2i 10.5,20.6,30.2'");
   parser.addOptionalParameter("j2e", "joined2extracted", cbica::Parameter::BOOLEAN, "0-1", "Axis to extract is always the final axis (axis '3' for a 4D image)", "The '-o' parameter can be used for output: '-o /path/to/extracted_'");
   parser.addOptionalParameter("e2j", "extracted2joined", cbica::Parameter::FLOAT, "0-10", "The spacing in the new direction", "Pass the folder containing all images in '-i'");
+  parser.addOptionalParameter("ejT", "extrac2joinToler", cbica::Parameter::FLOAT, "0-10", "The spacing tolerance in the original axes", "Only used of 'e2j' is passed", "Defaults to " + std::to_string(nifti2dicomSpacingTolerance));
   parser.addOptionalParameter("ls", "labelSimilarity", cbica::Parameter::FILE, "NIfTI Reference", "Calculate similarity measures for 2 label maps", "Pass the reference map after '-ls' and the comparison will be done with '-i'", "For images with more than 2 labels, individual label stats are also presented");
   parser.addOptionalParameter("lsb", "lSimilarityBrats", cbica::Parameter::FILE, "NIfTI Reference", "Calculate BraTS similarity measures for 2 brain labels", "Pass the reference map after '-lsb' and the comparison will be done with '-i'", "Assumed labels in image are '1,2,4' and missing labels will be populate with '0'");
   parser.addOptionalParameter("hd", "hausdorffDist", cbica::Parameter::FILE, "NIfTI Reference", "Calculate the Hausdorff Distance for the input image and", "the one passed after '-hd'");
@@ -1354,7 +1355,19 @@ int main(int argc, char** argv)
   else if (parser.isPresent("n2d"))
   {
     requestedAlgorithm = Nifti2Dicom;
-    parser.getParameterValue("n2d", targetImageFile); // in this case, it is the DICOM reference file
+    int tempPosition;
+    if (parser.compareParameter("n2d", tempPosition))
+    {
+      std::string tempValue;
+      if (argc > tempPosition + 1)
+      {
+        tempValue = argv[tempPosition + 1];
+        if (cbica::exists(tempValue))
+        {
+          targetImageFile = tempValue; // in this case, it is the DICOM reference file
+        }
+      }
+    }
     parser.getParameterValue("o", outputImageFile);
     if (parser.isPresent("ndD"))
     {
@@ -1647,6 +1660,10 @@ int main(int argc, char** argv)
       return EXIT_FAILURE;
     }
     parser.getParameterValue("e2j", imageStack2JoinSpacing);
+    if (parser.isPresent("ejT"))
+    {
+      parser.getParameterValue("ejT", nifti2dicomSpacingTolerance);
+    }
 
     auto imagesToJoin = cbica::filesInDirectory(inputImageFile);
     if (imagesToJoin.empty())

src/cbica_toolkit/src/cbicaITKSafeImageIO.h

@@ -743,6 +743,69 @@ namespace cbica
 
   This uses the dictionary created by the reference DICOM image
 
+  \param dicomImageReferenceDir Reference DICOM image
+  \param imageToWrite Pointer to processed image data which is to be written
+  \param outputDir File containing the image
+  \param nifti2dicomTolerance The direction tolerance for sanity check
+  \param outputPrefix The prefix for the output filenames
+  \return itk::Image of specified pixel and dimension type
+  */
+  template <typename ComputedImageType = ImageTypeFloat3D>
+  void WriteDicomImage(typename ComputedImageType::Pointer imageToWrite, const std::string &outputDir, const std::string& outputPrefix = "image")
+  {
+    using DicomImageType = itk::Image< int, ComputedImageType::ImageDimension >; // dicom currently doesn't support float properly
+
+    typedef itk::CastImageFilter<ComputedImageType, DicomImageType> CastFilterType;
+    typename CastFilterType::Pointer castFilter = CastFilterType::New();
+    castFilter->SetInput(imageToWrite);
+    castFilter->Update();
+
+    // initialize the image io
+    auto dicomIO = itk::GDCMImageIO::New();
+    dicomIO->SetComponentType(itk::ImageIOBase::IOComponentType::INT);
+
+    auto dictionary = dicomIO->GetMetaDataDictionary();
+
+    itk::EncapsulateMetaData<std::string>(dictionary, "0008|0008", "DERIVED\\SECONDARY"); // Image Type
+    itk::EncapsulateMetaData<std::string>(dictionary, "0008|103e", "Processed-CaPTk"); // # Series Description
+    itk::EncapsulateMetaData<std::string>(dictionary, "0008|0030", cbica::getCurrentLocalTimestamp()); // # Study Time - ensures unique string
+
+    dicomIO->SetMetaDataDictionary(dictionary);
+
+    auto seriesWriter = itk::ImageSeriesWriter< DicomImageType, itk::Image<typename DicomImageType::PixelType, 2> >::New(); // always write in 2D series
+
+    cbica::createDir(outputDir);
+
+    auto namesGenerator = itk::NumericSeriesFileNames::New();
+    //namesGenerator->SetUseSeriesDetails(false);
+    auto start = imageToWrite->GetLargestPossibleRegion().GetIndex();
+    auto size = imageToWrite->GetLargestPossibleRegion().GetSize();
+    namesGenerator->SetSeriesFormat((outputDir + "/" + outputPrefix + "%03d.dcm").c_str());
+    namesGenerator->SetStartIndex(start[2]);
+    namesGenerator->SetEndIndex(start[2] + size[2] - 1);
+    namesGenerator->SetIncrementIndex(1);
+
+    seriesWriter->SetInput(castFilter->GetOutput());
+    seriesWriter->SetFileNames(namesGenerator->GetFileNames());
+    seriesWriter->SetImageIO(dicomIO);
+    seriesWriter->SetMetaDataDictionary(dictionary);
+
+    try
+    {
+      seriesWriter->Write();
+    }
+    catch (itk::ExceptionObject &e)
+    {
+      std::cerr << "Error occurred while trying to write the image '" << outputDir << "': " << e.what() << "\n";
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  /**
+  \brief Write itk::ImageReader as DICOM to specified directory
+
+  This uses the dictionary created by the reference DICOM image
+
   \param dicomImageReferenceDir Reference DICOM image
   \param imageToWrite Pointer to processed image data which is to be written
   \param outputDir File containing the image
@@ -756,6 +819,11 @@ namespace cbica
     const std::string &outputDir, const float nifti2dicomTolerance = 0.0, const float nifti2dicomOriginTolerance = 0.0, const float nifti2dicomSpacingTolerance = 0.0,
     const std::string& outputPrefix = "image")
   {
+    if (dicomImageReferenceDir.empty()) // when no reference is present, switch to using default dicomIO
+    {
+      return WriteDicomImage< ComputedImageType >(imageToWrite, outputDir, outputPrefix);
+    }
+
     if (cbica::isDir(dicomImageReferenceDir))
     {
       using DicomImageType = itk::Image< int, ComputedImageType::ImageDimension >;

src/cbica_toolkit/src/cbicaITKUtilities.h

@@ -371,20 +371,20 @@ namespace cbica
     // perform an absolute check when either of the direction cosine location is zero
     if ((input_1 == 0) || (input_2 == 0))
     {
-      if (percentageDifference < threshold)
+      if (percentageDifference <= threshold)
       {
         percentageDifference = 0;
       }
       else
       {
         // otherwise, calculate the percentage based on the non-zero cosine
-        if (input_1 == 0)
+        if (input_2 != 0)
         {
-          percentageDifference /= input_1;
+          percentageDifference /= std::abs(input_2);
         }
-        else if (input_2 == 0)
+        else if (input_1 != 0)
         {
-          percentageDifference /= input_1;
+          percentageDifference /= std::abs(input_1);
         }
       }
     }
@@ -610,9 +610,10 @@ namespace cbica
 
   \param inputImage The vector of images from which the larger image is to be extracted
   \param newSpacing The spacing in the new dimension
+  \param spacingTolerance The spacing tolerance when processing individual images
   */
   template< class TInputImageType, class TOutputImageType >
-  typename TOutputImageType::Pointer GetJoinedImage(std::vector< typename TInputImageType::Pointer > &inputImages, double newSpacing = 1.0)
+  typename TOutputImageType::Pointer GetJoinedImage(std::vector< typename TInputImageType::Pointer > &inputImages, double newSpacing = 1.0, double spacingTolerance = 5.0)
   {
     if (TOutputImageType::ImageDimension - 1 != TInputImageType::ImageDimension)
     {
@@ -625,7 +626,7 @@ namespace cbica
 
     for (size_t N = 0; N < inputImages.size(); N++)
     {
-      if (!ImageSanityCheck< TInputImageType >(inputImages[0], inputImages[N]))
+      if (!ImageSanityCheck< TInputImageType >(inputImages[0], inputImages[N], 0, 0, spacingTolerance))
       {
         std::cerr << "Image Sanity check failed in index '" << N << "'\n";
         //return typename TOutputImageType::New();

src/cbica_toolkit/src/cbicaUtilities.cpp

@@ -2232,7 +2232,9 @@ namespace cbica
         return 0L;			/* Failed? */
 
     #endif
-  #endif
+#else
+    return 0;
+#endif
   }
 
   size_t getCurrentlyUsedMemory()

src/view/gui/SlicerManager.h

@@ -115,6 +115,11 @@ class SlicerManager : public QObject
     return mMask;
   }
 
+  void SetPathFileName(std::string file)
+  {
+    mPathFileName = file;
+  }
+
   //! Get the full file path of input image
   std::string GetPathFileName()
   {

src/view/gui/fMainWindow.cpp

@@ -5749,38 +5749,39 @@ void fMainWindow::openImages(QStringList files, bool callingFromCmd)
   /**** Image Loading ****/
 
   int i = 0, fileSizeCheck = files.size() + 1;
+  auto fileForSanityCheck = files[0].toStdString();
   if (mSlicerManagers.empty())
   {
+    std::string fileName = files[i].toStdString();
+    fileName = cbica::normPath(fileName);
+    updateProgress(i + 1, "Opening " + fileName, files.size());
+    //auto extension = cbica::getFilenameExtension(fileName);
+    //if (!extension.empty())
+    //{
+    //  std::transform(extension.begin(), extension.end(), extension.begin(), ::tolower);
+    //}
+    //if ((extension == ".dcm") || (extension == ".dicom") || (extension == "") ||
+    //  (extension == ".ima"))
+    if (cbica::IsDicom(fileName))
     {
-      std::string fileName = files[i].toStdString();
-      fileName = cbica::normPath(fileName);
-      updateProgress(i + 1, "Opening " + fileName, files.size());
-      //auto extension = cbica::getFilenameExtension(fileName);
-      //if (!extension.empty())
-      //{
-      //  std::transform(extension.begin(), extension.end(), extension.begin(), ::tolower);
-      //}
-      //if ((extension == ".dcm") || (extension == ".dicom") || (extension == "") ||
-      //  (extension == ".ima"))
-      if (cbica::IsDicom(fileName))
-      {
-        QDir d = QFileInfo(fileName.c_str()).absoluteDir();
-        QString fname = d.absolutePath();
-        dicomfilename = fileName;
-        this->openDicomImages(fname);
-      }
-      else
-      {
-        LoadSlicerImages(fileName, CAPTK::ImageExtension::NIfTI);
-      }
+      QDir d = QFileInfo(fileName.c_str()).absoluteDir();
+      QString fname = d.absolutePath();
+      dicomfilename = fileName;
+      this->openDicomImages(fname);
+    }
+    else
+    {
+      LoadSlicerImages(fileName, CAPTK::ImageExtension::NIfTI);
     }
     fileSizeCheck = 1;
   }
   else
   {
     fileSizeCheck = 0;
+    fileForSanityCheck = mSlicerManagers[0]->GetPathFileName();
   }
 
+  bool c = false;
   // basic sanity check
   if (files.size() > fileSizeCheck)
   {
@@ -5800,7 +5801,7 @@ void fMainWindow::openImages(QStringList files, bool callingFromCmd)
       {
         unsupportedExtension += fileName + "\n";
       }
-      else if (!cbica::ImageSanityCheck(files[0].toStdString(), files[i].toStdString()))
+      else if (!cbica::ImageSanityCheck(fileForSanityCheck, files[i].toStdString()))
       {
         erroredFiles += fileName + "\n";
       }
@@ -5822,7 +5823,7 @@ void fMainWindow::openImages(QStringList files, bool callingFromCmd)
     }
     if (!erroredFiles.empty())
     {
-      ShowErrorMessage("Extensions for the following files were not supported, CaPTk will try to load the rest:\n\n" + unsupportedExtension, this);
+      ShowErrorMessage("Sanity check (dimensions/size/origin/spacing/directions with original image) for the following files did not pass, CaPTk will try to load the rest:\n\n" + unsupportedExtension, this);
     }
 
     for (int i = 0; i < basicSanityChecksPassedFiles.size(); i++)
@@ -5916,6 +5917,8 @@ void fMainWindow::openDicomImages(QString dir)
   QDir d(dir);
   seriesDescValue = d.dirName().toStdString(); 
   imageManager->SetFilename(seriesDescValue);
+  auto filesInDir = cbica::filesInDirectory(dir.toStdString());
+  imageManager->SetPathFileName(filesInDir[0].c_str());
 
   QString id = QString(seriesDescValue.c_str()) + QString::number(mSlicerManagers.size() - 1);
   //