+ gdcm::Reader hreader;
+ hreader.SetFileName(m_ArgsInfo.inputModelFilename_arg);
+ hreader.Read();
+ gdcm::DataSet& ds = hreader.GetFile().GetDataSet();
+ int energyNumber, headNumber, rotationNumber; //Read the number of energy, of head and rotation steps
+ gdcm::Attribute<0x54,0x11> energyNumber_att;
+ energyNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
+ energyNumber = energyNumber_att.GetValue();
+
+ gdcm::Attribute<0x54,0x21> headNumber_att;
+ headNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet());
+ headNumber = headNumber_att.GetValue();
+
+ gdcm::Tag squenceTag(0x54,0x52); //rotation, read a sequence first
+ const gdcm::DataElement &sequence = hreader.GetFile().GetDataSet().GetDataElement( squenceTag );
+ gdcm::DataSet &sequenceDataSet = sequence.GetValueAsSQ()->GetItem(1).GetNestedDataSet();
+ gdcm::Tag rotationNumber_Tag(0x54,0x53);
+ const gdcm::DataElement &rotationNumber_Element = sequenceDataSet.GetDataElement( rotationNumber_Tag );
+ gdcm::Attribute<0x54,0x53> rotationNumber_att;
+ rotationNumber_att.SetFromDataElement(rotationNumber_Element);
+ rotationNumber = rotationNumber_att.GetValue();
+
+energyNumber = 3;
+headNumber = 4;
+rotationNumber = 6;
+ // Create a new mhd image with the dicom order slices
+ typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New();
+ mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion());
+ mhdCorrectOrder->Allocate();
+ unsigned int zAxis(0); //z value for the input mhd image
+ for (unsigned int energy = 0; energy < energyNumber; ++energy) {
+ for (unsigned int head = 0; head < headNumber; ++head) {
+ for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) {
+ std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl;
+
+ typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder
+ startIteratorIndexCorrectOrder[0] = 0;
+ startIteratorIndexCorrectOrder[1] = 0;
+ startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber;
+
+ typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd
+ startIteratorIndexOriginalOrder[0] = 0;
+ startIteratorIndexOriginalOrder[1] = 0;
+ startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber;
+
+ typename InputImageType::SizeType regionSizeIterator;
+ regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0];
+ regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1];
+ regionSizeIterator[2] = 1;
+
+ typename InputImageType::RegionType regionIteratorCorrectOrder;
+ regionIteratorCorrectOrder.SetSize(regionSizeIterator);
+ regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder);
+
+ typename InputImageType::RegionType regionIteratorOriginalOrder;
+ regionIteratorOriginalOrder.SetSize(regionSizeIterator);
+ regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder);
+
+ itk::ImageRegionIterator<InputImageType> CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder);
+ itk::ImageRegionIterator<InputImageType> OriginalOrderIterator(input,regionIteratorOriginalOrder);
+ while(!CorrectOrderIterator.IsAtEnd()) {
+ CorrectOrderIterator.Set(OriginalOrderIterator.Get());
+ ++CorrectOrderIterator;
+ ++OriginalOrderIterator;
+ }
+
+ ++zAxis;
+ }
+ }
+ }