#include #include #include #include #include #include #include namespace sitk = itk::simple; sitk::Transform ImageRegistrationExample1(sitk::Image fixed, sitk::Image moving, std::vector rotation, std::vector translation) { sitk::ImageRegistrationMethod R; R.SetMetricAsMeanSquares(); const double maxStep = 4.0; const double minStep = 0.01; const unsigned int numberOfIterations = 200; const double relaxationFactor = 0.5; R.SetOptimizerAsRegularStepGradientDescent(maxStep, minStep, numberOfIterations, relaxationFactor); R.SetInitialTransform(sitk::TranslationTransform(2, translation)); R.SetInterpolator(sitk::sitkLinear); sitk::Transform outTx = R.Execute(fixed, moving); return outTx; } int main() { //////GET INITIAL GUESSES (Ref starting points) std::ifstream initGuess; initGuess.open("initGuess.txt"); std::vector initialVals_x, initialVals_y; for (int i = 0; i < 100; i++) { int x, y; initGuess >> x; initGuess >> y; initialVals_x.push_back(x); initialVals_y.push_back(y); } for (int i = 0; i < initialVals_x.size(); i++) { std::cout << initialVals_x[i] << " " << initialVals_y[i] << std::endl; } for (int ii = 0; ii < 100; ii++) { sitk::Image moving = sitk::ReadImage("subset00_ref.tif", sitk::sitkFloat32); sitk::Image fixed = sitk::ReadImage("area00.tif", sitk::sitkFloat32); ////rotation matrix std::vector rot; rot.push_back(1.0); rot.push_back(0.0); rot.push_back(0.0); rot.push_back(1.0); //translation std::vector tr; double init_x = -1.0 * initialVals_x[ii]; double init_y = -1.0 * initialVals_y[ii]; tr.push_back(init_x); tr.push_back(init_y); sitk::Transform outTx = ImageRegistrationExample1(fixed, moving, rot, tr); sitk::WriteTransform(outTx, "TRANSFORM.txt"); //reverse parameters maps the moving image to the fixed image sitk::Transform rev = outTx.GetInverse(); std::vector params = rev.GetParameters(); sitk::ResampleImageFilter resample; resample.SetReferenceImage(fixed); resample.SetOutputOrigin(fixed.GetOrigin()); resample.SetOutputSpacing(fixed.GetSpacing()); resample.SetOutputDirection(fixed.GetDirection()); //resample.SetDefaultPixelValue(1); resample.SetSize(fixed.GetSize()); resample.SetTransform(outTx); sitk::Image output = resample.Execute(moving); sitk::RescaleIntensityImageFilter rescale; rescale.SetOutputMinimum(0); rescale.SetOutputMaximum(1); sitk::Image rescaled = rescale.Execute(output); sitk::WriteImage(rescaled, "rescaled.tif"); } }