[ITK-users] Loading a specifing region of an image from disk

[Pol Monsó Purtí]

2014-11-04 15:31 GMT+01:00 Bradley Lowekamp <blowekamp at mail.nih.gov>:

> Discussion inline below:
>
> On Nov 3, 2014, at 5:58 PM, Pol Monsó Purtí <lluna.nova at gmail.com> wrote:
>
> I've read the article early this afternoon and I found it very good and
> encouraging. Specially *"The other difference is that the reader fully
> supports streaming and only reads the required region from the file."*
>
> That is exactly what we need. I was just wondering wether it splitted the
> volume and took the chunks that intersected with the requested one or if it
> really just loaded the requested one. I replicated that code taking out the
> processing and just having a reader and a writer, but I didn't know how to
> test that indeed it was just loading that requested piece.
>
>
> For the ImageIO file formats I wrote which are just headers with binary
> data, they can read an arbitrary ImageIORegion and don't read any more.
> Practically however, its extremely in efficient to not read complete
> contiguous scanlines as there is to much random access and it slows down.
> For Tiff ( similar to some JPEG2000 formats), the image is encoded as
> strips or tiles. The underlying libtiff library can only easily decode
> these entire chunks, so the plan is that the requested region will be
> expanded by the image io to the union of the strips/tiles needed to read
> the requested region.
>
>
>
>> > 2.- Will the pipeline only load the requested region, including the
>> reader?
>>
>> The pipeline has a series of steps which include propagating a requested
>> region, and giving each filter the opportunity to enlarge this region. I
>> have frequently assemble pipelines which perform out-of-core processing as
>> you desire. However all filters in the pipeline must support streaming, and
>> not enlarging the requested region to the largest possible region.
>>
>>
> Yes, I've seen that not many support streaming. But in our case we would
> just need a reader that can read a specific region.
>
>
>> Additionally the ImageIO for the reader and writer must also support
>> streaming. There are several ImageIO's which support streaming reading
>> (MRC,MetaIO, VTK, Nift, JPEG2000, etc), I believe that only MRC, VTK and
>> MetaIO  fully support it for writing. So TIFF is currently not there.
>>
>
> Yeah, we were thinking on using HDF5, but maybe you have some idea on what
> would be best. That's another topic thought.
>
>
> I find the the header + raw image data are rather efficient for random
> access of the whole image. These files are archived with compression, and
> decompressed when needed.
>

Yes indeed.


>
>
>> Recently I have been working on improving the itkTIFFImageIO [2], and
>> have performance gains of upto 3-5X from the prior version. However, I
>> haven't gotten to adding streaming yet. It's is a back burner project for
>> me to add support for stream reading to it, so that may be coming shortly (
>> this poking helps ). However, adding stream and paste writing is not
>> planned.
>>
>>
> Good work! 3-5x! That might be of interest actually for something else we
> have. Can we pull your code?
>
>
> Its in ITK master please test :)
>

Will do.


>
>> For this case you should add in an ExtractImageFilter, or a
>> CropImageFilter.
>>
>>
> But if I use Extract or Crop, the reader is still loading the full image
> to disk, isn't it?
>
>
> Nope. The only the sub-region is requested to the reader, and if the
> reader supports streaming that region will be read.
>
> I'd suggest reading section 13.3 of the ITK Software guide a couple times.
> And then assemble the pipeline Reader->Extract, and printing the output of
> each filter for each step in Extract->UpdateOutputInformation(),
> Extract->PropagateRequestedRegion() and then Extract->UpdateOutputData().
> This will help you to under stand the step in the pipeline execution.
>
>
So do readers (hdf5 for example) transparently do streaming? Or do I have
to set the reader imageIO so that it know that it has to stream?

Section 13.3? There are 9 sections, which one do you refer to? I think I
get what you mean and I'll split up the process to undrstand it better as
you suggest.

Meanwhile, below I paste a piece of code more close to what we would like
to do. The code results in exception because the ImageIO starts at index
[0,0,0] while piece1 starts at 51,80,5 and therefore the Extract fails. I
didn't understand how to set the index on a ImageIO. Then I saw that in
your article you put the whole ImageIO and then you just select one slice
with coronalSlice (which would be my piece1/2). I tried that too, although
I wonder what's the difference with not specifying ImageIO and read
normally. setting imageio1/2 to the full image failes on the reader, saying
that it expected  24300094 bytes (that's the size of the image, so it's
trying to load all I guess), but read 411866 bytes.

How can I do this right? And how would the code change if I were to use
HDF5?

the code:

int main( int argc, char* argv[] )
{

   const unsigned int Dimension = 3;

  typedef unsigned char              PixelType;
  typedef itk::Image< PixelType, Dimension > ImageType;

  typedef itk::ImageFileReader< ImageType > ReaderFilterType;
  ReaderFilterType::Pointer reader1 = ReaderFilterType::New();
  reader1->SetFileName("../data/image.raw");

  //size is not known until we load the image
//  ImageType::RegionType largest =
reader1->GetOutput()->GetLargestPossibleRegion();
//  ImageType::SizeType size   = largest.GetSize();
//  ImageType::IndexType index = largest.GetIndex();

  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size[0] = 511;
  size[1] = 806;
  size[2] = 59;

  ImageType::RegionType largest(start,size);

  ReaderFilterType::Pointer reader2 = ReaderFilterType::New();
  reader2->SetFileName("../data/image.raw");

  ImageType::RegionType piece1;
  piece1.SetIndex(0,start[0] + 0.1*size[0]);
  piece1.SetIndex(1,start[1] + 0.1*size[1]);
  piece1.SetIndex(2,start[2] + 0.1*size[2]);

  piece1.SetSize(0,0.1*size[0]);
  piece1.SetSize(1,0.1*size[1]);
  piece1.SetSize(2,0.1*size[2]);

  ImageType::RegionType piece2;
  piece2.SetIndex(0,start[0] + 0.5*size[0]);
  piece2.SetIndex(1,start[1] + 0.5*size[1]);
  piece2.SetIndex(2,start[2] + 0.5*size[2]);

  piece2.SetSize(piece1.GetSize());

  std::cout << "largest: " << largest << std::endl;
  std::cout << "piece1: "  << piece1  << std::endl;
  std::cout << "piece2: "  << piece2  << std::endl;

  typedef itk::RawImageIO <PixelType, Dimension> ImageIOType;
  ImageIOType::Pointer imageio1 = ImageIOType::New();
  ImageIOType::Pointer imageio2 = ImageIOType::New();

  imageio1->SetDimensions(0,piece1.GetSize()[0]);
  imageio1->SetDimensions(1,piece1.GetSize()[1]);
  imageio1->SetDimensions(2,piece1.GetSize()[2]);


imageio1->SetHeaderSize(piece1.GetSize()[0]*piece1.GetSize()[1]*piece1.GetSize()[2]);

  imageio2->SetDimensions(0,piece2.GetSize()[0]);
  imageio2->SetDimensions(1,piece2.GetSize()[1]);
  imageio2->SetDimensions(2,piece2.GetSize()[2]);


imageio2->SetHeaderSize(piece2.GetSize()[0]*piece2.GetSize()[1]*piece2.GetSize()[2]);

  reader1->SetImageIO(imageio1);
  reader2->SetImageIO(imageio2);

  //this should now work
  reader1->UpdateOutputInformation();
  reader2->UpdateOutputInformation();

  std::cout << "reader largest " <<
reader1->GetOutput()->GetLargestPossibleRegion() << std::endl;

  typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractFilterType;
  ExtractFilterType::Pointer extract1 = ExtractFilterType::New();
  ExtractFilterType::Pointer extract2 = ExtractFilterType::New();
  extract1->SetExtractionRegion(piece1);
  extract2->SetExtractionRegion(piece2);
  extract1->SetInput(reader1->GetOutput());
  extract2->SetInput(reader2->GetOutput());
#if ITK_VERSION_MAJOR >= 4
  extract1->SetDirectionCollapseToIdentity();
  extract2->SetDirectionCollapseToIdentity();
#endif

  try
  {
      extract1->Update();
      extract2->Update();
  }
  catch( itk::ExceptionObject & error )
  {
      std::cerr << "Extract error: " << error << std::endl;
      return EXIT_FAILURE;
  }

  typedef itk::SubtractImageFilter <ImageType, ImageType >
     SubtractImageFilterType;

   SubtractImageFilterType::Pointer subtractFilter
     = SubtractImageFilterType::New ();
   subtractFilter->SetInput1(extract1->GetOutput());
   subtractFilter->SetInput2(extract2->GetOutput());


  typedef itk::ImageFileWriter< ImageType > WriterFilterType;
  WriterFilterType::Pointer writer = WriterFilterType::New();
  writer->SetFileName("image_out.raw");
  writer->SetInput(subtractFilter->GetOutput());

  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & error )
    {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
    }

  return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20141104/7d61f278/attachment.html>