Merge pull request #8 from yhuang43/dev

fix collapsed label handling (freesurfer/freesurfer commit db3bdc7)

gems/Testing/CMakeLists.txt

@@ -8,7 +8,10 @@ endif()
 add_executable(kvlAtlasMeshRasterizorTestGPU kvlAtlasMeshRasterizorTestGPU.cxx)
 target_link_libraries(kvlAtlasMeshRasterizorTestGPU kvlGEMSCommon)
 
-#
+# testCompressionLookupTable
+add_executable(testCompressionLookupTable testCompressionLookupTable.cxx)
+target_link_libraries(testCompressionLookupTable kvlGEMSCommon)
+
 add_executable(testdeformMeshPython testdeformMeshPython.cxx)
 target_link_libraries(testdeformMeshPython kvlGEMSCommon)
 

gems/Testing/testCompressionLookupTable.cxx

@@ -0,0 +1,62 @@
+#include "itkImageFileReader.h"
+#include "itkMGHImageIOFactory.h"
+#include "kvlCompressionLookupTable.h"
+
+int main(int argc, char **argv)
+{
+  // Add support for MGH file format to ITK. An alternative way to add this by default would be
+  // to edit ITK's itkImageIOFactory.cxx and explicitly adding it in the code there.
+  itk::ObjectFactoryBase::RegisterFactory( itk::MGHImageIOFactory::New() );
+
+  // Read the input images
+  typedef kvl::CompressionLookupTable::ImageType  LabelImageType;
+  std::vector< LabelImageType::ConstPointer >  labelImages;
+  for ( int argumentNumber = 1; argumentNumber < argc; argumentNumber++ )
+  {
+    std::cout << "Reading input image: " << argv[ argumentNumber ] << std::endl;
+    // Read the input image
+    typedef itk::ImageFileReader< LabelImageType >  ReaderType;
+    ReaderType::Pointer  reader = ReaderType::New();
+    reader->SetFileName( argv[ argumentNumber ] );
+    reader->Update();
+    LabelImageType::ConstPointer  labelImage = reader->GetOutput();
+
+    // Over-ride the spacing and origin since at this point we can't deal with that
+    const double spacing[] = { 1, 1, 1 };
+    const double origin[] = { 0, 0, 0 };
+    const_cast< LabelImageType* >( labelImage.GetPointer() )->SetSpacing( spacing );
+    const_cast< LabelImageType* >( labelImage.GetPointer() )->SetOrigin( origin );
+
+    // Remember this image
+    labelImages.push_back( labelImage );
+  }  
+
+  // Build a lookup table that maps the original intensities onto class numbers starting
+  // at 0 and densely packed
+  kvl::CompressionLookupTable::Pointer  lookupTable = kvl::CompressionLookupTable::New();
+  lookupTable->Construct( labelImages );
+  lookupTable->Write( "compressionLookupTable.txt" );
+
+  /*********** START OF SIMULATION ***********/
+  printf("\nClasses Contributed To Label Cost Calculations:\n");
+  // loop through each label,
+  // report which classes will contribute to cost/gradient/likelihood/color calculation
+  std::vector<LabelImageType::PixelType> labels = lookupTable->GetLabels();
+  std::vector<LabelImageType::PixelType>::const_iterator labelIt;
+  for (labelIt = labels.begin(); labelIt != labels.end(); labelIt++)
+  {
+    printf("label %5d <= ", *labelIt);
+    const std::vector< int >& classNumbers = lookupTable->GetClassNumbers(*labelIt);
+    for ( std::vector< int >::const_iterator classIt = classNumbers.begin(); 
+          classIt != classNumbers.end();
+          ++classIt )
+      printf(" %3d, ", *classIt);
+
+    printf("\n");
+  }
+
+  printf("\nTotal # of Labels: %d\n", labels.size());
+  printf("Total # of Classes: %d\n", lookupTable->GetNumberOfClasses());
+
+  exit(0);
+}

gems/kvlCompressionLookupTable.cxx

@@ -108,6 +108,8 @@ ::Construct( const std::vector< ImageType::ConstPointer >& images )
   m_LabelStringLookupTable.clear();
   m_ColorLookupTable.clear();
   m_NumberOfClasses = 0;
+
+  printf("\n");
 
   // Loop over all collapsed labels, if any, adding a new class for 
   // each real label that is encountered while also pointing the
@@ -128,25 +130,35 @@ ::Construct( const std::vector< ImageType::ConstPointer >& images )
         collapsedIt != collapsedLabels.end(); ++collapsedIt )
     {
     m_CompressionLookupTable[ collapsedIt->first ] = std::vector< int >();
+    printf("Collapsed Label %-5d : \n", collapsedIt->first);
+    // loop through the Collapsed Label list, add any labels not in m_CompressionLookupTable yet
     for ( std::vector< int >::const_iterator  labelIt = collapsedIt->second.begin();
           labelIt != collapsedIt->second.end(); ++labelIt )
       {
-      if ( m_CompressionLookupTable.find( *labelIt ) == m_CompressionLookupTable.end() )
-        {
-        std::cout << "Encountered new real label " << *labelIt << std::endl;
-
+        CompressionLookupTableType::const_iterator  compressionIt = m_CompressionLookupTable.find( *labelIt );	  
+	if ( compressionIt == m_CompressionLookupTable.end() )  //if ( m_CompressionLookupTable.find( *labelIt ) == m_CompressionLookupTable.end() )
+        {      
         const int  newClassNumber = m_NumberOfClasses;
         m_CompressionLookupTable[ *labelIt ] = std::vector< int >( 1, newClassNumber );
         m_CompressionLookupTable[ collapsedIt->first ].push_back( newClassNumber );
+	printf("\t Encountered new real label (collapsed label list): %-5d -> %-3d\n", *labelIt, newClassNumber);
+        //std::cout << "Encountered new real label (collapsedLabels.txt): " << *labelIt << ", assigned class " << newClassNumber << std::endl;	
         m_NumberOfClasses++;
         }
+	else
+	{
+	  // the label is in m_CompressionLookupTable already, find the class assigned to it
+          const int newClassNumber = compressionIt->second[ 0 ];
+	  m_CompressionLookupTable[ collapsedIt->first ].push_back( newClassNumber );
+	  printf("\t                            (collapsed label list): %-5d -> %-3d\n", *labelIt, newClassNumber);
+	}
 
       }  // End loop over all real labels belonging to a certain collapsed label
 
     } // End loop over all collapsed labels  
 
 
-  // Also loop over all pixels of all images, and create a new entry for every new intensity encountered
+  // Also loop over all pixels of all images, and create a new entry for any labels not in m_CompressionLookupTable yet
   for ( std::vector< ImageType::ConstPointer >::const_iterator it = images.begin();
          it != images.end(); ++it )
     {
@@ -157,10 +169,11 @@ ::Construct( const std::vector< ImageType::ConstPointer >& images )
       {
       if ( m_CompressionLookupTable.find( voxelIt.Get() ) == m_CompressionLookupTable.end() )
         {
-        std::cout << "Encountered new real label " << voxelIt.Get() << std::endl;
-
+        // label not in m_CompressionLookupTable yet, add it
         const int  newClassNumber = m_NumberOfClasses;
         m_CompressionLookupTable[ voxelIt.Get() ] = std::vector< int >( 1, newClassNumber );
+	printf("\t Encountered new real label (image) %-5d -> %-3d\n" , voxelIt.Get(), newClassNumber);
+        //std::cout << "Encountered new real label (image) " << voxelIt.Get() << ", assigned class " << newClassNumber << std::endl;	
         m_NumberOfClasses++;
         }  
       } // End loop over all pixels

gems/kvlCompressionLookupTable.h

@@ -49,6 +49,16 @@ public :
     {
     return m_CompressionLookupTable.find( label )->second;
     }
+
+  const std::vector<ImageType::PixelType> GetLabels() const
+    {
+      std::vector<ImageType::PixelType> labels;
+      CompressionLookupTableType::const_iterator it;
+      for (it = m_CompressionLookupTable.begin(); it != m_CompressionLookupTable.end(); ++it)
+        labels.push_back(it->first);
+
+      return labels;
+    }
 
   const ColorType&  GetColor( int classNumber ) const
     {

samseg/Samseg.py

@@ -502,15 +502,50 @@ def writeResults(self, biasFields, posteriors):
         volumeOfOneVoxel = np.abs(np.linalg.det(exampleImage.transform_matrix.as_numpy_array[:3, :3]))
         volumesInCubicMm = np.sum(posteriors, axis=0) * volumeOfOneVoxel
 
+        # Write intracranial volume
+        sbtiv = icv(zip(*[names, volumesInCubicMm]))
+        with open(os.path.join(self.savePath, 'sbtiv.stats'), 'w') as fid:
+            fid.write('# Measure Intra-Cranial, %.6f, mm^3\n' % sbtiv)
+
         # Write structural volumes
         with open(os.path.join(self.savePath, 'samseg.stats'), 'w') as fid:
+            fid.write('# Measure %s, %.6f, mm^3\n' % ('Intra-Cranial', sbtiv))
             for volume, name in zip(volumesInCubicMm, names):
                 fid.write('# Measure %s, %.6f, mm^3\n' % (name, volume))
 
-        # Write intracranial volume
-        sbtiv = icv(zip(*[names, volumesInCubicMm]))
-        with open(os.path.join(self.savePath, 'sbtiv.stats'), 'w') as fid:
-            fid.write('# Measure Intra-Cranial, %.6f, mm^3\n' % sbtiv)
+        # Write structural volumes in a csv
+        with open(os.path.join(self.savePath, 'samseg.csv'), 'w') as fid:
+            fid.write('ROI,volume_mm3,volume_ICV_x1000\n');
+            fid.write('%s,%.6f,1000\n' % ('Intra-Cranial', sbtiv))
+            for volume, name in zip(volumesInCubicMm, names):
+                fid.write('%s,%.6f,%.6f\n' % (name, volume,1000*volume/sbtiv))
+
+        # Write out a freesurfer-style stats file (good for use with asegstats2table)
+        with open(os.path.join(self.savePath, 'samseg.fs.stats'), 'w') as fid:
+            fid.write('# Measure EstimatedTotalIntraCranialVol, eTIV, Estimated Total Intracranial Volume, %0.6f, mm^3\n'%(sbtiv));
+            # Could add other measures here like total brain volume
+            fid.write('# ColHeaders  Index SegId NVoxels Volume_mm3 StructName\n');
+            k = 0;
+            voxsize = self.voxelSpacing[0]*self.voxelSpacing[1]*self.voxelSpacing[2];
+            # Sort them by seg index like with the aseg.stats. Exclude Unknown, but keep WM and Cortex
+            seglist = [];
+            k = 0;
+            for volume, name in zip(volumesInCubicMm, names):
+                if(name == 'Unknown'): 
+                    k = k+1;
+                    continue
+                idx = fslabels[k];
+                seglist.append([idx,volume,name]);
+                k = k+1;
+            seglist.sort()
+            k = 0;
+            for seg in seglist:
+                idx = seg[0];
+                volume = seg[1];
+                name = seg[2];
+                nvox = round(volume/voxsize,2);
+                fid.write('%3d %4d %7d %9.1f %s\n' % (k+1,idx,nvox,volume,name));
+                k = k+1;
 
         return volumesInCubicMm