code: def test_joint_pdf_gradients_dense(rng):
      # Compare the analytical and numerical (finite differences) gradient of
      # the joint distribution (i.e. derivatives of each histogram cell) w.r.t.
      # the transform parameters. Since the histograms are discrete partitions
      # of the image intensities, the finite difference approximation is
      # normally not very close to the analytical derivatives. Other sources of
      # error are the interpolation used when transforming the images and the
      # boundary intensities introduced when interpolating outside of the image
      # (i.e. some "zeros" are introduced at the boundary which affect the
      # numerical derivatives but is not taken into account by the analytical
      # derivatives). Thus, we need to relax the verification. Instead of
      # looking for the analytical and numerical gradients to be very close to
      # each other, we will verify that they approximately point in the same
      # direction by testing if the angle they form is close to zero.
      h = 1e-4
  
      # Make sure dictionary entries are processed in the same order regardless
      # of the platform. Otherwise any random numbers drawn within the loop
      # would make the test non-deterministic even if we fix the seed before
      # the loop. Right now, this test does not draw any samples, but we still
      # sort the entries to prevent future related failures.
      for ttype in sorted(factors):
          dim = ttype[1]
          if dim == 2:
              nslices = 1
              transform_method = vf.transform_2d_affine
          else:
              nslices = 45
              transform_method = vf.transform_3d_affine
  
          transform = regtransforms[ttype]
          factor = factors[ttype]
          theta = transform.get_identity_parameters()
  
          (
              static,
              moving,
              static_g2w,
              moving_g2w,
              smask,
              mmask,
              M,
          ) = setup_random_transform(
              transform, factor, nslices, 5.0, rng=rng
          )
          parzen_hist = ParzenJointHistogram(32)
          parzen_hist.setup(static, moving, smask, mmask)
  
          # Compute the gradient at theta with the implementation under test
          M = transform.param_to_matrix(theta)
          shape = np.array(static.shape, dtype=np.int32)
  
          moved = transform_method(moving.astype(np.float32), shape, M)
          moved = np.array(moved)
          parzen_hist.update_pdfs_dense(
              static.astype(np.float64), moved.astype(np.float64)
          )
          # Get the joint distribution evaluated at theta
          J0 = np.copy(parzen_hist.joint)
          grid_to_space = np.eye(dim + 1)
          spacing = np.ones(dim, dtype=np.float64)
          mgrad, inside = vf.gradient(
              moving.astype(np.float32),
              moving_g2w,
              spacing,
              shape,
              grid_to_space,
          )
          params = transform.get_identity_parameters()
          parzen_hist.update_gradient_dense(
              params,
              transform,
              static.astype(np.float64),
              moved.astype(np.float64),
              grid_to_space,
              mgrad,
              smask,
              mmask,
          )
          actual = np.copy(parzen_hist.joint_grad)
          # Now we have the gradient of the joint distribution w.r.t. the
          # transform parameters
  
          # Compute the gradient using finite differences
          n = transform.get_number_of_parameters()
          expected = np.empty_like(actual)
          for i in range(n):
              dtheta = theta.copy()
              dtheta[i] += h
              # Update the joint distribution with the transformed moving image
              M = transform.param_to_matrix(dtheta)
              shape = np.array(static.shape, dtype=np.int32)
              moved = transform_method(
                  moving.astype(np.float32), shape, M
              )
              moved = np.array(moved)
              parzen_hist.update_pdfs_dense(
                  static.astype(np.float64), moved.astype(np.float64)
              )
              J1 = np.copy(parzen_hist.joint)
              expected[..., i] = (J1 - J0) / h
  
          # Dot product and norms of gradients of each joint histogram cell
          # i.e. the derivatives of each cell w.r.t. all parameters
          P = (expected * actual).sum(2)
          enorms = np.sqrt((expected**2).sum(2))
          anorms = np.sqrt((actual**2).sum(2))
          prodnorms = enorms * anorms
          # Cosine of angle between the expected and actual gradients.
          # Exclude very small gradients
          P[prodnorms > 1e-6] /= prodnorms[prodnorms > 1e-6]
          P[prodnorms <= 1e-6] = 0
          # Verify that a large proportion of the gradients point almost in
          # the same direction. Disregard very small gradients
          mean_cosine = P[P != 0].mean()
          std_cosine = P[P != 0].std()
          assert mean_cosine > 0.9
          assert std_cosine < 0.25
  
imports: import numpy
  import scipy
  from functools import reduce
  from operator import mul
  from dipy.core.ndindex import ndindex
  from dipy.core.interpolation import interpolate_scalar_2d, interpolate_scalar_3d
  from dipy.data import get_fnames
  from dipy.align import vector_fields
  from dipy.align.transforms import regtransforms
  from dipy.align.parzenhist import ParzenJointHistogram, cubic_spline, cubic_spline_derivative, sample_domain_regular
  from numpy.testing import assert_array_equal, assert_array_almost_equal, assert_almost_equal, assert_equal, assert_raises
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 7 for calling a global function with example usage: test_joint_pdf_gradients_dense(rng) without return types