code: def test_grayscale_iter():
      nclasses = 4
      beta = np.float64(0.1)
      max_iter = 15
      background_noise = True
  
      image = create_image()
      com = ConstantObservationModel()
      icm = IteratedConditionalModes()
  
      mu, sigmasq = com.initialize_param_uniform(image, nclasses)
      neglogl = com.negloglikelihood(image, mu, sigmasq, nclasses)
      initial_segmentation = icm.initialize_maximum_likelihood(neglogl)
      npt.assert_(initial_segmentation.max() == nclasses - 1)
      npt.assert_(initial_segmentation.min() == 0)
  
      mu, sigmasq = com.seg_stats(image, initial_segmentation, nclasses)
      npt.assert_(mu[0] >= 0.0)
      npt.assert_(mu[1] >= 0.0)
      npt.assert_(mu[2] >= 0.0)
      npt.assert_(mu[3] >= 0.0)
      npt.assert_(sigmasq[0] >= 0.0)
      npt.assert_(sigmasq[1] >= 0.0)
      npt.assert_(sigmasq[2] >= 0.0)
      npt.assert_(sigmasq[3] >= 0.0)
  
      if background_noise:
          zero = np.zeros_like(image) + 0.001
          zero_noise = add_noise(zero, 10000, 1, noise_type="gaussian")
          image_gauss = np.where(image == 0, zero_noise, image)
      else:
          image_gauss = image
  
      final_segmentation = np.empty_like(image)
      seg_init = initial_segmentation
      energies = []
  
      for i in range(max_iter):
          PLN = icm.prob_neighborhood(
              initial_segmentation, beta, nclasses
          )
          npt.assert_(np.all((PLN >= 0) & (PLN <= 1.0)))
  
          if beta == 0.0:
              npt.assert_almost_equal(PLN[50, 50, 1, 0], 0.25, True)
              npt.assert_almost_equal(PLN[50, 50, 1, 1], 0.25, True)
              npt.assert_almost_equal(PLN[50, 50, 1, 2], 0.25, True)
              npt.assert_almost_equal(PLN[50, 50, 1, 3], 0.25, True)
              npt.assert_almost_equal(PLN[147, 129, 1, 0], 0.25, True)
              npt.assert_almost_equal(PLN[147, 129, 1, 1], 0.25, True)
              npt.assert_almost_equal(PLN[147, 129, 1, 2], 0.25, True)
              npt.assert_almost_equal(PLN[147, 129, 1, 3], 0.25, True)
              npt.assert_almost_equal(PLN[61, 152, 1, 0], 0.25, True)
              npt.assert_almost_equal(PLN[61, 152, 1, 1], 0.25, True)
              npt.assert_almost_equal(PLN[61, 152, 1, 2], 0.25, True)
              npt.assert_almost_equal(PLN[61, 152, 1, 3], 0.25, True)
              npt.assert_almost_equal(PLN[100, 100, 1, 0], 0.25, True)
              npt.assert_almost_equal(PLN[100, 100, 1, 1], 0.25, True)
              npt.assert_almost_equal(PLN[100, 100, 1, 2], 0.25, True)
              npt.assert_almost_equal(PLN[100, 100, 1, 3], 0.25, True)
  
          PLY = com.prob_image(image_gauss, nclasses, mu, sigmasq, PLN)
          npt.assert_(np.all((PLY >= 0) & (PLY <= 1.0)))
          npt.assert_(PLY[50, 50, 1, 0] > PLY[50, 50, 1, 1])
          npt.assert_(PLY[50, 50, 1, 0] > PLY[50, 50, 1, 2])
          npt.assert_(PLY[50, 50, 1, 0] > PLY[50, 50, 1, 3])
          npt.assert_(PLY[100, 100, 1, 3] > PLY[100, 100, 1, 0])
          npt.assert_(PLY[100, 100, 1, 3] > PLY[100, 100, 1, 1])
          npt.assert_(PLY[100, 100, 1, 3] > PLY[100, 100, 1, 2])
  
          mu_upd, sigmasq_upd = com.update_param(
              image_gauss, PLY, mu, nclasses
          )
          npt.assert_(mu_upd[0] >= 0.0)
          npt.assert_(mu_upd[1] >= 0.0)
          npt.assert_(mu_upd[2] >= 0.0)
          npt.assert_(mu_upd[3] >= 0.0)
          npt.assert_(sigmasq_upd[0] >= 0.0)
          npt.assert_(sigmasq_upd[1] >= 0.0)
          npt.assert_(sigmasq_upd[2] >= 0.0)
          npt.assert_(sigmasq_upd[3] >= 0.0)
  
          negll = com.negloglikelihood(
              image_gauss, mu_upd, sigmasq_upd, nclasses
          )
          npt.assert_(negll[50, 50, 1, 0] < negll[50, 50, 1, 1])
          npt.assert_(negll[50, 50, 1, 0] < negll[50, 50, 1, 2])
          npt.assert_(negll[50, 50, 1, 0] < negll[50, 50, 1, 3])
          npt.assert_(negll[100, 100, 1, 3] < negll[100, 100, 1, 0])
          npt.assert_(negll[100, 100, 1, 3] < negll[100, 100, 1, 1])
          npt.assert_(negll[100, 100, 1, 3] < negll[100, 100, 1, 2])
  
          final_segmentation, energy = icm.icm_ising(
              negll, beta, initial_segmentation
          )
          print(energy[energy > -np.inf].sum())
          energies.append(energy[energy > -np.inf].sum())
  
          initial_segmentation = final_segmentation
          mu = mu_upd
          sigmasq = sigmasq_upd
  
      npt.assert_(energies[-1] < energies[0])
  
      difference_map = np.abs(seg_init - final_segmentation)
      npt.assert_(np.abs(np.sum(difference_map)) != 0)
  
file_path: dipy/dipy/segment/tests/test_mrf.py
imports: import numpy
  import numpy.testing
  from dipy.data import get_fnames
  from dipy.sims.voxel import add_noise
  from dipy.segment.mrf import ConstantObservationModel, IteratedConditionalModes
  from dipy.segment.tissue import TissueClassifierHMRF
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 6 for calling a global function with example usage: test_grayscale_iter() without return types