Histogram of bias fields for a DWI

bias_hist.py

#!/usr/bin/env python

from matplotlib.backends.backend_pdf import PdfPages
from matplotlib import pyplot as plt
from conversion import read_bvals
import nibabel as nib
import sys
import numpy as np
from os.path import join as pjoin



def save_figure(temp, pp, bvals, i):
    count, bins= np.histogram(temp, bins=10)
    count=count/sum(count)

    plt.ylim([0,1])
    plt.bar(bins[:-1], height=count, width= bins[1]-bins[0])            
    plt.grid() 
    plt.title(f'Vol {i}, bvalue {bvals[i]}')
    plt.xlabel('Bias magnitude')
    plt.ylabel('Histogram converted to [0,1] range')
    plt.savefig(pp, format='pdf') 
    # plt.show() 
    plt.delaxes() 


def calc_hist():
    
    for bval in np.unique(bvals):
        print('bvalue', bval)        
        pp = PdfPages(outPrefix+ f'_bias_fields_b{bval}.pdf')

        for i in range(bias.shape[-1]):

            if bvals[i]!=bval:
                pass
                continue

            print('Vol',i)
            temp= data[...,i]
            save_figure(temp, pp, bvals, i)
                 
        pp.close()   
    
    print('bvalue all')
    pp = PdfPages(outPrefix+ '_bias_fields.pdf')

    for i in range(bias.shape[-1]):

        print('Vol',i)
        temp= data[...,i]
        save_figure(temp, pp, bvals, i)
             
    pp.close()   


if __name__=='__main__':
    bias= nib.load(sys.argv[1])
    inPrefix= sys.argv[1].split('_Bias.nii')[0]
    outPrefix= sys.argv[2]

    bvals= [int(b) for b in read_bvals(inPrefix+'.bval')]
    data= bias.get_fdata()
    calc_hist()

Usage ./bias_hist.py sub-XYZ*_Bias.nii.gz sub-XYZ