diff --git a/pynemaiqpet/__init__.py b/pynemaiqpet/__init__.py index 2e79b7b..19b4b72 100644 --- a/pynemaiqpet/__init__.py +++ b/pynemaiqpet/__init__.py @@ -3,4 +3,4 @@ try: __version__ = version("pynemaiqpet") except PackageNotFoundError: - __version__ = "unknown" \ No newline at end of file + __version__ = "unknown" diff --git a/pynemaiqpet/nema_smallanimal.py b/pynemaiqpet/nema_smallanimal.py index f77a7f2..6f8eadf 100644 --- a/pynemaiqpet/nema_smallanimal.py +++ b/pynemaiqpet/nema_smallanimal.py @@ -2,7 +2,7 @@ import numpy as np import pylab as py -import pymirc.metrics as pymr +import pymirc.metrics as pymr from pymirc.image_operations import kul_aff, aff_transform @@ -11,668 +11,750 @@ from scipy.ndimage.measurements import center_of_mass from scipy.ndimage.morphology import binary_fill_holes -from scipy.special import erf -from scipy.integrate import quad +from scipy.special import erf +from scipy.integrate import quad -from scipy.signal import argrelextrema +from scipy.signal import argrelextrema from scipy.optimize import minimize from lmfit import Model -#---------------------------------------------------------------------- + +# ---------------------------------------------------------------------- def cylinder_prof_integrand(eta, z, Z): - """ Integrand of the convolution of a disk convolved with a 2D Gaussian + """Integrand of the convolution of a disk convolved with a 2D Gaussian + + Parameters + ---------- + eta : float + integration variable - Parameters - ---------- - eta : float - integration variable + z : float + radial distance from center divided by sigma of Gaussian - z : float - radial distance from center divided by sigma of Gaussian + Z : float + disk radius divided by sigma of Gaussian - Z : float - disk radius divided by sigma of Gaussian + Returns + ------- + float + """ + return math.exp(-0.5 * eta**2) * ( + erf((math.sqrt(Z**2 - eta**2) - z) / math.sqrt(2)) + + erf((math.sqrt(Z**2 - eta**2) + z) / math.sqrt(2)) + ) - Returns - ------- - float - """ - return math.exp(-0.5*eta**2) * (erf( (math.sqrt(Z**2 - eta**2) - z)/math.sqrt(2) ) + - erf( (math.sqrt(Z**2 - eta**2) + z)/math.sqrt(2))) -#---------------------------------------------------------------------- +# ---------------------------------------------------------------------- def cylinder_prof(z, Z): - """ Normalized radial profile of a disk convolved with a 2D Gaussian - - Parameters - ---------- - z : float - radial distance from center divided by sigma of Gaussian - - Z : float - disk radius divided by sigma of Gaussian - - Returns - ------- - float - - Note - ---- - There is no analytic expression for the convolution. - The profile is numerically integrated using quad() from scipy.integrate - """ - return quad(cylinder_prof_integrand, 0, Z, args = (z,Z))[0] / math.sqrt(2*math.pi) - -#---------------------------------------------------------------------- -def cylinder_profile(r, S = 1., R = 1., fwhm = 1.): - """ Radial profile of a disk convolved with a 2D Gaussian - - Parameters - ---------- - r : 1D numpy float array - radial distance from center - - S : float - signal in the disk - - R : float - radius of the disk - - fwhm : float - FWHM of the Gaussian smoothing kernel - - Returns - ------- - 1D numpy float array - """ - sig = fwhm / 2.35 - - cp = np.frompyfunc(cylinder_prof, 2, 1) - - return S*cp(r/sig, R/sig).astype(float) - -#---------------------------------------------------------------------- -def fit_nema_2008_cylinder_profiles(vol, - voxsize, - Rrod_init = [2.5,2,1.5,1,0.5], - fwhm_init = 1.5, - S_init = 1, - fix_S = True, - fix_R = False, - fix_fwhm = False, - nrods = 4, - phantom = 'standard'): - """ Fit the radial profiles of the rods in a nema 2008 small animal PET phantom - - Parameters - ---------- - vol : 3D numpy float array - containing the image - - voxsize : 3 element 1D numpy array - containing the voxel size - - Rrod_init : list or 1D numpy array of floats, optional - containing the initial values of the rod radii - - S_init, fwhm_init: float, optional - initial values for the signal and the FWHM in the fit - - fix_S, fix_R, fix_fwhm : bool, optional - whether to keep the initial values of signal, radius and FWHM fixed during the fix - - nrods: int, optional - number of rods to fit - - phantom : string - phantom version ('standard' or 'mini') - - Returns - ------- - a list of lmfit fit results - - Note - ---- - - The axial direction should be the right most direction in the 3D numpy array. - The slices containing the rods are found automatically and summed. - In the summed image, all rods (disks) are segmented followed by a fit - of the radial profile. - """ - roi_vol = nema_2008_small_animal_pet_rois(vol, voxsize, phantom = phantom) - - rod_bbox = find_objects(roi_vol==4) - - # find the rods in the summed image - sum_img = vol[:,:,rod_bbox[0][2].start:rod_bbox[0][2].stop].mean(2) - - label_img, nlab = label(sum_img > 0.1*sum_img.max()) - labels = np.arange(1,nlab+1) - # sort the labels according to volume - npix = labeled_comprehension(sum_img, label_img, labels, len, int, 0) - sort_inds = npix.argsort()[::-1] - labels = labels[sort_inds] - npix = npix[sort_inds] - - #---------------------------------------------------------------------- - ncols = 2 - nrows = int(np.ceil(nrods/ncols)) - fig, ax = py.subplots(nrows,ncols,figsize = (12,7*nrows/2), sharey = True, sharex = True) - - retval = [] - - for irod in range(nrods): - rod_bbox = find_objects(label_img == labels[irod]) - - rod_bbox = [(slice(rod_bbox[0][0].start - 2,rod_bbox[0][0].stop + 2), - slice(rod_bbox[0][1].start - 2,rod_bbox[0][1].stop + 2))] - - rod_img = sum_img[rod_bbox[0]] - com = np.array(center_of_mass(rod_img)) - - x0 = (np.arange(rod_img.shape[0]) - com[0]) * voxsize[0] - x1 = (np.arange(rod_img.shape[1]) - com[1]) * voxsize[1] - - X0, X1 = np.meshgrid(x0, x1, indexing = 'ij') - RHO = np.sqrt(X0**2 + X1**2) - - rho = RHO.flatten() - signal = rod_img.flatten() - - # sort the values according to rho - sort_inds = rho.argsort() - rho = rho[sort_inds] - signal = signal[sort_inds] - - pmodel = Model(cylinder_profile) - params = pmodel.make_params(S = S_init, R = Rrod_init[irod], fwhm = fwhm_init) - - if fix_S: - params['S'].vary = False - if fix_R: - params['R'].vary = False - if fix_fwhm: - params['fwhm'].vary = False - - fitres = pmodel.fit(signal, r = rho, params = params) - retval.append(fitres) - fit_report = fitres.fit_report() - - iplot = np.unravel_index(irod, ax.shape) - ax[iplot].plot(rho,signal,'k.') - - rfit = np.linspace(0,rho.max(),100) - ax[iplot].plot(rfit,fitres.eval(r = rfit),'r-') - ax[iplot].text(0.99, 0.99, fit_report, fontsize = 6, transform = ax[iplot].transAxes, - verticalalignment='top', horizontalalignment = 'right', - backgroundcolor = 'white', bbox = {'pad':0, 'facecolor':'white','lw':0}) - ax[iplot].grid() - - for axx in ax[-1,:]: axx.set_xlabel('R (mm)') - for axx in ax[:,0]: axx.set_ylabel('signal') - - fig.tight_layout() - fig.show() - - return retval - -#---------------------------------------------------------------------- -def nema_2008_small_animal_pet_rois(vol, voxsize, lp_voxel = 'max', rod_th = 0.15, - phantom = 'standard'): - """ generate a label volume indicating the ROIs needed in the analysis of the - NEMA small animal PET IQ phantom - - Parameters - ---------- - vol : 3D numpy float array - containing the image - - voxsize : 3 element 1D numpy array - containing the voxel size - - lp_voxel: string, optional - method of how to compute the ROIs around the line profiles - in the summed images of the hot rods. - 'max' means the maximum voxels in the summed 2D image. - anything else means use all pixels that are within the rod radius - around the center of mass - - rod_th : float, optional - threshold to find the rod in the summed 2D image relative to the - mean of the big uniform region - - phantom : string - phantom version ('standard' or 'mini') - - Returns - ------- - a 3D integer numpy array - encoding the following ROIs: - 1 ... ROI of the big uniform region - 2 ... first cold insert - 3 ... second cold insert - 4 ... central line profile in 5mm rod - 5 ... central line profile in 4mm rod - 6 ... central line profile in 3mm rod - 7 ... central line profile in 2mm rod - 8 ... central line profile in 1mm rod - - Note - ---- - The rod ROIs in the summed 2D image are found by thresholding. - If the activity in the small rods is too low, they might be missed. - """ - roi_vol = np.zeros(vol.shape, dtype = np.uint) - - # calculate the summed z profile to place the ROIs - zprof = vol.sum(0).sum(0) - zprof_grad = np.gradient(zprof) - zprof_grad[np.abs(zprof_grad) < 0.13*np.abs(zprof_grad).max()] = 0 - - rising_edges = argrelextrema(zprof_grad, np.greater, order = 10)[0] - falling_edges = argrelextrema(zprof_grad, np.less, order = 10)[0] - - # if we only have 2 falling edges because the volume is cropped, we add the last slices as - # falling edge - - if falling_edges.shape[0] == 2: - falling_edges = np.concatenate([falling_edges,[vol.shape[2]]]) - - # define and analyze the big uniform ROI - uni_region_start_slice = rising_edges[1] - uni_region_end_slice = falling_edges[1] - uni_region_center_slice = 0.5*(uni_region_start_slice + uni_region_end_slice) - - uni_roi_start_slice = int(np.floor(uni_region_center_slice - 5./voxsize[2])) - uni_roi_end_slice = int(np.ceil(uni_region_center_slice + 5./voxsize[2])) - - uni_com = np.array(center_of_mass(vol[:,:,uni_roi_start_slice:(uni_roi_end_slice+1)])) - - x0 = (np.arange(vol.shape[0]) - uni_com[0]) * voxsize[0] - x1 = (np.arange(vol.shape[1]) - uni_com[1]) * voxsize[1] - x2 = (np.arange(vol.shape[2]) - uni_com[2]) * voxsize[2] - - X0,X1,X2 = np.meshgrid(x0,x1,x2,indexing='ij') - RHO = np.sqrt(X0**2 + X1**2) - - uni_mask = np.zeros(vol.shape, dtype = np.uint) - if phantom == 'standard': - uni_mask[RHO <= 11.25] = 1 - elif phantom == 'mini': - uni_mask[RHO <= 6.25] = 1 - uni_mask[:,:,:uni_roi_start_slice] = 0 - uni_mask[:,:,(uni_roi_end_slice+1):] = 0 - - uni_inds = np.where(uni_mask == 1) - roi_vol[uni_inds] = 1 - - # define and analyze the two cold ROIs - insert_region_start_slice = falling_edges[1] - insert_region_end_slice = falling_edges[2] - insert_region_center_slice = 0.5*(insert_region_start_slice + insert_region_end_slice) - - insert_roi_start_slice = int(np.floor(insert_region_center_slice - 3.75/voxsize[2])) - insert_roi_end_slice = int(np.ceil(insert_region_center_slice + 3.75/voxsize[2])) - - # sum the insert slices and subtract them from the max to find the two cold inserts - sum_insert_img = vol[:,:,insert_roi_start_slice:(insert_roi_end_slice+1)].mean(2) - - ref = np.percentile(sum_insert_img,99) - if phantom == 'standard': - insert_label_img, nlab_insert = label(sum_insert_img <= 0.5*ref) - elif phantom == 'mini': - # reset pixels outside the phantom, since inserts sometimes leak into background - tmp_inds = RHO[:,:,0] > 9 - sum_insert_img[tmp_inds] = ref - insert_label_img, nlab_insert = label(binary_erosion(sum_insert_img <= 0.5*ref)) - - # add backgroud low activity ROI to be compliant with standard phantom - insert_label_img[tmp_inds] = 3 - nlab_insert += 1 - - insert_labels = np.arange(1,nlab_insert+1) - # sort the labels according to volume - npix_insert = labeled_comprehension(sum_insert_img, insert_label_img, insert_labels, len, int, 0) - insert_sort_inds = npix_insert.argsort()[::-1] - insert_labels = insert_labels[insert_sort_inds] - npix_insert = npix_insert[insert_sort_inds] - - for i_insert in [1,2]: - tmp = insert_label_img.copy() - tmp[insert_label_img != insert_labels[i_insert]] = 0 - com_pixel = np.round(np.array(center_of_mass(tmp))) - - x0 = (np.arange(vol.shape[0]) - com_pixel[0]) * voxsize[0] - x1 = (np.arange(vol.shape[1]) - com_pixel[1]) * voxsize[1] - x2 = (np.arange(vol.shape[2])) * voxsize[2] - - X0,X1,X2 = np.meshgrid(x0,x1,x2,indexing='ij') - RHO = np.sqrt(X0**2 + X1**2) - - insert_mask = np.zeros(vol.shape, dtype = np.uint) - insert_mask[RHO <= 2] = 1 - insert_mask[:,:,:insert_roi_start_slice] = 0 - insert_mask[:,:,(insert_roi_end_slice+1):] = 0 - - insert_inds = np.where(insert_mask == 1) - roi_vol[insert_inds] = i_insert + 1 - - # find the rod z slices - rod_start_slice = falling_edges[0] - rod_end_slice = rising_edges[1] - rod_center = 0.5*(rod_start_slice + rod_end_slice) - - rod_roi_start_slice = int(np.floor(rod_center - 5./voxsize[2])) - rod_roi_end_slice = int(np.ceil(rod_center + 5./voxsize[2])) - - # sum the rod slices - sum_img = vol[:,:,rod_roi_start_slice:(rod_roi_end_slice+1)].mean(2) - - # label the summed image - label_img, nlab = label(sum_img > rod_th*sum_img.max()) - labels = np.arange(1,nlab+1) - - # sort the labels according to volume - npix = labeled_comprehension(sum_img, label_img, labels, len, int, 0) - sort_inds = npix.argsort()[::-1] - labels = labels[sort_inds] - npix = npix[sort_inds] - - # find the center for the line profiles - for i, lab in enumerate(labels): - - if lp_voxel == 'max': - rod_sum_img = sum_img.copy() - rod_sum_img[label_img != lab] = 0 - - central_pixel = np.unravel_index(rod_sum_img.argmax(),rod_sum_img.shape) - roi_vol[central_pixel[0],central_pixel[1],rod_roi_start_slice:(rod_roi_end_slice+1)] = i + 4 - else: - bbox = find_objects(label_img == lab)[0] - rod_sum_img = np.zeros(sum_img.shape) - rod_sum_img[bbox] = sum_img[bbox] + """Normalized radial profile of a disk convolved with a 2D Gaussian + + Parameters + ---------- + z : float + radial distance from center divided by sigma of Gaussian + + Z : float + disk radius divided by sigma of Gaussian + + Returns + ------- + float + + Note + ---- + There is no analytic expression for the convolution. + The profile is numerically integrated using quad() from scipy.integrate + """ + return quad(cylinder_prof_integrand, 0, Z, args=(z, Z))[0] / math.sqrt(2 * math.pi) + + +# ---------------------------------------------------------------------- +def cylinder_profile(r, S=1.0, R=1.0, fwhm=1.0): + """Radial profile of a disk convolved with a 2D Gaussian - central_pixel = np.round(np.array(center_of_mass(rod_sum_img))).astype(np.int) + Parameters + ---------- + r : 1D numpy float array + radial distance from center - x0 = (np.arange(rod_sum_img.shape[0]) - central_pixel[0]) * voxsize[0] - x1 = (np.arange(rod_sum_img.shape[1]) - central_pixel[1]) * voxsize[1] - X0, X1 = np.meshgrid(x0, x1, indexing = 'ij') + S : float + signal in the disk + + R : float + radius of the disk + + fwhm : float + FWHM of the Gaussian smoothing kernel + + Returns + ------- + 1D numpy float array + """ + sig = fwhm / 2.35 - dist_to_rod = np.sqrt(X0**2 + X1**2) + cp = np.frompyfunc(cylinder_prof, 2, 1) - for r in range(rod_roi_start_slice,(rod_roi_end_slice+1)): - roi_vol[...,r][dist_to_rod <= 0.5*(5-i)] = i + 4 + return S * cp(r / sig, R / sig).astype(float) - #------------------------------------------------------- - # if we only have 4 labels (rods), we find the last (smallest) one based on symmetries - if nlab == 4: - roi_img = roi_vol[...,rod_roi_start_slice] - com = center_of_mass(roi_vol == 1) - x0 = (np.arange(sum_img.shape[0]) - com[0]) * voxsize[0] - x1 = (np.arange(sum_img.shape[1]) - com[1]) * voxsize[1] - X0,X1 = np.meshgrid(x0, x1, indexing = 'ij') +# ---------------------------------------------------------------------- +def fit_nema_2008_cylinder_profiles( + vol, + voxsize, + Rrod_init=[2.5, 2, 1.5, 1, 0.5], + fwhm_init=1.5, + S_init=1, + fix_S=True, + fix_R=False, + fix_fwhm=False, + nrods=4, + phantom="standard", +): + """Fit the radial profiles of the rods in a nema 2008 small animal PET phantom + + Parameters + ---------- + vol : 3D numpy float array + containing the image + + voxsize : 3 element 1D numpy array + containing the voxel size + + Rrod_init : list or 1D numpy array of floats, optional + containing the initial values of the rod radii + + S_init, fwhm_init: float, optional + initial values for the signal and the FWHM in the fit + + fix_S, fix_R, fix_fwhm : bool, optional + whether to keep the initial values of signal, radius and FWHM fixed during the fix + + nrods: int, optional + number of rods to fit + + phantom : string + phantom version ('standard' or 'mini') + + Returns + ------- + a list of lmfit fit results + + Note + ---- + + The axial direction should be the right most direction in the 3D numpy array. + The slices containing the rods are found automatically and summed. + In the summed image, all rods (disks) are segmented followed by a fit + of the radial profile. + """ + roi_vol = nema_2008_small_animal_pet_rois(vol, voxsize, phantom=phantom) + + rod_bbox = find_objects(roi_vol == 4) + + # find the rods in the summed image + sum_img = vol[:, :, rod_bbox[0][2].start : rod_bbox[0][2].stop].mean(2) + + label_img, nlab = label(sum_img > 0.1 * sum_img.max()) + labels = np.arange(1, nlab + 1) + # sort the labels according to volume + npix = labeled_comprehension(sum_img, label_img, labels, len, int, 0) + sort_inds = npix.argsort()[::-1] + labels = labels[sort_inds] + npix = npix[sort_inds] + + # ---------------------------------------------------------------------- + ncols = 2 + nrows = int(np.ceil(nrods / ncols)) + fig, ax = py.subplots( + nrows, ncols, figsize=(12, 7 * nrows / 2), sharey=True, sharex=True + ) + + retval = [] + + for irod in range(nrods): + rod_bbox = find_objects(label_img == labels[irod]) + + rod_bbox = [ + ( + slice(rod_bbox[0][0].start - 2, rod_bbox[0][0].stop + 2), + slice(rod_bbox[0][1].start - 2, rod_bbox[0][1].stop + 2), + ) + ] + + rod_img = sum_img[rod_bbox[0]] + com = np.array(center_of_mass(rod_img)) + + x0 = (np.arange(rod_img.shape[0]) - com[0]) * voxsize[0] + x1 = (np.arange(rod_img.shape[1]) - com[1]) * voxsize[1] + + X0, X1 = np.meshgrid(x0, x1, indexing="ij") + RHO = np.sqrt(X0**2 + X1**2) + + rho = RHO.flatten() + signal = rod_img.flatten() + + # sort the values according to rho + sort_inds = rho.argsort() + rho = rho[sort_inds] + signal = signal[sort_inds] + + pmodel = Model(cylinder_profile) + params = pmodel.make_params(S=S_init, R=Rrod_init[irod], fwhm=fwhm_init) + + if fix_S: + params["S"].vary = False + if fix_R: + params["R"].vary = False + if fix_fwhm: + params["fwhm"].vary = False + + fitres = pmodel.fit(signal, r=rho, params=params) + retval.append(fitres) + fit_report = fitres.fit_report() + + iplot = np.unravel_index(irod, ax.shape) + ax[iplot].plot(rho, signal, "k.") + + rfit = np.linspace(0, rho.max(), 100) + ax[iplot].plot(rfit, fitres.eval(r=rfit), "r-") + ax[iplot].text( + 0.99, + 0.99, + fit_report, + fontsize=6, + transform=ax[iplot].transAxes, + verticalalignment="top", + horizontalalignment="right", + backgroundcolor="white", + bbox={"pad": 0, "facecolor": "white", "lw": 0}, + ) + ax[iplot].grid() + + for axx in ax[-1, :]: + axx.set_xlabel("R (mm)") + for axx in ax[:, 0]: + axx.set_ylabel("signal") + + fig.tight_layout() + fig.show() + + return retval + + +# ---------------------------------------------------------------------- +def nema_2008_small_animal_pet_rois( + vol, voxsize, lp_voxel="max", rod_th=0.15, phantom="standard" +): + """generate a label volume indicating the ROIs needed in the analysis of the + NEMA small animal PET IQ phantom + + Parameters + ---------- + vol : 3D numpy float array + containing the image + + voxsize : 3 element 1D numpy array + containing the voxel size + + lp_voxel: string, optional + method of how to compute the ROIs around the line profiles + in the summed images of the hot rods. + 'max' means the maximum voxels in the summed 2D image. + anything else means use all pixels that are within the rod radius + around the center of mass + + rod_th : float, optional + threshold to find the rod in the summed 2D image relative to the + mean of the big uniform region + + phantom : string + phantom version ('standard' or 'mini') + + Returns + ------- + a 3D integer numpy array + encoding the following ROIs: + 1 ... ROI of the big uniform region + 2 ... first cold insert + 3 ... second cold insert + 4 ... central line profile in 5mm rod + 5 ... central line profile in 4mm rod + 6 ... central line profile in 3mm rod + 7 ... central line profile in 2mm rod + 8 ... central line profile in 1mm rod + + Note + ---- + The rod ROIs in the summed 2D image are found by thresholding. + If the activity in the small rods is too low, they might be missed. + """ + roi_vol = np.zeros(vol.shape, dtype=np.uint) + + # calculate the summed z profile to place the ROIs + zprof = vol.sum(0).sum(0) + zprof_grad = np.gradient(zprof) + zprof_grad[np.abs(zprof_grad) < 0.13 * np.abs(zprof_grad).max()] = 0 + + rising_edges = argrelextrema(zprof_grad, np.greater, order=10)[0] + falling_edges = argrelextrema(zprof_grad, np.less, order=10)[0] + + # if we only have 2 falling edges because the volume is cropped, we add the last slices as + # falling edge + + if falling_edges.shape[0] == 2: + falling_edges = np.concatenate([falling_edges, [vol.shape[2]]]) + + # define and analyze the big uniform ROI + uni_region_start_slice = rising_edges[1] + uni_region_end_slice = falling_edges[1] + uni_region_center_slice = 0.5 * (uni_region_start_slice + uni_region_end_slice) + + uni_roi_start_slice = int(np.floor(uni_region_center_slice - 5.0 / voxsize[2])) + uni_roi_end_slice = int(np.ceil(uni_region_center_slice + 5.0 / voxsize[2])) + + uni_com = np.array( + center_of_mass(vol[:, :, uni_roi_start_slice : (uni_roi_end_slice + 1)]) + ) + + x0 = (np.arange(vol.shape[0]) - uni_com[0]) * voxsize[0] + x1 = (np.arange(vol.shape[1]) - uni_com[1]) * voxsize[1] + x2 = (np.arange(vol.shape[2]) - uni_com[2]) * voxsize[2] + + X0, X1, X2 = np.meshgrid(x0, x1, x2, indexing="ij") RHO = np.sqrt(X0**2 + X1**2) - PHI = np.arctan2(X1,X0) - rod_phis = np.array([PHI[roi_img == x][0] for x in np.arange(4,nlab+4)]) - PHI = ((PHI - rod_phis[3]) % (2*np.pi)) - np.pi - rod_phis = ((rod_phis - rod_phis[3]) % (2*np.pi)) - np.pi - - missing_phi = ((rod_phis[3] - rod_phis[2]) % (2*np.pi)) - np.pi - - mask = np.logical_and(np.abs(PHI - missing_phi) < 0.25, np.abs(RHO - 6.4) < 2) - - central_pixel = np.unravel_index(np.argmax(sum_img*mask), sum_img.shape) - if lp_voxel == 'max': - roi_vol[central_pixel[0],central_pixel[1],rod_roi_start_slice:(rod_roi_end_slice+1)] = 8 - else: - x0 = (np.arange(rod_sum_img.shape[0]) - central_pixel[0]) * voxsize[0] - x1 = (np.arange(rod_sum_img.shape[1]) - central_pixel[1]) * voxsize[1] - X0, X1 = np.meshgrid(x0, x1, indexing = 'ij') - - dist_to_rod = np.sqrt(X0**2 + X1**2) - - for r in range(rod_roi_start_slice,(rod_roi_end_slice+1)): - roi_vol[...,r][dist_to_rod <= 0.5] = 8 - - nlab += 1 - #------------------------------------------------------- - - - return roi_vol - -#-------------------------------------------------------------------- -def nema_2008_small_animal_iq_phantom(voxsize, shape, version = 'standard'): - """ generate a digital version of the upper part of the NEMA small animal PET - IQ phantom that can be used to align a NEMA scan - - Parameters - ---------- - voxsize : 3 element 1D numpy array - containing the voxel size - - shape: 3 element tuple of integers - shape of the volume - - version : string - phantom version ('standard' or 'mini') - - Returns - ------- - a 3D numpy array - """ - x0 = (np.arange(shape[0]) - 0.5*shape[0] - 0.5) * voxsize[0] - x1 = (np.arange(shape[1]) - 0.5*shape[1] - 0.5) * voxsize[1] - x2 = (np.arange(shape[2]) - 0.5*shape[2] - 0.5) * voxsize[2] - - X0,X1,X2 = np.meshgrid(x0,x1,x2,indexing='ij') - RHO = np.sqrt(X0**2 + X1**2) - - phantom = np.zeros(shape) - if version == 'standard': - phantom[RHO <= 30./2] = 1 - phantom[X2 < 0] = 0 - phantom[X2 > 33.] = 0 - - RHO1 = np.sqrt((X0 - 7.5)**2 + X1**2) - RHO2 = np.sqrt((X0 + 7.5)**2 + X1**2) - - #phantom[np.logical_and(RHO1 <= 9.2/2, X2 > 15)] = 0 - #phantom[np.logical_and(RHO2 <= 9.2/2, X2 > 15)] = 0 - phantom[np.logical_and(RHO1 <= 11/2, X2 > 15)] = 0 - phantom[np.logical_and(RHO2 <= 11/2, X2 > 15)] = 0 - elif version == 'mini': - phantom[RHO <= 20./2] = 1 - phantom[X2 < 0] = 0 - phantom[X2 > 34.] = 0 - - RHO1 = np.sqrt((X0 - 6)**2 + X1**2) - RHO2 = np.sqrt((X0 + 6)**2 + X1**2) - - phantom[np.logical_and(RHO1 <= 7.2/2, X2 > 16)] = 0 - phantom[np.logical_and(RHO2 <= 7.2/2, X2 > 16)] = 0 - - return phantom - -#-------------------------------------------------------------------- -def align_nema_2008_small_animal_iq_phantom(vol, voxsize, ftol = 1e-2, xtol = 1e-2, maxiter = 10, maxfev = 500, version = 'standard'): - """ align a reconstruction of the NEMA small animal PET IQ phantom to its digital version - - Parameters - ---------- - vol : 3D numpy float array - containing the image - - voxsize : 3 element 1D numpy array - containing the voxel size - - ftol, xtol, maxiter, maxfev : float / int - parameter for the optimizer used to minimze the cost function - - version : string - phantom version ('standard' or 'mini') - - Returns - ------- - a 3D numpy array - - Note - ---- - This routine can be useful to make sure that the rods in the NEMA scan are - parallel to the axial direction. - """ - phantom = nema_2008_small_animal_iq_phantom(voxsize, vol.shape, version = version) - phantom *= vol[vol>0.5*vol.max()].mean() - - reg_params = np.zeros(6) - - # registration of down sampled volumes - dsf = 3 - ds_aff = np.diag([dsf,dsf,dsf,1.]) - - phantom_ds = aff_transform(phantom, ds_aff, np.ceil(np.array(phantom.shape)/dsf).astype(int)) - - res = minimize(pymr.regis_cost_func, reg_params, - args = (phantom_ds, vol, True, True, lambda x,y: ((x-y)**2).mean(), ds_aff), - method = 'Powell', - options = {'ftol':ftol, 'xtol':xtol, 'disp':True, 'maxiter':maxiter, 'maxfev':maxfev}) - - reg_params = res.x.copy() - # we have to scale the translations by the down sample factor since they are in voxels - reg_params[:3] *= dsf - - res = minimize(pymr.regis_cost_func, reg_params, - args = (phantom, vol, True, True, lambda x,y: ((x-y)**2).mean()), - method = 'Powell', - options = {'ftol':ftol, 'xtol':xtol, 'disp':True, 'maxiter':maxiter, 'maxfev':maxfev}) - - regis_aff = kul_aff(res.x, origin = np.array(vol.shape)/2) - vol_aligned = aff_transform(vol, regis_aff, vol.shape) - - return vol_aligned - -#-------------------------------------------------------------------- + uni_mask = np.zeros(vol.shape, dtype=np.uint) + if phantom == "standard": + uni_mask[RHO <= 11.25] = 1 + elif phantom == "mini": + uni_mask[RHO <= 6.25] = 1 + uni_mask[:, :, :uni_roi_start_slice] = 0 + uni_mask[:, :, (uni_roi_end_slice + 1) :] = 0 + + uni_inds = np.where(uni_mask == 1) + roi_vol[uni_inds] = 1 + + # define and analyze the two cold ROIs + insert_region_start_slice = falling_edges[1] + insert_region_end_slice = falling_edges[2] + insert_region_center_slice = 0.5 * ( + insert_region_start_slice + insert_region_end_slice + ) + + insert_roi_start_slice = int( + np.floor(insert_region_center_slice - 3.75 / voxsize[2]) + ) + insert_roi_end_slice = int(np.ceil(insert_region_center_slice + 3.75 / voxsize[2])) + + # sum the insert slices and subtract them from the max to find the two cold inserts + sum_insert_img = vol[ + :, :, insert_roi_start_slice : (insert_roi_end_slice + 1) + ].mean(2) + + ref = np.percentile(sum_insert_img, 99) + if phantom == "standard": + insert_label_img, nlab_insert = label(sum_insert_img <= 0.5 * ref) + elif phantom == "mini": + # reset pixels outside the phantom, since inserts sometimes leak into background + tmp_inds = RHO[:, :, 0] > 9 + sum_insert_img[tmp_inds] = ref + insert_label_img, nlab_insert = label( + binary_erosion(sum_insert_img <= 0.5 * ref) + ) + + # add backgroud low activity ROI to be compliant with standard phantom + insert_label_img[tmp_inds] = 3 + nlab_insert += 1 + + insert_labels = np.arange(1, nlab_insert + 1) + # sort the labels according to volume + npix_insert = labeled_comprehension( + sum_insert_img, insert_label_img, insert_labels, len, int, 0 + ) + insert_sort_inds = npix_insert.argsort()[::-1] + insert_labels = insert_labels[insert_sort_inds] + npix_insert = npix_insert[insert_sort_inds] + + for i_insert in [1, 2]: + tmp = insert_label_img.copy() + tmp[insert_label_img != insert_labels[i_insert]] = 0 + com_pixel = np.round(np.array(center_of_mass(tmp))) + + x0 = (np.arange(vol.shape[0]) - com_pixel[0]) * voxsize[0] + x1 = (np.arange(vol.shape[1]) - com_pixel[1]) * voxsize[1] + x2 = (np.arange(vol.shape[2])) * voxsize[2] + + X0, X1, X2 = np.meshgrid(x0, x1, x2, indexing="ij") + RHO = np.sqrt(X0**2 + X1**2) + + insert_mask = np.zeros(vol.shape, dtype=np.uint) + insert_mask[RHO <= 2] = 1 + insert_mask[:, :, :insert_roi_start_slice] = 0 + insert_mask[:, :, (insert_roi_end_slice + 1) :] = 0 + + insert_inds = np.where(insert_mask == 1) + roi_vol[insert_inds] = i_insert + 1 + + # find the rod z slices + rod_start_slice = falling_edges[0] + rod_end_slice = rising_edges[1] + rod_center = 0.5 * (rod_start_slice + rod_end_slice) + + rod_roi_start_slice = int(np.floor(rod_center - 5.0 / voxsize[2])) + rod_roi_end_slice = int(np.ceil(rod_center + 5.0 / voxsize[2])) + + # sum the rod slices + sum_img = vol[:, :, rod_roi_start_slice : (rod_roi_end_slice + 1)].mean(2) + + # label the summed image + label_img, nlab = label(sum_img > rod_th * sum_img.max()) + labels = np.arange(1, nlab + 1) + + # sort the labels according to volume + npix = labeled_comprehension(sum_img, label_img, labels, len, int, 0) + sort_inds = npix.argsort()[::-1] + labels = labels[sort_inds] + npix = npix[sort_inds] + + # find the center for the line profiles + for i, lab in enumerate(labels): + + if lp_voxel == "max": + rod_sum_img = sum_img.copy() + rod_sum_img[label_img != lab] = 0 + + central_pixel = np.unravel_index(rod_sum_img.argmax(), rod_sum_img.shape) + roi_vol[ + central_pixel[0], + central_pixel[1], + rod_roi_start_slice : (rod_roi_end_slice + 1), + ] = ( + i + 4 + ) + else: + bbox = find_objects(label_img == lab)[0] + rod_sum_img = np.zeros(sum_img.shape) + rod_sum_img[bbox] = sum_img[bbox] + + central_pixel = np.round(np.array(center_of_mass(rod_sum_img))).astype( + np.int + ) + + x0 = (np.arange(rod_sum_img.shape[0]) - central_pixel[0]) * voxsize[0] + x1 = (np.arange(rod_sum_img.shape[1]) - central_pixel[1]) * voxsize[1] + X0, X1 = np.meshgrid(x0, x1, indexing="ij") + + dist_to_rod = np.sqrt(X0**2 + X1**2) + + for r in range(rod_roi_start_slice, (rod_roi_end_slice + 1)): + roi_vol[..., r][dist_to_rod <= 0.5 * (5 - i)] = i + 4 + + # ------------------------------------------------------- + # if we only have 4 labels (rods), we find the last (smallest) one based on symmetries + if nlab == 4: + roi_img = roi_vol[..., rod_roi_start_slice] + + com = center_of_mass(roi_vol == 1) + x0 = (np.arange(sum_img.shape[0]) - com[0]) * voxsize[0] + x1 = (np.arange(sum_img.shape[1]) - com[1]) * voxsize[1] + X0, X1 = np.meshgrid(x0, x1, indexing="ij") + RHO = np.sqrt(X0**2 + X1**2) + + PHI = np.arctan2(X1, X0) + rod_phis = np.array([PHI[roi_img == x][0] for x in np.arange(4, nlab + 4)]) + PHI = ((PHI - rod_phis[3]) % (2 * np.pi)) - np.pi + rod_phis = ((rod_phis - rod_phis[3]) % (2 * np.pi)) - np.pi + + missing_phi = ((rod_phis[3] - rod_phis[2]) % (2 * np.pi)) - np.pi + + mask = np.logical_and(np.abs(PHI - missing_phi) < 0.25, np.abs(RHO - 6.4) < 2) + + central_pixel = np.unravel_index(np.argmax(sum_img * mask), sum_img.shape) + if lp_voxel == "max": + roi_vol[ + central_pixel[0], + central_pixel[1], + rod_roi_start_slice : (rod_roi_end_slice + 1), + ] = 8 + else: + x0 = (np.arange(rod_sum_img.shape[0]) - central_pixel[0]) * voxsize[0] + x1 = (np.arange(rod_sum_img.shape[1]) - central_pixel[1]) * voxsize[1] + X0, X1 = np.meshgrid(x0, x1, indexing="ij") + + dist_to_rod = np.sqrt(X0**2 + X1**2) + + for r in range(rod_roi_start_slice, (rod_roi_end_slice + 1)): + roi_vol[..., r][dist_to_rod <= 0.5] = 8 + + nlab += 1 + # ------------------------------------------------------- + + return roi_vol + + +# -------------------------------------------------------------------- +def nema_2008_small_animal_iq_phantom(voxsize, shape, version="standard"): + """generate a digital version of the upper part of the NEMA small animal PET + IQ phantom that can be used to align a NEMA scan + + Parameters + ---------- + voxsize : 3 element 1D numpy array + containing the voxel size + + shape: 3 element tuple of integers + shape of the volume + + version : string + phantom version ('standard' or 'mini') + + Returns + ------- + a 3D numpy array + """ + x0 = (np.arange(shape[0]) - 0.5 * shape[0] - 0.5) * voxsize[0] + x1 = (np.arange(shape[1]) - 0.5 * shape[1] - 0.5) * voxsize[1] + x2 = (np.arange(shape[2]) - 0.5 * shape[2] - 0.5) * voxsize[2] + + X0, X1, X2 = np.meshgrid(x0, x1, x2, indexing="ij") + RHO = np.sqrt(X0**2 + X1**2) + + phantom = np.zeros(shape) + if version == "standard": + phantom[RHO <= 30.0 / 2] = 1 + phantom[X2 < 0] = 0 + phantom[X2 > 33.0] = 0 + + RHO1 = np.sqrt((X0 - 7.5) ** 2 + X1**2) + RHO2 = np.sqrt((X0 + 7.5) ** 2 + X1**2) + + # phantom[np.logical_and(RHO1 <= 9.2/2, X2 > 15)] = 0 + # phantom[np.logical_and(RHO2 <= 9.2/2, X2 > 15)] = 0 + phantom[np.logical_and(RHO1 <= 11 / 2, X2 > 15)] = 0 + phantom[np.logical_and(RHO2 <= 11 / 2, X2 > 15)] = 0 + elif version == "mini": + phantom[RHO <= 20.0 / 2] = 1 + phantom[X2 < 0] = 0 + phantom[X2 > 34.0] = 0 + + RHO1 = np.sqrt((X0 - 6) ** 2 + X1**2) + RHO2 = np.sqrt((X0 + 6) ** 2 + X1**2) + + phantom[np.logical_and(RHO1 <= 7.2 / 2, X2 > 16)] = 0 + phantom[np.logical_and(RHO2 <= 7.2 / 2, X2 > 16)] = 0 + + return phantom + + +# -------------------------------------------------------------------- +def align_nema_2008_small_animal_iq_phantom( + vol, voxsize, ftol=1e-2, xtol=1e-2, maxiter=10, maxfev=500, version="standard" +): + """align a reconstruction of the NEMA small animal PET IQ phantom to its digital version + + Parameters + ---------- + vol : 3D numpy float array + containing the image + + voxsize : 3 element 1D numpy array + containing the voxel size + + ftol, xtol, maxiter, maxfev : float / int + parameter for the optimizer used to minimze the cost function + + version : string + phantom version ('standard' or 'mini') + + Returns + ------- + a 3D numpy array + + Note + ---- + This routine can be useful to make sure that the rods in the NEMA scan are + parallel to the axial direction. + """ + phantom = nema_2008_small_animal_iq_phantom(voxsize, vol.shape, version=version) + phantom *= vol[vol > 0.5 * vol.max()].mean() + + reg_params = np.zeros(6) + + # registration of down sampled volumes + dsf = 3 + ds_aff = np.diag([dsf, dsf, dsf, 1.0]) + + phantom_ds = aff_transform( + phantom, ds_aff, np.ceil(np.array(phantom.shape) / dsf).astype(int) + ) + + res = minimize( + pymr.regis_cost_func, + reg_params, + args=(phantom_ds, vol, True, True, lambda x, y: ((x - y) ** 2).mean(), ds_aff), + method="Powell", + options={ + "ftol": ftol, + "xtol": xtol, + "disp": True, + "maxiter": maxiter, + "maxfev": maxfev, + }, + ) + + reg_params = res.x.copy() + # we have to scale the translations by the down sample factor since they are in voxels + reg_params[:3] *= dsf + + res = minimize( + pymr.regis_cost_func, + reg_params, + args=(phantom, vol, True, True, lambda x, y: ((x - y) ** 2).mean()), + method="Powell", + options={ + "ftol": ftol, + "xtol": xtol, + "disp": True, + "maxiter": maxiter, + "maxfev": maxfev, + }, + ) + + regis_aff = kul_aff(res.x, origin=np.array(vol.shape) / 2) + vol_aligned = aff_transform(vol, regis_aff, vol.shape) + + return vol_aligned + + +# -------------------------------------------------------------------- def nema_2008_small_animal_iq_phantom_report(vol, roi_vol): - """ generate the report for the NEMA 2008 small animal PET IQ phantom analysis - - Parameters - ---------- - vol : 3D numpy float array - containing the image - - roi_vol : 3D numpy float array - containing the ROI label image with following ROIs: - 1 ... ROI of the big uniform region - 2 ... first cold insert - 3 ... second cold insert - 4 ... central line profile in 5mm rod - 5 ... central line profile in 4mm rod - 6 ... central line profile in 3mm rod - 7 ... central line profile in 2mm rod - 8 ... central line profile in 1mm rod - - Returns - ------- - a 3D numpy array - - Note - ---- - The ROIs for the smaller rods are optional. - """ - np.set_printoptions(precision=3) - - # get the ROI values of the big uniform ROI - uni_values = vol[roi_vol == 1] - uni_mean = uni_values.mean() - uni_max = uni_values.max() - uni_min = uni_values.min() - uni_std = uni_values.std() - uni_perc_std = 100 * uni_std / uni_mean - - print("\nuniform ROI results") - print("------------------------------") - print("mean ...:", "%.3f" % uni_mean) - print("max ...:", "%.3f" % uni_max) - print("min ...:", "%.3f" % uni_min) - print("%std ...:", "%.3f" % uni_perc_std, "\n") - - # get the ROI values of the 2 cold inserts - insert_mean = np.zeros(2) - insert_std = np.zeros(2) - - insert1_values = vol[roi_vol == 2] - insert_mean[0] = insert1_values.mean() - insert_std[0] = insert1_values.std() - - insert2_values = vol[roi_vol == 3] - insert_mean[1] = insert2_values.mean() - insert_std[1] = insert2_values.std() - - insert_ratio = insert_mean / uni_mean - insert_perc_std = 100 * np.sqrt((insert_std/insert_mean)**2 + (uni_std/uni_mean)**2) - - print("\ncold insert results") - print("------------------------------") - print("spill over ratio ...:", insert_ratio) - print("%std ...:", insert_perc_std, "\n") - - # analyze the rod profiles - nrods = int(roi_vol.max() - 3) - lp_mean = np.zeros(nrods) - lp_std = np.zeros(nrods) - - # find the center for the line profiles - for i in range(nrods): - lp_values = vol[roi_vol == i + 4] - lp_mean[i] = lp_values.mean() - lp_std[i] = lp_values.std() - - lp_rc = lp_mean / uni_mean - lp_perc_std = 100 * np.sqrt((lp_std/lp_mean)**2 + (uni_std/uni_mean)**2) - - print("\nrod results") - print("------------------------------") - print("mean ...:", lp_mean) - print("recovery coeff...:", lp_rc) - print("%std ...:", lp_perc_std, "\n") - - - np.set_printoptions(precision=None) - - res = {'uniform_ROI_mean':uni_mean, - 'uniform_ROI_max':uni_max, - 'uniform_ROI_min':uni_min, - 'uniform_ROI_perc_std':uni_perc_std, - 'spill_over_ratio':insert_ratio, - 'spill_over_perc_std':insert_perc_std, - 'rod_mean':lp_mean, - 'rod_recovery_coeff':lp_rc, - 'rod_perc_std':lp_perc_std} - - return res - -#-------------------------------------------------------------------- + """generate the report for the NEMA 2008 small animal PET IQ phantom analysis + + Parameters + ---------- + vol : 3D numpy float array + containing the image + + roi_vol : 3D numpy float array + containing the ROI label image with following ROIs: + 1 ... ROI of the big uniform region + 2 ... first cold insert + 3 ... second cold insert + 4 ... central line profile in 5mm rod + 5 ... central line profile in 4mm rod + 6 ... central line profile in 3mm rod + 7 ... central line profile in 2mm rod + 8 ... central line profile in 1mm rod + + Returns + ------- + a 3D numpy array + + Note + ---- + The ROIs for the smaller rods are optional. + """ + np.set_printoptions(precision=3) + + # get the ROI values of the big uniform ROI + uni_values = vol[roi_vol == 1] + uni_mean = uni_values.mean() + uni_max = uni_values.max() + uni_min = uni_values.min() + uni_std = uni_values.std() + uni_perc_std = 100 * uni_std / uni_mean + + print("\nuniform ROI results") + print("------------------------------") + print("mean ...:", "%.3f" % uni_mean) + print("max ...:", "%.3f" % uni_max) + print("min ...:", "%.3f" % uni_min) + print("%std ...:", "%.3f" % uni_perc_std, "\n") + + # get the ROI values of the 2 cold inserts + insert_mean = np.zeros(2) + insert_std = np.zeros(2) + + insert1_values = vol[roi_vol == 2] + insert_mean[0] = insert1_values.mean() + insert_std[0] = insert1_values.std() + + insert2_values = vol[roi_vol == 3] + insert_mean[1] = insert2_values.mean() + insert_std[1] = insert2_values.std() + + insert_ratio = insert_mean / uni_mean + insert_perc_std = 100 * np.sqrt( + (insert_std / insert_mean) ** 2 + (uni_std / uni_mean) ** 2 + ) + + print("\ncold insert results") + print("------------------------------") + print("spill over ratio ...:", insert_ratio) + print("%std ...:", insert_perc_std, "\n") + + # analyze the rod profiles + nrods = int(roi_vol.max() - 3) + lp_mean = np.zeros(nrods) + lp_std = np.zeros(nrods) + + # find the center for the line profiles + for i in range(nrods): + lp_values = vol[roi_vol == i + 4] + lp_mean[i] = lp_values.mean() + lp_std[i] = lp_values.std() + + lp_rc = lp_mean / uni_mean + lp_perc_std = 100 * np.sqrt((lp_std / lp_mean) ** 2 + (uni_std / uni_mean) ** 2) + + print("\nrod results") + print("------------------------------") + print("mean ...:", lp_mean) + print("recovery coeff...:", lp_rc) + print("%std ...:", lp_perc_std, "\n") + + np.set_printoptions(precision=None) + + res = { + "uniform_ROI_mean": uni_mean, + "uniform_ROI_max": uni_max, + "uniform_ROI_min": uni_min, + "uniform_ROI_perc_std": uni_perc_std, + "spill_over_ratio": insert_ratio, + "spill_over_perc_std": insert_perc_std, + "rod_mean": lp_mean, + "rod_recovery_coeff": lp_rc, + "rod_perc_std": lp_perc_std, + } + + return res + + +# -------------------------------------------------------------------- def get_phantom_name(vol, voxsize): - """derive NEMA small animal phantom version from volume""" - zprof = vol.sum(0).sum(0) - zprof_grad = np.gradient(zprof) - zprof_grad[np.abs(zprof_grad) < 0.1*np.abs(zprof_grad).max()] = 0 - - rising_edges = argrelextrema(zprof_grad, np.greater, order = 10)[0] - falling_edges = argrelextrema(zprof_grad, np.less, order = 10)[0] + """derive NEMA small animal phantom version from volume""" + zprof = vol.sum(0).sum(0) + zprof_grad = np.gradient(zprof) + zprof_grad[np.abs(zprof_grad) < 0.1 * np.abs(zprof_grad).max()] = 0 + + rising_edges = argrelextrema(zprof_grad, np.greater, order=10)[0] + falling_edges = argrelextrema(zprof_grad, np.less, order=10)[0] - if falling_edges.shape[0] < 2: - falling_edges = np.concatenate([falling_edges,[vol.shape[2]]]) + if falling_edges.shape[0] < 2: + falling_edges = np.concatenate([falling_edges, [vol.shape[2]]]) - center_sl = (rising_edges[1] + falling_edges[1]) // 2 + center_sl = (rising_edges[1] + falling_edges[1]) // 2 - area = (vol[:,:,center_sl] > 0.5*np.percentile(vol[:,:,center_sl],99)).sum()*voxsize[0]*voxsize[1] + area = ( + (vol[:, :, center_sl] > 0.5 * np.percentile(vol[:, :, center_sl], 99)).sum() + * voxsize[0] + * voxsize[1] + ) - if area > 450: - phantom_name = 'standard' - else: - phantom_name = 'mini' + if area > 450: + phantom_name = "standard" + else: + phantom_name = "mini" - return phantom_name + return phantom_name diff --git a/pynemaiqpet/nema_wb.py b/pynemaiqpet/nema_wb.py index 97cf129..4cf3940 100644 --- a/pynemaiqpet/nema_wb.py +++ b/pynemaiqpet/nema_wb.py @@ -8,91 +8,101 @@ from scipy.ndimage.measurements import center_of_mass from scipy.ndimage.morphology import binary_fill_holes -from scipy.special import erf -from scipy.signal import argrelextrema, find_peaks_cwt - -from lmfit import Model - -def find_background_roi(vol, voxsize, Rcenter = 82, edge_margin = 28., showROI= False): - """ find the 2D background ROI for a NEMA sphere phantom - - Parameters - ---------- - - vol : 3d numpy array - containing the volume - - voxsize: 1d numpy array - containing the voxel size (mm) - - Rcenter: float (optional) - the radius (mm) of the sphere part that is not considered - default 82. - - edge_margin : float (optional) - margin (mm) to stay away from the boarder of the phantom - default 28. - - showROI : bool (optional) - whether to show 2D background ROI in a figure - default False - - Returns - ------- - - A tuple containing the indices of the background voxels. - """ - - # find the background value from a histogram analysis - h = np.histogram(vol[vol>0.01*vol.max()].flatten(),200) - bg = 0.5*(h[1][np.argmax(h[0])] + h[1][np.argmax(h[0]) + 1]) - - # get the axial slice with the maximum activity (spheres) - zprof = vol.sum((0,1)) - sphere_sl = np.argmax(zprof) - - sphere_2d_img = vol[...,sphere_sl] - - bg_mask = binary_fill_holes(median_filter(np.clip(sphere_2d_img,0,0.8*bg), size = 7) > 0.5*bg) - - # erode mask by ca. 2.8 cm to stay away from the boundary - nerode = int(edge_margin / voxsize[0]) - if (nerode % 2) == 0: nerode += 1 - - bg_mask = binary_erosion(bg_mask, np.ones((nerode,nerode))) - - # set center where spheres are to 0 - com = center_of_mass(binary_fill_holes(sphere_2d_img < 0.5*bg)) - - x = voxsize[0]*(np.arange(vol.shape[0]) - com[0]) - 5 - y = voxsize[1]*(np.arange(vol.shape[1]) - com[1]) - - X,Y = np.meshgrid(x,y) - R = np.sqrt(X**2 + Y**2) - - bg_mask[R<=Rcenter] = 0 - - # generate voxel indices for background voxels - tmp = np.zeros(vol.shape, dtype = np.int8) - tmp[...,sphere_sl] = bg_mask - bg_inds = np.where(tmp == 1) - - if showROI: - fig, ax = py.subplots(figsize = (5,5)) - ax.imshow(sphere_2d_img.T, cmap = py.cm.Greys, vmin = 0, vmax = np.percentile(sphere_2d_img,99.9)) - ax.contour(bg_mask.T) - ax.set_title(f'background ROI slice {sphere_sl}') - fig.tight_layout() - fig.show() +from scipy.special import erf +from scipy.signal import argrelextrema, find_peaks_cwt + +from lmfit import Model + + +def find_background_roi(vol, voxsize, Rcenter=82, edge_margin=28.0, showROI=False): + """find the 2D background ROI for a NEMA sphere phantom + + Parameters + ---------- + + vol : 3d numpy array + containing the volume + + voxsize: 1d numpy array + containing the voxel size (mm) + + Rcenter: float (optional) + the radius (mm) of the sphere part that is not considered - default 82. + + edge_margin : float (optional) + margin (mm) to stay away from the boarder of the phantom - default 28. + + showROI : bool (optional) + whether to show 2D background ROI in a figure - default False + + Returns + ------- + + A tuple containing the indices of the background voxels. + """ + + # find the background value from a histogram analysis + h = np.histogram(vol[vol > 0.01 * vol.max()].flatten(), 200) + bg = 0.5 * (h[1][np.argmax(h[0])] + h[1][np.argmax(h[0]) + 1]) + + # get the axial slice with the maximum activity (spheres) + zprof = vol.sum((0, 1)) + sphere_sl = np.argmax(zprof) + + sphere_2d_img = vol[..., sphere_sl] + + bg_mask = binary_fill_holes( + median_filter(np.clip(sphere_2d_img, 0, 0.8 * bg), size=7) > 0.5 * bg + ) + + # erode mask by ca. 2.8 cm to stay away from the boundary + nerode = int(edge_margin / voxsize[0]) + if (nerode % 2) == 0: + nerode += 1 + + bg_mask = binary_erosion(bg_mask, np.ones((nerode, nerode))) + + # set center where spheres are to 0 + com = center_of_mass(binary_fill_holes(sphere_2d_img < 0.5 * bg)) + + x = voxsize[0] * (np.arange(vol.shape[0]) - com[0]) - 5 + y = voxsize[1] * (np.arange(vol.shape[1]) - com[1]) + + X, Y = np.meshgrid(x, y) + R = np.sqrt(X**2 + Y**2) + + bg_mask[R <= Rcenter] = 0 + + # generate voxel indices for background voxels + tmp = np.zeros(vol.shape, dtype=np.int8) + tmp[..., sphere_sl] = bg_mask + bg_inds = np.where(tmp == 1) + + if showROI: + fig, ax = py.subplots(figsize=(5, 5)) + ax.imshow( + sphere_2d_img.T, + cmap=py.cm.Greys, + vmin=0, + vmax=np.percentile(sphere_2d_img, 99.9), + ) + ax.contour(bg_mask.T) + ax.set_title(f"background ROI slice {sphere_sl}") + fig.tight_layout() + fig.show() + + return bg_inds - return bg_inds -#-------------------------------------------------------------------------------------------------- +# -------------------------------------------------------------------------------------------------- -def gausssphere_profile(z = np.linspace(0, 2, 100), - Z = 0.8): - """ Radial profile of a sphere convolved with a 3D radial symmetric Gaussian + +def gausssphere_profile(z=np.linspace(0, 2, 100), Z=0.8): + """Radial profile of a sphere convolved with a 3D radial symmetric Gaussian Parameters ---------- - z : 1D numpy float array + z : 1D numpy float array normalized radial coordinate (r / (sqrt(2) * sigma)) Z : float @@ -104,32 +114,30 @@ def gausssphere_profile(z = np.linspace(0, 2, 100), """ sqrtpi = np.sqrt(np.pi) - + P = np.zeros_like(z) inds0 = np.argwhere(z == 0) inds1 = np.argwhere(z != 0) - - P[inds0] = erf(Z) - 2 * Z * np.exp(-Z**2) / sqrtpi - P[inds1] = ( 0.5 * (erf(z[inds1] + Z) - erf(z[inds1] - Z)) - - (0.5/sqrtpi) * ((np.exp(-(z[inds1] - Z)**2) - np.exp(-(z[inds1] + Z)**2)) / z[inds1])) - + P[inds0] = erf(Z) - 2 * Z * np.exp(-(Z**2)) / sqrtpi + + P[inds1] = 0.5 * (erf(z[inds1] + Z) - erf(z[inds1] - Z)) - (0.5 / sqrtpi) * ( + (np.exp(-((z[inds1] - Z) ** 2)) - np.exp(-((z[inds1] + Z) ** 2))) / z[inds1] + ) + return P -#-------------------------------------------------------------------------------------------------- -def glasssphere_profile(r, - R = 18.5, - FWHM = 5, - d = 1.5, - S = 10.0, - B = 1.0): - """ Radial profile of a hot sphere with cold glass wall in warm background +# -------------------------------------------------------------------------------------------------- + + +def glasssphere_profile(r, R=18.5, FWHM=5, d=1.5, S=10.0, B=1.0): + """Radial profile of a hot sphere with cold glass wall in warm background Parameters ---------- - r : 1D numpy float array + r : 1D numpy float array array with radial coordinates R : float, optional @@ -153,197 +161,202 @@ def glasssphere_profile(r, """ sqrt2 = np.sqrt(2) - sigma = FWHM / (2*np.sqrt(2*np.log(2))) - Z = R / (sigma*sqrt2) - w = d / (sigma*sqrt2) - z = r / (sigma*sqrt2) + sigma = FWHM / (2 * np.sqrt(2 * np.log(2))) + Z = R / (sigma * sqrt2) + w = d / (sigma * sqrt2) + z = r / (sigma * sqrt2) - P = S*gausssphere_profile(z, Z) - B*gausssphere_profile(z, Z + w) + B + P = S * gausssphere_profile(z, Z) - B * gausssphere_profile(z, Z + w) + B return P -#-------------------------------------------------------------------------------------------------- -def plot1dprofiles(vol, - voxsizes): - """ Plot profiles along the x, y and z axis through the center of mass of a sphere +# -------------------------------------------------------------------------------------------------- + + +def plot1dprofiles(vol, voxsizes): + """Plot profiles along the x, y and z axis through the center of mass of a sphere Parameters ---------- - vol : 3d numpy array + vol : 3d numpy array containing the volume - voxsizes : 3 component array + voxsizes : 3 component array with the voxel sizes - """ + """ # now we have to find the activity weighted center of gravity of the sphere # to do so we do a coarse delineation of the sphere (30% over bg) - bg = np.mean(vol[:,:,0]) - absth = relth*(vol.max() - bg) + bg + bg = np.mean(vol[:, :, 0]) + absth = relth * (vol.max() - bg) + bg - mask = np.zeros_like(vol, dtype = np.uint8) + mask = np.zeros_like(vol, dtype=np.uint8) mask[vol > absth] = 1 i0, i1, i2 = np.indices(vol.shape) - i0 = i0*voxsizes[0] - i1 = i1*voxsizes[1] - i2 = i2*voxsizes[2] + i0 = i0 * voxsizes[0] + i1 = i1 * voxsizes[1] + i2 = i2 * voxsizes[2] # calculate the maxmimum radius of the subvolumes # all voxels with a distance bigger than rmax will not be included in the fit - rmax = np.min((i0.max(),i1.max(),i2.max()))/2 + rmax = np.min((i0.max(), i1.max(), i2.max())) / 2 - # first try to get the center of mass via the coarse delineation - weights = vol[mask == 1] + # first try to get the center of mass via the coarse delineation + weights = vol[mask == 1] summedweights = np.sum(weights) - - c0 = np.sum(i0[mask == 1]*weights) / summedweights - c1 = np.sum(i1[mask == 1]*weights) / summedweights - c2 = np.sum(i2[mask == 1]*weights) / summedweights - r = np.sqrt((i0 - c0)**2 + (i1 - c1)**2 + (i2 - c2)**2) + c0 = np.sum(i0[mask == 1] * weights) / summedweights + c1 = np.sum(i1[mask == 1] * weights) / summedweights + c2 = np.sum(i2[mask == 1] * weights) / summedweights + + r = np.sqrt((i0 - c0) ** 2 + (i1 - c1) ** 2 + (i2 - c2) ** 2) # second try to get the center of mass - # use weights from a smoothed volume - sigmas = 4 / (2.355*voxsizes) - vol_sm = gaussian_filter(vol, sigma = sigmas) + # use weights from a smoothed volume + sigmas = 4 / (2.355 * voxsizes) + vol_sm = gaussian_filter(vol, sigma=sigmas) - weights = vol_sm[r <= rmax] + weights = vol_sm[r <= rmax] summedweights = np.sum(weights) - d0 = np.sum(i0[r <= rmax]*weights) / summedweights - d1 = np.sum(i1[r <= rmax]*weights) / summedweights - d2 = np.sum(i2[r <= rmax]*weights) / summedweights + d0 = np.sum(i0[r <= rmax] * weights) / summedweights + d1 = np.sum(i1[r <= rmax] * weights) / summedweights + d2 = np.sum(i2[r <= rmax] * weights) / summedweights - r = np.sqrt((i0 - d0)**2 + (i1 - d1)**2 + (i2 - d2)**2) + r = np.sqrt((i0 - d0) ** 2 + (i1 - d1) ** 2 + (i2 - d2) ** 2) if plot1dprofiles: - spherecenter = np.unravel_index(np.argmin(r),r.shape) + spherecenter = np.unravel_index(np.argmin(r), r.shape) prof0 = vol[:, spherecenter[1], spherecenter[2]] prof1 = vol[spherecenter[0], :, spherecenter[2]] prof2 = vol[spherecenter[0], spherecenter[1], :] dims = vol.shape - prof02 = vol.sum(axis = (1,2)) / (dims[1]*dims[2]) - prof12 = vol.sum(axis = (0,2)) / (dims[0]*dims[2]) - prof22 = vol.sum(axis = (0,1)) / (dims[0]*dims[1]) + prof02 = vol.sum(axis=(1, 2)) / (dims[1] * dims[2]) + prof12 = vol.sum(axis=(0, 2)) / (dims[0] * dims[2]) + prof22 = vol.sum(axis=(0, 1)) / (dims[0] * dims[1]) - c0 = sum(prof0*voxsizes[0]*np.arange(len(prof0))) / sum(prof0) - c1 = sum(prof1*voxsizes[1]*np.arange(len(prof1))) / sum(prof1) - c2 = sum(prof2*voxsizes[2]*np.arange(len(prof2))) / sum(prof2) + c0 = sum(prof0 * voxsizes[0] * np.arange(len(prof0))) / sum(prof0) + c1 = sum(prof1 * voxsizes[1] * np.arange(len(prof1))) / sum(prof1) + c2 = sum(prof2 * voxsizes[2] * np.arange(len(prof2))) / sum(prof2) - #bg = 0.5*(prof2[0:3].mean() + prof2[-3:].mean()) + # bg = 0.5*(prof2[0:3].mean() + prof2[-3:].mean()) - #fwhm0 = np.argwhere(prof0 - bg > 0.5*(prof0.max() - bg))[:,0].ptp()*voxsizes[0] - #fwhm1 = np.argwhere(prof1 - bg > 0.5*(prof1.max() - bg))[:,0].ptp()*voxsizes[1] - #fwhm2 = np.argwhere(prof2 - bg > 0.5*(prof2.max() - bg))[:,0].ptp()*voxsizes[2] + # fwhm0 = np.argwhere(prof0 - bg > 0.5*(prof0.max() - bg))[:,0].ptp()*voxsizes[0] + # fwhm1 = np.argwhere(prof1 - bg > 0.5*(prof1.max() - bg))[:,0].ptp()*voxsizes[1] + # fwhm2 = np.argwhere(prof2 - bg > 0.5*(prof2.max() - bg))[:,0].ptp()*voxsizes[2] fig1d, ax1d = py.subplots(1) - ax1d.plot(voxsizes[0]*np.arange(len(prof0)) - c0, prof0, label = 'x') - ax1d.plot(voxsizes[1]*np.arange(len(prof1)) - c1, prof1, label = 'y') - ax1d.plot(voxsizes[2]*np.arange(len(prof2)) - c2, prof2, label = 'z') + ax1d.plot(voxsizes[0] * np.arange(len(prof0)) - c0, prof0, label="x") + ax1d.plot(voxsizes[1] * np.arange(len(prof1)) - c1, prof1, label="y") + ax1d.plot(voxsizes[2] * np.arange(len(prof2)) - c2, prof2, label="z") py.legend() -#-------------------------------------------------------------------------------------------------- -def get_sphere_center(vol, - voxsizes, - relth = 0.25): - """ Get the center of gravity of a single hot sphere in a volume - +# -------------------------------------------------------------------------------------------------- + + +def get_sphere_center(vol, voxsizes, relth=0.25): + """Get the center of gravity of a single hot sphere in a volume + Parameters ---------- - vol : 3d numpy array + vol : 3d numpy array containing the volume - voxsizes : 3 component array + voxsizes : 3 component array with the voxel sizes relth : float, optional - the relative threshold (signal over background) for the first coarse + the relative threshold (signal over background) for the first coarse delination of the sphere - default 0.25 """ # now we have to find the activity weighted center of gravity of the sphere # to do so we do a coarse delineation of the sphere (30% over bg) - bg = np.mean(vol[:,:,0]) - absth = relth*(vol.max() - bg) + bg + bg = np.mean(vol[:, :, 0]) + absth = relth * (vol.max() - bg) + bg - mask = np.zeros_like(vol, dtype = np.uint8) + mask = np.zeros_like(vol, dtype=np.uint8) mask[vol > absth] = 1 i0, i1, i2 = np.indices(vol.shape) - i0 = i0*voxsizes[0] - i1 = i1*voxsizes[1] - i2 = i2*voxsizes[2] + i0 = i0 * voxsizes[0] + i1 = i1 * voxsizes[1] + i2 = i2 * voxsizes[2] # calculate the maxmimum radius of the subvolumes # all voxels with a distance bigger than rmax will not be included in the fit - rmax = np.min((i0.max(),i1.max(),i2.max()))/2 + rmax = np.min((i0.max(), i1.max(), i2.max())) / 2 - # first try to get the center of mass via the coarse delineation - weights = vol[mask == 1] + # first try to get the center of mass via the coarse delineation + weights = vol[mask == 1] summedweights = np.sum(weights) - - c0 = np.sum(i0[mask == 1]*weights) / summedweights - c1 = np.sum(i1[mask == 1]*weights) / summedweights - c2 = np.sum(i2[mask == 1]*weights) / summedweights - r = np.sqrt((i0 - c0)**2 + (i1 - c1)**2 + (i2 - c2)**2) + c0 = np.sum(i0[mask == 1] * weights) / summedweights + c1 = np.sum(i1[mask == 1] * weights) / summedweights + c2 = np.sum(i2[mask == 1] * weights) / summedweights + + r = np.sqrt((i0 - c0) ** 2 + (i1 - c1) ** 2 + (i2 - c2) ** 2) # second try to get the center of mass - # use weights from a smoothed volume - sigmas = 4 / (2.355*voxsizes) - vol_sm = gaussian_filter(vol, sigma = sigmas) + # use weights from a smoothed volume + sigmas = 4 / (2.355 * voxsizes) + vol_sm = gaussian_filter(vol, sigma=sigmas) - weights = vol_sm[r <= rmax] + weights = vol_sm[r <= rmax] summedweights = np.sum(weights) - d0 = np.sum(i0[r <= rmax]*weights) / summedweights - d1 = np.sum(i1[r <= rmax]*weights) / summedweights - d2 = np.sum(i2[r <= rmax]*weights) / summedweights + d0 = np.sum(i0[r <= rmax] * weights) / summedweights + d1 = np.sum(i1[r <= rmax] * weights) / summedweights + d2 = np.sum(i2[r <= rmax] * weights) / summedweights sphere_center = np.array([d0, d1, d2]) return sphere_center -#-------------------------------------------------------------------------------------------------- - -def fitspheresubvolume(vol, - voxsizes, - relth = 0.25, - Rfix = None, - FWHMfix = None, - dfix = None, - Sfix = None, - Bfix = None, - wm = 'dist', - cl = False, - sphere_center = None): + +# -------------------------------------------------------------------------------------------------- + + +def fitspheresubvolume( + vol, + voxsizes, + relth=0.25, + Rfix=None, + FWHMfix=None, + dfix=None, + Sfix=None, + Bfix=None, + wm="dist", + cl=False, + sphere_center=None, +): """Fit the radial sphere profile of a 3d volume containg 1 sphere Parameters ---------- - vol : 3d numpy array + vol : 3d numpy array containing the volume - voxsizes : 3 component array + voxsizes : 3 component array with the voxel sizes relth : float, optional - the relative threshold (signal over background) for the first coarse + the relative threshold (signal over background) for the first coarse delination of the sphere - dfix, Sfix, Bfix, Rfix : float, optional + dfix, Sfix, Bfix, Rfix : float, optional fixed values for the wall thickness, signal, background and radius - wm : string, optinal + wm : string, optinal the weighting method of the data (equal, dist, sqdist) cl : bool, optional bool whether to compute the confidence limits (this takes very long) - - sphere_center : 3 element np.array + + sphere_center : 3 element np.array containing the center of the spheres in mm this is the center of in voxel coordiantes multiplied by the voxel sizes @@ -351,69 +364,93 @@ def fitspheresubvolume(vol, ------- Dictionary with the fitresults (as returned by lmfit) - """ + """ + + if sphere_center is None: + sphere_center = get_sphere_center(vol, voxsizes, relth=relth) - if sphere_center is None: sphere_center = get_sphere_center(vol, voxsizes, relth = relth) - i0, i1, i2 = np.indices(vol.shape) - i0 = i0*voxsizes[0] - i1 = i1*voxsizes[1] - i2 = i2*voxsizes[2] + i0 = i0 * voxsizes[0] + i1 = i1 * voxsizes[1] + i2 = i2 * voxsizes[2] - rmax = np.min((i0.max(),i1.max(),i2.max()))/2 - r = np.sqrt((i0 - sphere_center[0])**2 + (i1 - sphere_center[1])**2 + (i2 - sphere_center[2])**2) + rmax = np.min((i0.max(), i1.max(), i2.max())) / 2 + r = np.sqrt( + (i0 - sphere_center[0]) ** 2 + + (i1 - sphere_center[1]) ** 2 + + (i2 - sphere_center[2]) ** 2 + ) data = vol[r <= rmax].flatten() rfit = r[r <= rmax].flatten() - if (Rfix == None): Rinit = 0.5*rmax - else: Rinit = Rfix - - if (FWHMfix == None): FWHMinit = 2*voxsizes[0] - else: FWHMinit = FWHMfix - - if (dfix == None): dinit = 0.15 - else: dinit = dfix - - if (Sfix == None): Sinit = data.max() - else: Sinit = Sfix - - if (Bfix == None): Binit = data.min() - else: Binit = Bfix + if Rfix == None: + Rinit = 0.5 * rmax + else: + Rinit = Rfix + + if FWHMfix == None: + FWHMinit = 2 * voxsizes[0] + else: + FWHMinit = FWHMfix + + if dfix == None: + dinit = 0.15 + else: + dinit = dfix + + if Sfix == None: + Sinit = data.max() + else: + Sinit = Sfix + + if Bfix == None: + Binit = data.min() + else: + Binit = Bfix # lets do the actual fit pmodel = Model(glasssphere_profile) - params = pmodel.make_params(R = Rinit, FWHM = FWHMinit, d = dinit, S = Sinit, B = Binit) + params = pmodel.make_params(R=Rinit, FWHM=FWHMinit, d=dinit, S=Sinit, B=Binit) # fix the parameters that should be fixed - if Rfix != None: params['R'].vary = False - if FWHMfix != None: params['FWHM'].vary = False - if dfix != None: params['d'].vary = False - if Sfix != None: params['S'].vary = False - if Bfix != None: params['B'].vary = False - - params['R'].min = 0 - params['FWHM'].min = 0 - params['d'].min = 0 - params['S'].min = 0 - params['B'].min = 0 - - if wm == 'equal' : weigths = np.ones_like(rfit) - elif wm == 'sqdist': weights = 1.0/(rfit**2) - else : weights = 1.0/rfit + if Rfix != None: + params["R"].vary = False + if FWHMfix != None: + params["FWHM"].vary = False + if dfix != None: + params["d"].vary = False + if Sfix != None: + params["S"].vary = False + if Bfix != None: + params["B"].vary = False + + params["R"].min = 0 + params["FWHM"].min = 0 + params["d"].min = 0 + params["S"].min = 0 + params["B"].min = 0 + + if wm == "equal": + weigths = np.ones_like(rfit) + elif wm == "sqdist": + weights = 1.0 / (rfit**2) + else: + weights = 1.0 / rfit weights[weights == np.inf] = 0 - fitres = pmodel.fit(data, r = rfit, params = params, weights = weights) + fitres = pmodel.fit(data, r=rfit, params=params, weights=weights) fitres.rdata = rfit - if cl: fitres.cls = fitres.conf_interval() + if cl: + fitres.cls = fitres.conf_interval() # calculate the a50 mean - fitres.a50th = fitres.values['B'] + 0.5*(vol.max() - fitres.values['B']) - fitres.mean_a50 = np.mean(data[data >= fitres.a50th]) + fitres.a50th = fitres.values["B"] + 0.5 * (vol.max() - fitres.values["B"]) + fitres.mean_a50 = np.mean(data[data >= fitres.a50th]) # calculate the mean - fitres.mean = np.mean(data[rfit <= fitres.values['R']]) + fitres.mean = np.mean(data[rfit <= fitres.values["R"]]) # calculate the max fitres.max = data.max() @@ -423,10 +460,12 @@ def fitspheresubvolume(vol, return fitres -#-------------------------------------------------------------------------------------------------- -def plotspherefit(fitres, ax = None, xlim = None, ylim = None, unit = 'mm', showres = True): - """Plot the results of a single sphere fit +# -------------------------------------------------------------------------------------------------- + + +def plotspherefit(fitres, ax=None, xlim=None, ylim=None, unit="mm", showres=True): + """Plot the results of a single sphere fit Parameters ---------- @@ -444,44 +483,68 @@ def plotspherefit(fitres, ax = None, xlim = None, ylim = None, unit = 'mm', show showres : bool, optional whether to add text about the fit results in the plot - """ - - rplot = np.linspace(0, fitres.rdata.max(),100) - - if ax == None: fig, ax = py.subplots(1) - if xlim == None: xlim = (0,rplot.max()) - - ax.plot(fitres.rdata, fitres.data, 'k.', ms = 2.5) - - ax.add_patch(patches.Rectangle((0, 0), fitres.values['R'], fitres.values['S'], - facecolor = 'lightgrey', edgecolor = 'None')) - x2 = fitres.values['R'] + fitres.values['d'] - dx2 = xlim[1] - x2 - ax.add_patch(patches.Rectangle((x2, 0), dx2, fitres.values['B'], - facecolor = 'lightgrey', edgecolor = 'None')) - - ax.plot(rplot, fitres.eval(r = rplot), 'r-') - ax.set_xlabel('R (' + unit + ')') - ax.set_ylabel('signal') + """ + + rplot = np.linspace(0, fitres.rdata.max(), 100) + + if ax == None: + fig, ax = py.subplots(1) + if xlim == None: + xlim = (0, rplot.max()) + + ax.plot(fitres.rdata, fitres.data, "k.", ms=2.5) + + ax.add_patch( + patches.Rectangle( + (0, 0), + fitres.values["R"], + fitres.values["S"], + facecolor="lightgrey", + edgecolor="None", + ) + ) + x2 = fitres.values["R"] + fitres.values["d"] + dx2 = xlim[1] - x2 + ax.add_patch( + patches.Rectangle( + (x2, 0), dx2, fitres.values["B"], facecolor="lightgrey", edgecolor="None" + ) + ) + + ax.plot(rplot, fitres.eval(r=rplot), "r-") + ax.set_xlabel("R (" + unit + ")") + ax.set_ylabel("signal") ax.set_xlim(xlim) - if ylim != None: ax.set_ylim(ylim) + if ylim != None: + ax.set_ylim(ylim) if showres: - ax.text(0.99, 0.99, fitres.fit_report(), fontsize = 6, transform = ax.transAxes, - verticalalignment='top', horizontalalignment = 'right') - -#-------------------------------------------------------------------------------------------------- - -def NEMASubvols(input_vol, - voxsizes, - relTh = 0.2, - minvol = 300, - margin = 9, - nbins = 100, - zignore = 38, - bgSignal = None): - """ Segment a complete NEMA PET volume with several hot sphere in different subvolumes containing + ax.text( + 0.99, + 0.99, + fitres.fit_report(), + fontsize=6, + transform=ax.transAxes, + verticalalignment="top", + horizontalalignment="right", + ) + + +# -------------------------------------------------------------------------------------------------- + + +def NEMASubvols( + input_vol, + voxsizes, + relTh=0.2, + minvol=300, + margin=9, + nbins=100, + zignore=38, + bgSignal=None, +): + """Segment a complete NEMA PET volume with several hot sphere in different subvolumes containing only one sphere Parameters @@ -489,7 +552,7 @@ def NEMASubvols(input_vol, input_vol : 3D numpy array the volume to be segmented - voxsizes a 1D numpy array + voxsizes a 1D numpy array containing the voxelsizes relTh : float, optional @@ -505,7 +568,7 @@ def NEMASubvols(input_vol, number of bins used in histogram for background determination zignore : float, optional - distance to edge of FOV that is ignored (same unit as voxelsize) + distance to edge of FOV that is ignored (same unit as voxelsize) bgSignal : float or None, optional the signal intensity of the background @@ -515,160 +578,179 @@ def NEMASubvols(input_vol, ------- list of slices to access the subvolumes from the original volume - """ + """ vol = input_vol.copy() xdim, ydim, zdim = vol.shape - + minvol = int(minvol / np.prod(voxsizes)) - + dx = int(np.ceil(margin / voxsizes[0])) dy = int(np.ceil(margin / voxsizes[1])) dz = int(np.ceil(margin / voxsizes[2])) nzignore = int(np.ceil(zignore / voxsizes[2])) - vol[:,:,:nzignore] = 0 - vol[:,:,-nzignore:] = 0 - + vol[:, :, :nzignore] = 0 + vol[:, :, -nzignore:] = 0 + # first do a quick search for the biggest sphere (noisy edge of FOV can spoil max value!) - histo = py.histogram(vol[vol > 0.01*vol.max()], nbins) - #bgSignal = histo[1][argrelextrema(histo[0], np.greater)[0][0]] + histo = py.histogram(vol[vol > 0.01 * vol.max()], nbins) + # bgSignal = histo[1][argrelextrema(histo[0], np.greater)[0][0]] if bgSignal is None: - bgSignal = histo[1][find_peaks_cwt(histo[0], np.arange(nbins/6,nbins))[0]] - thresh = bgSignal + relTh*(vol.max() - bgSignal) - - vol2 = np.zeros(vol.shape, dtype = int) + bgSignal = histo[1][find_peaks_cwt(histo[0], np.arange(nbins / 6, nbins))[0]] + thresh = bgSignal + relTh * (vol.max() - bgSignal) + + vol2 = np.zeros(vol.shape, dtype=int) vol2[vol > thresh] = 1 - + vol3, nrois = label(vol2) rois = np.arange(1, nrois + 1) roivols = labeled_comprehension(vol, vol3, rois, len, int, 0) - + i = 1 - - for roi in rois: - if(roivols[roi-1] < minvol): vol3[vol3 == roi] = 0 + + for roi in rois: + if roivols[roi - 1] < minvol: + vol3[vol3 == roi] = 0 else: vol3[vol3 == roi] = i i = i + 1 - - nspheres = vol3.max() + + nspheres = vol3.max() spherelabels = np.arange(1, nspheres + 1) - + bboxes = find_objects(vol3) - + nmaskvox = list() - slices = list() - + slices = list() + for bbox in bboxes: xstart = max(0, bbox[0].start - dx) - xstop = min(xdim, bbox[0].stop + dx + 1) - + xstop = min(xdim, bbox[0].stop + dx + 1) + ystart = max(0, bbox[1].start - dy) - ystop = min(xdim, bbox[1].stop + dy + 1) - + ystop = min(xdim, bbox[1].stop + dy + 1) + zstart = max(0, bbox[2].start - dz) - zstop = min(xdim, bbox[2].stop + dz + 1) - - slices.append((slice(xstart,xstop,None), slice(ystart,ystop,None), slice(zstart,zstop,None))) + zstop = min(xdim, bbox[2].stop + dz + 1) + + slices.append( + ( + slice(xstart, xstop, None), + slice(ystart, ystop, None), + slice(zstart, zstop, None), + ) + ) - nmaskvox.append((xstop-xstart)*(ystop-ystart)*(zstop-zstart)) + nmaskvox.append((xstop - xstart) * (ystop - ystart) * (zstop - zstart)) # sort subvols acc to number of voxel - slices = [ slices[kk] for kk in np.argsort(nmaskvox)[::-1] ] + slices = [slices[kk] for kk in np.argsort(nmaskvox)[::-1]] return slices -#-------------------------------------------------------------------------------------------------- - -def findNEMAROIs(vol, voxsizes, R = None, relth = 0.25, bgth = 0.5): - """ - image-based ROI definition in a NEMA IQ scan - - Arguments - --------- - - vol ... the volume to be segmented - - voxsizes ... a numpy array of voxelsizes in mm - - - Keyword arguments - ---------------- - - R ... (np.array or list) with the sphere radii in mm - default (None) means [18.5,14.,11.,8.5,6.5,5.] - - relth ... (float) relative threshold to find the spheres - above background (default 0.25) - - bgth ... (float) relative threshold to find homogenous - background around spheres - - Returns - ------- - - A volume containing labels for the ROIs: - 1 ... background - 2 - 7 ... NEMA spheres (large to small) - """ - - if R is None: R = np.array([18.5,14.,11.,8.5,6.5,5.]) - - slices = NEMASubvols(vol, voxsizes) - labelvol = np.zeros(vol.shape, dtype = np.uint8) - bgmask = np.zeros(vol.shape, dtype = np.uint8) - - for i in range(len(slices)): - subvol = vol[slices[i]] - sphere_center = get_sphere_center(subvol, voxsizes, relth = relth) - - i0, i1, i2 = np.indices(subvol.shape) - i0 = i0*voxsizes[0] - i1 = i1*voxsizes[1] - i2 = i2*voxsizes[2] - - r = np.sqrt((i0 - sphere_center[0])**2 + (i1 - sphere_center[1])**2 + (i2 - sphere_center[2])**2) - - if i == 0: bg = subvol[r > (R[i] + 7)].mean() - - submask = np.zeros(subvol.shape, dtype = np.uint8) - submask[r <= R[i]] = i + 2 - - labelvol[slices[i]] = submask - - # calculate the background - bgmask[vol > bgth*bg] = 1 - bgmask[labelvol >= 2] = 0 - - bgmask_eroded = binary_erosion(bgmask, iterations = int(15/voxsizes[0])) - labelvol[bgmask_eroded] = 1 - - return(labelvol) - -#-------------------------------------------------------------------------------------------------- - -def fit_WB_NEMA_sphere_profiles(vol, - voxsizes, - sm_fwhm = 0, - margin = 9.0, - dfix = 1.5, - Sfix = None, - Rfix = None, - FWHMfix = None, - wm = 'dist', - nmax_spheres = 6, - sameSignal = False, - showBGROI = False): - """ Analyse the sphere profiles of a NEMA scan + +# -------------------------------------------------------------------------------------------------- + + +def findNEMAROIs(vol, voxsizes, R=None, relth=0.25, bgth=0.5): + """ + image-based ROI definition in a NEMA IQ scan + + Arguments + --------- + + vol ... the volume to be segmented + + voxsizes ... a numpy array of voxelsizes in mm + + + Keyword arguments + ---------------- + + R ... (np.array or list) with the sphere radii in mm + default (None) means [18.5,14.,11.,8.5,6.5,5.] + + relth ... (float) relative threshold to find the spheres + above background (default 0.25) + + bgth ... (float) relative threshold to find homogenous + background around spheres + + Returns + ------- + + A volume containing labels for the ROIs: + 1 ... background + 2 - 7 ... NEMA spheres (large to small) + """ + + if R is None: + R = np.array([18.5, 14.0, 11.0, 8.5, 6.5, 5.0]) + + slices = NEMASubvols(vol, voxsizes) + labelvol = np.zeros(vol.shape, dtype=np.uint8) + bgmask = np.zeros(vol.shape, dtype=np.uint8) + + for i in range(len(slices)): + subvol = vol[slices[i]] + sphere_center = get_sphere_center(subvol, voxsizes, relth=relth) + + i0, i1, i2 = np.indices(subvol.shape) + i0 = i0 * voxsizes[0] + i1 = i1 * voxsizes[1] + i2 = i2 * voxsizes[2] + + r = np.sqrt( + (i0 - sphere_center[0]) ** 2 + + (i1 - sphere_center[1]) ** 2 + + (i2 - sphere_center[2]) ** 2 + ) + + if i == 0: + bg = subvol[r > (R[i] + 7)].mean() + + submask = np.zeros(subvol.shape, dtype=np.uint8) + submask[r <= R[i]] = i + 2 + + labelvol[slices[i]] = submask + + # calculate the background + bgmask[vol > bgth * bg] = 1 + bgmask[labelvol >= 2] = 0 + + bgmask_eroded = binary_erosion(bgmask, iterations=int(15 / voxsizes[0])) + labelvol[bgmask_eroded] = 1 + + return labelvol + + +# -------------------------------------------------------------------------------------------------- + + +def fit_WB_NEMA_sphere_profiles( + vol, + voxsizes, + sm_fwhm=0, + margin=9.0, + dfix=1.5, + Sfix=None, + Rfix=None, + FWHMfix=None, + wm="dist", + nmax_spheres=6, + sameSignal=False, + showBGROI=False, +): + """Analyse the sphere profiles of a NEMA scan Parameters ---------- vol : 3D numpy array the volume to be segmented - voxsizes : a 3 element numpy array + voxsizes : a 3 element numpy array of voxelsizes in mm sm_fwhm : float, optional @@ -693,143 +775,188 @@ def fit_WB_NEMA_sphere_profiles(vol, whether to forace all spheres to have the signal from the biggest sphere showBGROI : bool ,optional - whether to show the background ROI in a separate figure - default False + whether to show the background ROI in a separate figure - default False Returns ------- Dictionary containing the fitresults - """ + """ - unit = 'mm' + unit = "mm" if sm_fwhm > 0: - print('\nPost-smoothing with ' + str(sm_fwhm) + ' mm') - sigma = sm_fwhm / 2.355 - sigmas = sigma / voxsizes - vol = gaussian_filter(vol, sigma = sigmas) + print("\nPost-smoothing with " + str(sm_fwhm) + " mm") + sigma = sm_fwhm / 2.355 + sigmas = sigma / voxsizes + vol = gaussian_filter(vol, sigma=sigmas) # find the 2D background ROI - bg_inds = find_background_roi(vol, voxsizes, showROI = showBGROI) + bg_inds = find_background_roi(vol, voxsizes, showROI=showBGROI) bg_mean = vol[bg_inds].mean() - bg_cov = vol[bg_inds].std() / bg_mean + bg_cov = vol[bg_inds].std() / bg_mean - slices = NEMASubvols(vol, voxsizes, margin = margin, bgSignal = bg_mean) + slices = NEMASubvols(vol, voxsizes, margin=margin, bgSignal=bg_mean) subvols = list() - for iss, ss in enumerate(slices): - if iss < nmax_spheres: - subvols.append(vol[ss]) - - if Rfix == None: Rfix = [None] * len(subvols) - if len(Rfix) < len(subvols): Rfix = Rfix + [None] * (len(subvols) - len(Rfix)) + for iss, ss in enumerate(slices): + if iss < nmax_spheres: + subvols.append(vol[ss]) + + if Rfix == None: + Rfix = [None] * len(subvols) + if len(Rfix) < len(subvols): + Rfix = Rfix + [None] * (len(subvols) - len(Rfix)) # initial fit to get signal in the biggest sphere if (Sfix == None) and (sameSignal == True): - initfitres = fitspheresubvolume(subvols[0], voxsizes, dfix = dfix, Bfix = bg_mean, - FWHMfix = FWHMfix, Rfix = Rfix[0], wm = wm) - Sfix = initfitres.params['S'].value + initfitres = fitspheresubvolume( + subvols[0], + voxsizes, + dfix=dfix, + Bfix=bg_mean, + FWHMfix=FWHMfix, + Rfix=Rfix[0], + wm=wm, + ) + Sfix = initfitres.params["S"].value # fit of all spheres fitres = [] for i in range(len(subvols)): - fitres.append(fitspheresubvolume(subvols[i], voxsizes, dfix = dfix, Sfix = Sfix, - Bfix = bg_mean, FWHMfix = FWHMfix, Rfix = Rfix[i])) + fitres.append( + fitspheresubvolume( + subvols[i], + voxsizes, + dfix=dfix, + Sfix=Sfix, + Bfix=bg_mean, + FWHMfix=FWHMfix, + Rfix=Rfix[i], + ) + ) # summary of results - fwhms = np.array([x.values['FWHM'] for x in fitres]) - Rs = np.array([x.values['R'] for x in fitres]) - Bs = np.array([x.values['B'] for x in fitres]) - Ss = np.array([x.values['S'] for x in fitres]) - - sphere_mean_a50 = np.array([x.mean_a50 for x in fitres]) - sphere_mean = np.array([x.mean for x in fitres]) - sphere_max = np.array([x.max for x in fitres]) - - sphere_results = pd.DataFrame({'R':Rs,'FWHM':fwhms,'signal':Ss, 'mean_a50':sphere_mean_a50, - 'mean':sphere_mean, 'max':sphere_max, - 'background_mean': bg_mean, 'background_cov': bg_cov}, - index = np.arange(1,len(subvols)+1)) + fwhms = np.array([x.values["FWHM"] for x in fitres]) + Rs = np.array([x.values["R"] for x in fitres]) + Bs = np.array([x.values["B"] for x in fitres]) + Ss = np.array([x.values["S"] for x in fitres]) + + sphere_mean_a50 = np.array([x.mean_a50 for x in fitres]) + sphere_mean = np.array([x.mean for x in fitres]) + sphere_max = np.array([x.max for x in fitres]) + + sphere_results = pd.DataFrame( + { + "R": Rs, + "FWHM": fwhms, + "signal": Ss, + "mean_a50": sphere_mean_a50, + "mean": sphere_mean, + "max": sphere_max, + "background_mean": bg_mean, + "background_cov": bg_cov, + }, + index=np.arange(1, len(subvols) + 1), + ) return fitres, sphere_results -#------------------------------------------------------------------------------------------- -def show_WB_NEMA_profiles(fitres): - nsph = len(fitres) - - rmax = fitres[0].rdata.max() - fig, axes = py.subplots(2,3, figsize = (18,8.3)) - - ymax = 1.05*max([x.max for x in fitres]) - - for i in range(nsph): - plotspherefit(fitres[i], ylim = (0, ymax), - ax = axes[np.unravel_index(i,axes.shape)], xlim = (0,1.5*rmax)) - if nsph < 6: - for i in np.arange(nsph,6): - ax = axes[np.unravel_index(i,axes.shape)] - ax.set_axis_off() - - fig.tight_layout() - fig.show() - - return fig - -#------------------------------------------------------------------------------------------- -def show_WB_NEMA_recoveries(sphere_results, true_activity, - earlcolor = 'lightgreen', earlversion = 2): - - unit = 'mm' - - a50RCs = sphere_results['mean_a50'].values / true_activity - maxRCs = sphere_results['max'].values / true_activity - Rs = sphere_results['R'].values - - fig2, axes2 = py.subplots(1,2, figsize = (6,4.), sharex = True) - - # for the EARL limits see - # http://earl.eanm.org/cms/website.php?id=/en/projects/fdg_pet_ct_accreditation/accreditation_specifications.htm - - if earlversion == 1: - RCa50_min = np.array([0.76, 0.72, 0.63, 0.57, 0.44, 0.27]) - RCa50_max = np.array([0.89, 0.85, 0.78, 0.73, 0.60, 0.43]) - RCmax_min = np.array([0.95, 0.91, 0.83, 0.73, 0.59, 0.34]) - RCmax_max = np.array([1.16, 1.13, 1.09, 1.01, 0.85, 0.57]) - elif earlversion == 2: - RCa50_min = np.array([0.85, 0.82, 0.80, 0.76, 0.63, 0.39]) - RCa50_max = np.array([1.00, 0.97, 0.99, 0.97, 0.86, 0.61]) - RCmax_min = np.array([1.05, 1.01, 1.01, 1.00, 0.85, 0.52]) - RCmax_max = np.array([1.29, 1.26, 1.32, 1.38, 1.22, 0.88]) - - # add the EARL limits - if earlcolor != None: - for i in range(len(Rs)): - axes2[0].add_patch(patches.Rectangle((Rs[i] - 0.5, RCa50_min[i]), 1, - RCa50_max[i] - RCa50_min[i], - facecolor = earlcolor, edgecolor = 'None')) - axes2[0].set_title(f'EARL version {earlversion}') +# ------------------------------------------------------------------------------------------- +def show_WB_NEMA_profiles(fitres): + nsph = len(fitres) - axes2[0].plot(Rs, a50RCs, 'ko') - axes2[0].set_ylim(min(0.25,0.95*a50RCs.min()), max(1.02,1.05*a50RCs.max())) - axes2[0].set_xlabel('R (' + unit + ')') - axes2[0].set_ylabel('RC a50') + rmax = fitres[0].rdata.max() + fig, axes = py.subplots(2, 3, figsize=(18, 8.3)) - # add the EARL limits - if earlcolor != None: - for i in range(len(Rs)): - axes2[1].add_patch(patches.Rectangle((Rs[i] - 0.5, RCmax_min[i]), 1, - RCmax_max[i] - RCmax_min[i], - facecolor = earlcolor, edgecolor = 'None')) - axes2[1].set_title(f'EARL version {earlversion}') + ymax = 1.05 * max([x.max for x in fitres]) - axes2[1].plot(Rs, maxRCs, 'ko') - axes2[1].set_ylim(min(0.29,0.95*maxRCs.min()), max(1.4,1.05*maxRCs.max())) - axes2[1].set_xlabel('R (' + unit + ')') - axes2[1].set_ylabel('RC max') + for i in range(nsph): + plotspherefit( + fitres[i], + ylim=(0, ymax), + ax=axes[np.unravel_index(i, axes.shape)], + xlim=(0, 1.5 * rmax), + ) + if nsph < 6: + for i in np.arange(nsph, 6): + ax = axes[np.unravel_index(i, axes.shape)] + ax.set_axis_off() - fig2.tight_layout() - fig2.show() + fig.tight_layout() + fig.show() - return fig2 + return fig + + +# ------------------------------------------------------------------------------------------- +def show_WB_NEMA_recoveries( + sphere_results, true_activity, earlcolor="lightgreen", earlversion=2 +): + + unit = "mm" + + a50RCs = sphere_results["mean_a50"].values / true_activity + maxRCs = sphere_results["max"].values / true_activity + Rs = sphere_results["R"].values + + fig2, axes2 = py.subplots(1, 2, figsize=(6, 4.0), sharex=True) + + # for the EARL limits see + # http://earl.eanm.org/cms/website.php?id=/en/projects/fdg_pet_ct_accreditation/accreditation_specifications.htm + + if earlversion == 1: + RCa50_min = np.array([0.76, 0.72, 0.63, 0.57, 0.44, 0.27]) + RCa50_max = np.array([0.89, 0.85, 0.78, 0.73, 0.60, 0.43]) + RCmax_min = np.array([0.95, 0.91, 0.83, 0.73, 0.59, 0.34]) + RCmax_max = np.array([1.16, 1.13, 1.09, 1.01, 0.85, 0.57]) + elif earlversion == 2: + RCa50_min = np.array([0.85, 0.82, 0.80, 0.76, 0.63, 0.39]) + RCa50_max = np.array([1.00, 0.97, 0.99, 0.97, 0.86, 0.61]) + RCmax_min = np.array([1.05, 1.01, 1.01, 1.00, 0.85, 0.52]) + RCmax_max = np.array([1.29, 1.26, 1.32, 1.38, 1.22, 0.88]) + + # add the EARL limits + if earlcolor != None: + for i in range(len(Rs)): + axes2[0].add_patch( + patches.Rectangle( + (Rs[i] - 0.5, RCa50_min[i]), + 1, + RCa50_max[i] - RCa50_min[i], + facecolor=earlcolor, + edgecolor="None", + ) + ) + axes2[0].set_title(f"EARL version {earlversion}") + + axes2[0].plot(Rs, a50RCs, "ko") + axes2[0].set_ylim(min(0.25, 0.95 * a50RCs.min()), max(1.02, 1.05 * a50RCs.max())) + axes2[0].set_xlabel("R (" + unit + ")") + axes2[0].set_ylabel("RC a50") + + # add the EARL limits + if earlcolor != None: + for i in range(len(Rs)): + axes2[1].add_patch( + patches.Rectangle( + (Rs[i] - 0.5, RCmax_min[i]), + 1, + RCmax_max[i] - RCmax_min[i], + facecolor=earlcolor, + edgecolor="None", + ) + ) + axes2[1].set_title(f"EARL version {earlversion}") + + axes2[1].plot(Rs, maxRCs, "ko") + axes2[1].set_ylim(min(0.29, 0.95 * maxRCs.min()), max(1.4, 1.05 * maxRCs.max())) + axes2[1].set_xlabel("R (" + unit + ")") + axes2[1].set_ylabel("RC max") + + fig2.tight_layout() + fig2.show() + + return fig2 diff --git a/scripts/analyze_small_animal.py b/scripts/analyze_small_animal.py index 4c2e724..e36f6d7 100644 --- a/scripts/analyze_small_animal.py +++ b/scripts/analyze_small_animal.py @@ -5,126 +5,138 @@ import pynemaiqpet.nema_smallanimal as nsa import pandas as pd -from scipy.signal import argrelextrema +from scipy.signal import argrelextrema from scipy.ndimage import find_objects, center_of_mass import argparse from pathlib import Path -#----------------------------------------------------------------------------------------- + +# ----------------------------------------------------------------------------------------- # input parameter # master directory -mpath = Path.home() / 'Downloads/NEMA_Nifti_GS/' +mpath = Path.home() / "Downloads/NEMA_Nifti_GS/" # recursive search for nifti files to analyze -fnames = sorted(list(mpath.rglob('*Cropped*/*.nii'))) +fnames = sorted(list(mpath.rglob("*Cropped*/*.nii"))) # output directory for results -output_dir = Path('results') +output_dir = Path("results") # verbose output -verbose = True +verbose = True # method for ROI to calculate recovery for hot rods (max or mean) -lp_voxel = 'mean' -#----------------------------------------------------------------------------------------- +lp_voxel = "mean" +# ----------------------------------------------------------------------------------------- if not output_dir.exists(): - output_dir.mkdir(parents = True, exist_ok = True) + output_dir.mkdir(parents=True, exist_ok=True) df = pd.DataFrame() for nifti_file in fnames: - if verbose: - print(nifti_file) - - rel_path = nifti_file.relative_to(mpath) - - # read the PET volume from dicom - nii = nib.load(nifti_file) - nii = nib.as_closest_canonical(nii) - vol = np.flip(nii.get_fdata().squeeze(),(0,1)) - voxsize = nii.header['pixdim'][1:4] - - flipz = False - - # flip volume in z direction if necessary - if center_of_mass(vol)[2] < 0.5*vol.shape[2]: - vol = np.flip(vol,2) - flipz = True - - phantom_name = nsa.get_phantom_name(vol, voxsize) - if verbose: - print(phantom_name) - - # align the PET volume to "standard" space (a digitial version of the phantom) - vol_aligned = nsa.align_nema_2008_small_animal_iq_phantom(vol, voxsize, version = phantom_name) - - # generate the ROI label volume - roi_vol = nsa.nema_2008_small_animal_pet_rois(vol_aligned, voxsize, phantom = phantom_name, lp_voxel = lp_voxel) - - # generate the report - res = nsa.nema_2008_small_animal_iq_phantom_report(vol_aligned, roi_vol) - - # reformat the results such that we can append them to a single data frame - res_reformatted = {} - - for key, val in res.items(): - if isinstance(val,np.ndarray): - for ia, aval in enumerate(val): - res_reformatted[f'{key}_{ia+1}'] = aval - else: - res_reformatted[key] = val - - df = df.append(pd.DataFrame(res_reformatted, index = [rel_path])) - - #--------------------------------------------------------------------------------------------------- - # plots - - odir = output_dir / rel_path.parent - - odir.mkdir(parents=True, exist_ok=True) - - if flipz: - vol_aligned = np.flip(vol_aligned,2) - roi_vol = np.flip(roi_vol,2) - - # show the aligned volume and the ROI volume - vi = pv.ThreeAxisViewer([vol_aligned,vol_aligned], [None,roi_vol**0.1], voxsize = voxsize, ls = '') - fig_file = odir / (rel_path.stem + '_uniform_rois.png') - vi.fig.savefig(fig_file) - if verbose: - print(f'wrote {fig_file}') - plt.close(vi.fig) - - # show the summed rod planes - bbox_uni = find_objects(vol_aligned >= 0.5*np.percentile(vol_aligned,99)) - bbox_rods = find_objects(roi_vol == 4) - - rod_img = vol_aligned[bbox_uni[0][0],bbox_uni[0][1],bbox_rods[0][2]].mean(2).T - roi_img = roi_vol[bbox_uni[0][0],bbox_uni[0][1],bbox_rods[0][2]].mean(2).T - - fig, ax = plt.subplots(1,2, figsize = (10,5)) - ax[0].imshow(rod_img, cmap = plt.cm.Greys, interpolation='bilinear', vmax = res['uniform_ROI_mean']) - ax[0].contour(roi_img > 0, levels = 1, cmap = 'spring') - ax[1].imshow(rod_img**0.1, cmap = plt.cm.Greys, interpolation='bilinear') - ax[1].contour(roi_img > 0, levels = 1, cmap = 'spring') - - ax[0].set_title('rod image') - ax[1].set_title('rod image^0.1 (compressed contrast)') - - for axx in ax.ravel(): - axx.set_axis_off() - - fig.tight_layout() - fig.show() - fig_file = odir / (rel_path.stem + '_rod_rois.png') - fig.savefig(fig_file) - if verbose: - print(f'wrote {fig_file}') - plt.close(fig) - -#------------------------------------- + if verbose: + print(nifti_file) + + rel_path = nifti_file.relative_to(mpath) + + # read the PET volume from dicom + nii = nib.load(nifti_file) + nii = nib.as_closest_canonical(nii) + vol = np.flip(nii.get_fdata().squeeze(), (0, 1)) + voxsize = nii.header["pixdim"][1:4] + + flipz = False + + # flip volume in z direction if necessary + if center_of_mass(vol)[2] < 0.5 * vol.shape[2]: + vol = np.flip(vol, 2) + flipz = True + + phantom_name = nsa.get_phantom_name(vol, voxsize) + if verbose: + print(phantom_name) + + # align the PET volume to "standard" space (a digitial version of the phantom) + vol_aligned = nsa.align_nema_2008_small_animal_iq_phantom( + vol, voxsize, version=phantom_name + ) + + # generate the ROI label volume + roi_vol = nsa.nema_2008_small_animal_pet_rois( + vol_aligned, voxsize, phantom=phantom_name, lp_voxel=lp_voxel + ) + + # generate the report + res = nsa.nema_2008_small_animal_iq_phantom_report(vol_aligned, roi_vol) + + # reformat the results such that we can append them to a single data frame + res_reformatted = {} + + for key, val in res.items(): + if isinstance(val, np.ndarray): + for ia, aval in enumerate(val): + res_reformatted[f"{key}_{ia+1}"] = aval + else: + res_reformatted[key] = val + + df = df.append(pd.DataFrame(res_reformatted, index=[rel_path])) + + # --------------------------------------------------------------------------------------------------- + # plots + + odir = output_dir / rel_path.parent + + odir.mkdir(parents=True, exist_ok=True) + + if flipz: + vol_aligned = np.flip(vol_aligned, 2) + roi_vol = np.flip(roi_vol, 2) + + # show the aligned volume and the ROI volume + vi = pv.ThreeAxisViewer( + [vol_aligned, vol_aligned], [None, roi_vol**0.1], voxsize=voxsize, ls="" + ) + fig_file = odir / (rel_path.stem + "_uniform_rois.png") + vi.fig.savefig(fig_file) + if verbose: + print(f"wrote {fig_file}") + plt.close(vi.fig) + + # show the summed rod planes + bbox_uni = find_objects(vol_aligned >= 0.5 * np.percentile(vol_aligned, 99)) + bbox_rods = find_objects(roi_vol == 4) + + rod_img = vol_aligned[bbox_uni[0][0], bbox_uni[0][1], bbox_rods[0][2]].mean(2).T + roi_img = roi_vol[bbox_uni[0][0], bbox_uni[0][1], bbox_rods[0][2]].mean(2).T + + fig, ax = plt.subplots(1, 2, figsize=(10, 5)) + ax[0].imshow( + rod_img, + cmap=plt.cm.Greys, + interpolation="bilinear", + vmax=res["uniform_ROI_mean"], + ) + ax[0].contour(roi_img > 0, levels=1, cmap="spring") + ax[1].imshow(rod_img**0.1, cmap=plt.cm.Greys, interpolation="bilinear") + ax[1].contour(roi_img > 0, levels=1, cmap="spring") + + ax[0].set_title("rod image") + ax[1].set_title("rod image^0.1 (compressed contrast)") + + for axx in ax.ravel(): + axx.set_axis_off() + + fig.tight_layout() + fig.show() + fig_file = odir / (rel_path.stem + "_rod_rois.png") + fig.savefig(fig_file) + if verbose: + print(f"wrote {fig_file}") + plt.close(fig) + +# ------------------------------------- # write all results to a single csv file -csv_file = output_dir / 'results.csv' +csv_file = output_dir / "results.csv" df.to_csv(csv_file) if verbose: - print(f'wrote {csv_file}') + print(f"wrote {csv_file}") diff --git a/scripts/dmi_earl.py b/scripts/dmi_earl.py index 6bc114d..2e36c7e 100644 --- a/scripts/dmi_earl.py +++ b/scripts/dmi_earl.py @@ -17,93 +17,137 @@ from argparse import ArgumentParser -parser = ArgumentParser(description = 'Analyze GE PET NEMA in terms of recovery') -parser.add_argument('dcm_dir_pattern', help = 'dicom direcotry pattern') -parser.add_argument('--dcm_file_pattern', default = '*.dcm', help = 'dicom file pattern') -parser.add_argument('--injected_signal', default = None, type = float, - help = 'injected signal in the spheres. if not given, taken from biggest sph.') -parser.add_argument('--sm_fwhm_mm', default = 5, type = float, - help = 'FWHM of Gaussian kernel (mm) used for additional post-smoothing') -parser.add_argument('--force_smoothing', action = 'store_true', - help = 'force smoothing. even for already smoothed data') -parser.add_argument('--earl_version', default = 2, type = int, choices = [1,2], - help = 'EARL version to use for limits in recovery plots') +parser = ArgumentParser(description="Analyze GE PET NEMA in terms of recovery") +parser.add_argument("dcm_dir_pattern", help="dicom direcotry pattern") +parser.add_argument("--dcm_file_pattern", default="*.dcm", help="dicom file pattern") +parser.add_argument( + "--injected_signal", + default=None, + type=float, + help="injected signal in the spheres. if not given, taken from biggest sph.", +) +parser.add_argument( + "--sm_fwhm_mm", + default=5, + type=float, + help="FWHM of Gaussian kernel (mm) used for additional post-smoothing", +) +parser.add_argument( + "--force_smoothing", + action="store_true", + help="force smoothing. even for already smoothed data", +) +parser.add_argument( + "--earl_version", + default=2, + type=int, + choices=[1, 2], + help="EARL version to use for limits in recovery plots", +) args = parser.parse_args() # injected sphere activity concentration at start of acq. injected_signal = args.injected_signal -sm_fwhm_mm = args.sm_fwhm_mm +sm_fwhm_mm = args.sm_fwhm_mm force_smoothing = args.force_smoothing -earlversion = args.earl_version +earlversion = args.earl_version -dcm_dir_pattern = args.dcm_dir_pattern +dcm_dir_pattern = args.dcm_dir_pattern dcm_file_pattern = args.dcm_file_pattern -#------------------------------------------------------------------------------------------------ +# ------------------------------------------------------------------------------------------------ dcm_dirs = glob(dcm_dir_pattern) for dcm_dir in dcm_dirs: - print(os.path.basename(dcm_dir)) - - # load example data set included in package - dcm_pattern = os.path.join(dcm_dir, dcm_file_pattern) - dcm = pmf.DicomVolume(dcm_pattern) - vol = dcm.get_data() - voxsize = dcm.voxsize - - # FWHM of Gaussian kernel to apply before analysis - dcm_hdr = dcm.firstdcmheader - - # check if the data was already post-smoothed by analysing private GE tags - # for trans-axial and axial filter - if ((dcm_hdr[0x0009,0x10ba].value == 0) and (dcm_hdr[0x0009,0x10db].value == 0)) or force_smoothing: - sm_str = f'_{sm_fwhm_mm}mm_ps' - else: - sm_fwhm_mm = 0 - sm_str = '' - - # post smooth the image - if sm_fwhm_mm > 0: - print(f'Applying additional isotropic Gaussian post-smoothing of {sm_fwhm_mm}mm') - vol = gaussian_filter(vol, sigma = sm_fwhm_mm / (2.35*voxsize)) - - #------------------------------------------------------------------------------------------------- - # do a fit where we force all spheres to have the same signal (assuming that the activity - # concentration in all sphere was the same - # and we fix the radii of all spheres (using the values from the NEMA specs) - - # try also doing the fit without fixing the radii - fitres, sphere_results = nema.fit_WB_NEMA_sphere_profiles(vol, voxsize, sameSignal = True, - Rfix = [18.5, 14.0, 11.0, 8.5, 6.5, 5.], - showBGROI = False, - Sfix = injected_signal) - - print('fit with same signal and fixed radii') - print(sphere_results) - sphere_results.to_csv(os.path.join(dcm_dir,f'{dcm.firstdcmheader.SeriesDescription}{sm_str}.csv'.replace(' ','_'))) - - #------------------------------------------------------------------------------------------------- - # show the results - - fig = nema.show_WB_NEMA_profiles(fitres) - fig.savefig(os.path.join(dcm_dir,f'{dcm.firstdcmheader.SeriesDescription}{sm_str}_profiles_EARL_{earlversion}.png'.replace(' ','_'))) - - # plot the max and a50 recoveries - # the 2nd argument should be the true (expected) activity concentration in the spheres - # and also show the limits given by EARL (vesion 2) - - if injected_signal is None: - ref_signal = sphere_results.signal.values[0] - else: - ref_signal = injected_signal - - fig2 = nema.show_WB_NEMA_recoveries(sphere_results, ref_signal, - earlversion = earlversion) - fig2.suptitle(f'{dcm.firstdcmheader.SeriesDescription}{sm_str}') - fig2.tight_layout(pad = 2) - - fig2.savefig(os.path.join(dcm_dir,f'{dcm.firstdcmheader.SeriesDescription}{sm_str}_RCs_EARL_{earlversion}.png'.replace(' ','_'))) + print(os.path.basename(dcm_dir)) + + # load example data set included in package + dcm_pattern = os.path.join(dcm_dir, dcm_file_pattern) + dcm = pmf.DicomVolume(dcm_pattern) + vol = dcm.get_data() + voxsize = dcm.voxsize + + # FWHM of Gaussian kernel to apply before analysis + dcm_hdr = dcm.firstdcmheader + + # check if the data was already post-smoothed by analysing private GE tags + # for trans-axial and axial filter + if ( + (dcm_hdr[0x0009, 0x10BA].value == 0) and (dcm_hdr[0x0009, 0x10DB].value == 0) + ) or force_smoothing: + sm_str = f"_{sm_fwhm_mm}mm_ps" + else: + sm_fwhm_mm = 0 + sm_str = "" + + # post smooth the image + if sm_fwhm_mm > 0: + print( + f"Applying additional isotropic Gaussian post-smoothing of {sm_fwhm_mm}mm" + ) + vol = gaussian_filter(vol, sigma=sm_fwhm_mm / (2.35 * voxsize)) + + # ------------------------------------------------------------------------------------------------- + # do a fit where we force all spheres to have the same signal (assuming that the activity + # concentration in all sphere was the same + # and we fix the radii of all spheres (using the values from the NEMA specs) + + # try also doing the fit without fixing the radii + fitres, sphere_results = nema.fit_WB_NEMA_sphere_profiles( + vol, + voxsize, + sameSignal=True, + Rfix=[18.5, 14.0, 11.0, 8.5, 6.5, 5.0], + showBGROI=False, + Sfix=injected_signal, + ) + + print("fit with same signal and fixed radii") + print(sphere_results) + sphere_results.to_csv( + os.path.join( + dcm_dir, + f"{dcm.firstdcmheader.SeriesDescription}{sm_str}.csv".replace(" ", "_"), + ) + ) + + # ------------------------------------------------------------------------------------------------- + # show the results + + fig = nema.show_WB_NEMA_profiles(fitres) + fig.savefig( + os.path.join( + dcm_dir, + f"{dcm.firstdcmheader.SeriesDescription}{sm_str}_profiles_EARL_{earlversion}.png".replace( + " ", "_" + ), + ) + ) + + # plot the max and a50 recoveries + # the 2nd argument should be the true (expected) activity concentration in the spheres + # and also show the limits given by EARL (vesion 2) + + if injected_signal is None: + ref_signal = sphere_results.signal.values[0] + else: + ref_signal = injected_signal + + fig2 = nema.show_WB_NEMA_recoveries( + sphere_results, ref_signal, earlversion=earlversion + ) + fig2.suptitle(f"{dcm.firstdcmheader.SeriesDescription}{sm_str}") + fig2.tight_layout(pad=2) + + fig2.savefig( + os.path.join( + dcm_dir, + f"{dcm.firstdcmheader.SeriesDescription}{sm_str}_RCs_EARL_{earlversion}.png".replace( + " ", "_" + ), + ) + ) diff --git a/scripts/small_animal.py b/scripts/small_animal.py index 237d2e9..a8a5ef5 100644 --- a/scripts/small_animal.py +++ b/scripts/small_animal.py @@ -7,10 +7,14 @@ import argparse -parser = argparse.ArgumentParser(description='NEMA small animal IQ scan analyzer') -parser.add_argument('dcm_dir', help = 'absolute path of input dicom directory') -parser.add_argument('--phantom', help = 'phantom version', choices = ['standard','mini'], - default = 'standard') +parser = argparse.ArgumentParser(description="NEMA small animal IQ scan analyzer") +parser.add_argument("dcm_dir", help="absolute path of input dicom directory") +parser.add_argument( + "--phantom", + help="phantom version", + choices=["standard", "mini"], + default="standard", +) args = parser.parse_args() # read the PET volume from dicom @@ -18,16 +22,24 @@ vol = dcm.get_data() # align the PET volume to "standard" space (a digitial version of the phantom) -vol_aligned = nsa.align_nema_2008_small_animal_iq_phantom(vol, dcm.voxsize, version = args.phantom) +vol_aligned = nsa.align_nema_2008_small_animal_iq_phantom( + vol, dcm.voxsize, version=args.phantom +) # generate the ROI label volume -roi_vol = nsa.nema_2008_small_animal_pet_rois(vol_aligned, dcm.voxsize, phantom = args.phantom) +roi_vol = nsa.nema_2008_small_animal_pet_rois( + vol_aligned, dcm.voxsize, phantom=args.phantom +) # generate the report nsa.nema_2008_small_animal_iq_phantom_report(vol_aligned, roi_vol) # show the aligned volume and the ROI volume (cropped versions) -th = 0.3*np.percentile(vol_aligned, 99.9) +th = 0.3 * np.percentile(vol_aligned, 99.9) bbox = find_objects(vol_aligned > th)[0] -vi = pv.ThreeAxisViewer([vol_aligned[bbox],vol_aligned[bbox]], [None,roi_vol[bbox]**0.1], - voxsize = dcm.voxsize, ls = '') +vi = pv.ThreeAxisViewer( + [vol_aligned[bbox], vol_aligned[bbox]], + [None, roi_vol[bbox] ** 0.1], + voxsize=dcm.voxsize, + ls="", +)