docs: def iyang(imgIn, krnl, imgSeg, Cnt, itr=5):
      '''partial volume correction using iterative Yang method
      imgIn: input image which is blurred due to the PSF of the scanner
      krnl: shift invariant kernel of the PSF
      imgSeg: segmentation into regions starting with 0 (e.g., background) and then next integer numbers
      itr: number of iteration (default 5)
      '''
      dim = imgIn.shape
      m = np.int32(np.max(imgSeg))
      m_a = np.zeros(( m+1, itr ), dtype=np.float32)
  
      for jr in range(0,m+1): 
          m_a[jr, 0] = np.mean( imgIn[imgSeg==jr] )
  
      # init output image
      imgOut = np.copy(imgIn)
  
      # iterative Yang algorithm:
      for i in range(0, itr):
          if Cnt['VERBOSE']: print 'i> PVC Yang iteration =', i
          # piece-wise constant image
          imgPWC = imgOut
          imgPWC[imgPWC<0] = 0
          for jr in range(0,m+1):
              imgPWC[imgSeg==jr] = np.mean( imgPWC[imgSeg==jr] )
          
          #> blur the piece-wise constant image using either:
          #> (1) GPU convolution with a separable kernel (x,y,z), or
          #> (2) CPU, Python-based convolution
          if 'CCARCH' in Cnt and 'compute' in Cnt['CCARCH']:
  
              #> convert to dimensions of GPU processing [y,x,z]
              imin_d = np.transpose(imgPWC, (1, 2, 0))
              imout_d = np.zeros(imin_d.shape, dtype=np.float32)
              improc.convolve(imout_d, imin_d, krnl, Cnt)
              imgSmo = np.transpose(imout_d, (2,0,1))
          else:
              hxy = np.outer(krnl[1,:], krnl[2,:])
              hxyz = np.multiply.outer(krnl[0,:], hxy)
              imgSmo = ndi.convolve(imgPWC, hxyz, mode='constant', cval=0.)
          
          # correction factors
          imgCrr = np.ones(dim, dtype=np.float32)
          imgCrr[imgSmo>0] = imgPWC[imgSmo>0] / imgSmo[imgSmo>0]
          imgOut = imgIn * imgCrr;
          for jr in range(0,m+1):
              m_a[jr, i] = np.mean( imgOut[imgSeg==jr] )
  
      return imgOut, m_a