Deformable Image Registration (DIR) - cerr/CERR GitHub Wiki
Deformable Image Registration (DIR) is available in CERR. Three commonly available software packages are enabled for use: Plastimatch, Elastix & ANTs. Examples of implementation and usage are provided below.
Workflow
- register_scans - this routine accepts image input for moving and target images in multiple formats (planC, MHA, NIfTI) and generates the deformation mapping with the user-specified registration software and algorithm.
The resulting deformation mapping can be used to transform scans, label maps and RT dose into the fixed space. Or, pass the "inverseFlag" to use the deformation inverse to match fixed to moving space.
- warp_scan - apply the generated deformation to another scan in moving image space to match target space.
- warp_structures - deform a structure or label map to target space (nearest neighbor interpolation)
- warp_dose - warp RT dose distribution
Example
% Filename of base scan
baseCerrFile = '/path/to/base/basePlanCfile.mat';
% Filename of moving scan
movCerrFile = '/path/to/base/movPlanCfile.mat';
% Filename to save registered images
registeredFileName = '/path/to/base/registeredPlanCfile.mat';
% Plastimatch settings file
% ---- settings files that uses alignment of centers of scan volumes.
% Useful when initialTranslationXyzV is not passed to register_scans.m
% plmSettingsFile = fullfile(getCERRPath,'ImageRegistration',...
% 'plastimatch_command','bspline_register_cmd_dir.txt');
% ---- settings files that uses initial input transform constructed by
% passing initialTranslationXyzV to register_scans.m. Use this settings file
% in order to pass initial translation tranform defined by passing
% initialTranslationXyzV to register_scans.m
plmSettingsFile = fullfile(getCERRPath,'ImageRegistration',...
'plastimatch_command','plm_reg_settings_with_starting_transform.txt');
% temporary directory for cerr/plastimatch i/o
tmpDirPath = '/path/to/temp';
% index of scan to be used in base and moving planC's
baseScanNum = 1;
movScanNum = 1;
% Define structures to deform from moving to base scan
strsToDeformC = {'Struct1','Struct2'};
strCreationScanNum = baseScanNum;
% Load base image as planC
planC = loadPlanC(baseCerrFile,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(baseCerrFile,planC);
indexS = planC{end};
% Load moving image as planD
planD = loadPlanC(movCerrFile,tempdir);
planD = updatePlanFields(planD);
planD = quality_assure_planC(movCerrFile,planD);
indexSD = planD{end};
% NOTE - base and moving images can also be loaded from the same planC, if
% available
% Define registration settings
registration_tool = 'PLASTIMATCH';
baseMask3M = [];
movMask3M = [];
threshold_bone = -800;
inBspFile = '';
outBspFile = '';
algorithm = 'BSPLINE PLASTIMATCH';
% Generate landmarks
baseLandmarkStrC = {''};
movLandmarkStrC = {''};
% Register planD to planC
strBaseC = {planC{indexS.structures}.structureName};
strMovC = {planD{indexSD.structures}.structureName};
landmarkListM = [];
for iStr = 1:length(baseLandmarkStrC)
baseLandmarkStr = baseLandmarkStrC{iStr};
movLandmarkStr = movLandmarkStrC{iStr};
strIndV = getMatchingIndex(baseLandmarkStr,strBaseC,'exact');
assocScanNumV = getStructureAssociatedScan(strIndV,planC);
baseLandmarkInd = strIndV(assocScanNumV == baseScanNum);
strIndV = getMatchingIndex(movLandmarkStr,strMovC,'exact');
assocScanNumV = getStructureAssociatedScan(strIndV,planD);
movLandmarkInd = strIndV(assocScanNumV == movScanNum);
if length(baseLandmarkInd) ~= 1 || length(movLandmarkInd) ~= 1
continue;
end
[baseX,baseY,baseZ] = calcIsocenter(baseLandmarkInd,'COM',planC);
[movX,movY,movZ] = calcIsocenter(movLandmarkInd,'COM',planD);
[xBaseV,yBaseV,zBaseV] = getScanXYZVals(planC{indexS.scan}(baseScanNum));
[xMovV,yMovV,zMovV] = getScanXYZVals(planD{indexSD.scan}(movScanNum));
baseY = - baseY;
movY = - movY;
landmarkListM(iStr,:,1) = [iStr-1, [-baseX,-baseY,-baseZ]*10];
landmarkListM(iStr,:,2) = [iStr-1, [-movX,-movY,-movZ]*10];
end
% Initial transform
% Get center of scan volumes. This can be changed, for example,
% to get geometric centers to structures.
[xV,yV,zV] = getScanXYZVals(planC{indexS.scan}(baseScanNum));
xyzBaseV = [mean(xV), mean(yV), mean(zV)];
[xV,yV,zV] = getScanXYZVals(planD{indexSD.scan}(movScanNum));
xyzMovV = [mean(xV), mean(yV), mean(zV)];
initialTranslationXyzV = xyzMovV - xyzBaseV;
initialTranslationXyzV(1) = initialTranslationXyzV(1);
initialTranslationXyzV(2) = -initialTranslationXyzV(2); %planC to mha coords
initialTranslationXyzV(3) = -initialTranslationXyzV(3); %planC to mha coords
% Invoke Registration
planC = register_scans(planC, baseScanNum, planD, movScanNum,...
algorithm, registration_tool, tmpDirPath, baseMask3M,...
movMask3M, threshold_bone, plmSettingsFile, inBspFile,...
outBspFile, landmarkListM, initialTranslationXyzV);
% Warp structures from moving to the base scan (Note that scan is warped in register_scans.m
% Access the resulting defomration
deformS = planC{indexS.deform}(end);
% Get matching indices for structures to deform
movStructNumsV = [];
strNamC = {planD{indexSD.structures}.structureName};
for iStr = 1:length(strsToDeformC)
movStructNumsV = [movStructNumsV,...
getMatchingIndex(strsToDeformC{iStr},strNamC,'exact')];
end
% Invoke function to warp structures
planC = warp_structures(deformS,strCreationScanNum,movStructNumsV,planD,planC,tmpDirPath);
% Save planC
planC = save_planC(planC,[],'passed',registeredFileName);
Note: plmSettingsFile refers to file containing plastimatch settings. These can be adjusted as per the image types. Users can choose from the sample settings files or use them as templates to adjust settings.
Deformable Registration Software Links
Plastimatch https://plastimatch.org/
Elastix https://elastix.lumc.nl/