diff --git a/.github/workflows/integration.yml b/.github/workflows/integration.yml index 4b93632723..ce3598d17e 100644 --- a/.github/workflows/integration.yml +++ b/.github/workflows/integration.yml @@ -3,6 +3,9 @@ name: integration on: repository_dispatch: type: [integration-test-command] + push: + branches: + - add-spatial-resample jobs: integration-py3: @@ -32,7 +35,7 @@ jobs: - name: Install the dependencies run: | which python - python -m pip install --upgrade pip wheel + python -m pip install --upgrade pip==21.3.1 wheel # workaround monai##3752 python -m pip uninstall -y torch torchvision python -m pip install torch==1.10.1+cu111 torchvision==0.11.2+cu111 -f https://download.pytorch.org/whl/torch_stable.html python -m pip install -r requirements-dev.txt diff --git a/docs/source/transforms.rst b/docs/source/transforms.rst index 49ed4c9e6c..a8d04efaa5 100644 --- a/docs/source/transforms.rst +++ b/docs/source/transforms.rst @@ -557,6 +557,12 @@ Post-processing Spatial ^^^^^^^ +`SpatialResample` +""""""""""""""""" +.. autoclass:: SpatialResample + :members: + :special-members: __call__ + `Spacing` """"""""" .. image:: https://github.com/Project-MONAI/DocImages/raw/main/transforms/Spacing.png @@ -1398,6 +1404,12 @@ Post-processing (Dict) Spatial (Dict) ^^^^^^^^^^^^^^ +`SpatialResampled` +"""""""""""""""""" +.. autoclass:: SpatialResampled + :members: + :special-members: __call__ + `Spacingd` """""""""" .. image:: https://github.com/Project-MONAI/DocImages/raw/main/transforms/Spacingd.png diff --git a/monai/apps/deepgrow/transforms.py b/monai/apps/deepgrow/transforms.py index 310931d236..59dee17280 100644 --- a/monai/apps/deepgrow/transforms.py +++ b/monai/apps/deepgrow/transforms.py @@ -374,7 +374,7 @@ class SpatialCropForegroundd(MapTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `{key}_{meta_key_postfix}` to to fetch/store the meta data according + meta_key_postfix: if meta_keys is None, use `{key}_{meta_key_postfix}` to fetch/store the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -474,7 +474,7 @@ class AddGuidanceFromPointsd(Transform): for example, for data with key `image`, the metadata by default is in `image_meta_dict`. the meta data is a dictionary object which contains: filename, original_shape, etc. if None, will try to construct meta_keys by `{ref_image}_{meta_key_postfix}`. - meta_key_postfix: if meta_key is None, use `{ref_image}_{meta_key_postfix}` to to fetch the meta data according + meta_key_postfix: if meta_key is None, use `{ref_image}_{meta_key_postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -589,7 +589,7 @@ class SpatialCropGuidanced(MapTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -787,7 +787,7 @@ class RestoreLabeld(MapTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_key is None, use `key_{meta_key_postfix} to to fetch the meta data according + meta_key_postfix: if meta_key is None, use `key_{meta_key_postfix} to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -897,7 +897,7 @@ class Fetch2DSliced(MapTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: use `key_{meta_key_postfix}` to to fetch the meta data according to the key data, + meta_key_postfix: use `key_{meta_key_postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. diff --git a/monai/csrc/ext.cpp b/monai/csrc/ext.cpp index a2fa8bfc56..ac43e6fd3e 100644 --- a/monai/csrc/ext.cpp +++ b/monai/csrc/ext.cpp @@ -31,6 +31,7 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { py::enum_(m, "BoundType") .value("replicate", monai::BoundType::Replicate, "a a a | a b c d | d d d") .value("nearest", monai::BoundType::Replicate, "a a a | a b c d | d d d") + .value("border", monai::BoundType::Replicate, "a a a | a b c d | d d d") .value("dct1", monai::BoundType::DCT1, "d c b | a b c d | c b a") .value("mirror", monai::BoundType::DCT1, "d c b | a b c d | c b a") .value("dct2", monai::BoundType::DCT2, "c b a | a b c d | d c b") @@ -43,6 +44,7 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { .value("wrap", monai::BoundType::DFT, "b c d | a b c d | a b c") // .value("sliding", monai::BoundType::Sliding) .value("zero", monai::BoundType::Zero, "0 0 0 | a b c d | 0 0 0") + .value("zeros", monai::BoundType::Zero, "0 0 0 | a b c d | 0 0 0") .export_values(); // resample interpolation mode diff --git a/monai/csrc/resample/pushpull_cpu.cpp b/monai/csrc/resample/pushpull_cpu.cpp index d83557c6c3..4c9f6b9f9c 100644 --- a/monai/csrc/resample/pushpull_cpu.cpp +++ b/monai/csrc/resample/pushpull_cpu.cpp @@ -1527,10 +1527,10 @@ MONAI_NAMESPACE_DEVICE { // cpu iy0 = bound::index(bound1, iy0, src_Y); iz0 = bound::index(bound2, iz0, src_Z); - // Offsets into source volume offset_t o000, o100, o010, o001, o110, o011, o101, o111; if (do_pull || do_grad || do_sgrad) { + // Offsets into source volume o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; @@ -1539,18 +1539,20 @@ MONAI_NAMESPACE_DEVICE { // cpu o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; + } else { + // Offsets into 'push' volume + o000 = ix0 * out_sX + iy0 * out_sY + iz0 * out_sZ; + o100 = ix1 * out_sX + iy0 * out_sY + iz0 * out_sZ; + o010 = ix0 * out_sX + iy1 * out_sY + iz0 * out_sZ; + o001 = ix0 * out_sX + iy0 * out_sY + iz1 * out_sZ; + o110 = ix1 * out_sX + iy1 * out_sY + iz0 * out_sZ; + o011 = ix0 * out_sX + iy1 * out_sY + iz1 * out_sZ; + o101 = ix1 * out_sX + iy0 * out_sY + iz1 * out_sZ; + o111 = ix1 * out_sX + iy1 * out_sY + iz1 * out_sZ; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t gx = static_cast(0); scalar_t gy = static_cast(0); scalar_t gz = static_cast(0); @@ -1659,14 +1661,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t* out_ptr_NCXYZ = out_ptr + n * out_sN + w * out_sX + h * out_sY + d * out_sZ; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCXYZ += out_sC, src_ptr_NC += src_sC) { @@ -1678,14 +1672,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~ else if (do_sgrad) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t* out_ptr_NCXYZ = out_ptr + n * out_sN + w * out_sX + h * out_sY + d * out_sZ; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -1758,16 +1744,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o000 = ix0 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o100 = ix1 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o010 = ix0 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o001 = ix0 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o110 = ix1 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o011 = ix0 * out_sX + iy1 * out_sY + iz1 * out_sZ; - o101 = ix1 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o111 = ix1 * out_sX + iy1 * out_sY + iz1 * out_sZ; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o000, w000, s000); bound::add(out_ptr_N, o100, w100, s100); @@ -1822,21 +1798,23 @@ MONAI_NAMESPACE_DEVICE { // cpu ix0 = bound::index(bound0, ix0, src_X); iy0 = bound::index(bound1, iy0, src_Y); - // Offsets into source volume offset_t o00, o10, o01, o11; if (do_pull || do_grad || do_sgrad) { + // Offsets into source volume o00 = ix0 * src_sX + iy0 * src_sY; o10 = ix1 * src_sX + iy0 * src_sY; o01 = ix0 * src_sX + iy1 * src_sY; o11 = ix1 * src_sX + iy1 * src_sY; + } else { + // Offsets into 'push' volume + o00 = ix0 * out_sX + iy0 * out_sY; + o10 = ix1 * out_sX + iy0 * out_sY; + o01 = ix0 * out_sX + iy1 * out_sY; + o11 = ix1 * out_sX + iy1 * out_sY; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t gx = static_cast(0); scalar_t gy = static_cast(0); scalar_t* trgt_ptr_NCXY = trgt_ptr + n * trgt_sN + w * trgt_sX + h * trgt_sY; @@ -1895,10 +1873,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t* out_ptr_NCXY = out_ptr + n * out_sN + w * out_sX + h * out_sY; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCXY += out_sC, src_ptr_NC += src_sC) { @@ -1908,10 +1882,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_sgrad) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t* out_ptr_NCXY = out_ptr + n * out_sN + w * out_sX + h * out_sY; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -1926,11 +1896,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_push) { - // Offsets into 'push' volume - o00 = ix0 * out_sX + iy0 * out_sY; - o10 = ix1 * out_sX + iy0 * out_sY; - o01 = ix0 * out_sX + iy1 * out_sY; - o11 = ix1 * out_sX + iy1 * out_sY; scalar_t* trgt_ptr_NCXY = trgt_ptr + n * trgt_sN + w * trgt_sX + h * trgt_sY; scalar_t* out_ptr_NC = out_ptr + n * out_sN; if (trgt_K == 0) { @@ -1960,12 +1925,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o00 = ix0 * out_sX + iy0 * out_sY; - o10 = ix1 * out_sX + iy0 * out_sY; - o01 = ix0 * out_sX + iy1 * out_sY; - o11 = ix1 * out_sX + iy1 * out_sY; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o00, w00, s00); bound::add(out_ptr_N, o10, w10, s10); @@ -1996,20 +1955,21 @@ MONAI_NAMESPACE_DEVICE { // cpu ix1 = bound::index(bound0, ix0 + 1, src_X); ix0 = bound::index(bound0, ix0, src_X); - // Offsets into source volume offset_t o0, o1; if (do_pull || do_grad || do_sgrad) { + // Offsets into source volume o0 = ix0 * src_sX; o1 = ix1 * src_sX; + } else { + // Offsets into 'push' volume + o0 = ix0 * out_sX; + o1 = ix1 * out_sX; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { if (trgt_K == 0) { // backward w.r.t. push/pull - - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t gx = static_cast(0); scalar_t* trgt_ptr_NCX = trgt_ptr + n * trgt_sN + w * trgt_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -2037,8 +1997,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t* out_ptr_NCX = out_ptr + n * out_sN + w * out_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCX += out_sC, src_ptr_NC += src_sC) { @@ -2047,8 +2005,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_sgrad) { - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t* out_ptr_NCX = out_ptr + n * out_sN + w * out_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -2058,9 +2014,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_push) { - // Offsets into 'push' volume - o0 = ix0 * out_sX; - o1 = ix1 * out_sX; scalar_t* trgt_ptr_NCX = trgt_ptr + n * trgt_sN + w * trgt_sX; scalar_t* out_ptr_NC = out_ptr + n * out_sN; if (trgt_K == 0) { @@ -2081,10 +2034,6 @@ MONAI_NAMESPACE_DEVICE { // cpu } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o0 = ix0 * out_sX; - o1 = ix1 * out_sX; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o0, w0, s0); bound::add(out_ptr_N, o1, w1, s1); diff --git a/monai/csrc/resample/pushpull_cuda.cu b/monai/csrc/resample/pushpull_cuda.cu index 4a2d6c27ef..fe50154670 100644 --- a/monai/csrc/resample/pushpull_cuda.cu +++ b/monai/csrc/resample/pushpull_cuda.cu @@ -1491,10 +1491,10 @@ MONAI_NAMESPACE_DEVICE { // cuda iy0 = bound::index(bound1, iy0, src_Y); iz0 = bound::index(bound2, iz0, src_Z); - // Offsets into source volume offset_t o000, o100, o010, o001, o110, o011, o101, o111; if (do_pull || do_grad || do_sgrad) { + // Offsets into source volume o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; @@ -1503,18 +1503,20 @@ MONAI_NAMESPACE_DEVICE { // cuda o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; + } else { + // Offsets into 'push' volume + o000 = ix0 * out_sX + iy0 * out_sY + iz0 * out_sZ; + o100 = ix1 * out_sX + iy0 * out_sY + iz0 * out_sZ; + o010 = ix0 * out_sX + iy1 * out_sY + iz0 * out_sZ; + o001 = ix0 * out_sX + iy0 * out_sY + iz1 * out_sZ; + o110 = ix1 * out_sX + iy1 * out_sY + iz0 * out_sZ; + o011 = ix0 * out_sX + iy1 * out_sY + iz1 * out_sZ; + o101 = ix1 * out_sX + iy0 * out_sY + iz1 * out_sZ; + o111 = ix1 * out_sX + iy1 * out_sY + iz1 * out_sZ; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t gx = static_cast(0); scalar_t gy = static_cast(0); scalar_t gz = static_cast(0); @@ -1623,14 +1625,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t* out_ptr_NCXYZ = out_ptr + n * out_sN + w * out_sX + h * out_sY + d * out_sZ; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCXYZ += out_sC, src_ptr_NC += src_sC) { @@ -1642,14 +1636,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~ else if (do_sgrad) { - o000 = ix0 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o100 = ix1 * src_sX + iy0 * src_sY + iz0 * src_sZ; - o010 = ix0 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o001 = ix0 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o110 = ix1 * src_sX + iy1 * src_sY + iz0 * src_sZ; - o011 = ix0 * src_sX + iy1 * src_sY + iz1 * src_sZ; - o101 = ix1 * src_sX + iy0 * src_sY + iz1 * src_sZ; - o111 = ix1 * src_sX + iy1 * src_sY + iz1 * src_sZ; scalar_t* out_ptr_NCXYZ = out_ptr + n * out_sN + w * out_sX + h * out_sY + d * out_sZ; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -1672,15 +1658,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_push) { - // Offsets into 'push' volume - o000 = ix0 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o100 = ix1 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o010 = ix0 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o001 = ix0 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o110 = ix1 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o011 = ix0 * out_sX + iy1 * out_sY + iz1 * out_sZ; - o101 = ix1 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o111 = ix1 * out_sX + iy1 * out_sY + iz1 * out_sZ; scalar_t* trgt_ptr_NCXYZ = trgt_ptr + n * trgt_sN + w * trgt_sX + h * trgt_sY + d * trgt_sZ; scalar_t* out_ptr_NC = out_ptr + n * out_sN; if (trgt_K == 0) { @@ -1722,16 +1699,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o000 = ix0 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o100 = ix1 * out_sX + iy0 * out_sY + iz0 * out_sZ; - o010 = ix0 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o001 = ix0 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o110 = ix1 * out_sX + iy1 * out_sY + iz0 * out_sZ; - o011 = ix0 * out_sX + iy1 * out_sY + iz1 * out_sZ; - o101 = ix1 * out_sX + iy0 * out_sY + iz1 * out_sZ; - o111 = ix1 * out_sX + iy1 * out_sY + iz1 * out_sZ; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o000, w000, s000); bound::add(out_ptr_N, o100, w100, s100); @@ -1786,21 +1753,23 @@ MONAI_NAMESPACE_DEVICE { // cuda ix0 = bound::index(bound0, ix0, src_X); iy0 = bound::index(bound1, iy0, src_Y); - // Offsets into source volume offset_t o00, o10, o01, o11; if (do_pull || do_grad || do_sgrad) { + // Offsets into source volume o00 = ix0 * src_sX + iy0 * src_sY; o10 = ix1 * src_sX + iy0 * src_sY; o01 = ix0 * src_sX + iy1 * src_sY; o11 = ix1 * src_sX + iy1 * src_sY; + } else { + // Offsets into 'push' volume + o00 = ix0 * out_sX + iy0 * out_sY; + o10 = ix1 * out_sX + iy0 * out_sY; + o01 = ix0 * out_sX + iy1 * out_sY; + o11 = ix1 * out_sX + iy1 * out_sY; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t gx = static_cast(0); scalar_t gy = static_cast(0); scalar_t* trgt_ptr_NCXY = trgt_ptr + n * trgt_sN + w * trgt_sX + h * trgt_sY; @@ -1859,10 +1828,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t* out_ptr_NCXY = out_ptr + n * out_sN + w * out_sX + h * out_sY; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCXY += out_sC, src_ptr_NC += src_sC) { @@ -1872,10 +1837,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_sgrad) { - o00 = ix0 * src_sX + iy0 * src_sY; - o10 = ix1 * src_sX + iy0 * src_sY; - o01 = ix0 * src_sX + iy1 * src_sY; - o11 = ix1 * src_sX + iy1 * src_sY; scalar_t* out_ptr_NCXY = out_ptr + n * out_sN + w * out_sX + h * out_sY; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -1890,11 +1851,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_push) { - // Offsets into 'push' volume - o00 = ix0 * out_sX + iy0 * out_sY; - o10 = ix1 * out_sX + iy0 * out_sY; - o01 = ix0 * out_sX + iy1 * out_sY; - o11 = ix1 * out_sX + iy1 * out_sY; scalar_t* trgt_ptr_NCXY = trgt_ptr + n * trgt_sN + w * trgt_sX + h * trgt_sY; scalar_t* out_ptr_NC = out_ptr + n * out_sN; if (trgt_K == 0) { @@ -1924,12 +1880,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o00 = ix0 * out_sX + iy0 * out_sY; - o10 = ix1 * out_sX + iy0 * out_sY; - o01 = ix0 * out_sX + iy1 * out_sY; - o11 = ix1 * out_sX + iy1 * out_sY; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o00, w00, s00); bound::add(out_ptr_N, o10, w10, s10); @@ -1965,15 +1915,16 @@ MONAI_NAMESPACE_DEVICE { // cuda if (do_pull || do_grad || do_sgrad) { o0 = ix0 * src_sX; o1 = ix1 * src_sX; + } else { + // Offsets into 'push' volume + o0 = ix0 * out_sX; + o1 = ix1 * out_sX; } // ~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid gradient ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_grad) { if (trgt_K == 0) { // backward w.r.t. push/pull - - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t gx = static_cast(0); scalar_t* trgt_ptr_NCX = trgt_ptr + n * trgt_sN + w * trgt_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -2001,8 +1952,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Pull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (do_pull) { - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t* out_ptr_NCX = out_ptr + n * out_sN + w * out_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; for (offset_t c = 0; c < C; ++c, out_ptr_NCX += out_sC, src_ptr_NC += src_sC) { @@ -2011,8 +1960,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SGrad ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_sgrad) { - o0 = ix0 * src_sX; - o1 = ix1 * src_sX; scalar_t* out_ptr_NCX = out_ptr + n * out_sN + w * out_sX; scalar_t* src_ptr_NC = src_ptr + n * src_sN; @@ -2022,9 +1969,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_push) { - // Offsets into 'push' volume - o0 = ix0 * out_sX; - o1 = ix1 * out_sX; scalar_t* trgt_ptr_NCX = trgt_ptr + n * trgt_sN + w * trgt_sX; scalar_t* out_ptr_NC = out_ptr + n * out_sN; if (trgt_K == 0) { @@ -2045,10 +1989,6 @@ MONAI_NAMESPACE_DEVICE { // cuda } // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Push ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if (do_count) { - // Offsets into 'push' volume - o0 = ix0 * out_sX; - o1 = ix1 * out_sX; - scalar_t* out_ptr_N = out_ptr + n * out_sN; bound::add(out_ptr_N, o0, w0, s0); bound::add(out_ptr_N, o1, w1, s1); diff --git a/monai/data/nifti_writer.py b/monai/data/nifti_writer.py index 4e7a99f557..1ffae31183 100644 --- a/monai/data/nifti_writer.py +++ b/monai/data/nifti_writer.py @@ -18,6 +18,7 @@ from monai.config.type_definitions import NdarrayOrTensor from monai.data.utils import compute_shape_offset, to_affine_nd from monai.networks.layers import AffineTransform +from monai.transforms.utils_pytorch_numpy_unification import allclose from monai.utils import GridSampleMode, GridSamplePadMode, optional_import from monai.utils.type_conversion import convert_data_type @@ -98,10 +99,8 @@ def write_nifti( If None, use the data type of input data. output_dtype: data type for saving data. Defaults to ``np.float32``. """ - if isinstance(data, torch.Tensor): - data, *_ = convert_data_type(data, np.ndarray) - if isinstance(affine, torch.Tensor): - affine, *_ = convert_data_type(affine, np.ndarray) + data, *_ = convert_data_type(data, np.ndarray) + affine, *_ = convert_data_type(affine, np.ndarray) if not isinstance(data, np.ndarray): raise AssertionError("input data must be numpy array or torch tensor.") dtype = dtype or data.dtype @@ -111,12 +110,12 @@ def write_nifti( affine = to_affine_nd(sr, affine) # type: ignore if target_affine is None: - target_affine = affine - target_affine = to_affine_nd(sr, target_affine) + target_affine = affine # type: ignore + target_affine, *_ = convert_data_type(to_affine_nd(sr, target_affine), np.ndarray) # type: ignore - if np.allclose(affine, target_affine, atol=1e-3): + if allclose(affine, target_affine, atol=1e-3): # type: ignore # no affine changes, save (data, affine) - results_img = nib.Nifti1Image(data.astype(output_dtype, copy=False), to_affine_nd(3, target_affine)) + results_img = nib.Nifti1Image(data.astype(output_dtype, copy=False), to_affine_nd(3, target_affine)) # type: ignore nib.save(results_img, file_name) return @@ -127,7 +126,7 @@ def write_nifti( data_shape = data.shape data = nib.orientations.apply_orientation(data, ornt_transform) _affine = affine @ nib.orientations.inv_ornt_aff(ornt_transform, data_shape) - if np.allclose(_affine, target_affine, atol=1e-3) or not resample: + if allclose(_affine, target_affine, atol=1e-3) or not resample: # type: ignore results_img = nib.Nifti1Image(data.astype(output_dtype, copy=False), to_affine_nd(3, _affine)) # type: ignore nib.save(results_img, file_name) return @@ -138,8 +137,8 @@ def write_nifti( ) transform = np.linalg.inv(_affine) @ target_affine if output_spatial_shape is None: - output_spatial_shape, _ = compute_shape_offset(data.shape, _affine, target_affine) - output_spatial_shape_ = list(output_spatial_shape) if output_spatial_shape is not None else [] + output_spatial_shape, _ = compute_shape_offset(data.shape, _affine, target_affine) # type: ignore + output_spatial_shape_ = list(output_spatial_shape) if output_spatial_shape is not None else [] # type: ignore if data.ndim > 3: # multi channel, resampling each channel while len(output_spatial_shape_) < 3: output_spatial_shape_ = output_spatial_shape_ + [1] @@ -164,6 +163,6 @@ def write_nifti( ) data_np = data_torch.squeeze(0).squeeze(0).detach().cpu().numpy() - results_img = nib.Nifti1Image(data_np.astype(output_dtype, copy=False), to_affine_nd(3, target_affine)) + results_img = nib.Nifti1Image(data_np.astype(output_dtype, copy=False), to_affine_nd(3, target_affine)) # type: ignore nib.save(results_img, file_name) return diff --git a/monai/data/test_time_augmentation.py b/monai/data/test_time_augmentation.py index 13e9ebaad6..c75a61c56d 100644 --- a/monai/data/test_time_augmentation.py +++ b/monai/data/test_time_augmentation.py @@ -74,7 +74,7 @@ class TestTimeAugmentation: orig_meta_keys: the key of the meta data of original input data, will get the `affine`, `data_shape`, etc. the meta data is a dictionary object which contains: filename, original_shape, etc. if None, will try to construct meta_keys by `{orig_key}_{meta_key_postfix}`. - meta_key_postfix: use `key_{postfix}` to to fetch the meta data according to the key data, + meta_key_postfix: use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. diff --git a/monai/data/utils.py b/monai/data/utils.py index 79ef9bd7fb..ae784dc07f 100644 --- a/monai/data/utils.py +++ b/monai/data/utils.py @@ -27,12 +27,15 @@ import torch from torch.utils.data._utils.collate import default_collate -from monai.config.type_definitions import PathLike +from monai.config.type_definitions import NdarrayOrTensor, PathLike from monai.networks.layers.simplelayers import GaussianFilter from monai.utils import ( MAX_SEED, BlendMode, + Method, NumpyPadMode, + convert_data_type, + convert_to_dst_type, ensure_tuple, ensure_tuple_rep, ensure_tuple_size, @@ -42,7 +45,6 @@ look_up_option, optional_import, ) -from monai.utils.enums import Method pd, _ = optional_import("pandas") DataFrame, _ = optional_import("pandas", name="DataFrame") @@ -78,6 +80,7 @@ "no_collation", "convert_tables_to_dicts", "SUPPORTED_PICKLE_MOD", + "reorient_spatial_axes", ] # module to be used by `torch.save` @@ -658,7 +661,7 @@ def zoom_affine(affine: np.ndarray, scale: Union[np.ndarray, Sequence[float]], d def compute_shape_offset( - spatial_shape: Union[np.ndarray, Sequence[int]], in_affine: np.ndarray, out_affine: np.ndarray + spatial_shape: Union[np.ndarray, Sequence[int]], in_affine: NdarrayOrTensor, out_affine: NdarrayOrTensor ) -> Tuple[np.ndarray, np.ndarray]: """ Given input and output affine, compute appropriate shapes @@ -673,69 +676,101 @@ def compute_shape_offset( """ shape = np.array(spatial_shape, copy=True, dtype=float) sr = len(shape) - in_affine = to_affine_nd(sr, in_affine) - out_affine = to_affine_nd(sr, out_affine) + in_affine = convert_data_type(to_affine_nd(sr, in_affine), np.ndarray)[0] # type: ignore + out_affine = convert_data_type(to_affine_nd(sr, out_affine), np.ndarray)[0] # type: ignore in_coords = [(0.0, dim - 1.0) for dim in shape] - corners = np.asarray(np.meshgrid(*in_coords, indexing="ij")).reshape((len(shape), -1)) + corners: np.ndarray = np.asarray(np.meshgrid(*in_coords, indexing="ij")).reshape((len(shape), -1)) corners = np.concatenate((corners, np.ones_like(corners[:1]))) - corners = in_affine @ corners - corners_out = np.linalg.inv(out_affine) @ corners + corners = in_affine @ corners # type: ignore + try: + inv_mat = np.linalg.inv(out_affine) + except np.linalg.LinAlgError as e: + raise ValueError(f"Affine {out_affine} is not invertible") from e + corners_out = inv_mat @ corners corners_out = corners_out[:-1] / corners_out[-1] out_shape = np.round(corners_out.ptp(axis=1) + 1.0) - if np.allclose(nib.io_orientation(in_affine), nib.io_orientation(out_affine)): - # same orientation, get translate from the origin - offset = in_affine @ ([0] * sr + [1]) - offset = offset[:-1] / offset[-1] - else: - # different orientation, the min is the origin - corners = corners[:-1] / corners[-1] - offset = np.min(corners, 1) + mat = inv_mat[:-1, :-1] + k = 0 + for i in range(corners.shape[1]): + min_corner = np.min(mat @ corners[:-1, :] - mat @ corners[:-1, i : i + 1], 1) + if np.allclose(min_corner, 0.0, rtol=1e-3): + k = i + break + offset = corners[:-1, k] return out_shape.astype(int, copy=False), offset -def to_affine_nd(r: Union[np.ndarray, int], affine: np.ndarray) -> np.ndarray: +def to_affine_nd(r: Union[np.ndarray, int], affine: NdarrayOrTensor, dtype=np.float64) -> NdarrayOrTensor: """ Using elements from affine, to create a new affine matrix by assigning the rotation/zoom/scaling matrix and the translation vector. - when ``r`` is an integer, output is an (r+1)x(r+1) matrix, + When ``r`` is an integer, output is an (r+1)x(r+1) matrix, where the top left kxk elements are copied from ``affine``, the last column of the output affine is copied from ``affine``'s last column. `k` is determined by `min(r, len(affine) - 1)`. - when ``r`` is an affine matrix, the output has the same as ``r``, - the top left kxk elements are copied from ``affine``, + When ``r`` is an affine matrix, the output has the same shape as ``r``, + and the top left kxk elements are copied from ``affine``, the last column of the output affine is copied from ``affine``'s last column. `k` is determined by `min(len(r) - 1, len(affine) - 1)`. Args: r (int or matrix): number of spatial dimensions or an output affine to be filled. affine (matrix): 2D affine matrix + dtype: data type of the output array. Raises: ValueError: When ``affine`` dimensions is not 2. ValueError: When ``r`` is nonpositive. Returns: - an (r+1) x (r+1) matrix + an (r+1) x (r+1) matrix (tensor or ndarray depends on the input ``affine`` data type) """ - affine_np = np.array(affine, dtype=np.float64) + affine_np: np.ndarray + affine_np = convert_data_type(affine, output_type=np.ndarray, dtype=dtype, wrap_sequence=True)[0] # type: ignore + affine_np = affine_np.copy() if affine_np.ndim != 2: raise ValueError(f"affine must have 2 dimensions, got {affine_np.ndim}.") - new_affine = np.array(r, dtype=np.float64, copy=True) + new_affine = np.array(r, dtype=dtype, copy=True) if new_affine.ndim == 0: sr: int = int(new_affine.astype(np.uint)) if not np.isfinite(sr) or sr < 0: raise ValueError(f"r must be positive, got {sr}.") - new_affine = np.eye(sr + 1, dtype=np.float64) + new_affine = np.eye(sr + 1, dtype=dtype) d = max(min(len(new_affine) - 1, len(affine_np) - 1), 1) new_affine[:d, :d] = affine_np[:d, :d] if d > 1: new_affine[:d, -1] = affine_np[:d, -1] + new_affine, *_ = convert_to_dst_type(new_affine, affine, dtype=dtype) # type: ignore return new_affine +def reorient_spatial_axes( + data_shape: Sequence[int], init_affine: NdarrayOrTensor, target_affine: NdarrayOrTensor +) -> Tuple[np.ndarray, NdarrayOrTensor]: + """ + Given the input ``init_affine``, compute the orientation transform between + it and ``target_affine`` by rearranging/flipping the axes. + + Returns the orientation transform and the updated affine (tensor or ndarray + depends on the input ``affine`` data type). + Note that this function requires external module ``nibabel.orientations``. + """ + init_affine_, *_ = convert_data_type(init_affine, np.ndarray) + target_affine_, *_ = convert_data_type(target_affine, np.ndarray) + start_ornt = nib.orientations.io_orientation(init_affine_) + target_ornt = nib.orientations.io_orientation(target_affine_) + try: + ornt_transform = nib.orientations.ornt_transform(start_ornt, target_ornt) + except ValueError as e: + raise ValueError(f"The input affine {init_affine} and target affine {target_affine} are not compatible.") from e + new_affine = init_affine_ @ nib.orientations.inv_ornt_aff(ornt_transform, data_shape) + new_affine, *_ = convert_to_dst_type(new_affine, init_affine) + return ornt_transform, new_affine + + def create_file_basename( postfix: str, input_file_name: PathLike, diff --git a/monai/metrics/utils.py b/monai/metrics/utils.py index ccb6d93862..cf62cf9960 100644 --- a/monai/metrics/utils.py +++ b/monai/metrics/utils.py @@ -116,7 +116,7 @@ def get_mask_edges( The input images can be binary or labelfield images. If labelfield images are supplied, they are converted to binary images using `label_idx`. - `scipy`'s binary erosion is used to to calculate the edges of the binary + `scipy`'s binary erosion is used to calculate the edges of the binary labelfield. In order to improve the computing efficiency, before getting the edges, diff --git a/monai/networks/layers/spatial_transforms.py b/monai/networks/layers/spatial_transforms.py index 7aa3e110fc..56f13736b4 100644 --- a/monai/networks/layers/spatial_transforms.py +++ b/monai/networks/layers/spatial_transforms.py @@ -70,13 +70,13 @@ def grid_pull( `bound` can be an int, a string or a BoundType. Possible values are:: - - 0 or 'replicate' or 'nearest' or BoundType.replicate + - 0 or 'replicate' or 'nearest' or BoundType.replicate or 'border' - 1 or 'dct1' or 'mirror' or BoundType.dct1 - 2 or 'dct2' or 'reflect' or BoundType.dct2 - 3 or 'dst1' or 'antimirror' or BoundType.dst1 - 4 or 'dst2' or 'antireflect' or BoundType.dst2 - 5 or 'dft' or 'wrap' or BoundType.dft - - 7 or 'zero' or BoundType.zero + - 7 or 'zero' or 'zeros' or BoundType.zero A list of values can be provided, in the order [W, H, D], to specify dimension-specific boundary conditions. diff --git a/monai/transforms/__init__.py b/monai/transforms/__init__.py index 408edd6d56..26601de76b 100644 --- a/monai/transforms/__init__.py +++ b/monai/transforms/__init__.py @@ -306,6 +306,7 @@ Rotate, Rotate90, Spacing, + SpatialResample, Zoom, ) from .spatial.dictionary import ( @@ -360,6 +361,9 @@ Spacingd, SpacingD, SpacingDict, + SpatialResampled, + SpatialResampleD, + SpatialResampleDict, Zoomd, ZoomD, ZoomDict, @@ -552,6 +556,7 @@ zero_margins, ) from .utils_pytorch_numpy_unification import ( + allclose, any_np_pt, ascontiguousarray, clip, diff --git a/monai/transforms/croppad/dictionary.py b/monai/transforms/croppad/dictionary.py index d62d6fbdcb..c2033e3443 100644 --- a/monai/transforms/croppad/dictionary.py +++ b/monai/transforms/croppad/dictionary.py @@ -735,7 +735,7 @@ class RandSpatialCropSamplesd(Randomizable, MapTransform, InvertibleTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to add `patch_index` to the meta dict. allow_missing_keys: don't raise exception if key is missing. @@ -933,7 +933,7 @@ class RandWeightedCropd(Randomizable, MapTransform, InvertibleTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to add `patch_index` to the meta dict. allow_missing_keys: don't raise exception if key is missing. @@ -1075,7 +1075,7 @@ class RandCropByPosNegLabeld(Randomizable, MapTransform, InvertibleTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to add `patch_index` to the meta dict. allow_smaller: if `False`, an exception will be raised if the image is smaller than @@ -1279,7 +1279,7 @@ class RandCropByLabelClassesd(Randomizable, MapTransform, InvertibleTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to add `patch_index` to the meta dict. allow_smaller: if `False`, an exception will be raised if the image is smaller than diff --git a/monai/transforms/intensity/dictionary.py b/monai/transforms/intensity/dictionary.py index c66f768738..2d8f46ebd8 100644 --- a/monai/transforms/intensity/dictionary.py +++ b/monai/transforms/intensity/dictionary.py @@ -307,7 +307,7 @@ def __init__( the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to extract the factor value is `factor_key` is not None. allow_missing_keys: don't raise exception if key is missing. @@ -366,7 +366,7 @@ def __init__( the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to extract the factor value is `factor_key` is not None. prob: probability of rotating. diff --git a/monai/transforms/post/dictionary.py b/monai/transforms/post/dictionary.py index 96df434397..a7ffcb19bf 100644 --- a/monai/transforms/post/dictionary.py +++ b/monai/transforms/post/dictionary.py @@ -571,7 +571,7 @@ def __init__( it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `{orig_key}_{meta_key_postfix}`. meta data will also be inverted and stored in `meta_keys`. - meta_key_postfix: if `orig_meta_keys` is None, use `{orig_key}_{meta_key_postfix}` to to fetch the + meta_key_postfix: if `orig_meta_keys` is None, use `{orig_key}_{meta_key_postfix}` to fetch the meta data from dict, if `meta_keys` is None, use `{key}_{meta_key_postfix}`. default is `meta_dict`, the meta data is a dictionary object. For example, to handle orig_key `image`, read/write `affine` matrices from the diff --git a/monai/transforms/spatial/array.py b/monai/transforms/spatial/array.py index 60e4f564c8..206847df17 100644 --- a/monai/transforms/spatial/array.py +++ b/monai/transforms/spatial/array.py @@ -20,9 +20,9 @@ from monai.config import USE_COMPILED, DtypeLike from monai.config.type_definitions import NdarrayOrTensor -from monai.data.utils import compute_shape_offset, to_affine_nd, zoom_affine +from monai.data.utils import compute_shape_offset, reorient_spatial_axes, to_affine_nd, zoom_affine from monai.networks.layers import AffineTransform, GaussianFilter, grid_pull -from monai.networks.utils import meshgrid_ij +from monai.networks.utils import meshgrid_ij, normalize_transform from monai.transforms.croppad.array import CenterSpatialCrop, Pad from monai.transforms.transform import Randomizable, RandomizableTransform, ThreadUnsafe, Transform from monai.transforms.utils import ( @@ -34,6 +34,7 @@ create_translate, map_spatial_axes, ) +from monai.transforms.utils_pytorch_numpy_unification import allclose, moveaxis from monai.utils import ( GridSampleMode, GridSamplePadMode, @@ -46,15 +47,19 @@ fall_back_tuple, issequenceiterable, optional_import, + pytorch_after, ) from monai.utils.deprecate_utils import deprecated_arg from monai.utils.enums import TransformBackends from monai.utils.module import look_up_option from monai.utils.type_conversion import convert_data_type, convert_to_dst_type -nib, _ = optional_import("nibabel") +AFFINE_TOL = 1e-3 + +nib, has_nib = optional_import("nibabel") __all__ = [ + "SpatialResample", "Spacing", "Orientation", "Flip", @@ -82,12 +87,184 @@ RandRange = Optional[Union[Sequence[Union[Tuple[float, float], float]], float]] +class SpatialResample(Transform): + """ + Resample input image from the orientation/spacing defined by ``src_affine`` affine matrix into + the ones specified by ``dst_affine`` affine matrix. + + Internally this transform computes the affine transform matrix from ``src_affine`` to ``dst_affine``, + by ``xform = linalg.solve(src_affine, dst_affine)``, and call ``monai.transforms.Affine`` with ``xform``. + """ + + backend = [TransformBackends.TORCH] + + def __init__( + self, + mode: Union[GridSampleMode, str] = GridSampleMode.BILINEAR, + padding_mode: Union[GridSamplePadMode, str] = GridSamplePadMode.BORDER, + align_corners: bool = False, + dtype: DtypeLike = np.float64, + ): + """ + Args: + mode: {``"bilinear"``, ``"nearest"``} + Interpolation mode to calculate output values. Defaults to ``"bilinear"``. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull + padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} + Padding mode for outside grid values. Defaults to ``"border"``. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + align_corners: Geometrically, we consider the pixels of the input as squares rather than points. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + If ``None``, use the data type of input data. To be compatible with other modules, + the output data type is always ``np.float32``. + """ + self.mode = mode + self.padding_mode = padding_mode + self.align_corners = align_corners + self.dtype = dtype + + def __call__( + self, + img: NdarrayOrTensor, + src_affine: Optional[NdarrayOrTensor] = None, + dst_affine: Optional[NdarrayOrTensor] = None, + spatial_size: Optional[Union[Sequence[int], int]] = None, + mode: Union[GridSampleMode, str, None] = GridSampleMode.BILINEAR, + padding_mode: Union[GridSamplePadMode, str, None] = GridSamplePadMode.BORDER, + align_corners: Optional[bool] = False, + dtype: DtypeLike = None, + ) -> Tuple[NdarrayOrTensor, NdarrayOrTensor]: + """ + Args: + img: input image to be resampled. It currently supports channel-first arrays with + at most three spatial dimensions. + src_affine: source affine matrix. Defaults to ``None``, which means the identity matrix. + the shape should be `(r+1, r+1)` where `r` is the spatial rank of ``img``. + dst_affine: destination affine matrix. Defaults to ``None``, which means the same as `src_affine`. + the shape should be `(r+1, r+1)` where `r` is the spatial rank of ``img``. + when `dst` is None, the input will be returned without resampling, but the data type + will be `float32`. + spatial_size: output image spatial size. + if `spatial_size` and `self.spatial_size` are not defined, + the transform will compute a spatial size automatically containing the previous field of view. + if `spatial_size` is ``-1`` are the transform will use the corresponding input img size. + mode: {``"bilinear"``, ``"nearest"``} + Interpolation mode to calculate output values. Defaults to ``"bilinear"``. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull + padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} + Padding mode for outside grid values. Defaults to ``"border"``. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + align_corners: Geometrically, we consider the pixels of the input as squares rather than points. + See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + dtype: data type for resampling computation. Defaults to ``self.dtype`` or + ``np.float64`` (for best precision). If ``None``, use the data type of input data. + To be compatible with other modules, the output data type is always `float32`. + + The spatial rank is determined by the smallest among ``img.ndim -1``, ``len(src_affine) - 1``, and ``3``. + + When both ``monai.config.USE_COMPILED`` and ``align_corners`` are set to ``True``, + MONAI's resampling implementation will be used. + """ + if src_affine is None: + src_affine = np.eye(4, dtype=np.float64) + spatial_rank = min(len(img.shape) - 1, src_affine.shape[0] - 1, 3) + src_affine = to_affine_nd(spatial_rank, src_affine) + dst_affine = to_affine_nd(spatial_rank, dst_affine) if dst_affine is not None else src_affine + dst_affine, *_ = convert_to_dst_type(dst_affine, dst_affine, dtype=torch.float32) + + if allclose(src_affine, dst_affine, atol=AFFINE_TOL): + # no significant change, return original image + output_data, *_ = convert_to_dst_type(img, img, dtype=torch.float32) + return output_data, dst_affine + + if has_nib and isinstance(img, np.ndarray): + spatial_ornt, dst_r = reorient_spatial_axes(img.shape[1 : spatial_rank + 1], src_affine, dst_affine) + if allclose(dst_r, dst_affine, atol=AFFINE_TOL): + # simple reorientation achieves the desired affine + spatial_ornt[:, 0] += 1 + spatial_ornt = np.concatenate([np.array([[0, 1]]), spatial_ornt]) + img_ = nib.orientations.apply_orientation(img, spatial_ornt) + output_data, *_ = convert_to_dst_type(img_, img, dtype=torch.float32) + return output_data, dst_affine + + try: + src_affine, *_ = convert_to_dst_type(src_affine, dst_affine) + if isinstance(src_affine, np.ndarray): + xform = np.linalg.solve(src_affine, dst_affine) + else: + xform = ( + torch.linalg.solve(src_affine, dst_affine) + if pytorch_after(1, 8, 0) + else torch.solve(dst_affine, src_affine).solution # type: ignore + ) + except (np.linalg.LinAlgError, RuntimeError) as e: + raise ValueError(f"src affine is not invertible: {src_affine}") from e + xform = to_affine_nd(spatial_rank, xform) # type: ignore + # no resampling if it's identity transform + if allclose(xform, np.diag(np.ones(len(xform))), atol=AFFINE_TOL): + output_data, *_ = convert_to_dst_type(img, img, dtype=torch.float32) + return output_data, dst_affine + + _dtype = dtype or self.dtype or img.dtype + spatial_size = ensure_tuple(spatial_size) + in_spatial_size = list(img.shape[1 : spatial_rank + 1]) + if spatial_size[0] == -1: # if the spatial_size == -1 + spatial_size = in_spatial_size + elif spatial_size[0] is None: + spatial_size, _ = compute_shape_offset(in_spatial_size, src_affine, dst_affine) # type: ignore + spatial_size = spatial_size[:spatial_rank] + chns, additional_dims = img.shape[0], img.shape[spatial_rank + 1 :] # beyond three spatial dims + # resample + img_ = convert_data_type(img, torch.Tensor, dtype=_dtype)[0] # type: ignore + xform = convert_to_dst_type(xform, img_)[0] # type: ignore + align_corners = self.align_corners if align_corners is None else align_corners + mode = mode or self.mode + padding_mode = padding_mode or self.padding_mode + if additional_dims: + xform_shape = [-1] + in_spatial_size + img_ = img_.reshape(xform_shape) + if align_corners: + _t_r = torch.diag(torch.ones(len(xform), dtype=xform.dtype, device=xform.device)) # type: ignore + for idx, d_dst in enumerate(spatial_size[:spatial_rank]): + _t_r[idx, -1] = (max(d_dst, 2) - 1.0) / 2.0 + xform = xform @ _t_r + if not USE_COMPILED: + _t_l = normalize_transform(in_spatial_size, xform.device, xform.dtype, align_corners=True) # type: ignore + xform = _t_l @ xform # type: ignore + affine_xform = Affine( + affine=xform, spatial_size=spatial_size, norm_coords=False, image_only=True, dtype=_dtype + ) + output_data = affine_xform(img_, mode=mode, padding_mode=padding_mode) + else: + affine_xform = AffineTransform( + normalized=False, + mode=mode, # type: ignore + padding_mode=padding_mode, # type: ignore + align_corners=align_corners, + reverse_indexing=True, + ) + output_data = affine_xform(img_.unsqueeze(0), theta=xform, spatial_size=spatial_size).squeeze(0) + if additional_dims: + full_shape = (chns, *spatial_size, *additional_dims) + output_data = output_data.reshape(full_shape) + # output dtype float + output_data, *_ = convert_to_dst_type(output_data, img, dtype=torch.float32) + return output_data, dst_affine + + class Spacing(Transform): """ Resample input image into the specified `pixdim`. """ - backend = [TransformBackends.TORCH] + backend = SpatialResample.backend def __init__( self, @@ -122,6 +299,9 @@ def __init__( mode: {``"bilinear"``, ``"nearest"``} Interpolation mode to calculate output values. Defaults to ``"bilinear"``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``"border"``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html @@ -135,12 +315,15 @@ def __init__( """ self.pixdim = np.array(ensure_tuple(pixdim), dtype=np.float64) self.diagonal = diagonal - self.mode: GridSampleMode = look_up_option(mode, GridSampleMode) - self.padding_mode: GridSamplePadMode = look_up_option(padding_mode, GridSamplePadMode) - self.align_corners = align_corners - self.dtype = dtype self.image_only = image_only + self.sp_resample = SpatialResample( + mode=look_up_option(mode, GridSampleMode), + padding_mode=look_up_option(padding_mode, GridSamplePadMode), + align_corners=align_corners, + dtype=dtype, + ) + def __call__( self, data_array: NdarrayOrTensor, @@ -149,7 +332,7 @@ def __call__( padding_mode: Optional[Union[GridSamplePadMode, str]] = None, align_corners: Optional[bool] = None, dtype: DtypeLike = None, - output_spatial_shape: Optional[np.ndarray] = None, + output_spatial_shape: Optional[Union[Sequence[int], int]] = None, ) -> Union[NdarrayOrTensor, Tuple[NdarrayOrTensor, NdarrayOrTensor, NdarrayOrTensor]]: """ Args: @@ -158,6 +341,9 @@ def __call__( mode: {``"bilinear"``, ``"nearest"``} Interpolation mode to calculate output values. Defaults to ``self.mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``self.padding_mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html @@ -178,7 +364,6 @@ def __call__( data_array (resampled into `self.pixdim`), original affine, current affine. """ - _dtype = dtype or self.dtype or data_array.dtype sr = int(data_array.ndim - 1) if sr <= 0: raise ValueError("data_array must have at least one spatial dimension.") @@ -188,7 +373,7 @@ def __call__( affine_ = np.eye(sr + 1, dtype=np.float64) else: affine_np, *_ = convert_data_type(affine, np.ndarray) # type: ignore - affine_ = to_affine_nd(sr, affine_np) + affine_ = to_affine_nd(sr, affine_np) # type: ignore out_d = self.pixdim[:sr] if out_d.size < sr: @@ -198,32 +383,16 @@ def __call__( new_affine = zoom_affine(affine_, out_d, diagonal=self.diagonal) output_shape, offset = compute_shape_offset(data_array.shape[1:], affine_, new_affine) new_affine[:sr, -1] = offset[:sr] - transform = np.linalg.inv(affine_) @ new_affine - # adapt to the actual rank - transform = to_affine_nd(sr, transform) - - # no resampling if it's identity transform - if np.allclose(transform, np.diag(np.ones(len(transform))), atol=1e-3): - output_data = data_array - else: - # resample - affine_xform = AffineTransform( - normalized=False, - mode=look_up_option(mode or self.mode, GridSampleMode), - padding_mode=look_up_option(padding_mode or self.padding_mode, GridSamplePadMode), - align_corners=self.align_corners if align_corners is None else align_corners, - reverse_indexing=True, - ) - data_array_t: torch.Tensor - data_array_t, *_ = convert_data_type(data_array, torch.Tensor, dtype=_dtype) # type: ignore - output_data = affine_xform( - # AffineTransform requires a batch dim - data_array_t.unsqueeze(0), - convert_data_type(transform, torch.Tensor, data_array_t.device, dtype=_dtype)[0], - spatial_size=output_shape if output_spatial_shape is None else output_spatial_shape, - ).squeeze(0) - - output_data, *_ = convert_to_dst_type(output_data, data_array, dtype=torch.float32) + output_data, new_affine = self.sp_resample( + data_array, + src_affine=affine, + dst_affine=new_affine, + spatial_size=list(output_shape) if output_spatial_shape is None else output_spatial_shape, + mode=mode, + padding_mode=padding_mode, + align_corners=align_corners, + dtype=dtype, + ) new_affine = to_affine_nd(affine_np, new_affine) # type: ignore new_affine, *_ = convert_to_dst_type(src=new_affine, dst=affine, dtype=torch.float32) @@ -297,13 +466,14 @@ def __call__( sr = len(spatial_shape) if sr <= 0: raise ValueError("data_array must have at least one spatial dimension.") + affine_: np.ndarray if affine is None: # default to identity affine_np = affine = np.eye(sr + 1, dtype=np.float64) affine_ = np.eye(sr + 1, dtype=np.float64) else: affine_np, *_ = convert_data_type(affine, np.ndarray) # type: ignore - affine_ = to_affine_nd(sr, affine_np) + affine_ = to_affine_nd(sr, affine_np) # type: ignore src = nib.io_orientation(affine_) if self.as_closest_canonical: @@ -475,7 +645,7 @@ class Rotate(Transform, ThreadUnsafe): See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html align_corners: Defaults to False. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html - dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + dtype: data type for resampling computation. Defaults to ``np.float32``. If None, use the data type of input data. To be compatible with other modules, the output data type is always ``np.float32``. """ @@ -489,7 +659,7 @@ def __init__( mode: Union[GridSampleMode, str] = GridSampleMode.BILINEAR, padding_mode: Union[GridSamplePadMode, str] = GridSamplePadMode.BORDER, align_corners: bool = False, - dtype: Union[DtypeLike, torch.dtype] = np.float64, + dtype: Union[DtypeLike, torch.dtype] = np.float32, ) -> None: self.angle = angle self.keep_size = keep_size @@ -784,7 +954,7 @@ class RandRotate(RandomizableTransform): See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html align_corners: Defaults to False. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html - dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + dtype: data type for resampling computation. Defaults to ``np.float32``. If None, use the data type of input data. To be compatible with other modules, the output data type is always ``np.float32``. """ @@ -801,7 +971,7 @@ def __init__( mode: Union[GridSampleMode, str] = GridSampleMode.BILINEAR, padding_mode: Union[GridSamplePadMode, str] = GridSamplePadMode.BORDER, align_corners: bool = False, - dtype: Union[DtypeLike, torch.dtype] = np.float64, + dtype: Union[DtypeLike, torch.dtype] = np.float32, ) -> None: RandomizableTransform.__init__(self, prob) self.range_x = ensure_tuple(range_x) @@ -1085,6 +1255,9 @@ class AffineGrid(Transform): pixel/voxel relative to the center of the input image. Defaults to no translation. scale_params: scale factor for every spatial dims. a tuple of 2 floats for 2D, a tuple of 3 floats for 3D. Defaults to `1.0`. + dtype: data type for the grid computation. Defaults to ``np.float32``. + If ``None``, use the data type of input data (if `grid` is provided). + device: device on which the tensor will be allocated, if a new grid is generated. affine: If applied, ignore the params (`rotate_params`, etc.) and use the supplied matrix. Should be square with each side = num of image spatial dimensions + 1. @@ -1105,6 +1278,7 @@ def __init__( scale_params: Optional[Union[Sequence[float], float]] = None, as_tensor_output: bool = True, device: Optional[torch.device] = None, + dtype: DtypeLike = np.float32, affine: Optional[NdarrayOrTensor] = None, ) -> None: self.rotate_params = rotate_params @@ -1112,6 +1286,7 @@ def __init__( self.translate_params = translate_params self.scale_params = scale_params self.device = device + self.dtype = dtype self.affine = affine def __call__( @@ -1130,12 +1305,10 @@ def __call__( ValueError: When ``grid=None`` and ``spatial_size=None``. Incompatible values. """ - if grid is None: - if spatial_size is not None: - grid = create_grid(spatial_size, device=self.device, backend="torch") - else: + if grid is None: # create grid from spatial_size + if spatial_size is None: raise ValueError("Incompatible values: grid=None and spatial_size=None.") - + grid = create_grid(spatial_size, device=self.device, backend="torch", dtype=self.dtype) _b = TransformBackends.TORCH if isinstance(grid, torch.Tensor) else TransformBackends.NUMPY _device = grid.device if isinstance(grid, torch.Tensor) else self.device affine: NdarrayOrTensor @@ -1157,7 +1330,7 @@ def __call__( else: affine = self.affine - grid, *_ = convert_data_type(grid, torch.Tensor, device=_device, dtype=float) + grid, *_ = convert_data_type(grid, torch.Tensor, device=_device, dtype=self.dtype or grid.dtype) affine, *_ = convert_to_dst_type(affine, grid) grid = (affine @ grid.reshape((grid.shape[0], -1))).reshape([-1] + list(grid.shape[1:])) @@ -1339,7 +1512,9 @@ def __init__( mode: Union[GridSampleMode, str] = GridSampleMode.BILINEAR, padding_mode: Union[GridSamplePadMode, str] = GridSamplePadMode.BORDER, as_tensor_output: bool = True, + norm_coords: bool = True, device: Optional[torch.device] = None, + dtype: DtypeLike = np.float64, ) -> None: """ computes output image using values from `img`, locations from `grid` using pytorch. @@ -1352,7 +1527,17 @@ def __init__( padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``"border"``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull + norm_coords: whether to normalize the coordinates from `[-(size-1)/2, (size-1)/2]` to + `[0, size - 1]` (for ``monai/csrc`` implementation) or + `[-1, 1]` (for torch ``grid_sample`` implementation) to be compatible with the underlying + resampling API. device: device on which the tensor will be allocated. + dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + If ``None``, use the data type of input data. To be compatible with other modules, + the output data type is always `float32`. .. deprecated:: 0.6.0 ``as_tensor_output`` is deprecated. @@ -1360,7 +1545,9 @@ def __init__( """ self.mode: GridSampleMode = look_up_option(mode, GridSampleMode) self.padding_mode: GridSamplePadMode = look_up_option(padding_mode, GridSamplePadMode) + self.norm_coords = norm_coords self.device = device + self.dtype = dtype def __call__( self, @@ -1368,55 +1555,67 @@ def __call__( grid: Optional[NdarrayOrTensor] = None, mode: Optional[Union[GridSampleMode, str]] = None, padding_mode: Optional[Union[GridSamplePadMode, str]] = None, + dtype: DtypeLike = None, ) -> NdarrayOrTensor: """ Args: img: shape must be (num_channels, H, W[, D]). grid: shape must be (3, H, W) for 2D or (4, H, W, D) for 3D. + if ``norm_coords`` is True, the grid values must be in `[-(size-1)/2, (size-1)/2]`. + if ``USE_COMPILED=True`` and ``norm_coords=False``, grid values must be in `[0, size-1]`. + if ``USE_COMPILED=False`` and ``norm_coords=False``, grid values must be in `[-1, 1]`. mode: {``"bilinear"``, ``"nearest"``} Interpolation mode to calculate output values. Defaults to ``self.mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``self.padding_mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + dtype: data type for resampling computation. Defaults to ``self.dtype``. + To be compatible with other modules, the output data type is always `float32`. + + See also: + :py:const:`monai.config.USE_COMPILED` """ if grid is None: raise ValueError("Unknown grid.") _device = img.device if isinstance(img, torch.Tensor) else self.device img_t: torch.Tensor grid_t: torch.Tensor - img_t, *_ = convert_data_type(img, torch.Tensor, device=_device, dtype=torch.float32) # type: ignore - grid_t, *_ = convert_to_dst_type(grid, img_t) # type: ignore + _dtype = dtype or self.dtype or img.dtype + img_t, *_ = convert_data_type(img, torch.Tensor, device=_device, dtype=_dtype) # type: ignore + grid_t = convert_to_dst_type(grid, img_t)[0] # type: ignore + if grid_t is grid: # copy if needed (convert_data_type converts to contiguous) + grid_t = grid_t.clone(memory_format=torch.contiguous_format) + sr = min(len(img_t.shape[1:]), 3) if USE_COMPILED: - for i, dim in enumerate(img_t.shape[1:]): - grid_t[i] += (dim - 1.0) / 2.0 - grid_t = grid_t[:-1] / grid_t[-1:] - grid_t = grid_t.permute(list(range(grid_t.ndimension()))[1:] + [0]) - _padding_mode = look_up_option( - self.padding_mode if padding_mode is None else padding_mode, GridSamplePadMode - ).value - if _padding_mode == "zeros": - bound = 7 - elif _padding_mode == "border": - bound = 0 + if self.norm_coords: + for i, dim in enumerate(img_t.shape[1 : 1 + sr]): + grid_t[i] = (max(dim, 2) / 2.0 - 0.5 + grid_t[i]) / grid_t[-1:] + grid_t = moveaxis(grid_t[:sr], 0, -1) # type: ignore + _padding_mode = self.padding_mode if padding_mode is None else padding_mode + _padding_mode = _padding_mode.value if isinstance(_padding_mode, GridSamplePadMode) else _padding_mode + bound = 1 if _padding_mode == "reflection" else _padding_mode + _interp_mode = self.mode if mode is None else mode + _interp_mode = _interp_mode.value if isinstance(_interp_mode, GridSampleMode) else _interp_mode + if _interp_mode == "bicubic": + interp = 3 + elif _interp_mode == "bilinear": + interp = 1 else: - bound = 1 - _interp_mode = look_up_option(self.mode if mode is None else mode, GridSampleMode).value + interp = _interp_mode # type: ignore out = grid_pull( - img_t.unsqueeze(0), - grid_t.unsqueeze(0), - bound=bound, - extrapolate=True, - interpolation=1 if _interp_mode == "bilinear" else _interp_mode, + img_t.unsqueeze(0), grid_t.unsqueeze(0), bound=bound, extrapolate=True, interpolation=interp )[0] else: - for i, dim in enumerate(img_t.shape[1:]): - grid_t[i] = 2.0 * grid_t[i] / (dim - 1.0) - grid_t = grid_t[:-1] / grid_t[-1:] - index_ordering: List[int] = list(range(img_t.ndimension() - 2, -1, -1)) - grid_t = grid_t[index_ordering] - grid_t = grid_t.permute(list(range(grid_t.ndimension()))[1:] + [0]) + if self.norm_coords: + for i, dim in enumerate(img_t.shape[1 : 1 + sr]): + grid_t[i] = 2.0 / (max(2, dim) - 1.0) * grid_t[i] / grid_t[-1:] + index_ordering: List[int] = list(range(sr - 1, -1, -1)) + grid_t = moveaxis(grid_t[index_ordering], 0, -1) # type: ignore out = torch.nn.functional.grid_sample( img_t.unsqueeze(0), grid_t.unsqueeze(0), @@ -1425,7 +1624,7 @@ def __call__( align_corners=True, )[0] out_val: NdarrayOrTensor - out_val, *_ = convert_to_dst_type(out, dst=img, dtype=out.dtype) + out_val, *_ = convert_to_dst_type(out, dst=img, dtype=np.float32) return out_val @@ -1449,8 +1648,10 @@ def __init__( spatial_size: Optional[Union[Sequence[int], int]] = None, mode: Union[GridSampleMode, str] = GridSampleMode.BILINEAR, padding_mode: Union[GridSamplePadMode, str] = GridSamplePadMode.REFLECTION, + norm_coords: bool = True, as_tensor_output: bool = True, device: Optional[torch.device] = None, + dtype: DtypeLike = np.float32, image_only: bool = False, ) -> None: """ @@ -1485,10 +1686,21 @@ def __init__( mode: {``"bilinear"``, ``"nearest"``} Interpolation mode to calculate output values. Defaults to ``"bilinear"``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``"reflection"``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + norm_coords: whether to normalize the coordinates from `[-(size-1)/2, (size-1)/2]` to + `[0, size - 1]` or `[-1, 1]` to be compatible with the underlying resampling API. + If the coordinates are generated by ``monai.transforms.utils.create_grid`` + and the ``affine`` doesn't include the normalization, this argument should be set to ``True``. + If the output `self.affine_grid` is already normalized, this argument should be set to ``False``. device: device on which the tensor will be allocated. + dtype: data type for resampling computation. Defaults to ``np.float32``. + If ``None``, use the data type of input data. To be compatible with other modules, + the output data type is always `float32`. image_only: if True return only the image volume, otherwise return (image, affine). .. deprecated:: 0.6.0 @@ -1501,10 +1713,11 @@ def __init__( translate_params=translate_params, scale_params=scale_params, affine=affine, + dtype=dtype, device=device, ) self.image_only = image_only - self.resampler = Resample(device=device) + self.resampler = Resample(norm_coords=norm_coords, device=device, dtype=dtype) self.spatial_size = spatial_size self.mode: GridSampleMode = look_up_option(mode, GridSampleMode) self.padding_mode: GridSamplePadMode = look_up_option(padding_mode, GridSamplePadMode) @@ -1527,6 +1740,9 @@ def __call__( mode: {``"bilinear"``, ``"nearest"``} Interpolation mode to calculate output values. Defaults to ``self.mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html + When `USE_COMPILED` is `True`, this argument uses + ``"nearest"``, ``"bilinear"``, ``"bicubic"`` to indicate 0, 1, 3 order interpolations. + See also: https://docs.monai.io/en/stable/networks.html#grid-pull padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} Padding mode for outside grid values. Defaults to ``self.padding_mode``. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html diff --git a/monai/transforms/spatial/dictionary.py b/monai/transforms/spatial/dictionary.py index 963015f420..b319bb44b1 100644 --- a/monai/transforms/spatial/dictionary.py +++ b/monai/transforms/spatial/dictionary.py @@ -46,6 +46,7 @@ Rotate, Rotate90, Spacing, + SpatialResample, Zoom, ) from monai.transforms.transform import MapTransform, RandomizableTransform @@ -68,6 +69,7 @@ nib, _ = optional_import("nibabel") __all__ = [ + "SpatialResampled", "Spacingd", "Orientationd", "Rotate90d", @@ -86,6 +88,8 @@ "RandRotated", "Zoomd", "RandZoomd", + "SpatialResampleD", + "SpatialResampleDict", "SpacingD", "SpacingDict", "OrientationD", @@ -131,6 +135,160 @@ DEFAULT_POST_FIX = PostFix.meta() +class SpatialResampled(MapTransform, InvertibleTransform): + """ + Dictionary-based wrapper of :py:class:`monai.transforms.SpatialResample`. + + This transform assumes the ``data`` dictionary has a key for the input + data's metadata and contains ``src_affine`` and ``dst_affine`` required by + `SpatialResample`. The key is formed by ``key_{meta_key_postfix}``. The + transform will swap ``src_affine`` and ``dst_affine`` affine (with potential data type + changes) in the dictionary so that ``src_affine`` always refers to the current + status of affine. + + See also: + :py:class:`monai.transforms.SpatialResample` + """ + + backend = SpatialResample.backend + + def __init__( + self, + keys: KeysCollection, + mode: GridSampleModeSequence = GridSampleMode.BILINEAR, + padding_mode: GridSamplePadModeSequence = GridSamplePadMode.BORDER, + align_corners: Union[Sequence[bool], bool] = False, + dtype: Optional[Union[Sequence[DtypeLike], DtypeLike]] = np.float64, + meta_keys: Optional[KeysCollection] = None, + meta_key_postfix: str = DEFAULT_POST_FIX, + meta_src_keys: Optional[KeysCollection] = "src_affine", + meta_dst_keys: Optional[KeysCollection] = "dst_affine", + allow_missing_keys: bool = False, + ) -> None: + """ + Args: + keys: keys of the corresponding items to be transformed. + mode: {``"bilinear"``, ``"nearest"``} + Interpolation mode to calculate output values. Defaults to ``"bilinear"``. + See also: https://pytorch.org/docs/stable/nn.functional.html#grid-sample + It also can be a sequence of string, each element corresponds to a key in ``keys``. + padding_mode: {``"zeros"``, ``"border"``, ``"reflection"``} + Padding mode for outside grid values. Defaults to ``"border"``. + See also: https://pytorch.org/docs/stable/nn.functional.html#grid-sample + It also can be a sequence of string, each element corresponds to a key in ``keys``. + align_corners: Geometrically, we consider the pixels of the input as squares rather than points. + See also: https://pytorch.org/docs/stable/nn.functional.html#grid-sample + It also can be a sequence of bool, each element corresponds to a key in ``keys``. + dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + If None, use the data type of input data. To be compatible with other modules, + the output data type is always ``np.float32``. + It also can be a sequence of dtypes, each element corresponds to a key in ``keys``. + meta_keys: explicitly indicate the key of the corresponding meta data dictionary. + for example, for data with key `image`, the metadata by default is in `image_meta_dict`. + the meta data is a dictionary object which contains: filename, affine, original_shape, etc. + it can be a sequence of string, map to the `keys`. + if None, will try to construct meta_keys by `key_{meta_key_postfix}`. + meta_key_postfix: if meta_keys=None, use `key_{postfix}` to fetch the meta data according + to the key data, default is `meta_dict`, the meta data is a dictionary object. + For example, to handle key `image`, read/write affine matrices from the + metadata `image_meta_dict` dictionary's `affine` field. + meta_src_keys: the key of the corresponding ``src_affine`` in the metadata dictionary. + meta_dst_keys: the key of the corresponding ``dst_affine`` in the metadata dictionary. + allow_missing_keys: don't raise exception if key is missing. + """ + super().__init__(keys, allow_missing_keys) + self.sp_transform = SpatialResample() + self.mode = ensure_tuple_rep(mode, len(self.keys)) + self.padding_mode = ensure_tuple_rep(padding_mode, len(self.keys)) + self.align_corners = ensure_tuple_rep(align_corners, len(self.keys)) + self.dtype = ensure_tuple_rep(dtype, len(self.keys)) + self.meta_keys = ensure_tuple_rep(None, len(self.keys)) if meta_keys is None else ensure_tuple(meta_keys) + if len(self.keys) != len(self.meta_keys): + raise ValueError("meta_keys should have the same length as keys.") + self.meta_key_postfix = ensure_tuple_rep(meta_key_postfix, len(self.keys)) + self.meta_src_keys = ensure_tuple_rep(meta_src_keys, len(self.keys)) + self.meta_dst_keys = ensure_tuple_rep(meta_dst_keys, len(self.keys)) + + def __call__( + self, data: Mapping[Union[Hashable, str], Dict[str, NdarrayOrTensor]] + ) -> Dict[Hashable, NdarrayOrTensor]: + d: Dict = dict(data) + for (key, mode, padding_mode, align_corners, dtype, *metakeyinfo) in self.key_iterator( + d, + self.mode, + self.padding_mode, + self.align_corners, + self.dtype, + self.meta_keys, + self.meta_key_postfix, + self.meta_src_keys, + self.meta_dst_keys, + ): + meta_key, meta_key_postfix, meta_src_key, meta_dst_key = metakeyinfo + meta_key = meta_key or f"{key}_{meta_key_postfix}" + # create metadata if necessary + if meta_key not in d: + d[meta_key] = {meta_src_key: None, meta_dst_key: None} + meta_data = d[meta_key] + original_spatial_shape = d[key].shape[1:] + d[key], meta_data[meta_dst_key] = self.sp_transform( # write dst affine because the dtype might change + img=d[key], + src_affine=meta_data[meta_src_key], + dst_affine=meta_data[meta_dst_key], + spatial_size=None, # None means shape auto inferred + mode=mode, + padding_mode=padding_mode, + align_corners=align_corners, + dtype=dtype, + ) + meta_data[meta_dst_key], meta_data[meta_src_key] = meta_data[meta_src_key], meta_data[meta_dst_key] + self.push_transform( + d, + key, + extra_info={ + "meta_key": meta_key, + "meta_src_key": meta_src_key, + "meta_dst_key": meta_dst_key, + "mode": mode.value if isinstance(mode, Enum) else mode, + "padding_mode": padding_mode.value if isinstance(padding_mode, Enum) else padding_mode, + "align_corners": align_corners if align_corners is not None else TraceKeys.NONE, + }, + orig_size=original_spatial_shape, + ) + return d + + def inverse(self, data: Mapping[Hashable, NdarrayOrTensor]) -> Dict[Hashable, NdarrayOrTensor]: + d = deepcopy(dict(data)) + for key, dtype in self.key_iterator(d, self.dtype): + transform = self.get_most_recent_transform(d, key) + # Create inverse transform + meta_data = d[transform[TraceKeys.EXTRA_INFO]["meta_key"]] + src_key = transform[TraceKeys.EXTRA_INFO]["meta_src_key"] + dst_key = transform[TraceKeys.EXTRA_INFO]["meta_dst_key"] + src_affine = meta_data[src_key] # type: ignore + dst_affine = meta_data[dst_key] # type: ignore + mode = transform[TraceKeys.EXTRA_INFO]["mode"] + padding_mode = transform[TraceKeys.EXTRA_INFO]["padding_mode"] + align_corners = transform[TraceKeys.EXTRA_INFO]["align_corners"] + orig_size = transform[TraceKeys.ORIG_SIZE] + inverse_transform = SpatialResample() + # Apply inverse + d[key], dst_affine = inverse_transform( + img=d[key], + src_affine=src_affine, + dst_affine=dst_affine, + mode=mode, + padding_mode=padding_mode, + align_corners=False if align_corners == TraceKeys.NONE else align_corners, + dtype=dtype, + spatial_size=orig_size, + ) + meta_data[src_key], meta_data[dst_key] = dst_affine, meta_data[src_key] # type: ignore + # Remove the applied transform + self.pop_transform(d, key) + return d + + class Spacingd(MapTransform, InvertibleTransform): """ Dictionary-based wrapper of :py:class:`monai.transforms.Spacing`. @@ -201,7 +359,7 @@ def __init__( the meta data is a dictionary object which contains: filename, affine, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys=None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys=None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -335,7 +493,7 @@ def __init__( the meta data is a dictionary object which contains: filename, affine, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. For example, to handle key `image`, read/write affine matrices from the metadata `image_meta_dict` dictionary's `affine` field. @@ -607,6 +765,7 @@ def __init__( padding_mode: GridSamplePadModeSequence = GridSamplePadMode.REFLECTION, as_tensor_output: bool = True, device: Optional[torch.device] = None, + dtype: Union[DtypeLike, torch.dtype] = np.float32, allow_missing_keys: bool = False, ) -> None: """ @@ -646,6 +805,9 @@ def __init__( See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html It also can be a sequence of string, each element corresponds to a key in ``keys``. device: device on which the tensor will be allocated. + dtype: data type for resampling computation. Defaults to ``np.float32``. + If ``None``, use the data type of input data. To be compatible with other modules, + the output data type is always `float32`. allow_missing_keys: don't raise exception if key is missing. See also: @@ -665,6 +827,7 @@ def __init__( affine=affine, spatial_size=spatial_size, device=device, + dtype=dtype, ) self.mode = ensure_tuple_rep(mode, len(self.keys)) self.padding_mode = ensure_tuple_rep(padding_mode, len(self.keys)) @@ -1320,7 +1483,7 @@ class Rotated(MapTransform, InvertibleTransform): align_corners: Defaults to False. See also: https://pytorch.org/docs/stable/generated/torch.nn.functional.grid_sample.html It also can be a sequence of bool, each element corresponds to a key in ``keys``. - dtype: data type for resampling computation. Defaults to ``np.float64`` for best precision. + dtype: data type for resampling computation. Defaults to ``np.float32``. If None, use the data type of input data. To be compatible with other modules, the output data type is always ``np.float32``. It also can be a sequence of dtype or None, each element corresponds to a key in ``keys``. @@ -1337,7 +1500,7 @@ def __init__( mode: GridSampleModeSequence = GridSampleMode.BILINEAR, padding_mode: GridSamplePadModeSequence = GridSamplePadMode.BORDER, align_corners: Union[Sequence[bool], bool] = False, - dtype: Union[Sequence[Union[DtypeLike, torch.dtype]], Union[DtypeLike, torch.dtype]] = np.float64, + dtype: Union[Sequence[Union[DtypeLike, torch.dtype]], DtypeLike, torch.dtype] = np.float32, allow_missing_keys: bool = False, ) -> None: super().__init__(keys, allow_missing_keys) @@ -1451,7 +1614,7 @@ def __init__( mode: GridSampleModeSequence = GridSampleMode.BILINEAR, padding_mode: GridSamplePadModeSequence = GridSamplePadMode.BORDER, align_corners: Union[Sequence[bool], bool] = False, - dtype: Union[Sequence[Union[DtypeLike, torch.dtype]], Union[DtypeLike, torch.dtype]] = np.float64, + dtype: Union[Sequence[Union[DtypeLike, torch.dtype]], DtypeLike, torch.dtype] = np.float32, allow_missing_keys: bool = False, ) -> None: MapTransform.__init__(self, keys, allow_missing_keys) @@ -1870,6 +2033,7 @@ def __call__(self, data: Mapping[Hashable, NdarrayOrTensor]) -> Dict[Hashable, N return d +SpatialResampleD = SpatialResampleDict = SpatialResampled SpacingD = SpacingDict = Spacingd OrientationD = OrientationDict = Orientationd Rotate90D = Rotate90Dict = Rotate90d diff --git a/monai/transforms/utility/dictionary.py b/monai/transforms/utility/dictionary.py index e61aa5f512..96e2ea6361 100644 --- a/monai/transforms/utility/dictionary.py +++ b/monai/transforms/utility/dictionary.py @@ -1449,7 +1449,7 @@ class IntensityStatsd(MapTransform): the meta data is a dictionary object which contains: filename, original_shape, etc. it can be a sequence of string, map to the `keys`. if None, will try to construct meta_keys by `key_{meta_key_postfix}`. - meta_key_postfix: if meta_keys is None, use `key_{postfix}` to to fetch the meta data according + meta_key_postfix: if meta_keys is None, use `key_{postfix}` to fetch the meta data according to the key data, default is `meta_dict`, the meta data is a dictionary object. used to store the computed statistics to the meta dict. allow_missing_keys: don't raise exception if key is missing. diff --git a/monai/transforms/utils.py b/monai/transforms/utils.py index a53a58f520..a58d80c6f1 100644 --- a/monai/transforms/utils.py +++ b/monai/transforms/utils.py @@ -48,6 +48,7 @@ ensure_tuple_rep, ensure_tuple_size, fall_back_tuple, + get_equivalent_dtype, issequenceiterable, look_up_option, min_version, @@ -569,7 +570,7 @@ def create_grid( spatial_size: Sequence[int], spacing: Optional[Sequence[float]] = None, homogeneous: bool = True, - dtype=float, + dtype: Union[DtypeLike, torch.dtype] = float, device: Optional[torch.device] = None, backend=TransformBackends.NUMPY, ): @@ -580,16 +581,17 @@ def create_grid( spatial_size: spatial size of the grid. spacing: same len as ``spatial_size``, defaults to 1.0 (dense grid). homogeneous: whether to make homogeneous coordinates. - dtype: output grid data type. + dtype: output grid data type, defaults to `float`. device: device to compute and store the output (when the backend is "torch"). backend: APIs to use, ``numpy`` or ``torch``. """ _backend = look_up_option(backend, TransformBackends) + _dtype = dtype or float if _backend == TransformBackends.NUMPY: - return _create_grid_numpy(spatial_size, spacing, homogeneous, dtype) + return _create_grid_numpy(spatial_size, spacing, homogeneous, _dtype) if _backend == TransformBackends.TORCH: - return _create_grid_torch(spatial_size, spacing, homogeneous, dtype, device) + return _create_grid_torch(spatial_size, spacing, homogeneous, _dtype, device) raise ValueError(f"backend {backend} is not supported") @@ -597,14 +599,14 @@ def _create_grid_numpy( spatial_size: Sequence[int], spacing: Optional[Sequence[float]] = None, homogeneous: bool = True, - dtype: DtypeLike = float, + dtype: Union[DtypeLike, torch.dtype] = float, ): """ compute a `spatial_size` mesh with the numpy API. """ spacing = spacing or tuple(1.0 for _ in spatial_size) ranges = [np.linspace(-(d - 1.0) / 2.0 * s, (d - 1.0) / 2.0 * s, int(d)) for d, s in zip(spatial_size, spacing)] - coords = np.asarray(np.meshgrid(*ranges, indexing="ij"), dtype=dtype) + coords = np.asarray(np.meshgrid(*ranges, indexing="ij"), dtype=get_equivalent_dtype(dtype, np.ndarray)) if not homogeneous: return coords return np.concatenate([coords, np.ones_like(coords[:1])]) @@ -622,7 +624,13 @@ def _create_grid_torch( """ spacing = spacing or tuple(1.0 for _ in spatial_size) ranges = [ - torch.linspace(-(d - 1.0) / 2.0 * s, (d - 1.0) / 2.0 * s, int(d), device=device, dtype=dtype) + torch.linspace( + -(d - 1.0) / 2.0 * s, + (d - 1.0) / 2.0 * s, + int(d), + device=device, + dtype=get_equivalent_dtype(dtype, torch.Tensor), + ) for d, s in zip(spatial_size, spacing) ] coords = meshgrid_ij(*ranges) diff --git a/monai/transforms/utils_pytorch_numpy_unification.py b/monai/transforms/utils_pytorch_numpy_unification.py index 99d348b5df..4f93f4a982 100644 --- a/monai/transforms/utils_pytorch_numpy_unification.py +++ b/monai/transforms/utils_pytorch_numpy_unification.py @@ -19,6 +19,7 @@ from monai.utils.type_conversion import convert_data_type, convert_to_dst_type __all__ = [ + "allclose", "moveaxis", "in1d", "clip", @@ -43,6 +44,14 @@ ] +def allclose(a: NdarrayOrTensor, b: NdarrayOrTensor, rtol=1e-5, atol=1e-8, equal_nan=False) -> bool: + """`np.allclose` with equivalent implementation for torch.""" + b, *_ = convert_to_dst_type(b, a) + if isinstance(a, np.ndarray): + return np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan) + return torch.allclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan) # type: ignore + + def moveaxis(x: NdarrayOrTensor, src: Union[int, Sequence[int]], dst: Union[int, Sequence[int]]) -> NdarrayOrTensor: """`moveaxis` for pytorch and numpy, using `permute` for pytorch version < 1.7""" if isinstance(x, torch.Tensor): @@ -358,7 +367,7 @@ def isnan(x: NdarrayOrTensor) -> NdarrayOrTensor: Args: x: array/tensor - dim: dimension along which to stack + """ if isinstance(x, np.ndarray): return np.isnan(x) diff --git a/monai/utils/module.py b/monai/utils/module.py index c0fc10a7c0..8813301d8e 100644 --- a/monai/utils/module.py +++ b/monai/utils/module.py @@ -19,7 +19,7 @@ from pkgutil import walk_packages from re import match from types import FunctionType -from typing import Any, Callable, Collection, Hashable, Iterable, List, Mapping, Tuple, cast +from typing import Any, Callable, Collection, Hashable, Iterable, List, Mapping, Tuple, Union, cast import torch @@ -44,7 +44,7 @@ ] -def look_up_option(opt_str, supported: Collection, default="no_default"): +def look_up_option(opt_str, supported: Union[Collection, enum.EnumMeta], default="no_default"): """ Look up the option in the supported collection and return the matched item. Raise a value error possibly with a guess of the closest match. diff --git a/tests/test_affined.py b/tests/test_affined.py index 355833b858..665c93d23f 100644 --- a/tests/test_affined.py +++ b/tests/test_affined.py @@ -28,6 +28,13 @@ p(np.arange(9).reshape(1, 3, 3)), ] ) + TESTS.append( + [ + dict(keys="img", padding_mode="zeros", spatial_size=(-1, 0), device=device, dtype=None), + {"img": p(np.arange(9, dtype=float).reshape((1, 3, 3)))}, + p(np.arange(9).reshape(1, 3, 3)), + ] + ) TESTS.append( [ dict(keys="img", padding_mode="zeros", device=device), diff --git a/tests/test_global_mutual_information_loss.py b/tests/test_global_mutual_information_loss.py index a919396e3e..6aa5455352 100644 --- a/tests/test_global_mutual_information_loss.py +++ b/tests/test_global_mutual_information_loss.py @@ -16,7 +16,7 @@ from monai import transforms from monai.losses.image_dissimilarity import GlobalMutualInformationLoss -from tests.utils import SkipIfBeforePyTorchVersion, download_url_or_skip_test +from tests.utils import SkipIfBeforePyTorchVersion, download_url_or_skip_test, skip_if_quick device = "cuda" if torch.cuda.is_available() else "cpu" @@ -51,6 +51,7 @@ } +@skip_if_quick class TestGlobalMutualInformationLoss(unittest.TestCase): def setUp(self): download_url_or_skip_test(FILE_URL, FILE_PATH) @@ -100,6 +101,8 @@ def transformation(translate_params=(0.0, 0.0, 0.0), rotate_params=(0.0, 0.0, 0. result = loss_fn(a2, a1).detach().cpu().numpy() np.testing.assert_allclose(result, expected_value, rtol=1e-3, atol=5e-3) + +class TestGlobalMutualInformationLossIll(unittest.TestCase): def test_ill_shape(self): loss = GlobalMutualInformationLoss() with self.assertRaisesRegex(ValueError, ""): diff --git a/tests/test_integration_classification_2d.py b/tests/test_integration_classification_2d.py index 3572678b64..0e79a26ea7 100644 --- a/tests/test_integration_classification_2d.py +++ b/tests/test_integration_classification_2d.py @@ -68,7 +68,7 @@ def run_training_test(root_dir, train_x, train_y, val_x, val_y, device="cuda:0", AddChannel(), Transpose(indices=[0, 2, 1]), ScaleIntensity(), - RandRotate(range_x=np.pi / 12, prob=0.5, keep_size=True), + RandRotate(range_x=np.pi / 12, prob=0.5, keep_size=True, dtype=np.float64), RandFlip(spatial_axis=0, prob=0.5), RandZoom(min_zoom=0.9, max_zoom=1.1, prob=0.5), ToTensor(), diff --git a/tests/test_integration_fast_train.py b/tests/test_integration_fast_train.py index 3522a57342..51b2ac1d3f 100644 --- a/tests/test_integration_fast_train.py +++ b/tests/test_integration_fast_train.py @@ -123,6 +123,7 @@ def test_train_timing(self): range_x=np.pi / 4, mode=("bilinear", "nearest"), align_corners=True, + dtype=np.float64, ), RandAffined(keys=["image", "label"], prob=0.5, rotate_range=np.pi / 2, mode=("bilinear", "nearest")), RandGaussianNoised(keys="image", prob=0.5), diff --git a/tests/test_inverse.py b/tests/test_inverse.py index 4455009658..8696d974ba 100644 --- a/tests/test_inverse.py +++ b/tests/test_inverse.py @@ -209,10 +209,15 @@ TESTS.append(("RandZoom 3d", "3D", 9e-2, RandZoomd(KEYS, 1, [0.5, 0.6, 0.9], [1.1, 1, 1.05], keep_size=True))) -TESTS.append(("RandRotated, prob 0", "2D", 0, RandRotated(KEYS, prob=0))) +TESTS.append(("RandRotated, prob 0", "2D", 0, RandRotated(KEYS, prob=0, dtype=np.float64))) TESTS.append( - ("Rotated 2d", "2D", 8e-2, Rotated(KEYS, random.uniform(np.pi / 6, np.pi), keep_size=True, align_corners=False)) + ( + "Rotated 2d", + "2D", + 8e-2, + Rotated(KEYS, random.uniform(np.pi / 6, np.pi), keep_size=True, align_corners=False, dtype=np.float64), + ) ) TESTS.append( @@ -220,7 +225,7 @@ "Rotated 3d", "3D", 1e-1, - Rotated(KEYS, [random.uniform(np.pi / 6, np.pi) for _ in range(3)], True), # type: ignore + Rotated(KEYS, [random.uniform(np.pi / 6, np.pi) for _ in range(3)], True, dtype=np.float64), # type: ignore ) ) @@ -229,7 +234,7 @@ "RandRotated 3d", "3D", 1e-1, - RandRotated(KEYS, *[random.uniform(np.pi / 6, np.pi) for _ in range(3)], 1), # type: ignore + RandRotated(KEYS, *[random.uniform(np.pi / 6, np.pi) for _ in range(3)], 1, dtype=np.float64), # type: ignore ) ) diff --git a/tests/test_inverse_collation.py b/tests/test_inverse_collation.py index 3c293070bb..4e8c6b58cc 100644 --- a/tests/test_inverse_collation.py +++ b/tests/test_inverse_collation.py @@ -50,7 +50,7 @@ RandAxisFlipd(keys=KEYS, prob=0.5), Compose([RandRotate90d(keys=KEYS, spatial_axes=(1, 2)), ToTensord(keys=KEYS)]), RandZoomd(keys=KEYS, prob=0.5, min_zoom=0.5, max_zoom=1.1, keep_size=True), - RandRotated(keys=KEYS, prob=0.5, range_x=np.pi), + RandRotated(keys=KEYS, prob=0.5, range_x=np.pi, dtype=np.float64), RandAffined( keys=KEYS, prob=0.5, rotate_range=np.pi, device=torch.device("cuda" if torch.cuda.is_available() else "cpu") ), @@ -65,7 +65,7 @@ RandAxisFlipd(keys=KEYS, prob=0.5), Compose([RandRotate90d(keys=KEYS, prob=0.5, spatial_axes=(0, 1)), ToTensord(keys=KEYS)]), RandZoomd(keys=KEYS, prob=0.5, min_zoom=0.5, max_zoom=1.1, keep_size=True), - RandRotated(keys=KEYS, prob=0.5, range_x=np.pi), + RandRotated(keys=KEYS, prob=0.5, range_x=np.pi, dtype=np.float64), RandAffined( keys=KEYS, prob=0.5, rotate_range=np.pi, device=torch.device("cuda" if torch.cuda.is_available() else "cpu") ), diff --git a/tests/test_invertd.py b/tests/test_invertd.py index f28587eb6b..64c26c4012 100644 --- a/tests/test_invertd.py +++ b/tests/test_invertd.py @@ -58,7 +58,7 @@ def test_invert(self): RandAxisFlipd(KEYS, prob=0.5), RandRotate90d(KEYS, spatial_axes=(1, 2)), RandZoomd(KEYS, prob=0.5, min_zoom=0.5, max_zoom=1.1, keep_size=True), - RandRotated(KEYS, prob=0.5, range_x=np.pi, mode="bilinear", align_corners=True), + RandRotated(KEYS, prob=0.5, range_x=np.pi, mode="bilinear", align_corners=True, dtype=np.float64), RandAffined(KEYS, prob=0.5, rotate_range=np.pi, mode="nearest"), ResizeWithPadOrCropd(KEYS, 100), # test EnsureTensor for complicated dict data and invert it diff --git a/tests/test_pad_collation.py b/tests/test_pad_collation.py index a070fc760f..530e5f86a3 100644 --- a/tests/test_pad_collation.py +++ b/tests/test_pad_collation.py @@ -42,12 +42,12 @@ PadListDataCollate(method="end", mode="constant"), ]: TESTS.append((dict, pad_collate, RandSpatialCropd("image", roi_size=[8, 7], random_size=True))) - TESTS.append((dict, pad_collate, RandRotated("image", prob=1, range_x=np.pi, keep_size=False))) + TESTS.append((dict, pad_collate, RandRotated("image", prob=1, range_x=np.pi, keep_size=False, dtype=np.float64))) TESTS.append((dict, pad_collate, RandZoomd("image", prob=1, min_zoom=1.1, max_zoom=2.0, keep_size=False))) TESTS.append((dict, pad_collate, Compose([RandRotate90d("image", prob=1, max_k=2), ToTensord("image")]))) TESTS.append((list, pad_collate, RandSpatialCrop(roi_size=[8, 7], random_size=True))) - TESTS.append((list, pad_collate, RandRotate(prob=1, range_x=np.pi, keep_size=False))) + TESTS.append((list, pad_collate, RandRotate(prob=1, range_x=np.pi, keep_size=False, dtype=np.float64))) TESTS.append((list, pad_collate, RandZoom(prob=1, min_zoom=1.1, max_zoom=2.0, keep_size=False))) TESTS.append((list, pad_collate, Compose([RandRotate90(prob=1, max_k=2), ToTensor()]))) diff --git a/tests/test_rand_elastic_3d.py b/tests/test_rand_elastic_3d.py index d4eed3753f..39ce779cb0 100644 --- a/tests/test_rand_elastic_3d.py +++ b/tests/test_rand_elastic_3d.py @@ -88,7 +88,7 @@ def test_rand_3d_elastic(self, input_param, input_data, expected_val): g = Rand3DElastic(**input_param) g.set_random_state(123) result = g(**input_data) - assert_allclose(result, expected_val, rtol=1e-4, atol=1e-4) + assert_allclose(result, expected_val, rtol=1e-1, atol=1e-1) if __name__ == "__main__": diff --git a/tests/test_rand_elasticd_3d.py b/tests/test_rand_elasticd_3d.py index 84bc765ab0..c78ed1f42e 100644 --- a/tests/test_rand_elasticd_3d.py +++ b/tests/test_rand_elasticd_3d.py @@ -145,7 +145,7 @@ def test_rand_3d_elasticd(self, input_param, input_data, expected_val): for key in res: result = res[key] expected = expected_val[key] if isinstance(expected_val, dict) else expected_val - assert_allclose(result, expected, rtol=1e-4, atol=1e-4) + assert_allclose(result, expected, rtol=1e-2, atol=1e-2) if __name__ == "__main__": diff --git a/tests/test_rand_rotate.py b/tests/test_rand_rotate.py index b453b01884..7a85fce23b 100644 --- a/tests/test_rand_rotate.py +++ b/tests/test_rand_rotate.py @@ -71,6 +71,7 @@ def test_correct_results(self, im_type, degrees, keep_size, mode, padding_mode, mode=mode, padding_mode=padding_mode, align_corners=align_corners, + dtype=np.float64, ) rotate_fn.set_random_state(243) rotated = rotate_fn(im_type(self.imt[0])) @@ -104,6 +105,7 @@ def test_correct_results(self, im_type, x, y, z, keep_size, mode, padding_mode, mode=mode, padding_mode=padding_mode, align_corners=align_corners, + dtype=np.float64, ) rotate_fn.set_random_state(243) rotated = rotate_fn(im_type(self.imt[0])) diff --git a/tests/test_rand_rotated.py b/tests/test_rand_rotated.py index 23314720c1..464b37d925 100644 --- a/tests/test_rand_rotated.py +++ b/tests/test_rand_rotated.py @@ -116,6 +116,7 @@ def test_correct_results(self, im_type, degrees, keep_size, mode, padding_mode, mode=mode, padding_mode=padding_mode, align_corners=align_corners, + dtype=np.float64, ) rotate_fn.set_random_state(243) rotated = rotate_fn({"img": im_type(self.imt[0]), "seg": im_type(self.segn[0])}) @@ -151,6 +152,7 @@ def test_correct_shapes(self, im_type, x, y, z, keep_size, mode, padding_mode, a mode=mode, padding_mode=padding_mode, align_corners=align_corners, + dtype=np.float64, ) rotate_fn.set_random_state(243) rotated = rotate_fn({"img": im_type(self.imt[0]), "seg": im_type(self.segn[0])}) diff --git a/tests/test_rotate.py b/tests/test_rotate.py index 42947e7f72..01842f6d73 100644 --- a/tests/test_rotate.py +++ b/tests/test_rotate.py @@ -46,7 +46,7 @@ class TestRotate2D(NumpyImageTestCase2D): @parameterized.expand(TEST_CASES_2D) def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, align_corners): - rotate_fn = Rotate(angle, keep_size, mode, padding_mode, align_corners) + rotate_fn = Rotate(angle, keep_size, mode, padding_mode, align_corners, dtype=np.float64) rotated = rotate_fn(im_type(self.imt[0])) if keep_size: np.testing.assert_allclose(self.imt[0].shape, rotated.shape) @@ -74,7 +74,7 @@ def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, al class TestRotate3D(NumpyImageTestCase3D): @parameterized.expand(TEST_CASES_3D) def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, align_corners): - rotate_fn = Rotate([angle, 0, 0], keep_size, mode, padding_mode, align_corners) + rotate_fn = Rotate([angle, 0, 0], keep_size, mode, padding_mode, align_corners, dtype=np.float64) rotated = rotate_fn(im_type(self.imt[0])) if keep_size: np.testing.assert_allclose(self.imt[0].shape, rotated.shape) @@ -100,7 +100,7 @@ def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, al @parameterized.expand(TEST_CASES_SHAPE_3D) def test_correct_shape(self, im_type, angle, mode, padding_mode, align_corners): - rotate_fn = Rotate(angle, True, align_corners=align_corners) + rotate_fn = Rotate(angle, True, align_corners=align_corners, dtype=np.float64) rotated = rotate_fn(im_type(self.imt[0]), mode=mode, padding_mode=padding_mode) np.testing.assert_allclose(self.imt[0].shape, rotated.shape) diff --git a/tests/test_rotated.py b/tests/test_rotated.py index 1b759cfef5..43b5a68f61 100644 --- a/tests/test_rotated.py +++ b/tests/test_rotated.py @@ -40,7 +40,9 @@ class TestRotated2D(NumpyImageTestCase2D): @parameterized.expand(TEST_CASES_2D) def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, align_corners): - rotate_fn = Rotated(("img", "seg"), angle, keep_size, (mode, "nearest"), padding_mode, align_corners) + rotate_fn = Rotated( + ("img", "seg"), angle, keep_size, (mode, "nearest"), padding_mode, align_corners, dtype=np.float64 + ) rotated = rotate_fn({"img": im_type(self.imt[0]), "seg": im_type(self.segn[0])}) if keep_size: np.testing.assert_allclose(self.imt[0].shape, rotated["img"].shape) @@ -69,7 +71,9 @@ def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, al class TestRotated3D(NumpyImageTestCase3D): @parameterized.expand(TEST_CASES_3D) def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, align_corners): - rotate_fn = Rotated(("img", "seg"), [0, angle, 0], keep_size, (mode, "nearest"), padding_mode, align_corners) + rotate_fn = Rotated( + ("img", "seg"), [0, angle, 0], keep_size, (mode, "nearest"), padding_mode, align_corners, dtype=np.float64 + ) rotated = rotate_fn({"img": im_type(self.imt[0]), "seg": im_type(self.segn[0])}) if keep_size: np.testing.assert_allclose(self.imt[0].shape, rotated["img"].shape) @@ -98,7 +102,9 @@ def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, al class TestRotated3DXY(NumpyImageTestCase3D): @parameterized.expand(TEST_CASES_3D) def test_correct_results(self, im_type, angle, keep_size, mode, padding_mode, align_corners): - rotate_fn = Rotated(("img", "seg"), [0, 0, angle], keep_size, (mode, "nearest"), padding_mode, align_corners) + rotate_fn = Rotated( + ("img", "seg"), [0, 0, angle], keep_size, (mode, "nearest"), padding_mode, align_corners, dtype=np.float64 + ) rotated = rotate_fn({"img": im_type(self.imt[0]), "seg": im_type(self.segn[0])}) if keep_size: np.testing.assert_allclose(self.imt[0].shape, rotated["img"].shape) diff --git a/tests/test_spacing.py b/tests/test_spacing.py index ebff25712d..0dd10a54f3 100644 --- a/tests/test_spacing.py +++ b/tests/test_spacing.py @@ -78,7 +78,7 @@ TESTS.append( [ p, - {"pixdim": (1.0, 1.0)}, + {"pixdim": (1.0, 1.0), "align_corners": True}, np.arange(24).reshape((2, 3, 4)), # data {}, np.array( @@ -192,6 +192,15 @@ np.array([[[[1.0, 1.0, 1.0]], [[1.0, 1.0, 1.0]]]]), ] ) + TESTS.append( # 5D input + [ + p, + {"pixdim": [-1, -1, 0.5], "padding_mode": "zeros", "dtype": float, "align_corners": True}, + np.ones((1, 2, 2, 2, 1)), # data + {"affine": np.eye(4)}, + np.ones((1, 2, 2, 3, 1)), + ] + ) class TestSpacingCase(unittest.TestCase): @@ -203,8 +212,8 @@ def test_spacing(self, in_type, init_param, img, data_param, expected_output): self.assertEqual(_img.device, output_data.device) output_data = output_data.cpu() - np.testing.assert_allclose(output_data, expected_output, atol=1e-3, rtol=1e-3) - sr = len(output_data.shape) - 1 + np.testing.assert_allclose(output_data, expected_output, atol=1e-1, rtol=1e-1) + sr = min(len(output_data.shape) - 1, 3) if isinstance(init_param["pixdim"], float): init_param["pixdim"] = [init_param["pixdim"]] * sr init_pixdim = ensure_tuple(init_param["pixdim"]) diff --git a/tests/test_spatial_resample.py b/tests/test_spatial_resample.py new file mode 100644 index 0000000000..8bd23bdb2c --- /dev/null +++ b/tests/test_spatial_resample.py @@ -0,0 +1,146 @@ +# Copyright (c) MONAI Consortium +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import itertools +import unittest + +import numpy as np +from parameterized import parameterized + +from monai.config import USE_COMPILED +from monai.transforms import SpatialResample +from tests.utils import TEST_NDARRAYS, assert_allclose + +TESTS = [] + +for ind, dst in enumerate( + [ + np.asarray([[1.0, 0.0, 0.0], [0.0, -1.0, 1.0], [0.0, 0.0, 1.0]]), # flip the second + np.asarray([[-1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), # flip the first + ] +): + for p in TEST_NDARRAYS: + for p_src in TEST_NDARRAYS: + for align in (False, True): + for interp_mode in ("nearest", "bilinear"): + TESTS.append( + [ + {}, # default no params + np.arange(4).reshape((1, 2, 2)) + 1.0, # data + { + "src_affine": p_src(np.eye(3)), + "dst_affine": p(dst), + "dtype": np.float32, + "align_corners": align, + "mode": interp_mode, + "padding_mode": "zeros", + }, + np.array([[[2.0, 1.0], [4.0, 3.0]]]) if ind == 0 else np.array([[[3.0, 4.0], [1.0, 2.0]]]), + ] + ) + +for ind, dst in enumerate( + [ + np.asarray([[1.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + np.asarray([[-1.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + ] +): + for p_src in TEST_NDARRAYS: + for align in (True, False): + if align and USE_COMPILED: + interp = ("nearest", "bilinear", 0, 1) # type: ignore + else: + interp = ("nearest", "bilinear") # type: ignore + for interp_mode in interp: # type: ignore + for padding_mode in ("zeros", "border", "reflection"): + TESTS.append( + [ + {}, # default no params + np.arange(12).reshape((1, 2, 2, 3)) + 1.0, # data + { + "src_affine": p_src(np.eye(4)), + "dst_affine": p_src(dst), + "dtype": np.float64, + "align_corners": align, + "mode": interp_mode, + "padding_mode": padding_mode, + }, + np.array([[[[4.0, 5.0, 6.0], [1.0, 2.0, 3.0]], [[10.0, 11.0, 12.0], [7.0, 8.0, 9.0]]]]) + if ind == 0 + else np.array( + [[[[7.0, 8.0, 9.0], [10.0, 11.0, 12.0]], [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]]] + ), + ] + ) + + +class TestSpatialResample(unittest.TestCase): + @parameterized.expand(itertools.product(TEST_NDARRAYS, TESTS)) + def test_flips(self, p_type, args): + init_param, img, data_param, expected_output = args + _img = p_type(img) + _expected_output = p_type(expected_output) + output_data, output_dst = SpatialResample(**init_param)(img=_img, **data_param) + assert_allclose(output_data, _expected_output, rtol=1e-2, atol=1e-2) + expected_dst = ( + data_param.get("dst_affine") if data_param.get("dst_affine") is not None else data_param.get("src_affine") + ) + assert_allclose(output_dst, expected_dst, type_test=False, rtol=1e-2, atol=1e-2) + + @parameterized.expand(itertools.product([True, False], TEST_NDARRAYS)) + def test_4d_5d(self, is_5d, p_type): + new_shape = (1, 2, 2, 3, 1, 1) if is_5d else (1, 2, 2, 3, 1) + img = np.arange(12).reshape(new_shape) + img = np.tile(img, (1, 1, 1, 1, 2, 2) if is_5d else (1, 1, 1, 1, 2)) + _img = p_type(img) + dst = np.asarray([[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, -1.0, 1.5], [0.0, 0.0, 0.0, 1.0]]) + output_data, output_dst = SpatialResample(dtype=np.float32)( + img=_img, src_affine=p_type(np.eye(4)), dst_affine=dst + ) + expected_data = ( + np.asarray( + [ + [ + [[[0.0, 0.0], [0.0, 1.0]], [[0.5, 0.0], [1.5, 1.0]], [[1.0, 2.0], [2.0, 2.0]]], + [[[3.0, 3.0], [3.0, 4.0]], [[3.5, 3.0], [4.5, 4.0]], [[4.0, 5.0], [5.0, 5.0]]], + ], + [ + [[[6.0, 6.0], [6.0, 7.0]], [[6.5, 6.0], [7.5, 7.0]], [[7.0, 8.0], [8.0, 8.0]]], + [[[9.0, 9.0], [9.0, 10.0]], [[9.5, 9.0], [10.5, 10.0]], [[10.0, 11.0], [11.0, 11.0]]], + ], + ], + dtype=np.float32, + ) + if is_5d + else np.asarray( + [ + [[[0.5, 0.0], [0.0, 2.0], [1.5, 1.0]], [[3.5, 3.0], [3.0, 5.0], [4.5, 4.0]]], + [[[6.5, 6.0], [6.0, 8.0], [7.5, 7.0]], [[9.5, 9.0], [9.0, 11.0], [10.5, 10.0]]], + ], + dtype=np.float32, + ) + ) + assert_allclose(output_data, p_type(expected_data[None]), rtol=1e-2, atol=1e-2) + assert_allclose(output_dst, dst, type_test=False, rtol=1e-2, atol=1e-2) + + def test_ill_affine(self): + img = np.arange(12).reshape(1, 2, 2, 3) + ill_affine = np.asarray( + [[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, -1.0, 1.5], [0.0, 0.0, 0.0, 1.0]] + ) + with self.assertRaises(ValueError): + SpatialResample()(img=img, src_affine=np.eye(4), dst_affine=ill_affine) + with self.assertRaises(ValueError): + SpatialResample()(img=img, src_affine=ill_affine, dst_affine=np.eye(3)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_spatial_resampled.py b/tests/test_spatial_resampled.py new file mode 100644 index 0000000000..564635731c --- /dev/null +++ b/tests/test_spatial_resampled.py @@ -0,0 +1,113 @@ +# Copyright (c) MONAI Consortium +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import itertools +import unittest + +import numpy as np +from parameterized import parameterized + +from monai.config import USE_COMPILED +from monai.transforms import SpatialResampleD +from tests.utils import TEST_NDARRAYS, assert_allclose + +TESTS = [] + +for ind, dst in enumerate( + [ + np.asarray([[1.0, 0.0, 0.0], [0.0, -1.0, 1.0], [0.0, 0.0, 1.0]]), # flip the second + np.asarray([[-1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), # flip the first + ] +): + for p in TEST_NDARRAYS: + for p_src in TEST_NDARRAYS: + for align in (False, True): + for interp_mode in ("nearest", "bilinear"): + TESTS.append( + [ + {}, # default no params + np.arange(4).reshape((1, 2, 2)) + 1.0, # data + { + "src": p_src(np.eye(3)), + "dst": p(dst), + "dtype": np.float32, + "align_corners": align, + "mode": interp_mode, + "padding_mode": "zeros", + }, + np.array([[[2.0, 1.0], [4.0, 3.0]]]) if ind == 0 else np.array([[[3.0, 4.0], [1.0, 2.0]]]), + ] + ) + +for ind, dst in enumerate( + [ + np.asarray([[1.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + np.asarray([[-1.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + ] +): + for p_src in TEST_NDARRAYS: + for align in (True, False): + if align and USE_COMPILED: + interp = ("nearest", "bilinear", 0, 1) # type: ignore + else: + interp = ("nearest", "bilinear") # type: ignore + for interp_mode in interp: # type: ignore + for padding_mode in ("zeros", "border", "reflection"): + TESTS.append( + [ + {}, # default no params + np.arange(12).reshape((1, 2, 2, 3)) + 1.0, # data + { + "src": p_src(np.eye(4)), + "dst": p_src(dst), + "dtype": np.float64, + "align_corners": align, + "mode": interp_mode, + "padding_mode": padding_mode, + }, + np.array([[[[4.0, 5.0, 6.0], [1.0, 2.0, 3.0]], [[10.0, 11.0, 12.0], [7.0, 8.0, 9.0]]]]) + if ind == 0 + else np.array( + [[[[7.0, 8.0, 9.0], [10.0, 11.0, 12.0]], [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]]] + ), + ] + ) + + +class TestSpatialResample(unittest.TestCase): + @parameterized.expand(itertools.product(TEST_NDARRAYS, TESTS)) + def test_flips_inverse(self, p_type, args): + _, img, data_param, expected_output = args + _img = p_type(img) + _expected_output = p_type(expected_output) + input_dict = {"img": _img, "img_meta_dict": {"src": data_param.get("src"), "dst": data_param.get("dst")}} + xform = SpatialResampleD( + keys="img", + meta_src_keys="src", + meta_dst_keys="dst", + mode=data_param.get("mode"), + padding_mode=data_param.get("padding_mode"), + align_corners=data_param.get("align_corners"), + ) + output_data = xform(input_dict) + assert_allclose(output_data["img"], _expected_output, rtol=1e-2, atol=1e-2) + assert_allclose( + output_data["img_meta_dict"]["src"], data_param.get("dst"), type_test=False, rtol=1e-2, atol=1e-2 + ) + + inverted = xform.inverse(output_data) + self.assertEqual(inverted["img_transforms"], []) # no further invert after inverting + assert_allclose(inverted["img_meta_dict"]["src"], data_param.get("src"), type_test=False, rtol=1e-2, atol=1e-2) + assert_allclose(inverted["img"], _img, rtol=1e-2, atol=1e-2) + + +if __name__ == "__main__": + unittest.main()