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
def iyang ( imgIn , krnl , imgSeg , Cnt , itr = 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 ] ) imgOut = np . copy ( imgIn ) for i in range ( 0 , itr ) : if Cnt [ 'VERBOSE' ] : print 'i> PVC Yang iteration =' , i imgPWC = imgOut imgPWC [ imgPWC < 0 ] = 0 for jr in range ( 0 , m + 1 ) : imgPWC [ imgSeg == jr ] = np . mean ( imgPWC [ imgSeg == jr ] ) if 'CCARCH' in Cnt and 'compute' in Cnt [ 'CCARCH' ] : 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. ) 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
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)