<div dir="ltr">Hi,<div><br></div><div>I'm trying to do 3D registration using the Quaternion Classes and for some reason I can't..</div><div><br></div><div>I get an error and I think it is because of the OptimizerScales.. I'm unable to find an example that shows the use of Quaternions transforms so I'm a little lost.</div><div><br></div><div>The error I get is at this line of code within itkQuaternionRigidTransformGradientDescentOptimizer.cxx:</div><div><br></div><div>const unsigned int spaceDimension =  m_CostFunction->GetNumberOfParameters();<br></div><div><br></div><div>#include "itkAffineTransform.h"<br></div><div><div>#include "itkQuaternionRigidTransform.h"</div><div>#include "itkQuaternionRigidTransformGradientDescentOptimizer.h"</div><div>#include "itkCenteredTransformInitializer.h"</div><div>#include "itkMultiResolutionImageRegistrationMethod.h"</div><div>#include "itkMattesMutualInformationImageToImageMetric.h"</div><div>#include "itkRegularStepGradientDescentOptimizer.h"</div><div>#include "itkRecursiveMultiResolutionPyramidImageFilter.h"</div><div>#include "itkImage.h"</div><div>#include "itkImageFileReader.h"</div><div>#include "itkImageFileWriter.h"</div><div>#include "itkResampleImageFilter.h"</div><div>#include "itkCastImageFilter.h"</div><div>#include "itkNormalizeImageFilter.h"</div><div>#include "itkCheckerBoardImageFilter.h"</div><div><br></div><div>#include "itkCommand.h"</div><div>class CommandIterationUpdate : public itk::Command</div><div>{</div><div>public:</div><div><span class="" style="white-space:pre">       </span>typedef  CommandIterationUpdate   Self;</div><div><span class="" style="white-space:pre">  </span>typedef  itk::Command             Superclass;</div><div><span class="" style="white-space:pre">       </span>typedef  itk::SmartPointer<Self>  Pointer;</div><div><span class="" style="white-space:pre"> </span>itkNewMacro(Self);</div><div>protected:</div><div><span class="" style="white-space:pre">        </span>CommandIterationUpdate() : m_CumulativeIterationIndex(0) {};</div><div>public:</div><div><span class="" style="white-space:pre"> </span>typedef   itk::RegularStepGradientDescentOptimizer  OptimizerType;</div><div><span class="" style="white-space:pre">       </span>typedef   const OptimizerType *                     OptimizerPointer;</div><div><span class="" style="white-space:pre">   </span>void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE</div><div><span class="" style="white-space:pre"> </span>{</div><div><span class="" style="white-space:pre">          </span>Execute((const itk::Object *)caller, event);</div><div><span class="" style="white-space:pre">       </span>}</div><div><span class="" style="white-space:pre">  </span>void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE</div><div><span class="" style="white-space:pre">  </span>{</div><div><span class="" style="white-space:pre">          </span>OptimizerPointer optimizer = static_cast<OptimizerPointer>(object);</div><div><span class="" style="white-space:pre">          </span>if (!(itk::IterationEvent().CheckEvent(&event)))</div><div><span class="" style="white-space:pre">               </span>{</div><div><span class="" style="white-space:pre">                  </span>return;</div><div><span class="" style="white-space:pre">            </span>}</div><div><span class="" style="white-space:pre">          </span>std::cout << optimizer->GetCurrentIteration() << "   ";</div><div><span class="" style="white-space:pre">          </span>std::cout << optimizer->GetValue() << "   ";</div><div><span class="" style="white-space:pre">             </span>std::cout << optimizer->GetCurrentPosition() << "  " <<</div><div><span class="" style="white-space:pre">                    </span>m_CumulativeIterationIndex++ << std::endl;</div><div><span class="" style="white-space:pre">   </span>}</div><div>private:</div><div><span class="" style="white-space:pre">   </span>unsigned int m_CumulativeIterationIndex;</div><div>};</div><div><br></div><div>template <typename TRegistration></div><div>class RegistrationInterfaceCommand : public itk::Command</div><div>{</div><div>public:</div><div><span class="" style="white-space:pre">      </span>typedef  RegistrationInterfaceCommand   Self;</div><div><span class="" style="white-space:pre">    </span>typedef  itk::Command                   Superclass;</div><div><span class="" style="white-space:pre">      </span>typedef  itk::SmartPointer<Self>        Pointer;</div><div><span class="" style="white-space:pre">        </span>itkNewMacro(Self);</div><div>protected:</div><div><span class="" style="white-space:pre">        </span>RegistrationInterfaceCommand() {};</div><div>public:</div><div><span class="" style="white-space:pre">   </span>typedef   TRegistration                              RegistrationType;</div><div><span class="" style="white-space:pre">     </span>typedef   RegistrationType *                         RegistrationPointer;</div><div><span class="" style="white-space:pre">     </span>typedef   itk::RegularStepGradientDescentOptimizer   OptimizerType;</div><div><span class="" style="white-space:pre">      </span>typedef   OptimizerType *                            OptimizerPointer;</div><div><span class="" style="white-space:pre">      </span>void Execute(itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE</div><div><span class="" style="white-space:pre">        </span>{</div><div><span class="" style="white-space:pre">          </span>if (!(itk::IterationEvent().CheckEvent(&event)))</div><div><span class="" style="white-space:pre">               </span>{</div><div><span class="" style="white-space:pre">                  </span>return;</div><div><span class="" style="white-space:pre">            </span>}</div><div><span class="" style="white-space:pre">          </span>RegistrationPointer registration = static_cast<RegistrationPointer>(object);</div><div><span class="" style="white-space:pre">         </span>OptimizerPointer optimizer =</div><div><span class="" style="white-space:pre">                       </span>static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());</div><div><span class="" style="white-space:pre">            </span>std::cout << "-------------------------------------" << std::endl;</div><div><span class="" style="white-space:pre">           </span>std::cout << "MultiResolution Level : "</div><div><span class="" style="white-space:pre">                    </span><< registration->GetCurrentLevel() << std::endl;</div><div><span class="" style="white-space:pre">            </span>std::cout << std::endl;</div><div><span class="" style="white-space:pre">              </span>if (registration->GetCurrentLevel() == 0)</div><div><span class="" style="white-space:pre">               </span>{</div><div><span class="" style="white-space:pre">                  </span>optimizer->SetMaximumStepLength(1.00);</div><div><span class="" style="white-space:pre">                  </span>optimizer->SetMinimumStepLength(0.01);</div><div><span class="" style="white-space:pre">          </span>}</div><div><span class="" style="white-space:pre">          </span>else</div><div><span class="" style="white-space:pre">               </span>{</div><div><span class="" style="white-space:pre">                  </span>optimizer->SetMaximumStepLength(optimizer->GetMaximumStepLength() / 4.0);</div><div><span class="" style="white-space:pre">                    </span>optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() / 10.0);</div><div><span class="" style="white-space:pre">           </span>}</div><div><span class="" style="white-space:pre">  </span>}</div><div><span class="" style="white-space:pre">  </span>void Execute(const itk::Object *, const itk::EventObject &) ITK_OVERRIDE</div><div><span class="" style="white-space:pre">       </span>{ return; }</div><div>};</div><div>int main(int argc, char *argv[])</div><div>{</div><div><span class="" style="white-space:pre">        </span>if (argc < 4)</div><div><span class="" style="white-space:pre">   </span>{</div><div><span class="" style="white-space:pre">          </span>std::cerr << "Missing Parameters " << std::endl;</div><div><span class="" style="white-space:pre">             </span>std::cerr << "Usage: " << argv[0];</div><div><span class="" style="white-space:pre">           </span>std::cerr << " imagenfija  imagenflotante ";</div><div><span class="" style="white-space:pre">              </span>std::cerr << " salida";</div><div><span class="" style="white-space:pre">            </span>return EXIT_FAILURE;</div><div><span class="" style="white-space:pre">       </span>}</div><div><span class="" style="white-space:pre">  </span>const    unsigned int    Dimension = 3;</div><div><span class="" style="white-space:pre">        </span>typedef  signed short  PixelType;</div><div><span class="" style="white-space:pre">        </span>typedef itk::Image< PixelType, Dimension >  FixedImageType;</div><div><span class="" style="white-space:pre"> </span>typedef itk::Image< PixelType, Dimension >  MovingImageType;</div><div><span class="" style="white-space:pre">        </span>typedef   float                                    InternalPixelType;</div><div><span class="" style="white-space:pre">   </span>typedef itk::Image< InternalPixelType, Dimension > InternalImageType;</div><div><br></div><div><span class="" style="white-space:pre">       </span>typedef itk::QuaternionRigidTransform< double> TransformType;</div><div><br></div><div><span class="" style="white-space:pre">       </span>typedef itk::QuaternionRigidTransformGradientDescentOptimizer       OptimizerType;</div><div><span class="" style="white-space:pre">      </span>typedef itk::LinearInterpolateImageFunction<</div><div><span class="" style="white-space:pre">            </span>InternalImageType,</div><div><span class="" style="white-space:pre">         </span>double             > InterpolatorType;</div><div><span class="" style="white-space:pre">    </span>typedef itk::MattesMutualInformationImageToImageMetric<</div><div><span class="" style="white-space:pre">         </span>InternalImageType,</div><div><span class="" style="white-space:pre">         </span>InternalImageType >    MetricType;</div><div><span class="" style="white-space:pre">    </span>typedef OptimizerType::ScalesType       OptimizerScalesType;</div><div><span class="" style="white-space:pre">    </span>typedef itk::MultiResolutionImageRegistrationMethod<</div><div><span class="" style="white-space:pre">            </span>InternalImageType,</div><div><span class="" style="white-space:pre">         </span>InternalImageType    > RegistrationType;</div><div><span class="" style="white-space:pre">      </span>typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;</div><div><span class="" style="white-space:pre">       </span>typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</div><div><span class="" style="white-space:pre">       </span>typedef itk::NormalizeImageFilter<FixedImageType, InternalImageType> NormalizeFilterType;</div><div><span class="" style="white-space:pre">    </span>typedef itk::DiscreteGaussianImageFilter<InternalImageType, InternalImageType> GaussianFilterType;</div><div><span class="" style="white-space:pre">   </span>typedef itk::CenteredTransformInitializer<TransformType, FixedImageType, MovingImageType >  TransformInitializerType;</div><div><br></div><div><span class="" style="white-space:pre">      </span>OptimizerType::Pointer      optimizer = OptimizerType::New();</div><div><span class="" style="white-space:pre">   </span>InterpolatorType::Pointer   interpolator = InterpolatorType::New();</div><div><span class="" style="white-space:pre">       </span>RegistrationType::Pointer   registration = RegistrationType::New();</div><div><span class="" style="white-space:pre">       </span>MetricType::Pointer         metric = MetricType::New();</div><div><span class="" style="white-space:pre">        </span>TransformType::Pointer   transform = TransformType::New();</div><div><span class="" style="white-space:pre">        </span>FixedImageReaderType::Pointer  fixedImageReader = FixedImageReaderType::New();</div><div><span class="" style="white-space:pre">    </span>MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</div><div><span class="" style="white-space:pre">   </span>NormalizeFilterType::Pointer fixedNormalizer = NormalizeFilterType::New();</div><div><span class="" style="white-space:pre"> </span>NormalizeFilterType::Pointer movingNormalizer = NormalizeFilterType::New();</div><div><span class="" style="white-space:pre">        </span>GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();</div><div><span class="" style="white-space:pre">     </span>GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();</div><div><span class="" style="white-space:pre">    </span>TransformInitializerType::Pointer initializer = TransformInitializerType::New();</div><div><br></div><div><span class="" style="white-space:pre">  </span>fixedImageReader->SetFileName(argv[1]);</div><div><span class="" style="white-space:pre"> </span>movingImageReader->SetFileName(argv[2]);</div><div><br></div><div><span class="" style="white-space:pre">       </span>registration->SetOptimizer(optimizer);</div><div><span class="" style="white-space:pre">  </span>registration->SetInterpolator(interpolator);</div><div><span class="" style="white-space:pre">    </span>registration->SetMetric(metric);</div><div><span class="" style="white-space:pre">        </span>registration->SetTransform(transform);</div><div><br></div><div><span class="" style="white-space:pre"> </span>fixedNormalizer->SetInput(fixedImageReader->GetOutput());</div><div><span class="" style="white-space:pre">    </span>movingNormalizer->SetInput(movingImageReader->GetOutput());</div><div><br></div><div><span class="" style="white-space:pre"> </span>fixedSmoother->SetVariance(2.0);</div><div><span class="" style="white-space:pre">        </span>movingSmoother->SetVariance(2.0);</div><div><br></div><div><span class="" style="white-space:pre">      </span>fixedSmoother->SetInput(fixedNormalizer->GetOutput());</div><div><span class="" style="white-space:pre">       </span>movingSmoother->SetInput(movingNormalizer->GetOutput());</div><div><br></div><div><span class="" style="white-space:pre">    </span>registration->SetFixedImage(fixedSmoother->GetOutput());</div><div><span class="" style="white-space:pre">     </span>registration->SetMovingImage(movingSmoother->GetOutput());</div><div><br></div><div><span class="" style="white-space:pre">  </span>fixedNormalizer->Update();</div><div><span class="" style="white-space:pre">      </span>registration->SetFixedImageRegion(fixedNormalizer->GetOutput()->GetBufferedRegion());</div><div><br></div><div><span class="" style="white-space:pre">    </span>initializer->SetTransform(transform);</div><div><span class="" style="white-space:pre">   </span>initializer->SetFixedImage(fixedImageReader->GetOutput());</div><div><span class="" style="white-space:pre">   </span>initializer->SetMovingImage(movingImageReader->GetOutput());</div><div><span class="" style="white-space:pre"> </span>initializer->MomentsOn();</div><div><span class="" style="white-space:pre">       </span>initializer->InitializeTransform();</div><div><span class="" style="white-space:pre">     </span>//transform->SetIdentity(); //Esto generó Too Many Samples outside mapping..</div><div><br></div><div><span class="" style="white-space:pre">  </span>registration->SetInitialTransformParameters(transform->GetParameters());</div><div><br></div><div><br></div><div><span class="" style="white-space:pre">   </span>OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());</div><div><br></div><div><span class="" style="white-space:pre">       </span>if (Dimension == 3)</div><div><span class="" style="white-space:pre">        </span>{</div><div><span class="" style="white-space:pre">          </span>optimizerScales[0] = 1.0;</div><div><span class="" style="white-space:pre">          </span>optimizerScales[1] = 1.0;</div><div><span class="" style="white-space:pre">          </span>optimizerScales[2] = 1.0;</div><div><span class="" style="white-space:pre">          </span>optimizerScales[3] = 1.0;</div><div><span class="" style="white-space:pre">          </span>optimizerScales[4] = 0.0001;</div><div><span class="" style="white-space:pre">               </span>optimizerScales[5] = 0.0001;</div><div><span class="" style="white-space:pre">               </span>optimizerScales[6] = 0.0001;</div><div><span class="" style="white-space:pre">               </span>optimizerScales[7] = 0.0001;</div><div><span class="" style="white-space:pre">       </span>}</div><div><br></div><div><span class="" style="white-space:pre"> </span>optimizer->SetScales(optimizerScales);</div><div><span class="" style="white-space:pre">  </span></div><div><span class="" style="white-space:pre">   </span>optimizer->SetNumberOfIterations(20);</div><div><span class="" style="white-space:pre">   </span>optimizer->MaximizeOn();</div><div><br></div><div><span class="" style="white-space:pre">       </span>metric->SetNumberOfSpatialSamples(50);</div><div><br></div><div><span class="" style="white-space:pre"> </span>CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();</div><div><span class="" style="white-space:pre">  </span>optimizer->AddObserver(itk::IterationEvent(), observer);</div><div><br></div><div><br></div><div><span class="" style="white-space:pre">      </span>typedef RegistrationInterfaceCommand<RegistrationType> CommandType;</div><div><span class="" style="white-space:pre">  </span>CommandType::Pointer command = CommandType::New();</div><div><span class="" style="white-space:pre"> </span>registration->AddObserver(itk::IterationEvent(), command);</div><div><br></div><div><span class="" style="white-space:pre">     </span>try</div><div><span class="" style="white-space:pre">        </span>{</div><div><span class="" style="white-space:pre">          </span>registration->Update();</div><div><span class="" style="white-space:pre">         </span>std::cout << "Optimizer stop condition: "</div><div><span class="" style="white-space:pre">                  </span><< registration->GetOptimizer()->GetStopConditionDescription()</div><div><span class="" style="white-space:pre">                 </span><< std::endl;</div><div><span class="" style="white-space:pre">        </span>}</div><div><span class="" style="white-space:pre">  </span>catch (itk::ExceptionObject & err)</div><div><span class="" style="white-space:pre">     </span>{</div><div><span class="" style="white-space:pre">          </span>std::cout << "ExceptionObject caught !" << std::endl;</div><div><span class="" style="white-space:pre">                </span>std::cout << err << std::endl;</div><div><span class="" style="white-space:pre">         </span>return EXIT_FAILURE;</div><div><span class="" style="white-space:pre">       </span>}</div><div><span class="" style="white-space:pre">  </span>std::cout << "Optimizer Stopping Condition = "</div><div><span class="" style="white-space:pre">             </span><< optimizer->GetStopCondition() << std::endl;</div><div><span class="" style="white-space:pre">      </span>typedef RegistrationType::ParametersType ParametersType;</div><div><span class="" style="white-space:pre">   </span>ParametersType finalParameters = registration->GetLastTransformParameters();</div><div><span class="" style="white-space:pre">    </span>double TranslationAlongX = finalParameters[4];</div><div><span class="" style="white-space:pre">     </span>double TranslationAlongY = finalParameters[5];</div><div><span class="" style="white-space:pre">     </span>unsigned int numberOfIterations = optimizer->GetCurrentIteration();</div><div><span class="" style="white-space:pre">     </span>double bestValue = optimizer->GetValue();</div><div><br></div><div><span class="" style="white-space:pre">      </span>std::cout << "Result = " << std::endl;</div><div><span class="" style="white-space:pre">       </span>std::cout << " Translation X = " << TranslationAlongX << std::endl;</div><div><span class="" style="white-space:pre">    </span>std::cout << " Translation Y = " << TranslationAlongY << std::endl;</div><div><span class="" style="white-space:pre">    </span>std::cout << " Iterations    = " << numberOfIterations << std::endl;</div><div><span class="" style="white-space:pre"> </span>std::cout << " Metric value  = " << bestValue << std::endl;</div><div><br></div><div><br></div><div><span class="" style="white-space:pre"> </span>typedef itk::ResampleImageFilter<</div><div><span class="" style="white-space:pre">               </span>MovingImageType,</div><div><span class="" style="white-space:pre">           </span>FixedImageType >    ResampleFilterType;</div><div><span class="" style="white-space:pre">       </span>TransformType::Pointer finalTransform = TransformType::New();</div><div><span class="" style="white-space:pre">      </span>finalTransform->SetParameters(finalParameters);</div><div><span class="" style="white-space:pre"> </span>finalTransform->SetFixedParameters(transform->GetFixedParameters());</div><div><span class="" style="white-space:pre"> </span>ResampleFilterType::Pointer resample = ResampleFilterType::New();</div><div><span class="" style="white-space:pre">  </span>resample->SetTransform(finalTransform);</div><div><span class="" style="white-space:pre"> </span>resample->SetInput(movingImageReader->GetOutput());</div><div><span class="" style="white-space:pre">  </span>FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</div><div><span class="" style="white-space:pre">     </span>PixelType backgroundGrayLevel = 0;</div><div><br></div><div><span class="" style="white-space:pre">        </span>resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());</div><div><span class="" style="white-space:pre"> </span>resample->SetOutputOrigin(fixedImage->GetOrigin());</div><div><span class="" style="white-space:pre">  </span>resample->SetOutputSpacing(fixedImage->GetSpacing());</div><div><span class="" style="white-space:pre">        </span>resample->SetOutputDirection(fixedImage->GetDirection());</div><div><span class="" style="white-space:pre">    </span>resample->SetDefaultPixelValue(backgroundGrayLevel);</div><div><span class="" style="white-space:pre">    </span>typedef  signed short                           OutputPixelType;</div><div><span class="" style="white-space:pre">     </span>typedef itk::Image< OutputPixelType, Dimension > OutputImageType;</div><div><span class="" style="white-space:pre">    </span>typedef itk::CastImageFilter<</div><div><span class="" style="white-space:pre">           </span>FixedImageType,</div><div><span class="" style="white-space:pre">            </span>OutputImageType >          CastFilterType;</div><div><span class="" style="white-space:pre"> </span>typedef itk::ImageFileWriter< OutputImageType >  WriterType;</div><div><span class="" style="white-space:pre">        </span>WriterType::Pointer      writer = WriterType::New();</div><div><span class="" style="white-space:pre">    </span>CastFilterType::Pointer  caster = CastFilterType::New();</div><div><span class="" style="white-space:pre">  </span>writer->SetFileName(argv[3]);</div><div><span class="" style="white-space:pre">   </span>caster->SetInput(resample->GetOutput());</div><div><span class="" style="white-space:pre">     </span>writer->SetInput(caster->GetOutput());</div><div><span class="" style="white-space:pre">       </span>writer->Update();</div><div><br></div><div><span class="" style="white-space:pre">      </span>typedef itk::CheckerBoardImageFilter< FixedImageType > CheckerBoardFilterType;</div><div><span class="" style="white-space:pre">       </span>CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();</div><div><span class="" style="white-space:pre">   </span>checker->SetInput1(fixedImage);</div><div><span class="" style="white-space:pre"> </span>checker->SetInput2(resample->GetOutput());</div><div><span class="" style="white-space:pre">   </span>caster->SetInput(checker->GetOutput());</div><div><span class="" style="white-space:pre">      </span>writer->SetInput(caster->GetOutput());</div><div><span class="" style="white-space:pre">       </span>resample->SetDefaultPixelValue(0);</div><div><br></div><div><span class="" style="white-space:pre">     </span>resample->SetTransform(finalTransform);</div><div><span class="" style="white-space:pre"> </span>if (argc > 4)</div><div><span class="" style="white-space:pre">   </span>{</div><div><span class="" style="white-space:pre">          </span>writer->SetFileName(argv[11]);</div><div><span class="" style="white-space:pre">          </span>writer->Update();</div><div><span class="" style="white-space:pre">       </span>}</div><div><span class="" style="white-space:pre">  </span>return EXIT_SUCCESS;</div><div>}</div></div><div><br></div></div>