1 #include <tclap/CmdLine.h> 4 #include <itkTransformFileReader.h> 5 #include <itkAddImageFilter.h> 6 #include <itkImageRegionIterator.h> 9 int main(
int ac,
const char** av)
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";
16 TCLAP::CmdLine cmd(descriptionMessage,
' ',ANIMA_VERSION);
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);
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);
31 catch (TCLAP::ArgException& e)
33 std::cerr <<
"Error: " << e.error() <<
"for argument " << e.argId() << std::endl;
38 typedef itk::Image <double, 3> ImageType;
39 ImageType::Pointer geometryImage = anima::readImage <ImageType> (geomArg.getValue());
41 typedef itk::MatrixOffsetTransformBase <double,3> MatrixTransformType;
42 typedef MatrixTransformType::Pointer MatrixTransformPointer;
43 typedef itk::ContinuousIndex <double,3> ContinuousIndexType;
44 ContinuousIndexType lowerCornerBoundingBox, upperCornerBoundingBox;
46 for (
unsigned int i = 0;i < ImageType::ImageDimension;++i)
48 lowerCornerBoundingBox[i] = geometryImage->GetLargestPossibleRegion().GetIndex()[i];
49 upperCornerBoundingBox[i] = lowerCornerBoundingBox[i] + geometryImage->GetLargestPossibleRegion().GetSize()[i];
52 std::vector <std::string> inputImages = inArg.getValue();
53 std::vector <std::string> inputTrsfs = trArg.getValue();
56 for (
unsigned int k = 0;k < inputImages.size();++k)
58 ImageType::Pointer input = anima::readImage <ImageType> (inputImages[k]);
60 ImageType::IndexType lowerCornerIndex = input->GetLargestPossibleRegion().GetIndex();
61 ImageType::IndexType upperCornerIndex = lowerCornerIndex + input->GetLargestPossibleRegion().GetSize();
63 MatrixTransformPointer trsf = MatrixTransformType::New();
66 if (inputTrsfs.size() != 0)
68 itk::TransformFileReader::Pointer reader = itk::TransformFileReader::New();
69 reader->SetFileName(inputTrsfs[k]);
72 const itk::TransformFileReader::TransformListType *trsfList = reader->GetTransformList();
73 itk::TransformFileReader::TransformListType::const_iterator tr_it = trsfList->begin();
75 trsf = dynamic_cast <MatrixTransformType *> ((*tr_it).GetPointer());
76 MatrixTransformPointer tmpInvert = MatrixTransformType::New();
77 trsf->GetInverse(tmpInvert);
81 unsigned int numCorners = 1 << ImageType::ImageDimension;
82 for (
unsigned int i = 0;i < numCorners;++i)
84 ImageType::IndexType corner;
85 unsigned int upper = i;
87 for(
unsigned int dim = 0; dim < ImageType::ImageDimension; ++dim)
90 corner[dim] = upperCornerIndex[dim];
92 corner[dim] = lowerCornerIndex[dim];
97 ImageType::PointType tmpPoint;
98 input->TransformIndexToPhysicalPoint(corner,tmpPoint);
99 tmpPoint = trsf->TransformPoint(tmpPoint);
101 itk::ContinuousIndex <double,3> tmpRes;
102 geometryImage->TransformPhysicalPointToContinuousIndex(tmpPoint,tmpRes);
103 for (
unsigned int j = 0;j < ImageType::ImageDimension;++j)
105 if (tmpRes[j] < lowerCornerBoundingBox[j])
106 lowerCornerBoundingBox[j] = tmpRes[j];
108 if (tmpRes[j] > upperCornerBoundingBox[j])
109 upperCornerBoundingBox[j] = tmpRes[j];
115 ImageType::Pointer outputImage = ImageType::New();
116 outputImage->Initialize();
118 for (
unsigned int i = 0;i < ImageType::ImageDimension;++i)
120 lowerCornerBoundingBox[i] = std::floor(lowerCornerBoundingBox[i]);
121 upperCornerBoundingBox[i] = std::ceil(upperCornerBoundingBox[i]);
124 ImageType::PointType origin;
125 geometryImage->TransformContinuousIndexToPhysicalPoint(lowerCornerBoundingBox,origin);
126 outputImage->SetOrigin(origin);
127 outputImage->SetSpacing(geometryImage->GetSpacing());
128 outputImage->SetDirection(geometryImage->GetDirection());
130 ImageType::RegionType outputRegion;
131 for (
unsigned int i = 0;i < ImageType::ImageDimension;++i)
133 int lowerIndex = lowerCornerBoundingBox[i];
134 int upperIndex = upperCornerBoundingBox[i];
137 outputRegion.SetIndex(i,lowerIndex);
139 outputRegion.SetIndex(i,0);
141 outputRegion.SetSize(i,upperIndex - lowerIndex);
144 outputImage->SetRegions(outputRegion);
145 outputImage->Allocate();
146 outputImage->FillBuffer(0.0);
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);
161 typedef itk::AddImageFilter <ImageType,ImageType,ImageType> AddFilterType;
162 typedef itk::ImageRegionIterator <ImageType> ImageIteratorType;
163 typedef itk::ImageRegionIterator <MaskImageType> MaskImageIteratorType;
165 for (
unsigned int i = 0;i < inputImages.size();++i)
167 std::cout <<
"Mosaicing image " << i+1 <<
": " << inputImages[i];
168 if (inputTrsfs.size() > 0)
169 std::cout <<
" with transform " << inputTrsfs[i] << std::endl;
171 std::cout << std::endl;
173 ImageType::Pointer input = anima::readImage <ImageType> (inputImages[i]);
175 MatrixTransformPointer trsf = MatrixTransformType::New();
178 if (inputTrsfs.size() != 0)
180 itk::TransformFileReader::Pointer reader = itk::TransformFileReader::New();
181 reader->SetFileName(inputTrsfs[i]);
184 const itk::TransformFileReader::TransformListType *trsfList = reader->GetTransformList();
185 itk::TransformFileReader::TransformListType::const_iterator tr_it = trsfList->begin();
187 trsf = dynamic_cast <MatrixTransformType *> ((*tr_it).GetPointer());
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();
200 AddFilterType::Pointer addFilter = AddFilterType::New();
201 addFilter->SetInput1(outputImage);
202 addFilter->SetInput2(imageResampler->GetOutput());
203 addFilter->SetNumberOfWorkUnits(nbpArg.getValue());
206 outputImage = addFilter->GetOutput();
207 outputImage->DisconnectPipeline();
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());
216 tmpMask->FillBuffer(1);
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();
228 MaskImageIteratorType outputMaskItr(maskImage,maskImage->GetLargestPossibleRegion());
229 ImageIteratorType trsfMaskItr(maskImResampler->GetOutput(),maskImage->GetLargestPossibleRegion());
231 while (!outputMaskItr.IsAtEnd())
233 if (trsfMaskItr.Get() > 0)
234 outputMaskItr.Set(outputMaskItr.Get() + 1);
242 ImageIteratorType outItr(outputImage,outputImage->GetLargestPossibleRegion());
243 MaskImageIteratorType maskItr(maskImage,maskImage->GetLargestPossibleRegion());
244 while (!outItr.IsAtEnd())
246 if (maskItr.Get() != 0)
247 outItr.Set(outItr.Get() / maskItr.Get());
255 anima::writeImage <ImageType> (outArg.getValue(),outputImage);
257 if (outMaskArg.getValue() !=
"")
258 anima::writeImage <MaskImageType> (outMaskArg.getValue(),maskImage);
int main(int ac, const char **av)