diff --git a/amypet/align_brkdyn_ct.py b/amypet/align_brkdyn_ct.py index ddda0199..0be7bde2 100644 --- a/amypet/align_brkdyn_ct.py +++ b/amypet/align_brkdyn_ct.py @@ -39,6 +39,7 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F used as reference instead. ''' + # > output path from the input dictionary opth = niidat['outpath'].parent @@ -55,11 +56,16 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F outdct = np.load(fout, allow_pickle=True) outdct = outdct.item() return outdct + else: + # > output dictionary + outdct = {} # > number of frames for break dynamic acquisitions A and B nfrmA = len(niidat['series'][0]) nfrmB = len(niidat['series'][1]) + outdct['nfrm_t'] = (nfrmA, nfrmB) + # > sort out the correct index for PET reference frame if refpetidx is not None: if any([refpetidx[0]<-nfrmA, refpetidx[1]<-nfrmB, refpetidx[0]>=nfrmA, refpetidx[1]>=nfrmB]): @@ -72,6 +78,38 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F refpetidx[1] += nfrmB + # > DECAY CORRECTION (if requested) + if Cnt['align']['decay_corr']: + + # > get the start time of each series for decay correction if requested + ts = [sri['time'][0] for sri in niidat['descr']] + # > index the earliest ref time + i_tref = np.argmin(ts) + idxs = list(range(len(ts))) + idxs.pop(i_tref) + i_tsrs = idxs + # > time difference + td = [ts[i] - ts[i_tref] for i in i_tsrs] + if len(td) > 1: + raise ValueError( + 'currently only one dynamic break is allowed - detected more than one') + else: + td = td[0] + + # > what tracer / radionuclide is used? + istp = 'F18' * (niidat['tracer'] in Cnt['tracer']['f18']) + 'C11' * (niidat['tracer'] in Cnt['tracer']['c11']) + + # > decay constant using half-life + lmbd = np.log(2) / nimpa.resources.riLUT[istp]['thalf'] + + # > decay correction factor + dcycrr = np.exp(lmbd * td) + else: + dcycrr = 1. + + outdct['dcycrr'] = dcycrr + + # > reference images based on CT fctref = [None, None] fpetref = [None, None] @@ -235,9 +273,22 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F for fi, f in enumerate(algn_frm[1]['faligned']): print(f) - fcopy = algnFpth/Path(f).name - shutil.copyfile(f, fcopy) - faligned.append(fcopy) + + # > correction for decay + if Cnt['align']['decay_corr']: + fdcrr = algnFpth/(Path(f).name.split('.nii')[0]+'_dcrr.nii.gz') + imd = nimpa.getnii(f, output='all') + nimpa.array2nii( + imd['im']*dcycrr, + imd['affine'], + fdcrr, + trnsp=imd['transpose'], + flip=imd['flip']) + faligned.append(fdcrr) + else: + fcopy = algnFpth/Path(f).name + shutil.copyfile(f, fcopy) + faligned.append(fcopy) faff = affsF/(Path(f).name.split('.nii')[0]+'_affine.txt') shutil.copyfile(algn_frm[1]['affines'][fi], faff) @@ -287,11 +338,23 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F for fi, f in enumerate(algn_frm[i_acq]['faligned']): print(f) - #fnii_org.append(niidat['series'][i_acq][keys2[fi]]['fnii']) - fcopy = algnFpth/Path(f).name - shutil.copyfile(f, fcopy) - faligned.append(fcopy) + # > correction for decay + if Cnt['align']['decay_corr']: + fdcrr = algnFpth/(Path(f).name.split('.nii')[0]+'_dcrr.nii.gz') + imd = nimpa.getnii(f, output='all') + nimpa.array2nii( + imd['im']*dcycrr, + imd['affine'], + fdcrr, + trnsp=imd['transpose'], + flip=imd['flip']) + faligned.append(fdcrr) + else: + fcopy = algnFpth/Path(f).name + shutil.copyfile(f, fcopy) + faligned.append(fcopy) + faff = affsF/(Path(f).name.split('.nii')[0]+'_affine.txt') shutil.copyfile(algn_frm[i_acq]['affines'][fi], faff) @@ -335,9 +398,6 @@ def align_break_petct(niidat, cts, Cnt, qcpth=None, refpetidx=None, use_stored=F #-------------------------------------------- - - # > output dictionary - outdct = {} outdct['align_acq'] = algn_frm outdct['ctref'] = fctref outdct['petref'] = fpetref diff --git a/amypet/params.toml b/amypet/params.toml index 04644f3d..2743fde2 100644 --- a/amypet/params.toml +++ b/amypet/params.toml @@ -5,6 +5,7 @@ pib = ["pib"] flute = ["flt", "flut", "flute", "flutemetamol"] fbb = ["fbb", "florbetaben"] fbp = ["fbp", "florbetapir"] +mk6240 = ["mk6240", "mk-6240"] [regpars] fwhm_t1_mni = 3 @@ -55,5 +56,5 @@ fbb = [5400, 6600, 1200] fbp = [3000, 3600, 600] [tracer] -f18 = ['fbb', 'fbp', 'flute', 'mk-6240'] +f18 = ['fbb', 'fbp', 'flute', 'mk6240'] c11 = ['pib'] \ No newline at end of file