code: def test_compute_ssd_demons_step_2d(rng):
      r"""
      Compares the output of the demons step in 2d against an analytical
      step. The fixed image is given by $F(x) = \frac{1}{2}||x - c_f||^2$, the
      moving image is given by $G(x) = \frac{1}{2}||x - c_g||^2$,
      $x, c_f, c_g \in R^{2}$
  
      References
      ----------
      [Vercauteren09] Vercauteren, T., Pennec, X., Perchant, A., & Ayache, N.
                      (2009). Diffeomorphic demons: efficient non-parametric
                      image registration. NeuroImage, 45(1 Suppl), S61-72.
                      doi:10.1016/j.neuroimage.2008.10.040
      """
      # Select arbitrary images' shape (same shape for both images)
      sh = (20, 10)
  
      # Select arbitrary centers
      c_f = np.asarray(sh) / 2
      c_g = c_f + 0.5
  
      # Compute the identity vector field I(x) = x in R^2
      x_0 = np.asarray(range(sh[0]))
      x_1 = np.asarray(range(sh[1]))
      X = np.ndarray(sh + (2,), dtype=np.float64)
      O = np.ones(sh)
      X[..., 0] = x_0[:, None] * O
      X[..., 1] = x_1[None, :] * O
  
      # Compute the gradient fields of F and G
      grad_F = X - c_f
      grad_G = X - c_g
  
      Fnoise = (
          rng.random(np.size(grad_F)).reshape(grad_F.shape)
          * grad_F.max()
          * 0.1
      )
      Fnoise = Fnoise.astype(floating)
      grad_F += Fnoise
  
      Gnoise = (
          rng.random(np.size(grad_G)).reshape(grad_G.shape)
          * grad_G.max()
          * 0.1
      )
      Gnoise = Gnoise.astype(floating)
      grad_G += Gnoise
  
      # The squared norm of grad_G to be used later
      sq_norm_grad_G = np.sum(grad_G**2, -1)
  
      # Compute F and G
      F = 0.5 * np.sum(grad_F**2, -1)
      G = 0.5 * sq_norm_grad_G
  
      Fnoise = rng.random(np.size(F)).reshape(F.shape) * F.max() * 0.1
      Fnoise = Fnoise.astype(floating)
      F += Fnoise
  
      Gnoise = rng.random(np.size(G)).reshape(G.shape) * G.max() * 0.1
      Gnoise = Gnoise.astype(floating)
      G += Gnoise
  
      delta_field = np.array(G - F, dtype=floating)
  
      # Select some pixels to force gradient = 0 and F=G
      random_labels = rng.integers(0, 2, sh[0] * sh[1])
      random_labels = random_labels.reshape(sh)
  
      F[random_labels == 0] = G[random_labels == 0]
      delta_field[random_labels == 0] = 0
      grad_G[random_labels == 0, ...] = 0
      sq_norm_grad_G[random_labels == 0, ...] = 0
  
      # Set arbitrary values for $\sigma_i$ (eq. 4 in [Vercauteren09])
      # The original Demons algorithm used simply |F(x) - G(x)| as an
      # estimator, so let's use it as well
      sigma_i_sq = (F - G) ** 2
  
      # Now select arbitrary parameters for $\sigma_x$ (eq 4 in [Vercauteren09])
      for sigma_x_sq in [0.01, 1.5, 4.2]:
          # Directly compute the demons step according to eq. 4 in
          # [Vercauteren09]
          num = (sigma_x_sq * (F - G))[random_labels == 1]
          den = (sigma_x_sq * sq_norm_grad_G + sigma_i_sq)[
              random_labels == 1
          ]
          # This is $J^{P}$ in eq. 4 [Vercauteren09]
          expected = -1 * np.array(grad_G)
          expected[random_labels == 1, 0] *= num / den
          expected[random_labels == 1, 1] *= num / den
          expected[random_labels == 0, ...] = 0
  
          # Now compute it using the implementation under test
          actual = np.empty_like(expected, dtype=floating)
  
          ssd.compute_ssd_demons_step_2d(
              delta_field,
              np.array(grad_G, dtype=floating),
              sigma_x_sq,
              actual,
          )
  
          assert_array_almost_equal(actual, expected)
  
file_path: dipy/dipy/align/tests/test_sumsqdiff.py
imports: import numpy
  from dipy.align import floating
  from dipy.align import sumsqdiff
  from numpy.testing import assert_equal, assert_almost_equal, assert_array_almost_equal, assert_allclose
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 9 for calling a global function with example usage: test_compute_ssd_demons_step_2d(rng) without return types