"""
Benchmark for DROP.
DROP is a approach for image registration and motion estimation based on Markov Random Fields.
.. seealso:: https://www.mrf-registration.net
Related Publication:
1. Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods
Authors: Ben Glocker, Aristeidis Sotiras, Nikos Komodakis, Nikos Paragios
Published in: Annual Review of Biomedical Engineering, Vol. 12, 2011, pp. 219-244
**Installation for Linux**
1. Download executable according your operation system,
https://www.mrf-registration.net/deformable/index.html
2. Copy/extract executables and libraries to you favourite destination
3. Install all missing libraries such as QT4 with OpenGL support
4. Test calling the executable `./dropreg2d` which should return something like::
Usage: dropreg2d <source> <target> <result> <paramfile> [mask]
Usage
-----
To see the explanation of particular parameters see the User Manual
http://www.mrf-registration.net/download/drop_user_guide_V1.05.pdf
Sample run::
mkdir ./results
python bm_experiments/bm_DROP.py \
-t ./data-images/pairs-imgs-lnds_histol.csv \
-d ./data-images \
-o ./results \
-DROP $HOME/Applications/DROP/dropreg2d \
--path_config ./configs/DROP.txt \
--visual --unique
.. note:: experiments was tested on Ubuntu (Linux) based OS system
.. note:: to check whether you have all needed libraries on Linux use `ldd dropreg2d`,
see: https://askubuntu.com/a/709271/863070
AND set path to the `libdroplib.so` as `export LD_LIBRARY_PATH=/lib:/usr/lib:/usr/local/lib`,
see: https://unix.stackexchange.com/a/67783 ; https://stackoverflow.com/a/49660575/4521646
.. note:: This method is not optimized nor suitable for large images, so all used images
are first scaled to be 2000x2000 pixels and then the registration is performed.
After registration is resulting image scaled back. The landmarks are scalded accordingly.
Glocker, Ben, et al. "Deformable medical image registration: setting the state of the art
with discrete methods." Annual review of biomedical engineering 13 (2011): 219-244.
https://pdfs.semanticscholar.org/0c5f/277d357e3667b18b5420fd660f221c935fcc.pdf
Copyright (C) 2017-2019 Jiri Borovec <jiri.borovec@fel.cvut.cz>
"""
import glob
import logging
import os
import shutil
import sys
import time
import numpy as np
import SimpleITK as sitk
sys.path += [os.path.abspath('.'), os.path.abspath('..')] # Add path to root
from birl.benchmark import ImRegBenchmark
from birl.utilities.data_io import (
convert_image_from_mhd,
convert_image_to_mhd,
image_sizes,
load_landmarks,
save_landmarks,
)
from bm_experiments import bm_comp_perform
class BmDROP(ImRegBenchmark):
""" Benchmark for DROP
no run test while this method requires manual installation of DROP
For the app installation details, see module details.
.. note:: DROP requires gray scale images in MHD format where pixel values
are in range (0, 255) of uint8.
Example
-------
>>> from birl.utilities.data_io import create_folder, update_path
>>> path_out = create_folder('temp_results')
>>> path_csv = os.path.join(update_path('data-images'), 'pairs-imgs-lnds_mix.csv')
>>> params = {'path_table': path_csv,
... 'path_out': path_out,
... 'nb_workers': 2,
... 'unique': False,
... 'visual': True,
... 'exec_DROP': 'dropreg2d',
... 'path_config': os.path.join(update_path('configs'), 'DROP.txt')}
>>> benchmark = BmDROP(params)
>>> benchmark.run() # doctest: +SKIP
>>> import shutil
>>> shutil.rmtree(path_out, ignore_errors=True)
"""
#: required experiment parameters
REQUIRED_PARAMS = ImRegBenchmark.REQUIRED_PARAMS + ['exec_DROP', 'path_config']
#: maximal image size (diagonal in pixels) recommended for DROP registration
MAX_IMAGE_DIAGONAL = int(np.sqrt(2e3**2 + 2e3**2))
#: time need for image conversion and optional scaling
COL_TIME_CONVERT = 'conversion time [s]'
def _prepare(self):
logging.info('-> copy configuration...')
self._copy_config_to_expt('path_config')
def _prepare_img_registration(self, item):
""" converting the input images to gra-scale and MHD format
:param dict item: dictionary with registration params
:return dict: the same or updated registration info
"""
logging.debug('.. converting images to MHD')
path_im_ref, path_im_move, _, _ = self._get_paths(item)
path_reg_dir = self._get_path_reg_dir(item)
diags = [image_sizes(p_img)[1] for p_img in (path_im_ref, path_im_move)]
item['scaling'] = max(1, max(diags) / float(self.MAX_IMAGE_DIAGONAL))
t_start = time.time()
for path_img, col in [(path_im_ref, self.COL_IMAGE_REF), (path_im_move, self.COL_IMAGE_MOVE)]:
item[col + self.COL_IMAGE_EXT_TEMP] = convert_image_to_mhd(
path_img, path_out_dir=path_reg_dir, overwrite=False, to_gray=True, scaling=item.get('scaling', 1.)
)
item[self.COL_TIME_CONVERT] = time.time() - t_start
return item
def _generate_regist_command(self, item):
""" generate the registration command
:param dict item: dictionary with registration params
:return str|list(str): the execution commands
"""
logging.debug('.. prepare DROP registration command')
path_im_ref, path_im_move, _, _ = self._get_paths(item)
path_dir = self._get_path_reg_dir(item)
# NOTE: for some reason " in the command makes it crash in Run mode
# somehow it works fine even with " in Debug mode
command = ' '.join([
self.params['exec_DROP'],
path_im_move,
path_im_ref,
os.path.join(path_dir, 'output'),
self.params['path_config'],
])
return command
def _extract_warped_image_landmarks(self, item):
""" get registration results - warped registered images and landmarks
:param dict item: dictionary with registration params
:return dict: paths to warped images/landmarks
"""
path_reg_dir = self._get_path_reg_dir(item)
_, path_im_move, path_lnds_ref, _ = self._get_paths(item)
# convert MHD image
path_img_ = convert_image_from_mhd(os.path.join(path_reg_dir, 'output.mhd'), scaling=item.get('scaling', 1.))
img_name, _ = os.path.splitext(os.path.basename(path_im_move))
_, img_ext = os.path.splitext(os.path.basename(path_img_))
path_img_warp = path_img_.replace('output' + img_ext, img_name + img_ext)
shutil.move(path_img_, path_img_warp)
# load transform and warp landmarks
# lnds_move = load_landmarks(path_lnds_move)
lnds_ref = load_landmarks(path_lnds_ref)
lnds_name = os.path.basename(path_lnds_ref)
path_lnds_warp = os.path.join(path_reg_dir, lnds_name)
if lnds_ref is None:
raise ValueError('missing landmarks to be transformed "%s"' % lnds_name)
# down-scale landmarks if defined
lnds_ref = lnds_ref / item.get('scaling', 1.)
# extract deformation
path_deform_x = os.path.join(path_reg_dir, 'output_x.mhd')
path_deform_y = os.path.join(path_reg_dir, 'output_y.mhd')
try:
shift = self.extract_landmarks_shift_from_mhd(path_deform_x, path_deform_y, lnds_ref)
except Exception:
logging.exception(path_reg_dir)
shift = np.zeros(lnds_ref.shape)
# lnds_warp = lnds_move - shift
lnds_warp = lnds_ref + shift
# upscale landmarks if defined
lnds_warp = lnds_warp * item.get('scaling', 1.)
save_landmarks(path_lnds_warp, lnds_warp)
# return formatted results
return {
self.COL_IMAGE_MOVE_WARP: path_img_warp,
self.COL_POINTS_REF_WARP: path_lnds_warp,
}
def _clear_after_registration(self, item, patterns=('output*', '*.mhd', '*.raw')):
""" clean unnecessarily files after the registration
:param dict item: dictionary with registration information
:param list(str) patterns: string patterns of file names
:return dict: the same or updated registration info
"""
logging.debug('.. cleaning after registration experiment, remove `output`')
path_reg_dir = self._get_path_reg_dir(item)
for ptn in patterns:
for p_file in glob.glob(os.path.join(path_reg_dir, ptn)):
os.remove(p_file)
return item
@staticmethod
def extend_parse(arg_parser):
""" extent the basic arg parses by some extra required parameters
:return object:
"""
# SEE: https://docs.python.org/3/library/argparse.html
arg_parser.add_argument(
'-DROP', '--exec_DROP', type=str, required=True, help='path to DROP executable, use `dropreg2d`'
)
arg_parser.add_argument(
'-cfg', '--path_config', type=str, required=True, help='parameters for DROP registration'
)
return arg_parser
@staticmethod
def extract_landmarks_shift_from_mhd(path_deform_x, path_deform_y, lnds):
""" given pair of deformation fields and landmark positions get shift
:param str path_deform_x: path to deformation field in X axis
:param str path_deform_y: path to deformation field in Y axis
:param ndarray lnds: landmarks
:return ndarray: shift for each landmarks
"""
# define function for parsing particular shift from MHD
def __parse_shift(path_deform_, lnds):
if not os.path.isfile(path_deform_):
raise FileNotFoundError('missing deformation: %s' % path_deform_)
deform_ = sitk.GetArrayFromImage(sitk.ReadImage(path_deform_))
if deform_ is None:
raise ValueError('loaded deformation is Empty - %s' % path_deform_)
lnds_max = np.max(lnds, axis=0)[::-1]
if not all(ln < dim for ln, dim in zip(lnds_max, deform_.shape)):
raise ValueError(
'landmarks max %s is larger then (exceeded) deformation shape %s' %
(lnds_max.tolist(), deform_.shape)
)
shift_ = deform_[lnds[:, 1], lnds[:, 0]]
return shift_
lnds = np.array(np.round(lnds), dtype=int)
# get shift in both axis
shift_x = __parse_shift(path_deform_x, lnds)
shift_y = __parse_shift(path_deform_y, lnds)
# concatenate
shift = np.array([shift_x, shift_y]).T
return shift
# RUN by given parameters
if __name__ == "__main__":
logging.basicConfig(level=logging.INFO)
logging.info(__doc__)
arg_params, path_expt = BmDROP.main()
if arg_params.get('run_comp_benchmark', False):
bm_comp_perform.main(path_expt)
