WIP: adding ground truth creator. Closes InsightSoftwareConsortium#108.

test/CMakeLists.txt

@@ -8,6 +8,7 @@ set(MontageTests
   itkMontagePCMTestFiles.cxx
   itkMontageGenericTests.cxx
   itkMontageTest2D.cxx
+  itkMontageTruthCreator.cxx
   )
 
 CreateTestDriver(Montage "${Montage-Test_LIBRARIES}" "${MontageTests}")

test/itkMontageTest2D.cxx

@@ -25,7 +25,7 @@ int itkMontageTest2D(int argc, char* argv[])
 {
   if ( argc < 4 )
     {
-    std::cerr << "Usage: " << argv[0] << " <directoryWtihInputData> <montageTSV> <mockTSV>";
+    std::cerr << "Usage: " << argv[0] << " <directoryWithInputData> <montageTSV> <mockTSV>";
     std::cerr << " [ varyPaddingMethods peakMethod loadIntoMemory streamSubdivisions doPairs";
     std::cerr << " writeTransforms allowDrift positionTolerance writeImage ]" << std::endl;
     return EXIT_FAILURE;

test/itkMontageTruthCreator.cxx

@@ -0,0 +1,212 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  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 "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkResampleImageFilter.h"
+#include "itkRGBPixel.h"
+#include "itkTileMontage.h"
+
+template< unsigned Dimension >
+struct Tile
+{
+  using PointType = itk::ContinuousIndex< double, Dimension >;
+  PointType Position; // x, y... coordinates
+  std::string FileName;
+};
+
+// calculate nD-index given linear index and montage size
+template< unsigned Dimension >
+itk::Size< Dimension >
+linearToNDindex( itk::SizeValueType linearIndex, const std::vector< unsigned >& montageSize )
+{
+  assert( montageSize.size() == Dimension );
+  itk::Size< Dimension > ind;
+  itk::SizeValueType stride = 1u;
+  for ( unsigned d = 0; d < Dimension; d++ )
+    {
+    stride *= montageSize[d];
+    ind[d] = linearIndex % stride;
+    linearIndex /= stride;
+    }
+  return ind;
+}
+
+template< typename PixelType, unsigned Dimension >
+int
+CreateGroundTruth2D( char* inFilename, const std::vector< unsigned >& montageSize, const std::vector< double >& overlap,
+                     const std::string& outDir )
+{
+  using ImageType = itk::Image< PixelType, Dimension >;
+  using ReaderType = itk::ImageFileReader< ImageType >;
+  typename ReaderType::Pointer reader = ReaderType::New();
+  reader->SetFileName( inFilename );
+  reader->Update();
+  typename ImageType::Pointer inImage = reader->GetOutput();
+  typename ImageType::SizeType size = inImage->GetLargestPossibleRegion().GetSize();
+  typename ImageType::SizeType tileSize;
+
+  // calculate tile sizes in index space
+  std::vector< double > tileSizes( Dimension );
+  std::vector< double > slopeFactor{ 0.3, 0.2, 0.1 }; // consistent variation as fraction of overlap
+  std::vector< double > slope( Dimension );
+  std::vector< double > intercept( Dimension );
+  itk::SizeValueType totalTiles = 1;
+  for ( unsigned d = 0; d < Dimension; d++ )
+    {
+    tileSizes[d] = size[d] / ( montageSize[d] + overlap[d] - montageSize[d] * overlap[d] );
+    tileSize[d] = std::ceil( tileSizes[d] );
+    double slopeMax = tileSizes[d] * overlap[d] / montageSize[d];
+    slope[d] = slopeFactor[0] * slopeMax;
+    intercept[d] = -0.5 * slopeFactor[0] * tileSizes[d] * overlap[d];
+    totalTiles *= montageSize[d];
+    }
+
+  using ScalarImage = itk::Image< typename itk::NumericTraits< PixelType >::ValueType, Dimension >;
+  using MontageType = itk::TileMontage< ScalarImage, double >;
+  std::ofstream f( outDir + "TileConfiguration.txt" );
+  f << "# Tile coordinates are in index space, not physical space\n";
+  f << "dim = " << Dimension << "\n\n";
+  std::ofstream fr( outDir + "TileConfiguration.registered.txt" );
+  fr << "# Tile coordinates are in index space, not physical space\n";
+  fr << "dim = " << Dimension << "\n\n";
+
+  // calculate tile positions
+  std::vector< Tile< Dimension > > tiles( totalTiles );
+  for ( itk::SizeValueType linearIndex = 0; linearIndex < totalTiles; linearIndex++ )
+    {
+    typename MontageType::TileIndexType ind = linearToNDindex< Dimension >( linearIndex, montageSize );
+    tiles[linearIndex].FileName = "tile_" + std::to_string( linearIndex ) + ".nrrd";
+    f << tiles[linearIndex].FileName << ";;(";
+    fr << tiles[linearIndex].FileName << ";;(";
+    for ( unsigned d = 0; d < Dimension; d++ )
+      {
+      tiles[linearIndex].Position[d] = ind[d] * ( ( 1 - overlap[d] ) * tileSizes[d] + slope[d] ) + intercept[d];
+      if ( d > 0 )
+        {
+        f << ", ";
+        fr << ", ";
+        }
+      f << ind[d] * ( 1 - overlap[d] ) * tileSizes[d];
+      fr << tiles[linearIndex].Position[d];
+      }
+    f << ')' << std::endl;
+    fr << ')' << std::endl;
+    }
+  f.close();
+  fr.close();
+
+  // now resample the tiles
+  using ResampleFilterType = itk::ResampleImageFilter< ImageType, ImageType >;
+  typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
+  resampleFilter->SetReferenceImage( inImage );
+  resampleFilter->SetInput( inImage );
+  resampleFilter->SetSize( tileSize );
+  // resampleFilter->SetTransform( scaleTransform ); // identity
+  // resampleFilter->SetInterpolator( interpolator );
+  using WriterType = itk::ImageFileWriter< ImageType >;
+  typename WriterType::Pointer writer = WriterType::New();
+  writer->SetInput( resampleFilter->GetOutput() );
+
+  for ( itk::SizeValueType linearIndex = 0; linearIndex < totalTiles; linearIndex++ )
+    {
+    typename ImageType::PointType p;
+    inImage->TransformContinuousIndexToPhysicalPoint( tiles[linearIndex].Position, p );
+    resampleFilter->SetOutputOrigin( p );
+    std::cout << "Resampling " << tiles[linearIndex].FileName << std::flush;
+    resampleFilter->Update();
+
+    std::cout << " Writing " << std::flush;
+    writer->SetFileName( outDir + tiles[linearIndex].FileName );
+    writer->Update();
+    std::cout << " Done" << std::endl;
+    }
+
+
+  return EXIT_SUCCESS;
+}
+
+template< typename PixelType >
+int CreateGroundTruth( int argc, char* argv[], unsigned dimension )
+{
+  std::string outputPath;
+  if ( argc > 2 + 2 * int( dimension ) )
+    {
+    outputPath = argv[2 + 2 * dimension];
+    if ( outputPath.back() != '/' && outputPath.back() != '\\' )
+      {
+      outputPath += '/';
+      }
+    }
+
+  std::vector< unsigned > montageSize( dimension );
+  std::vector< double > overlap( dimension );
+  for ( unsigned i = 0; i < dimension; i++ )
+    {
+    montageSize[i] = std::stoul( argv[2 + i] );
+    overlap[i] = std::stod( argv[2 + dimension + i] ) / 100.0;
+    }
+
+  if ( dimension == 2 )
+    {
+    return CreateGroundTruth2D< PixelType, 2 >( argv[1], montageSize, overlap, outputPath );
+    }
+  else if ( dimension == 3 )
+    {
+    return CreateGroundTruth2D< PixelType, 3 >( argv[1], montageSize, overlap, outputPath );
+    }
+  else
+    {
+    std::cerr << "Only dimensions 2 and 3 are supported!" << std::endl;
+    return EXIT_FAILURE;
+    }
+}
+
+int itkMontageTruthCreator( int argc, char* argv[] )
+{
+  if ( argc < 6 )
+    {
+    std::cerr << "Usage: " << argv[0] << " <inputImage> <montageSize> <overlap%>";
+    std::cerr << " [outputDirectory]" << std::endl;
+    std::cerr << "Example:" << argv[0] << " slice2D.png 8 6 15.0 10.0 ./outGT/" << std::endl;
+    return EXIT_FAILURE;
+    }
+
+  itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO( argv[1], itk::ImageIOFactory::ReadMode );
+  imageIO->SetFileName( argv[1] );
+  imageIO->ReadImageInformation();
+  unsigned dim = imageIO->GetNumberOfDimensions();
+  const itk::ImageIOBase::IOPixelType pixelType = imageIO->GetPixelType();
+
+  if ( dim == 3 && argc < 8 )
+    {
+    std::cerr << "Usage: " << argv[0] << " <inputImage> <montageSize> <overlap%>";
+    std::cerr << " [outputDirectory]" << std::endl;
+    std::cerr << "Example:" << argv[0] << " volume3D.nrrd 8 6 4 15.0 10.0 5.0 ./outGT/" << std::endl;
+    return EXIT_FAILURE;
+    }
+
+  if ( pixelType == itk::ImageIOBase::IOPixelType::RGB )
+    {
+    return CreateGroundTruth< itk::RGBPixel< unsigned char > >( argc, argv, dim );
+    }
+  else
+    {
+    return CreateGroundTruth< unsigned short >( argc, argv, dim );
+    }
+}