diff --git a/Modules/Core/Common/wrapping/CMakeLists.txt b/Modules/Core/Common/wrapping/CMakeLists.txt
index 7f764672096..dedc4880f7d 100644
--- a/Modules/Core/Common/wrapping/CMakeLists.txt
+++ b/Modules/Core/Common/wrapping/CMakeLists.txt
@@ -40,6 +40,8 @@ set(WRAPPER_SUBMODULE_ORDER
   itkSymmetricSecondRankTensor
   itkImage
   itkVectorImage
+  itkSpecialCoordinatesImage
+  itkPhasedArray3DSpecialCoordinatesImage
   itkNeighborhood
 #  SwigExtras
   itkImageSourceCommon
diff --git a/Modules/Core/Common/wrapping/itkImageSource.wrap b/Modules/Core/Common/wrapping/itkImageSource.wrap
index 4debe5c1bc0..7dedfc8b85b 100644
--- a/Modules/Core/Common/wrapping/itkImageSource.wrap
+++ b/Modules/Core/Common/wrapping/itkImageSource.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ImageSource" POINTER)
   # Force uchar and ulong image sources for saving in 8 bits and watershed filter
   UNIQUE(image_types "UC;RGBUC;RGBAUC;SI;UI;UL;ULL;${ITKM_IT};${WRAP_ITK_ALL_TYPES}")
@@ -24,4 +25,12 @@ itk_wrap_class("itk::ImageSource" POINTER)
     endforeach()
   endif()
 
+  # PhasedArray3DSpecialCoordinatesImage
+  foreach(real_type ${WRAP_ITK_SCALAR})
+    itk_wrap_template("PA3DSCI${ITKM_${real_type}}" "itk::PhasedArray3DSpecialCoordinatesImage<${ITKT_${real_type}}>")
+  endforeach()
+  foreach(complex_type ${WRAP_ITK_COMPLEX_REAL})
+    itk_wrap_template("PA3DSCI${ITKM_${complex_type}}" "itk::PhasedArray3DSpecialCoordinatesImage<${ITKT_${complex_type}}>")
+  endforeach()
+
 itk_end_wrap_class()
diff --git a/Modules/Core/Common/wrapping/itkImageToImageFilterB.wrap b/Modules/Core/Common/wrapping/itkImageToImageFilterB.wrap
index 04c8e3f3db8..26d743be5d0 100644
--- a/Modules/Core/Common/wrapping/itkImageToImageFilterB.wrap
+++ b/Modules/Core/Common/wrapping/itkImageToImageFilterB.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ImageToImageFilter" POINTER)
 
   # scalar <-> RGB
@@ -71,4 +72,31 @@ itk_wrap_class("itk::ImageToImageFilter" POINTER)
     endforeach()
   endif()
 
+  # PhasedArray3DSpecialCoordinatesImage
+  foreach(t ${WRAP_ITK_SCALAR})
+    foreach(ut ${WRAP_ITK_INT})
+      itk_wrap_template("I${ITKM_${t}}3PA3DSCI${ITKM_${ut}}"
+        "itk::Image< ${ITKT_${t}}, 3 >, itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${ut}} >")
+    endforeach()
+    foreach(t2 ${WRAP_ITK_SCALAR})
+      itk_wrap_template("PA3DSCI${ITKM_${t}}PA3DSCI${ITKM_${t2}}"
+        "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t2}} >")
+      itk_wrap_template("PA3DSCI${ITKM_${t}}I${ITKM_${t2}}3"
+        "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, itk::Image< ${ITKT_${t2}}, 3 >")
+    endforeach()
+  endforeach()
+
+  foreach(rt ${WRAP_ITK_REAL})
+    foreach(ct ${WRAP_ITK_COMPLEX_REAL})
+      itk_wrap_template("PA3DSCI${ITKM_${rt}}PA3DSCI${ITKM_${ct}}"
+      "itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${rt}} >, itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${ct}} >")
+      itk_wrap_template("PA3DSCI${ITKM_${ct}}PA3DSCI${ITKM_${rt}}"
+      "itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${ct}} >, itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${rt}} >")
+    endforeach()
+  endforeach()
+
+  foreach(t3 ${WRAP_ITK_COMPLEX_REAL})
+    itk_wrap_template("PA3DSCI${ITKM_${t3}}PA3DSCI${ITKM_${t3}}"
+      "itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${t3}} >, itk::PhasedArray3DSpecialCoordinatesImage<  ${ITKT_${t3}} >")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/Common/wrapping/itkInPlaceImageFilterB.wrap b/Modules/Core/Common/wrapping/itkInPlaceImageFilterB.wrap
index 3e2773c3047..a57a987e07e 100644
--- a/Modules/Core/Common/wrapping/itkInPlaceImageFilterB.wrap
+++ b/Modules/Core/Common/wrapping/itkInPlaceImageFilterB.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::InPlaceImageFilter" POINTER)
   # VectorImage <-> scalar
   UNIQUE(to_types "UC;${WRAP_ITK_SCALAR}")
@@ -38,4 +39,12 @@ itk_wrap_class("itk::InPlaceImageFilter" POINTER)
     endforeach()
   endforeach()
 
+  # PhasedArray3DSpecialCoordinatesImage
+  foreach(t ${WRAP_ITK_SCALAR})
+    foreach(ut ${WRAP_ITK_INT})
+      itk_wrap_template("I${ITKM_${t}}3PA3DSCI${ITKM_${ut}}"
+        "itk::Image< ${ITKT_${t}}, 3 >, itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${ut}} >")
+    endforeach()
+  endforeach()
+
 itk_end_wrap_class()
diff --git a/Modules/Core/Common/wrapping/itkPhasedArray3DSpecialCoordinatesImage.wrap b/Modules/Core/Common/wrapping/itkPhasedArray3DSpecialCoordinatesImage.wrap
new file mode 100644
index 00000000000..981ff044e24
--- /dev/null
+++ b/Modules/Core/Common/wrapping/itkPhasedArray3DSpecialCoordinatesImage.wrap
@@ -0,0 +1,49 @@
+itk_wrap_include("complex")
+
+# Explicitly override template method wrappings so that implicit
+# scalar type is always `double` for greatest precision.
+# Adds wrapping overrides to `itkPhasedArray3DSpecialCoordinatesImage_ext.i` configured with
+# CMake for input to SWIG wrapping generation.
+# See `DECL_PYTHON_IMAGEBASE_CLASS` definition in `ITK/Wrapping/Generators/Python/PyBase/pyBase.i`
+# for precedent.
+string(APPEND ITK_WRAP_PYTHON_SWIG_EXT "
+%inline %{
+#include \"itkContinuousIndexSwigInterface.h\"
+%}
+
+%define DECL_PYTHON_PHASEDARRAY3DSPECIALCOORDINATESIMAGE_CLASS(swig_name, template_params)
+
+  %extend swig_name {
+      itkIndex##template_params TransformPhysicalPointToIndex(const itkPointD##template_params & point ) {
+          return self->TransformPhysicalPointToIndex<double>( point );
+        }
+
+      itkContinuousIndexD##template_params TransformPhysicalPointToContinuousIndex(const itkPointD##template_params & point ) {
+          return self->TransformPhysicalPointToContinuousIndex<double>( point );
+      }
+
+      itkPointD##template_params TransformContinuousIndexToPhysicalPoint(const itkContinuousIndexD##template_params & idx ) {
+          return self->TransformContinuousIndexToPhysicalPoint<double>( idx );
+      }
+
+      itkPointD##template_params TransformIndexToPhysicalPoint(const itkIndex##template_params & idx ) {
+          return self->TransformIndexToPhysicalPoint<double>( idx );
+      }
+  }
+
+%enddef
+")
+
+# Wrap class for real and complex pixel types. Dimension is always 3.
+itk_wrap_class("itk::PhasedArray3DSpecialCoordinatesImage" POINTER_WITH_SUPERCLASS)
+  foreach(t ${WRAP_ITK_SCALAR})
+    itk_wrap_template("${ITKM_${t}}" "${ITKT_${t}}")
+    string(APPEND ITK_WRAP_PYTHON_SWIG_EXT "DECL_PYTHON_PHASEDARRAY3DSPECIALCOORDINATESIMAGE_CLASS(${WRAPPER_SWIG_NAME}${ITKM_${t}}, 3)\n")
+  endforeach()
+  foreach(t3 ${WRAP_ITK_COMPLEX_REAL})
+    itk_wrap_template("${ITKM_${t3}}" "${ITKT_${t3}}")
+    string(APPEND ITK_WRAP_PYTHON_SWIG_EXT "DECL_PYTHON_PHASEDARRAY3DSPECIALCOORDINATESIMAGE_CLASS(${WRAPPER_SWIG_NAME}${ITKM_${t3}}, 3)\n")
+  endforeach()
+itk_end_wrap_class()
+
+# Then wrap consuming filters in itkPhasedArray3DSpecialCoordinatesImageFilters.wrap
diff --git a/Modules/Core/Common/wrapping/itkSpecialCoordinatesImage.wrap b/Modules/Core/Common/wrapping/itkSpecialCoordinatesImage.wrap
new file mode 100644
index 00000000000..962e4697d1c
--- /dev/null
+++ b/Modules/Core/Common/wrapping/itkSpecialCoordinatesImage.wrap
@@ -0,0 +1,13 @@
+itk_wrap_include("complex")
+
+# Wrap underlying class hierarchy
+itk_wrap_class("itk::SpecialCoordinatesImage" POINTER)
+  foreach(d ${ITK_WRAP_IMAGE_DIMS})
+    foreach(t ${WRAP_ITK_SCALAR})
+      itk_wrap_template("${ITKM_${t}}${d}" "${ITKT_${t}}, ${d}")
+    endforeach()
+    foreach(t3 ${WRAP_ITK_COMPLEX_REAL})
+      itk_wrap_template("${ITKM_${t3}}${d}" "${ITKT_${t3}}, ${d}")
+    endforeach()
+  endforeach()
+itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkGaussianInterpolateImageFunction.wrap b/Modules/Core/ImageFunction/wrapping/itkGaussianInterpolateImageFunction.wrap
index 0128182f826..14a2542e29f 100644
--- a/Modules/Core/ImageFunction/wrapping/itkGaussianInterpolateImageFunction.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkGaussianInterpolateImageFunction.wrap
@@ -1,7 +1,17 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::GaussianInterpolateImageFunction" POINTER)
   foreach(d ${ITK_WRAP_IMAGE_DIMS})
     foreach(t ${WRAP_ITK_SCALAR})
       itk_wrap_template("${ITKM_I${t}${d}}${ITKM_D}" "${ITKT_I${t}${d}},${ITKT_D}")
     endforeach()
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}${ITKM_D}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${ITKT_D}")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkImageFunctionBase.wrap b/Modules/Core/ImageFunction/wrapping/itkImageFunctionBase.wrap
index ad7f06a9245..2b3188d42a0 100644
--- a/Modules/Core/ImageFunction/wrapping/itkImageFunctionBase.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkImageFunctionBase.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ImageFunction" POINTER)
   foreach(d ${ITK_WRAP_IMAGE_DIMS})
     # UC is required for InterpolateImageFunction
@@ -31,4 +32,13 @@ itk_wrap_class("itk::ImageFunction" POINTER)
     endforeach()
 
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}${ITKM_D}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${ITKT_D}, ${ITKT_D}")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkInterpolateImageFunction.wrap b/Modules/Core/ImageFunction/wrapping/itkInterpolateImageFunction.wrap
index 82a58e0c5f9..c2e785c27b1 100644
--- a/Modules/Core/ImageFunction/wrapping/itkInterpolateImageFunction.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkInterpolateImageFunction.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::InterpolateImageFunction" POINTER)
   UNIQUE(types "${WRAP_ITK_SCALAR};UC")
   foreach(d ${ITK_WRAP_IMAGE_DIMS})
@@ -15,4 +16,12 @@ itk_wrap_class("itk::InterpolateImageFunction" POINTER)
       itk_wrap_template("${ITKM_I${v}${d}}${ITKM_D}" "${ITKT_I${v}${d}},${ITKT_D}")
     endforeach()
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}${ITKM_D}" "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${ITKT_D}")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkLinearInterpolateImageFunction.wrap b/Modules/Core/ImageFunction/wrapping/itkLinearInterpolateImageFunction.wrap
index 10460f6dc6f..8479eaf0cd9 100644
--- a/Modules/Core/ImageFunction/wrapping/itkLinearInterpolateImageFunction.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkLinearInterpolateImageFunction.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::LinearInterpolateImageFunction" POINTER)
   foreach(d ${ITK_WRAP_IMAGE_DIMS})
     foreach(t ${WRAP_ITK_SCALAR})
@@ -13,4 +14,13 @@ itk_wrap_class("itk::LinearInterpolateImageFunction" POINTER)
       itk_wrap_template("${ITKM_I${v}${d}}${ITKM_D}" "${ITKT_I${v}${d}},${ITKT_D}")
     endforeach()
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}${ITKM_D}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${ITKT_D}")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkNearestNeighborInterpolateImageFunction.wrap b/Modules/Core/ImageFunction/wrapping/itkNearestNeighborInterpolateImageFunction.wrap
index c19d0192622..7fc7d4141bf 100644
--- a/Modules/Core/ImageFunction/wrapping/itkNearestNeighborInterpolateImageFunction.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkNearestNeighborInterpolateImageFunction.wrap
@@ -1,7 +1,17 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::NearestNeighborInterpolateImageFunction" POINTER)
   foreach(d ${ITK_WRAP_IMAGE_DIMS})
     foreach(t ${WRAP_ITK_SCALAR})
       itk_wrap_template("${ITKM_I${t}${d}}${ITKM_D}" "${ITKT_I${t}${d}},${ITKT_D}")
     endforeach()
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}${ITKM_D}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${ITKT_D}")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Core/ImageFunction/wrapping/itkWindowedSincInterpolateImageFunction.wrap b/Modules/Core/ImageFunction/wrapping/itkWindowedSincInterpolateImageFunction.wrap
index 14f06b15a94..e3c79750efb 100644
--- a/Modules/Core/ImageFunction/wrapping/itkWindowedSincInterpolateImageFunction.wrap
+++ b/Modules/Core/ImageFunction/wrapping/itkWindowedSincInterpolateImageFunction.wrap
@@ -1,5 +1,6 @@
 set(WRAPPER_AUTO_INCLUDE_HEADERS OFF)
 itk_wrap_include("itkWindowedSincInterpolateImageFunction.h")
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 
 set(window_functions "Hamming" "Cosine" "Welch" "Lanczos")
 set(radii 2 3)
@@ -21,4 +22,16 @@ itk_wrap_class("itk::WindowedSincInterpolateImageFunction" POINTER)
       endforeach()
     endforeach()
   endforeach()
+
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    foreach(r ${radii}) # radius
+      foreach(function ${window_functions})
+        itk_wrap_template("PA3DSCI${ITKM_${t}}${r}${function}" "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, ${r}, itk::Function::${function}WindowFunction< ${r} >")
+      endforeach()
+    endforeach()
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Filtering/ImageFilterBase/wrapping/itkCastImageFilter.wrap b/Modules/Filtering/ImageFilterBase/wrapping/itkCastImageFilter.wrap
index d3a705d408f..fbfd3ac07f4 100644
--- a/Modules/Filtering/ImageFilterBase/wrapping/itkCastImageFilter.wrap
+++ b/Modules/Filtering/ImageFilterBase/wrapping/itkCastImageFilter.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::CastImageFilter" POINTER_WITH_SUPERCLASS)
   # Create cast filters between all scalar types. Also force that cast-to-uchar
   # filters are created for all scalar types, and unsigned long for
@@ -80,4 +81,14 @@ itk_wrap_class("itk::CastImageFilter" POINTER_WITH_SUPERCLASS)
   if(NOT ITK_WRAP_unsigned_long_long AND NOT ${ITKM_IT} STREQUAL "ULL")
     itk_wrap_image_filter_combinations("ULL" "${WRAP_ITK_INT}")
   endif()
+
+  # PhasedArray3DSpecialCoordinatesImage
+  foreach(t1 ${WRAP_ITK_SCALAR})
+    foreach(t2 ${WRAP_ITK_SCALAR})
+      itk_wrap_template("PA3DSCI${ITKM_${t1}}PA3DSCI${ITKM_${t2}}"
+        "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t1}} >, itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t2}} >")
+      itk_wrap_template("PA3DSCI${ITKM_${t1}}I${ITKM_${t2}}3"
+        "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t1}} >, itk::Image< ${ITKT_${t2}}, 3 >")
+    endforeach()
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/Filtering/ImageGrid/test/Input/ScanConvertPhasedArray3DTestInput.mha.sha512 b/Modules/Filtering/ImageGrid/test/Input/ScanConvertPhasedArray3DTestInput.mha.sha512
new file mode 100644
index 00000000000..a95f1d6697b
--- /dev/null
+++ b/Modules/Filtering/ImageGrid/test/Input/ScanConvertPhasedArray3DTestInput.mha.sha512
@@ -0,0 +1 @@
+7d5e6c2b070107a279277c48ead400bab15edd7aaa1b1b2238a8437da1435c41f75c3b2323fdba17f6415d10ea37992ff053f52a5f6d7a0937f78cfaf1f3687a
diff --git a/Modules/Filtering/ImageGrid/wrapping/itkResampleImageFilter.wrap b/Modules/Filtering/ImageGrid/wrapping/itkResampleImageFilter.wrap
index abb6a7faa19..5cd9102a8f3 100644
--- a/Modules/Filtering/ImageGrid/wrapping/itkResampleImageFilter.wrap
+++ b/Modules/Filtering/ImageGrid/wrapping/itkResampleImageFilter.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ResampleImageFilter" POINTER)
   itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2)
   itk_wrap_image_filter("${WRAP_ITK_VECTOR_REAL}" 2)
@@ -11,4 +12,13 @@ itk_wrap_class("itk::ResampleImageFilter" POINTER)
     endforeach()
   endforeach()
 
+  set(PA3DSCI_types "F")
+  if(ITK_WRAP_unsigned_char)
+    list(APPEND PA3DSCI_types "UC")
+  endif()
+  foreach(t ${PA3DSCI_types})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}I${ITKM_${t}}3"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >, itk::Image< ${ITKT_${t}}, 3 >")
+  endforeach()
+
 itk_end_wrap_class()
diff --git a/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha.sha512 b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha.sha512
new file mode 100644
index 00000000000..0a2d3b4a03f
--- /dev/null
+++ b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha.sha512
@@ -0,0 +1 @@
+8fc8129454b6ec93fb5c133d9835f8a49df0ab14584d44b3454e1b36dee60ef92e29bb9e8c3b9a1340129877fc0ec81d7beb5022eaf3bcbbed8fa3ebab8e6c74
diff --git a/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha.sha512 b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha.sha512
new file mode 100644
index 00000000000..d5d6f9622c3
--- /dev/null
+++ b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha.sha512
@@ -0,0 +1 @@
+cf4bac01a2692100f1c2b2445f3468b26c754e380f4d0119b466c079639f309ad18f775106f341a72c6126a7a636d55c494ed40edf441de262418ff65b2a792e
diff --git a/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DTestOutput.mha.sha512 b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DTestOutput.mha.sha512
new file mode 100644
index 00000000000..d545bf6f3dc
--- /dev/null
+++ b/Modules/Filtering/ImageGrid/wrapping/test/Baseline/ScanConvertPhasedArray3DTestOutput.mha.sha512
@@ -0,0 +1 @@
+6e589887b660f79513d1c09479cc70bd815e532e013fcb3bb185a21c641e417210d6a8286199a265b8e6b9159254e6237108539d0c5eb865a7bedb03d64cbf50
diff --git a/Modules/Filtering/ImageGrid/wrapping/test/CMakeLists.txt b/Modules/Filtering/ImageGrid/wrapping/test/CMakeLists.txt
index 74ca426e136..d8f0f2882b1 100644
--- a/Modules/Filtering/ImageGrid/wrapping/test/CMakeLists.txt
+++ b/Modules/Filtering/ImageGrid/wrapping/test/CMakeLists.txt
@@ -46,7 +46,7 @@ endif()
 
 # let's make sure 3D uchar images are wrapped
 list(FIND ITK_WRAP_IMAGE_DIMS 3 wrap_3_index)
-if(ITK_WRAP_PYTHON AND ITK_WRAP_unsigned_char AND wrap_2_index GREATER -1)
+if(ITK_WRAP_PYTHON AND ITK_WRAP_unsigned_char AND wrap_3_index GREATER -1)
   itk_python_add_test(NAME OrientImageFilterTest
     TEST_DRIVER_ARGS
       --compare DATA{Baseline/OrientImageFilterTest.mha}
@@ -55,4 +55,56 @@ if(ITK_WRAP_PYTHON AND ITK_WRAP_unsigned_char AND wrap_2_index GREATER -1)
       DATA{${ITK_DATA_ROOT}/Baseline/Filtering/SubsampleVolume.mha}
       ${ITK_TEST_OUTPUT_DIR}/OrientImageFilterTest.mha
     )
+
+  itk_python_expression_add_test(NAME PythonInstantiateGaussianInterpolateImageFunctionPA3DSCI
+    EXPRESSION "iFunc = itk.GaussianInterpolateImageFunction[itk.PhasedArray3DSpecialCoordinatesImage[itk.F], itk.D].New()"
+    )
+
+  itk_python_add_test(NAME PythonPhasedArray3DResampleNearestTest
+    TEST_DRIVER_ARGS
+    --compareIntensityTolerance 1
+    --compare
+      DATA{Baseline/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha}
+      ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha
+    COMMAND PhasedArray3DResampleTest.py
+      -i DATA{${test_input_dir}/ScanConvertPhasedArray3DTestInput.mha}
+      -o ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DITKNearestNeighborTestOutput.mha
+      --azimuth-angular-separation 0.0872665
+      --elevation-angular-separation 0.0174533
+      --radius-sample-size 0.2
+      --first-sample-distance 8.0
+      --interpolation-type nearest
+    )
+
+  itk_python_add_test(NAME PythonPhasedArray3DResampleLinearTest
+    TEST_DRIVER_ARGS
+    --compareIntensityTolerance 1
+    --compare
+      DATA{Baseline/ScanConvertPhasedArray3DTestOutput.mha}
+      ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DTestOutput.mha
+    COMMAND PhasedArray3DResampleTest.py
+      -i DATA{${test_input_dir}/ScanConvertPhasedArray3DTestInput.mha}
+      -o ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DTestOutput.mha
+      --azimuth-angular-separation 0.0872665
+      --elevation-angular-separation 0.0174533
+      --radius-sample-size 0.2
+      --first-sample-distance 8.0
+      # --interpolation-type linear  # linear is the default
+    )
+
+  itk_python_add_test(NAME PythonPhasedArray3DResampleSincTest
+    TEST_DRIVER_ARGS
+    --compareIntensityTolerance 1
+    --compare
+      DATA{Baseline/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha}
+      ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha
+    COMMAND PhasedArray3DResampleTest.py
+      -i DATA{${test_input_dir}/ScanConvertPhasedArray3DTestInput.mha}
+      -o ${ITK_TEST_OUTPUT_DIR}/ScanConvertPhasedArray3DITKWindowedSincTestOutput.mha
+      --azimuth-angular-separation 0.0872665
+      --elevation-angular-separation 0.0174533
+      --radius-sample-size 0.2
+      --first-sample-distance 8.0
+      --interpolation-type sinc
+    )
 endif()
diff --git a/Modules/Filtering/ImageGrid/wrapping/test/PhasedArray3DResampleTest.py b/Modules/Filtering/ImageGrid/wrapping/test/PhasedArray3DResampleTest.py
new file mode 100644
index 00000000000..70c559851e4
--- /dev/null
+++ b/Modules/Filtering/ImageGrid/wrapping/test/PhasedArray3DResampleTest.py
@@ -0,0 +1,75 @@
+#==========================================================================
+#
+#   Copyright NumFOCUS
+#
+#   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
+#
+#          https://www.apache.org/licenses/LICENSE-2.0.txt
+#
+#   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 math
+import argparse
+
+import itk
+
+parser = argparse.ArgumentParser(description="Estimate three back-scatter coefficients.")
+parser.add_argument("-i", "--input-image", required=True)
+parser.add_argument("-o", "--output-image", required=True)
+parser.add_argument("-a", "--azimuth-angular-separation", type=float, required=True)
+parser.add_argument("-e", "--elevation-angular-separation", type=float, required=True)
+parser.add_argument("-r", "--radius-sample-size", type=float, required=True)
+parser.add_argument("-f", "--first-sample-distance", type=float, required=True)
+parser.add_argument("-t", "--interpolation-type", choices=["nearest", "linear", "sinc"], default="linear")
+args = parser.parse_args()
+
+itk.auto_progress(2)
+
+pixel_type = itk.UC
+dimension = 3
+image_type = itk.PhasedArray3DSpecialCoordinatesImage[pixel_type]
+reader = itk.ImageFileReader[image_type].New()
+reader.SetFileName(args.input_image)
+reader.Update()
+image = reader.GetOutput()
+image.DisconnectPipeline()
+
+image.SetAzimuthAngularSeparation(args.azimuth_angular_separation)
+image.SetElevationAngularSeparation(args.elevation_angular_separation)
+image.SetFirstSampleDistance(args.first_sample_distance)
+image.SetRadiusSampleSize(args.radius_sample_size)
+print(image)
+
+print("Verify resampling works with PhasedArray3DSpecialCoordinatesImage input")
+output_size = [128] * dimension
+output_spacing = [0.2] * dimension
+origin0 = -1 * output_spacing[0] * output_size[0] / 2.0
+origin1 = -1 * output_spacing[1] * output_size[1] / 2.0
+origin2 = 0.0;
+output_origin = [origin0, origin1, origin2]
+
+if args.interpolation_type == "nearest":
+  interpolator = itk.NearestNeighborInterpolateImageFunction.New(image)
+elif args.interpolation_type == "linear":
+  interpolator = itk.LinearInterpolateImageFunction.New(image)
+elif args.interpolation_type == "sinc":
+  window_type = itk.LanczosWindowFunction[dimension]
+  interpolator = itk.WindowedSincInterpolateImageFunction[type(image), dimension, window_type].New()
+
+result = itk.resample_image_filter(
+  image,
+  size=output_size,
+  output_spacing=output_spacing,
+  output_origin=output_origin,
+  interpolator=interpolator,
+  )
+itk.imwrite(result, args.output_image, compression=True)
+print(f"Result image written to: {args.output_image}")
diff --git a/Modules/Filtering/ImageIntensity/wrapping/itkRescaleIntensityImageFilter.wrap b/Modules/Filtering/ImageIntensity/wrapping/itkRescaleIntensityImageFilter.wrap
index ce6ef2ad350..d759ac3a974 100644
--- a/Modules/Filtering/ImageIntensity/wrapping/itkRescaleIntensityImageFilter.wrap
+++ b/Modules/Filtering/ImageIntensity/wrapping/itkRescaleIntensityImageFilter.wrap
@@ -1,6 +1,15 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::RescaleIntensityImageFilter" POINTER_WITH_SUPERCLASS)
   # Create rescale filters from every scalar type to every scalar type. Also force
   # filters from every scalar to uchar, to allow for saving as 8-bit images.
   UNIQUE(to_types "UC;${WRAP_ITK_SCALAR}")
   itk_wrap_image_filter_combinations("${WRAP_ITK_SCALAR}" "${to_types}")
+
+  # PhasedArray3DSpecialCoordinatesImage
+  foreach(t ${WRAP_ITK_SCALAR})
+    foreach(ut ${WRAP_ITK_INT})
+      itk_wrap_template("I${ITKM_${t}}3PA3DSCI${ITKM_${ut}}"
+        "itk::Image< ${ITKT_${t}}, 3 >, itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${ut}} >")
+    endforeach()
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/IO/ImageBase/wrapping/itkImageFileReader.wrap b/Modules/IO/ImageBase/wrapping/itkImageFileReader.wrap
index c8e2fab32dd..6c581c92b0b 100644
--- a/Modules/IO/ImageBase/wrapping/itkImageFileReader.wrap
+++ b/Modules/IO/ImageBase/wrapping/itkImageFileReader.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ImageFileReader" POINTER)
   # Force uchar image IO
   UNIQUE(image_types "ULL;SI;UI;UC;${WRAP_ITK_ALL_TYPES}")
@@ -9,4 +10,9 @@ itk_wrap_class("itk::ImageFileReader" POINTER)
       itk_wrap_template("${ITKM_VI${type}${d}}"  "${ITKT_VI${type}${d}}")
     endforeach()
   endforeach()
+
+  foreach(t ${WRAP_ITK_SCALAR})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Modules/IO/ImageBase/wrapping/itkImageFileWriter.wrap b/Modules/IO/ImageBase/wrapping/itkImageFileWriter.wrap
index 9511e5354d2..6ee18a6c248 100644
--- a/Modules/IO/ImageBase/wrapping/itkImageFileWriter.wrap
+++ b/Modules/IO/ImageBase/wrapping/itkImageFileWriter.wrap
@@ -1,3 +1,4 @@
+itk_wrap_include("itkPhasedArray3DSpecialCoordinatesImage.h")
 itk_wrap_class("itk::ImageFileWriter" POINTER)
   # Force uchar image IO
   UNIQUE(image_types "ULL;UI;SI;UC;RGBUC;${WRAP_ITK_ALL_TYPES}")
@@ -9,4 +10,9 @@ itk_wrap_class("itk::ImageFileWriter" POINTER)
       itk_wrap_template("${ITKM_VI${type}${d}}"  "${ITKT_VI${type}${d}}")
     endforeach()
   endforeach()
+
+  foreach(t ${WRAP_ITK_SCALAR})
+    itk_wrap_template("PA3DSCI${ITKM_${t}}"
+      "itk::PhasedArray3DSpecialCoordinatesImage< ${ITKT_${t}} >")
+  endforeach()
 itk_end_wrap_class()
diff --git a/Wrapping/Generators/Python/PyBase/pyBase.i b/Wrapping/Generators/Python/PyBase/pyBase.i
index ebe4bdafc60..431f4632d43 100644
--- a/Wrapping/Generators/Python/PyBase/pyBase.i
+++ b/Wrapping/Generators/Python/PyBase/pyBase.i
@@ -484,15 +484,11 @@ str = str
         }
 
         itkPointD##template_params TransformContinuousIndexToPhysicalPoint( itkContinuousIndexD##template_params & idx ) {
-            itkPointD##template_params point;
-            self->TransformContinuousIndexToPhysicalPoint<double>( idx, point );
-            return point;
+            return self->TransformContinuousIndexToPhysicalPoint<double>( idx );
         }
 
         itkPointD##template_params TransformIndexToPhysicalPoint( itkIndex##template_params & idx ) {
-            itkPointD##template_params point;
-            self->TransformIndexToPhysicalPoint<double>( idx, point );
-            return point;
+            return self->TransformIndexToPhysicalPoint<double>( idx );
         }
 
         %pythoncode %{