code: def test_parzen_densities(rng):
      # Test the computation of the joint intensity distribution
      # using a dense and a sparse set of values
      nbins = 32
      nr = 30
      nc = 35
      ns = 20
      nvals = 50
  
      for dim in [2, 3]:
          if dim == 2:
              shape = (nr, nc)
              static, moving = create_random_image_pair(
                  shape, nvals, rng
              )
          else:
              shape = (ns, nr, nc)
              static, moving = create_random_image_pair(
                  shape, nvals, rng
              )
  
          # Initialize
          parzen_hist = ParzenJointHistogram(nbins)
          parzen_hist.setup(static, moving)
          # Get distributions computed by dense sampling
          parzen_hist.update_pdfs_dense(static, moving)
          actual_joint_dense = parzen_hist.joint
          actual_mmarginal_dense = parzen_hist.mmarginal
          actual_smarginal_dense = parzen_hist.smarginal
  
          # Get distributions computed by sparse sampling
          sval = static.reshape(-1)
          mval = moving.reshape(-1)
          parzen_hist.update_pdfs_sparse(sval, mval)
          actual_joint_sparse = parzen_hist.joint
          actual_mmarginal_sparse = parzen_hist.mmarginal
          actual_smarginal_sparse = parzen_hist.smarginal
  
          # Compute the expected joint distribution with dense sampling
          expected_joint_dense = np.zeros(shape=(nbins, nbins))
          for index in ndindex(shape):
              sv = parzen_hist.bin_normalize_static(static[index])
              mv = parzen_hist.bin_normalize_moving(moving[index])
              sbin = parzen_hist.bin_index(sv)
              # The spline is centered at mv, will evaluate for all row
              spline_arg = np.array([i - mv for i in range(nbins)])
              contribution = cubic_spline(spline_arg)
              expected_joint_dense[sbin, :] += contribution
  
          # Compute the expected joint distribution with sparse sampling
          expected_joint_sparse = np.zeros(shape=(nbins, nbins))
          for index in range(sval.shape[0]):
              sv = parzen_hist.bin_normalize_static(sval[index])
              mv = parzen_hist.bin_normalize_moving(mval[index])
              sbin = parzen_hist.bin_index(sv)
              # The spline is centered at mv, will evaluate for all row
              spline_arg = np.array([i - mv for i in range(nbins)])
              contribution = cubic_spline(spline_arg)
              expected_joint_sparse[sbin, :] += contribution
  
          # Verify joint distributions
          expected_joint_dense /= expected_joint_dense.sum()
          expected_joint_sparse /= expected_joint_sparse.sum()
          assert_array_almost_equal(
              actual_joint_dense, expected_joint_dense
          )
          assert_array_almost_equal(
              actual_joint_sparse, expected_joint_sparse
          )
  
          # Verify moving marginals
          expected_mmarginal_dense = expected_joint_dense.sum(0)
          expected_mmarginal_dense /= expected_mmarginal_dense.sum()
          expected_mmarginal_sparse = expected_joint_sparse.sum(0)
          expected_mmarginal_sparse /= expected_mmarginal_sparse.sum()
          assert_array_almost_equal(
              actual_mmarginal_dense, expected_mmarginal_dense
          )
          assert_array_almost_equal(
              actual_mmarginal_sparse, expected_mmarginal_sparse
          )
  
          # Verify static marginals
          expected_smarginal_dense = expected_joint_dense.sum(1)
          expected_smarginal_dense /= expected_smarginal_dense.sum()
          expected_smarginal_sparse = expected_joint_sparse.sum(1)
          expected_smarginal_sparse /= expected_smarginal_sparse.sum()
          assert_array_almost_equal(
              actual_smarginal_dense, expected_smarginal_dense
          )
          assert_array_almost_equal(
              actual_smarginal_sparse, expected_smarginal_sparse
          )
  
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 5 for calling a global function with example usage: test_parzen_densities(rng) without return types