ANIMA  4.0
animaImageMosaicing.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
4 #include <itkTransformFileReader.h>
5 #include <itkAddImageFilter.h>
6 #include <itkImageRegionIterator.h>
8 
9 int main(int ac, const char** av)
10 {
11  std::string descriptionMessage = "Resampler tool for stitching several image into one.\n"
12  "Applies linear transforms associated to the input images, then\n"
13  "generates a larger image containing the mosaic of the transformed input images.\n"
14  "INRIA / IRISA - VisAGeS/Empenn Team";
15 
16  TCLAP::CmdLine cmd(descriptionMessage, ' ',ANIMA_VERSION);
17 
18  TCLAP::MultiArg<std::string> inArg("i","input","Input image (should be used several times)",true,"input image",cmd);
19  TCLAP::MultiArg<std::string> trArg("t","trsf","Transformation (one linear transform for each input)",false,"transformation",cmd);
20  TCLAP::ValueArg<std::string> outArg("o","output","Output mosaic image",true,"","output mosaic image",cmd);
21  TCLAP::ValueArg<std::string> outMaskArg("O","out-mask","Output mosaic mask",false,"","output mosaic mask",cmd);
22  TCLAP::ValueArg<std::string> geomArg("g","geometry","Geometry image",true,"","geometry image",cmd);
23 
24  TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","Number of threads to run on (default : all cores)",
25  false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
26 
27  try
28  {
29  cmd.parse(ac,av);
30  }
31  catch (TCLAP::ArgException& e)
32  {
33  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
34  return EXIT_FAILURE;
35  }
36 
37  // Find out the type of the image in file
38  typedef itk::Image <double, 3> ImageType;
39  ImageType::Pointer geometryImage = anima::readImage <ImageType> (geomArg.getValue());
40 
41  typedef itk::MatrixOffsetTransformBase <double,3> MatrixTransformType;
42  typedef MatrixTransformType::Pointer MatrixTransformPointer;
43  typedef itk::ContinuousIndex <double,3> ContinuousIndexType;
44  ContinuousIndexType lowerCornerBoundingBox, upperCornerBoundingBox;
45 
46  for (unsigned int i = 0;i < ImageType::ImageDimension;++i)
47  {
48  lowerCornerBoundingBox[i] = geometryImage->GetLargestPossibleRegion().GetIndex()[i];
49  upperCornerBoundingBox[i] = lowerCornerBoundingBox[i] + geometryImage->GetLargestPossibleRegion().GetSize()[i];
50  }
51 
52  std::vector <std::string> inputImages = inArg.getValue();
53  std::vector <std::string> inputTrsfs = trArg.getValue();
54 
55  // Explore images and transforms to get the boundaries of the future image
56  for (unsigned int k = 0;k < inputImages.size();++k)
57  {
58  ImageType::Pointer input = anima::readImage <ImageType> (inputImages[k]);
59 
60  ImageType::IndexType lowerCornerIndex = input->GetLargestPossibleRegion().GetIndex();
61  ImageType::IndexType upperCornerIndex = lowerCornerIndex + input->GetLargestPossibleRegion().GetSize();
62 
63  MatrixTransformPointer trsf = MatrixTransformType::New();
64  trsf->SetIdentity();
65 
66  if (inputTrsfs.size() != 0)
67  {
68  itk::TransformFileReader::Pointer reader = itk::TransformFileReader::New();
69  reader->SetFileName(inputTrsfs[k]);
70  reader->Update();
71 
72  const itk::TransformFileReader::TransformListType *trsfList = reader->GetTransformList();
73  itk::TransformFileReader::TransformListType::const_iterator tr_it = trsfList->begin();
74 
75  trsf = dynamic_cast <MatrixTransformType *> ((*tr_it).GetPointer());
76  MatrixTransformPointer tmpInvert = MatrixTransformType::New();
77  trsf->GetInverse(tmpInvert);
78  trsf = tmpInvert;
79  }
80 
81  unsigned int numCorners = 1 << ImageType::ImageDimension;
82  for (unsigned int i = 0;i < numCorners;++i)
83  {
84  ImageType::IndexType corner;
85  unsigned int upper = i; // each bit indicates upper/lower on dim
86 
87  for(unsigned int dim = 0; dim < ImageType::ImageDimension; ++dim)
88  {
89  if (upper & 1)
90  corner[dim] = upperCornerIndex[dim];
91  else
92  corner[dim] = lowerCornerIndex[dim];
93 
94  upper >>= 1;
95  }
96 
97  ImageType::PointType tmpPoint;
98  input->TransformIndexToPhysicalPoint(corner,tmpPoint);
99  tmpPoint = trsf->TransformPoint(tmpPoint);
100 
101  itk::ContinuousIndex <double,3> tmpRes;
102  geometryImage->TransformPhysicalPointToContinuousIndex(tmpPoint,tmpRes);
103  for (unsigned int j = 0;j < ImageType::ImageDimension;++j)
104  {
105  if (tmpRes[j] < lowerCornerBoundingBox[j])
106  lowerCornerBoundingBox[j] = tmpRes[j];
107 
108  if (tmpRes[j] > upperCornerBoundingBox[j])
109  upperCornerBoundingBox[j] = tmpRes[j];
110  }
111  }
112  }
113 
114  // Now define output geometry
115  ImageType::Pointer outputImage = ImageType::New();
116  outputImage->Initialize();
117 
118  for (unsigned int i = 0;i < ImageType::ImageDimension;++i)
119  {
120  lowerCornerBoundingBox[i] = std::floor(lowerCornerBoundingBox[i]);
121  upperCornerBoundingBox[i] = std::ceil(upperCornerBoundingBox[i]);
122  }
123 
124  ImageType::PointType origin;
125  geometryImage->TransformContinuousIndexToPhysicalPoint(lowerCornerBoundingBox,origin);
126  outputImage->SetOrigin(origin);
127  outputImage->SetSpacing(geometryImage->GetSpacing());
128  outputImage->SetDirection(geometryImage->GetDirection());
129 
130  ImageType::RegionType outputRegion;
131  for (unsigned int i = 0;i < ImageType::ImageDimension;++i)
132  {
133  int lowerIndex = lowerCornerBoundingBox[i];
134  int upperIndex = upperCornerBoundingBox[i];
135 
136  if (lowerIndex >= 0)
137  outputRegion.SetIndex(i,lowerIndex);
138  else
139  outputRegion.SetIndex(i,0);
140 
141  outputRegion.SetSize(i,upperIndex - lowerIndex);
142  }
143 
144  outputImage->SetRegions(outputRegion);
145  outputImage->Allocate();
146  outputImage->FillBuffer(0.0);
147 
148  typedef itk::Image <unsigned short,3> MaskImageType;
149  MaskImageType::Pointer maskImage = MaskImageType::New();
150  maskImage->Initialize();
151  maskImage->SetOrigin(origin);
152  maskImage->SetSpacing(geometryImage->GetSpacing());
153  maskImage->SetDirection(geometryImage->GetDirection());
154  maskImage->SetRegions(outputRegion);
155  maskImage->Allocate();
156  maskImage->FillBuffer(0);
157 
158  // Now finally resample images in new geometry and average
159  typedef anima::ResampleImageFilter <ImageType, ImageType> ResampleFilterType;
160  typedef anima::ResampleImageFilter <MaskImageType, ImageType> MaskResampleFilterType;
161  typedef itk::AddImageFilter <ImageType,ImageType,ImageType> AddFilterType;
162  typedef itk::ImageRegionIterator <ImageType> ImageIteratorType;
163  typedef itk::ImageRegionIterator <MaskImageType> MaskImageIteratorType;
164 
165  for (unsigned int i = 0;i < inputImages.size();++i)
166  {
167  std::cout << "Mosaicing image " << i+1 << ": " << inputImages[i];
168  if (inputTrsfs.size() > 0)
169  std::cout << " with transform " << inputTrsfs[i] << std::endl;
170  else
171  std::cout << std::endl;
172 
173  ImageType::Pointer input = anima::readImage <ImageType> (inputImages[i]);
174 
175  MatrixTransformPointer trsf = MatrixTransformType::New();
176  trsf->SetIdentity();
177 
178  if (inputTrsfs.size() != 0)
179  {
180  itk::TransformFileReader::Pointer reader = itk::TransformFileReader::New();
181  reader->SetFileName(inputTrsfs[i]);
182  reader->Update();
183 
184  const itk::TransformFileReader::TransformListType *trsfList = reader->GetTransformList();
185  itk::TransformFileReader::TransformListType::const_iterator tr_it = trsfList->begin();
186 
187  trsf = dynamic_cast <MatrixTransformType *> ((*tr_it).GetPointer());
188  }
189 
190  ResampleFilterType::Pointer imageResampler = ResampleFilterType::New();
191  imageResampler->SetTransform(trsf);
192  imageResampler->SetSize(outputRegion.GetSize());
193  imageResampler->SetOutputOrigin(origin);
194  imageResampler->SetOutputSpacing(outputImage->GetSpacing());
195  imageResampler->SetOutputDirection(outputImage->GetDirection());
196  imageResampler->SetInput(input);
197  imageResampler->SetNumberOfWorkUnits(nbpArg.getValue());
198  imageResampler->Update();
199 
200  AddFilterType::Pointer addFilter = AddFilterType::New();
201  addFilter->SetInput1(outputImage);
202  addFilter->SetInput2(imageResampler->GetOutput());
203  addFilter->SetNumberOfWorkUnits(nbpArg.getValue());
204  addFilter->Update();
205 
206  outputImage = addFilter->GetOutput();
207  outputImage->DisconnectPipeline();
208 
209  MaskImageType::Pointer tmpMask = MaskImageType::New();
210  tmpMask->Initialize();
211  tmpMask->SetOrigin(input->GetOrigin());
212  tmpMask->SetSpacing(input->GetSpacing());
213  tmpMask->SetDirection(input->GetDirection());
214  tmpMask->SetRegions(input->GetLargestPossibleRegion());
215  tmpMask->Allocate();
216  tmpMask->FillBuffer(1);
217 
218  MaskResampleFilterType::Pointer maskImResampler = MaskResampleFilterType::New();
219  maskImResampler->SetTransform(trsf);
220  maskImResampler->SetSize(outputRegion.GetSize());
221  maskImResampler->SetOutputOrigin(origin);
222  maskImResampler->SetOutputSpacing(outputImage->GetSpacing());
223  maskImResampler->SetOutputDirection(outputImage->GetDirection());
224  maskImResampler->SetInput(tmpMask);
225  maskImResampler->SetNumberOfWorkUnits(nbpArg.getValue());
226  maskImResampler->Update();
227 
228  MaskImageIteratorType outputMaskItr(maskImage,maskImage->GetLargestPossibleRegion());
229  ImageIteratorType trsfMaskItr(maskImResampler->GetOutput(),maskImage->GetLargestPossibleRegion());
230 
231  while (!outputMaskItr.IsAtEnd())
232  {
233  if (trsfMaskItr.Get() > 0)
234  outputMaskItr.Set(outputMaskItr.Get() + 1);
235 
236  ++trsfMaskItr;
237  ++outputMaskItr;
238  }
239  }
240 
241  // Finally divide by the number of images at each pixel
242  ImageIteratorType outItr(outputImage,outputImage->GetLargestPossibleRegion());
243  MaskImageIteratorType maskItr(maskImage,maskImage->GetLargestPossibleRegion());
244  while (!outItr.IsAtEnd())
245  {
246  if (maskItr.Get() != 0)
247  outItr.Set(outItr.Get() / maskItr.Get());
248  else
249  outItr.Set(0);
250 
251  ++outItr;
252  ++maskItr;
253  }
254 
255  anima::writeImage <ImageType> (outArg.getValue(),outputImage);
256 
257  if (outMaskArg.getValue() != "")
258  anima::writeImage <MaskImageType> (outMaskArg.getValue(),maskImage);
259 
260  return EXIT_SUCCESS;
261 }
int main(int ac, const char **av)