From 07b4889171f8efaac364dc3ed1828a5780b53a28 Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Thu, 7 Nov 2024 11:14:57 -0500 Subject: [PATCH 1/8] initial commit for snr branch --- scanbuddy/proc/__init__.py | 273 +++++++++++++++++++++++++++++++++++++ scanbuddy/proc/snr.py | 253 ++++++++++++++++++++++++++++++++++ scanbuddy/view/dash.py | 47 ++++++- scripts/start.py | 4 + setup.py | 3 +- 5 files changed, 578 insertions(+), 2 deletions(-) create mode 100644 scanbuddy/proc/snr.py diff --git a/scanbuddy/proc/__init__.py b/scanbuddy/proc/__init__.py index 487f88c..15dcc78 100644 --- a/scanbuddy/proc/__init__.py +++ b/scanbuddy/proc/__init__.py @@ -1,3 +1,5 @@ +import sys +import pdb import time import math import json @@ -16,6 +18,8 @@ def __init__(self): def reset(self): self._instances = SortedDict() + self._slice_means = SortedDict() + pub.sendMessage('plot_snr', snr_metric=str(0.0)) logger.debug('received message to reset') def listener(self, ds, path): @@ -24,16 +28,48 @@ def listener(self, ds, path): 'path': path, 'volreg': None } + self._slice_means[key] = { + 'path': path, + 'slice_means': None, + 'mask_threshold': None + } logger.debug('current state of instances') logger.debug(json.dumps(self._instances, default=list, indent=2)) tasks = self.check_volreg(key) + ### edits need to be made here + snr_tasks = self.check_snr(key) logger.debug('publishing message to volreg topic with the following tasks') logger.debug(json.dumps(tasks, indent=2)) pub.sendMessage('volreg', tasks=tasks) logger.debug(f'publishing message to params topic') pub.sendMessage('params', ds=ds) + logger.debug(f'publishing message to snr topic') + logger.debug(f'snr task sorted dict: {snr_tasks}') + pub.sendMessage('snr', tasks=snr_tasks) + logger.debug('after snr calculation') + + logger.debug(json.dumps(self._instances, indent=2)) + if key == 5: + self._fdata_array = self._slice_means[key]['slice_means'] + #self._slice_means_array = self._slice_means[key]['slice_means'] + #self._mask_array = self._slice_means[key]['mask'] + elif key > 5: + self._fdata_array = np.concatenate((self._fdata_array, self._slice_means[key]['slice_means']), axis=3) + #self._slice_means_array = np.column_stack((self._slice_means_array, self._slice_means[key]['slice_means'])) + #self._mask_array = np.column_stack((self._mask_array, self._slice_means[key]['mask'])) + #if len(self._instances) % self._snr_interval == 0: + if key == 10: + snr_metric = self.calc_snr() + if key > 53: + #logger.info(f'shape of slice_means_array: {self._slice_means_array.shape}') + logger.info(f'shape of fdata_array: {self._fdata_array.shape}') + logger.info('publishing message to plot_snr topic') + snr_metric = round(self.calc_snr(), 2) + logger.info(f'running snr metric: {snr_metric}') + pub.sendMessage('plot_snr', snr_metric=snr_metric) + logger.debug(f'after volreg') logger.debug(json.dumps(self._instances, indent=2)) project = ds.get('StudyDescription', '[STUDY]') @@ -70,3 +106,240 @@ def check_volreg(self, key): return tasks + def calc_snr(self): + + slice_intensity_means, slice_voxel_counts, data = self.get_mean_slice_intensitites() + + slice_count = slice_intensity_means.shape[0] + volume_count = slice_intensity_means.shape[1] + + slice_weighted_mean_mean = 0 + slice_weighted_stdev_mean = 0 + slice_weighted_snr_mean = 0 + slice_weighted_max_mean = 0 + slice_weighted_min_mean = 0 + outlier_count = 0 + total_voxel_count = 0 + + for slice_idx in range(slice_count): + slice_data = slice_intensity_means[slice_idx] + slice_voxel_count = slice_voxel_counts[slice_idx] + slice_mean = slice_data.mean() + slice_stdev = slice_data.std(ddof=1) + slice_snr = slice_mean / slice_stdev + + slice_weighted_mean_mean += (slice_mean * slice_voxel_count) + slice_weighted_stdev_mean += (slice_stdev * slice_voxel_count) + slice_weighted_snr_mean += (slice_snr * slice_voxel_count) + + total_voxel_count += slice_voxel_count + + return slice_weighted_snr_mean / total_voxel_count + + + def get_mean_slice_intensitites(self): + + data = self.generate_mask() + mask = np.ma.getmask(data) + dim_x, dim_y, dim_z, dim_t = data.shape + + slice_intensity_means = np.zeros( (dim_z,dim_t) ) + slice_voxel_counts = np.zeros( (dim_z), dtype='uint32' ) + slice_size = dim_x * dim_y + + for slice_idx in range(dim_z): + slice_voxel_counts[slice_idx] = slice_size - mask[:,:,slice_idx,0].sum() + + for volume_idx in range(dim_t): + for slice_idx in range(dim_z): + slice_data = data[:,:,slice_idx,volume_idx] + slice_intensity_means[slice_idx,volume_idx] = slice_data.mean() + + return slice_intensity_means, slice_voxel_counts, data + + + def generate_mask(self): + + mean_data = np.mean(self._fdata_array,axis=3) + + pdb.set_trace() + + numpy_3d_mask = np.zeros(mean_data.shape, dtype=bool) + + to_mask = (mean_data <= self._slice_means[5]['mask_threshold']) + + mask_lower_count = int(to_mask.sum()) + + numpy_3d_mask = numpy_3d_mask | to_mask + + numpy_4d_mask = np.zeros(self._fdata_array.shape, dtype=bool) + + numpy_4d_mask[numpy_3d_mask] = 1 + + masked_data = np.ma.masked_array(self._fdata_array,mask=numpy_4d_mask) + + return masked_data + + + + + + #pdb.set_trace() + + ''' + + if self._slice_means_array.shape[0] != self._mask_array.shape[2]: + raise ValueError("The number of slices in mean_voxel_intensities must equal the number of slices in binary_mask.") + + + slice_weighted_mean_mean = 0.0 + slice_weighted_stdev_mean = 0.0 + total_voxel_count = 0 + + for slice_idx in range(self._slice_means_array.shape[0]): + slice_voxel_count = np.sum(self._mask_array[:,:,slice_idx] == False) + slice_data = self._slice_means_array[slice_idx, :] + + if slice_voxel_count > 0: + slice_mean = np.mean(slice_data) + slice_stdev = np.std(slice_data, ddof=1) + + if np.isfinite(slice_mean) and np.isfinite(slice_stdev): + slice_weighted_mean_mean += (slice_mean * slice_voxel_count) + slice_weighted_stdev_mean += (slice_stdev * slice_voxel_count) + total_voxel_count += slice_voxel_count + + if total_voxel_count == 0: + raise ValueError("Total voxel count is zero. Check your binary mask.") + + slice_weighted_mean = slice_weighted_mean_mean / total_voxel_count + slice_weighted_stdev = slice_weighted_stdev_mean / total_voxel_count + + slice_weighted_snr = slice_weighted_mean / slice_weighted_stdev if slice_weighted_stdev != 0 else 0 + + return slice_weighted_snr + ''' + + ''' + slice_wmean_sum = 0 + slice_wstd_sum = 0 + slice_wsnr_sum = 0 + total_masked_voxels = 0 + total_voxel_count = 0 + slice_summary = [] + + slice_voxel_counts = self.get_slice_voxel_count() + + for slice_idx in range(self._slice_means_array.shape[0]): + + volume_means = self._slice_means_array[slice_idx, :] + # compute the raw slice mean across time + mean = volume_means.mean() + # compute the raw slice std across time + std = volume_means.std(ddof=1) + if std == 0: + snr = 0 + else: + snr = mean / std + + # store raw statistics within the slice summary + #slice_summary.append(SliceRow(i + 1, mean, std, snr)) + # count the number of un-masked voxels for the i'th slice + + #num_masked_voxels = (self._mask_array[:,:,slice_idx] == True).sum() + + slice_voxel_count = slice_voxel_counts[slice_idx] + + slice_wsnr_sum += snr * slice_voxel_count + + total_voxel_count += slice_voxel_count + + # keep a running sum of the number of un-masked voxels + + #total_masked_voxels += num_masked_voxels + + # keep running sums of the weighted mean, std, and snr + #slice_wmean_sum += mean * num_masked_voxels + #slice_wstd_sum += std * num_masked_voxels + + #slice_wsnr_sum += snr * num_masked_voxels + + # compute the weighted slice-based mean, standard deviation, and snr + #wm = slice_wmean_sum / total_masked_voxels + #ws = slice_wstd_sum / total_masked_voxels + wsnr = slice_wsnr_sum / total_voxel_count + + logger.debug(wsnr) + + return wsnr + ''' + + def get_slice_voxel_count(self, lower_threshold_to_zero=None): + """ + Compute the voxel counts per slice based on a masking array and optional thresholding. + """ + + # Apply mask to the data with given mask array + #masked_data = np.ma.masked_array(self._slice_means_array[:, i], mask=self._mask_array[]) + + # Apply threshold if specified + #if lower_threshold_to_zero is not None: + # masked_data = np.ma.where(masked_data > lower_threshold_to_zero, masked_data, 0) + + # Number of slices (z-dimension) + num_slices = self._slice_means_array.shape[0] + + # Initialize the array to hold voxel counts for each slice + slice_voxel_counts = np.zeros(num_slices, dtype=int) + + # Calculate voxels for each slice in the z-dimension + for slice_idx in range(num_slices): + #pdb.set_trace() + #masked_data = np.ma.masked_array(self._slice_means_array[slice_idx,:], mask=self._mask_array[:,:,slice_idx]) + # Counting unmasked (valid) voxels in each slice + slice_voxel_counts[slice_idx] = np.sum(~self._mask_array[slice_idx]) + + return slice_voxel_counts + + + def check_snr(self, key): + tasks = list() + current = self._slice_means[key] + + current_idx = self._slice_means.bisect_left(key) + + try: + value = self._slice_means.values()[current_idx] + tasks.append(value) + except IndexError: + pass + + return tasks + + ## get the index thats snr_interval (10, for example) to the left of the current + ''' + try: + left_index = current_idx - int(self._snr_interval) + for idx in range(left_index, current_idx): + value = self._instances.values()[idx] + logger.debug(f'adding dicom {value["path"]} to snr calculation task') + tasks.append(value) + except IndexError: + pass + + return tasks + + # extract the latest snr value calculation + + if len(self._instances) == self._snr_interval: + prev_snr = 0.0 + else: + for key in reversed(self._instances.keys()): + if self._instances[key]['snr'] is None: + continue + else: + prev_snr = float(self._instances[key]['snr']) + logger.info(f'found previously calculated snr at index {key}') + return tasks, prev_snr + ''' + diff --git a/scanbuddy/proc/snr.py b/scanbuddy/proc/snr.py new file mode 100644 index 0000000..bdf81b1 --- /dev/null +++ b/scanbuddy/proc/snr.py @@ -0,0 +1,253 @@ +import os +import sys +import pdb +import glob +import json +import time +import shutil +import random +import logging +import pydicom +import subprocess +import numpy as np +import nilearn.image +import nibabel as nib +from pubsub import pub +import collections as c +from pathlib import Path + +logger = logging.getLogger(__name__) + +class SNR: + def __init__(self): + pub.subscribe(self.listener, 'snr') + + + def listener(self, tasks): + logger.info('received tasks for snr calculation') + self.snr_tasks = tasks + #self._prev_snr = prev_snr + + self.run() + + def run(self): + self.get_num_tasks() + + start = time.time() + + dcm = self.read_dicoms(self._num_tasks-1) + + logger.info(f'calculating slice means for volume {int(dcm.InstanceNumber)}') + + snr_nii = self.run_dcm2niix(dcm) + + #slice_means, mask = self.calculate_snr(dcm, snr_nii) + + data_array, mask_threshold = self.calculate_snr(dcm, snr_nii) + + self.insert_snr(data_array, self.snr_tasks[0], mask_threshold) + + self.clean_dir() + + elapsed = time.time() - start + + logger.info(f'snr processing took {elapsed} seconds') + + + def calculate_snr(self, first_dcm, snr_nii): + ''' + dcm 1 has already had an snr calculation run for it, + so that value will be used as part of the running average. + Calculate the snr for dcm2 and then add that to dcm 1 and + take the average + ''' + + mask_threshold = self.get_mask_threshold(first_dcm) + + mask, data_array = self.generate_mask(mask_threshold, snr_nii) + + #current_slice_means = self.slice_based_stats(data_array, mask) + + return data_array, mask_threshold + + ''' + if self._prev_snr == 0.0: + running_snr = round(current_snr, 2) + + else: + running_snr = float(round((current_snr + self._prev_snr) / 2, 2)) + logger.info(f'value of calculated snr: {current_snr}') + logger.info(f'value of previously calculated snr: {self._prev_snr}') + logger.info(f'value of running snr: {running_snr}') + + return running_snr + ''' + + def insert_snr(self, slice_means, task, mask): + logger.debug(f'state of tasks when inserting {self.snr_tasks}') + x, y, z = slice_means.shape + data_array_4d = slice_means.reshape(x, y, z, 1) + task['slice_means'] = data_array_4d + task['mask_threshold'] = mask + + + def run_dcm2niix(self, dicom): + + self._outdir = Path(os.path.dirname(self.snr_tasks[0]['path']), 'tmp_snr_dir') + os.makedirs(self._outdir, exist_ok=True) + + for idx, dicom in enumerate(self.snr_tasks): + shutil.copy(self.snr_tasks[idx]['path'], self._outdir) + + dcm2niix_cmd = [ + 'dcm2niix', + '-b', 'y', + '-s', 'y' + '-f', f'snr_calc', + '-o', self._outdir, + self._outdir + ] + + output = subprocess.check_output(dcm2niix_cmd, stderr=subprocess.STDOUT) + + logger.debug(f'dcm2niix output: {output}') + + nii_file = self.find_nii(self._outdir) + + return nii_file + + + def find_nii(self, directory): + for file in os.listdir(directory): + if f'tmp_snr_dir' in file and file.endswith('.nii'): + return Path(directory, file) + + def generate_mask(self, mask_threshold, nii): + nii_array = nib.load(nii).get_fdata() + binary_mask = nii_array <= mask_threshold + return binary_mask, nii_array + ''' + keeping this code for posterity. maybe one day we will + want a pre-masked array + + thresholded_nii = nilearn.image.threshold_img(nii, mask_threshold) + output_file_path = Path(os.path.dirname(nii), 'thresholded_img.nii') + nib.save(thresholded_nii, output_file_path) + masked_data_array = nib.load(output_file_path).get_fdata() + return (masked_data_array) + ''' + + def get_mask_threshold(self, ds): + bits_stored = ds.get('BitsStored', None) + receive_coil = self.find_coil(ds) + + if bits_stored == 12: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 150.0') + return 150.0 + if bits_stored == 16: + if receive_coil in ['Head_32']: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 1500.0') + return 1500.0 + if receive_coil in ['Head_64', 'HeadNeck_64']: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 3000.0') + return 3000.0 + raise MaskThresholdError(f'unexpected bits stored "{bits_stored}" + receive coil "{receive_coil}"') + + def slice_based_stats(self, dat, mask): + ''' + Compute slice-based statistics + + :param dat: Data matrix + :type dat: numpy.array + :param mask: Binary mask + :type mask: numpy.array + :returns: Tuple of ['summary', 'mean', 'std', 'snr'] + :rtype: collections.namedtuple + ''' + # get the length of each dimension into named variables + _,_,z = dat.shape + # create slice means for every volume independently + slice_means =np.zeros(z) + volume = np.ma.array(dat[:,:,:], mask=mask[:,:,:], fill_value=0.0) + for j in range(z): + # Get the j'th slice + slice_ij = volume[:, :, j] + # Compute the mean of the slice, automatically handling mask + x = slice_ij.mean() + # Store the mean; if x is masked, it will return the fill_value by default + slice_means[j] = x if not np.ma.is_masked(x) else 0.0 + + return slice_means + + # we are breaking this function into two parts for now. The rest of this + # this function will be in the proc __init__.py file + ''' + + # compute the raw and weighted means, stds, and snrs + slice_wmean_sum = 0 + slice_wstd_sum = 0 + slice_wsnr_sum = 0 + total_masked_voxels = 0 + slice_summary = [] + + for i,volume_means in iter(slice_means.items()): + # cast volume means to a numpy array to make life easier + volume_means = np.array(volume_means) + # compute the raw slice mean across time + mean = volume_means.mean() + # compute the raw slice std across time + std = volume_means.std(ddof=1) + if std == 0: + snr = 0 + else: + snr = mean / std + + # store raw statistics within the slice summary + #slice_summary.append(SliceRow(i + 1, mean, std, snr)) + # count the number of un-masked voxels for the i'th slice + num_masked_voxels = (mask[:,:,i] == False).sum() + # keep a running sum of the number of un-masked voxels + total_masked_voxels += num_masked_voxels + # keep running sums of the weighted mean, std, and snr + slice_wmean_sum += volume_means.mean() * num_masked_voxels + slice_wstd_sum += volume_means.std(ddof=1) * num_masked_voxels + slice_wsnr_sum += snr * num_masked_voxels + # compute the weighted slice-based mean, standard deviation, and snr + wm = slice_wmean_sum / total_masked_voxels + ws = slice_wstd_sum / total_masked_voxels + wsnr = slice_wsnr_sum / total_masked_voxels + + return wsnr + + ''' + + def read_dicoms(self, last_idx): + logger.debug(f'state of tasks when reading dicom: {self.snr_tasks}') + dcm1 = self.snr_tasks[0]['path'] + #dcm2 = self.snr_tasks[last_idx]['path'] + + ds1 = pydicom.dcmread(dcm1, force=True, stop_before_pixels=True) + #ds2 = pydicom.dcmread(dcm2, force=True, stop_before_pixels=True) + + return ds1 + + def find_coil(self, ds): + seq = ds[(0x5200, 0x9229)][0] + seq = seq[(0x0018, 0x9042)][0] + return seq[(0x0018, 0x1250)].value + + def clean_dir(self): + for file in glob.glob(f'{self._outdir}/*.json'): + os.remove(file) + for file in glob.glob(f'{self._outdir}/*.nii'): + os.remove(file) + for file in glob.glob(f'{self._outdir}/*.dcm'): + os.remove(file) + + def get_num_tasks(self): + self._num_tasks = len(self.snr_tasks) + + +class MaskThresholdError(Exception): + pass + diff --git a/scanbuddy/view/dash.py b/scanbuddy/view/dash.py index c316c24..6afa1f0 100644 --- a/scanbuddy/view/dash.py +++ b/scanbuddy/view/dash.py @@ -34,6 +34,7 @@ def __init__(self, host='127.0.0.1', port=8080, config=None, debug=False): self._subtitle = 'Ready' self._num_warnings = 0 self._instances = dict() + self._current_snr = 0.0 self._redis_client = redis.StrictRedis( host=self._config.find_one('$.broker.host', default='127.0.0.1'), port=self._config.find_one('$.broker.port', default=6379), @@ -43,6 +44,7 @@ def __init__(self, host='127.0.0.1', port=8080, config=None, debug=False): self.init_page() self.init_callbacks() pub.subscribe(self.listener, 'plot') + pub.subscribe(self.plot_snr, 'plot_snr') def init_app(self): self._app = Dash( @@ -276,6 +278,36 @@ def init_page(self): ], style={"margin": "0px"} ), + dbc.Row( + [ + dbc.Col( + "SNR", + width=8, + style={ + "borderRight": "1px solid black", + "textAlign": "center", + "display": "flex", + "alignItems": "center", + "justifyContent": "flex-end", + "paddingRight": "5px", + "fontSize": "1.5vw", + "borderBottom": "1px solid black" + } + ), + dbc.Col( + id='snr', + children="0", + width=4, + style={ + "borderBottom": "1px solid black", + "textAlign": "center", + "padding": "1rem", + "fontSize": "1.25vw" + } + ) + ], + style={"margin": "0px"} + ), ], style={"border": "1px solid black", "padding": "0"} ) @@ -393,6 +425,7 @@ def init_callbacks(self): Output('movements-05mm', 'children'), Output('movements-1mm', 'children'), Output('max-abs-motion', 'children'), + Output('snr', 'children'), Input('plot-interval-component', 'n_intervals'), )(self.update_metrics) @@ -429,7 +462,15 @@ def update_metrics(self, n): else: max_abs_motion = 0 - return str(num_vols), str(movements_05mm), str(movements_1mm), str(max_abs_motion) + snr = self.get_snr() + + return str(num_vols), str(movements_05mm), str(movements_1mm), str(max_abs_motion), str(snr) + + def get_snr(self): + if not self._current_snr: + return 0.0 + else: + return self._current_snr def get_subtitle(self): return self._subtitle @@ -643,6 +684,7 @@ def todataframe(self): arr = list() for i,instance in enumerate(self._instances.values(), start=1): volreg = instance['volreg'] + #snr = instance.get('snr', '0.0') if volreg: arr.append([i] + volreg) df = pd.DataFrame(arr, columns=['N', 'roll', 'pitch', 'yaw', 'x', 'y', 'z']) @@ -659,5 +701,8 @@ def listener(self, instances, subtitle_string): self._instances = instances self._subtitle = subtitle_string + def plot_snr(self, snr_metric): + self._current_snr = snr_metric + class AuthError(Exception): pass diff --git a/scripts/start.py b/scripts/start.py index 72ef336..5566999 100755 --- a/scripts/start.py +++ b/scripts/start.py @@ -10,6 +10,7 @@ from scanbuddy.proc import Processor from scanbuddy.proc.volreg import VolReg from scanbuddy.proc.params import Params +from scanbuddy.proc.snr import SNR from scanbuddy.view.dash import View from scanbuddy.broker.redis import MessageBroker from scanbuddy.config import Config @@ -25,6 +26,8 @@ def main(): parser.add_argument('--port', type=int, default=8080) parser.add_argument('--folder', type=Path, required=True) parser.add_argument('-v', '--verbose', action='store_true') + parser.add_argument('--snr-interval', default=10, + help='Every N volumes snr should be calculated') args = parser.parse_args() config = Config(args.config) @@ -37,6 +40,7 @@ def main(): config=config ) volreg = VolReg(mock=args.mock) + snr = SNR() view = View( host=args.host, port=args.port, diff --git a/setup.py b/setup.py index 9217fb7..72eae9a 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,8 @@ 'redis', 'pyyaml', 'jsonpath-ng', - 'dash_auth' + 'dash_auth', + 'nibabel' ] about = dict() From cf45af2cc60e368336b10f0a784d4aff81a3d32e Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Mon, 11 Nov 2024 12:37:39 -0500 Subject: [PATCH 2/8] init empty numpy array --- scanbuddy/proc/__init__.py | 50 +++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 17 deletions(-) diff --git a/scanbuddy/proc/__init__.py b/scanbuddy/proc/__init__.py index 15dcc78..5488f75 100644 --- a/scanbuddy/proc/__init__.py +++ b/scanbuddy/proc/__init__.py @@ -4,6 +4,7 @@ import math import json import logging +import threading import numpy as np from pubsub import pub from sortedcontainers import SortedDict @@ -36,6 +37,8 @@ def listener(self, ds, path): logger.debug('current state of instances') logger.debug(json.dumps(self._instances, default=list, indent=2)) + num_vols = ds[(0x0020, 0x0105)].value + tasks = self.check_volreg(key) ### edits need to be made here snr_tasks = self.check_snr(key) @@ -51,24 +54,33 @@ def listener(self, ds, path): logger.debug(json.dumps(self._instances, indent=2)) - if key == 5: - self._fdata_array = self._slice_means[key]['slice_means'] - #self._slice_means_array = self._slice_means[key]['slice_means'] - #self._mask_array = self._slice_means[key]['mask'] - elif key > 5: - self._fdata_array = np.concatenate((self._fdata_array, self._slice_means[key]['slice_means']), axis=3) - #self._slice_means_array = np.column_stack((self._slice_means_array, self._slice_means[key]['slice_means'])) - #self._mask_array = np.column_stack((self._mask_array, self._slice_means[key]['mask'])) - #if len(self._instances) % self._snr_interval == 0: - if key == 10: - snr_metric = self.calc_snr() + if key == 1: + x, y, z, _ = self._slice_means[key]['slice_means'].shape + + #self._fdata_array = np.zeros((x, y, z, num_vols)) + self._fdata_array = np.empty((x, y, z, num_vols)) + + logger.info(f'shape of zeros: {self._fdata_array.shape}') + logger.info(f'shape of first slice means: {self._slice_means[key]['slice_means'].shape}') + time.sleep(10) + + if key >= 5: + insert_position = key - 5 + #self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'] + self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'].squeeze() + + #elif key > 5: + #self._fdata_array = np.concatenate((self._fdata_array, self._slice_means[key]['slice_means']), axis=3) if key > 53: #logger.info(f'shape of slice_means_array: {self._slice_means_array.shape}') logger.info(f'shape of fdata_array: {self._fdata_array.shape}') logger.info('publishing message to plot_snr topic') - snr_metric = round(self.calc_snr(), 2) - logger.info(f'running snr metric: {snr_metric}') - pub.sendMessage('plot_snr', snr_metric=snr_metric) + + snr_thread = threading.Thread(target=self.calculate_and_publish_snr) + snr_thread.start() + #snr_metric = round(self.calc_snr(), 2) + #logger.info(f'running snr metric: {snr_metric}') + #pub.sendMessage('plot_snr', snr_metric=snr_metric) logger.debug(f'after volreg') logger.debug(json.dumps(self._instances, indent=2)) @@ -79,6 +91,11 @@ def listener(self, ds, path): subtitle_string = f'{project} • {session} • {scandesc} • {scannum}' pub.sendMessage('plot', instances=self._instances, subtitle_string=subtitle_string) + def calculate_and_publish_snr(self): + snr_metric = round(self.calc_snr(), 2) + logger.info(f'running snr metric: {snr_metric}') + pub.sendMessage('plot_snr', snr_metric=snr_metric) + def check_volreg(self, key): tasks = list() current = self._instances[key] @@ -134,6 +151,7 @@ def calc_snr(self): total_voxel_count += slice_voxel_count + #return slice_weighted_mean_mean / total_voxel_count return slice_weighted_snr_mean / total_voxel_count @@ -160,9 +178,7 @@ def get_mean_slice_intensitites(self): def generate_mask(self): - mean_data = np.mean(self._fdata_array,axis=3) - - pdb.set_trace() + mean_data = np.mean(self._fdata_array, axis=3) numpy_3d_mask = np.zeros(mean_data.shape, dtype=bool) From 6733fb292cfa50a407e4f9f63e49cdcbf44d38f8 Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Wed, 27 Nov 2024 14:08:04 -0500 Subject: [PATCH 3/8] threading using python v 3.13 (no gil) --- scanbuddy/proc/__init__.py | 362 ++++++++++++++++++------------------- scanbuddy/proc/snr.py | 204 ++------------------- scanbuddy/proc/volreg.py | 18 +- 3 files changed, 209 insertions(+), 375 deletions(-) diff --git a/scanbuddy/proc/__init__.py b/scanbuddy/proc/__init__.py index 5488f75..5cf4a6d 100644 --- a/scanbuddy/proc/__init__.py +++ b/scanbuddy/proc/__init__.py @@ -1,9 +1,12 @@ +import os import sys import pdb import time import math import json +import shutil import logging +import datetime import threading import numpy as np from pubsub import pub @@ -27,17 +30,19 @@ def listener(self, ds, path): key = int(ds.InstanceNumber) self._instances[key] = { 'path': path, - 'volreg': None + 'volreg': None, + 'nii_path': None } self._slice_means[key] = { 'path': path, 'slice_means': None, - 'mask_threshold': None + 'mask_threshold': None, + 'mask': None } logger.debug('current state of instances') logger.debug(json.dumps(self._instances, default=list, indent=2)) - num_vols = ds[(0x0020, 0x0105)].value + tasks = self.check_volreg(key) ### edits need to be made here @@ -47,54 +52,86 @@ def listener(self, ds, path): pub.sendMessage('volreg', tasks=tasks) logger.debug(f'publishing message to params topic') pub.sendMessage('params', ds=ds) - logger.debug(f'publishing message to snr topic') + logger.debug(f'publishing message to snr_fdata topic') logger.debug(f'snr task sorted dict: {snr_tasks}') - pub.sendMessage('snr', tasks=snr_tasks) + pub.sendMessage('snr', nii_path=self._instances[key]['nii_path'], tasks=snr_tasks) logger.debug('after snr calculation') logger.debug(json.dumps(self._instances, indent=2)) - if key == 1: - x, y, z, _ = self._slice_means[key]['slice_means'].shape - #self._fdata_array = np.zeros((x, y, z, num_vols)) - self._fdata_array = np.empty((x, y, z, num_vols)) + logger.debug(f'after volreg') + logger.debug(json.dumps(self._instances, indent=2)) + project = ds.get('StudyDescription', '[STUDY]') + session = ds.get('PatientID', '[PATIENT]') + scandesc = ds.get('SeriesDescription', '[SERIES]') + scannum = ds.get('SeriesNumber', '[NUMBER]') + subtitle_string = f'{project} • {session} • {scandesc} • {scannum}' + pub.sendMessage('plot', instances=self._instances, subtitle_string=subtitle_string) + + if key < 5: + self._num_vols = ds[(0x0020, 0x0105)].value + self._mask_threshold, self._decrement = self.get_mask_threshold(ds) + x, y, self._z, _ = self._slice_means[key]['slice_means'].shape + + self._fdata_array = np.empty((x, y, self._z, self._num_vols)) + #self._slice_intensity_means = np.zeros( (self._z, math.ceil(self._num_vols // 2)) ) + self._slice_intensity_means = np.zeros( (self._z, self._num_vols) ) + + logger.info(f'shape of zeros: {self._fdata_array.shape}') logger.info(f'shape of first slice means: {self._slice_means[key]['slice_means'].shape}') - time.sleep(10) if key >= 5: insert_position = key - 5 - #self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'] self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'].squeeze() - #elif key > 5: - #self._fdata_array = np.concatenate((self._fdata_array, self._slice_means[key]['slice_means']), axis=3) - if key > 53: - #logger.info(f'shape of slice_means_array: {self._slice_means_array.shape}') - logger.info(f'shape of fdata_array: {self._fdata_array.shape}') - logger.info('publishing message to plot_snr topic') + #if key == 54: + # logger.info('launching calculate and publish snr thread') + + # snr_thread = threading.Thread(target=self.calculate_and_publish_snr, args=(key,)) + # snr_thread.start() + - snr_thread = threading.Thread(target=self.calculate_and_publish_snr) + if key > 53 and (key % 3 == 0) and key < self._num_vols: + logger.info('launching calculate and publish snr thread') + + snr_thread = threading.Thread(target=self.calculate_and_publish_snr, args=(key,)) snr_thread.start() - #snr_metric = round(self.calc_snr(), 2) + #snr_metric = round(self.calc_snr(key), 2) #logger.info(f'running snr metric: {snr_metric}') - #pub.sendMessage('plot_snr', snr_metric=snr_metric) + #pub.sendMessage('plot_snr', snr_metric=snr_metric) + + if key == self._num_vols: + time.sleep(2) + data_path = os.path.dirname(self._instances[key]['path']) + logger.info(f'removing dicom dir: {data_path}') + shutil.rmtree(data_path) + + #if key == self._num_vols: + # logger.info('RUNNING FINAL SNR CALCULATION') + # snr_metric = round(self.calc_snr(key), 2) + # logger.info(f'final snr metric: {snr_metric}') + # pub.sendMessage('plot_snr', snr_metric=snr_metric) - logger.debug(f'after volreg') - logger.debug(json.dumps(self._instances, indent=2)) - project = ds.get('StudyDescription', '[STUDY]') - session = ds.get('PatientID', '[PATIENT]') - scandesc = ds.get('SeriesDescription', '[SERIES]') - scannum = ds.get('SeriesNumber', '[NUMBER]') - subtitle_string = f'{project} • {session} • {scandesc} • {scannum}' - pub.sendMessage('plot', instances=self._instances, subtitle_string=subtitle_string) - def calculate_and_publish_snr(self): - snr_metric = round(self.calc_snr(), 2) + + def calculate_and_publish_snr(self, key): + start = time.time() + snr_metric = round(self.calc_snr(key), 2) + elapsed = time.time() - start + #self._plot_dict[self._key] = elapsed + logger.info(f'snr calculation took {elapsed} seconds') logger.info(f'running snr metric: {snr_metric}') - pub.sendMessage('plot_snr', snr_metric=snr_metric) + if np.isnan(snr_metric): + logger.info(f'snr is a nan, decrementing mask threshold by {self._decrement}') + self._mask_threshold = self._mask_threshold - self._decrement + logger.info(f'new threshold: {self._mask_threshold}') + #self._numpy_4d_mask = np.zeros(self._fdata_array.shape, dtype=bool) + self._slice_intensity_means = np.zeros( (self._z, self._num_vols) ) + else: + pub.sendMessage('plot_snr', snr_metric=snr_metric) def check_volreg(self, key): tasks = list() @@ -123,12 +160,16 @@ def check_volreg(self, key): return tasks - def calc_snr(self): + def calc_snr(self, key): + + slice_intensity_means, slice_voxel_counts, data = self.get_mean_slice_intensitites(key) - slice_intensity_means, slice_voxel_counts, data = self.get_mean_slice_intensitites() + non_zero_columns = ~np.all(slice_intensity_means == 0, axis=0) - slice_count = slice_intensity_means.shape[0] - volume_count = slice_intensity_means.shape[1] + slice_intensity_means_2 = slice_intensity_means[:, non_zero_columns] + + slice_count = slice_intensity_means_2.shape[0] + volume_count = slice_intensity_means_2.shape[1] slice_weighted_mean_mean = 0 slice_weighted_stdev_mean = 0 @@ -139,7 +180,7 @@ def calc_snr(self): total_voxel_count = 0 for slice_idx in range(slice_count): - slice_data = slice_intensity_means[slice_idx] + slice_data = slice_intensity_means_2[slice_idx] slice_voxel_count = slice_voxel_counts[slice_idx] slice_mean = slice_data.mean() slice_stdev = slice_data.std(ddof=1) @@ -151,171 +192,149 @@ def calc_snr(self): total_voxel_count += slice_voxel_count - #return slice_weighted_mean_mean / total_voxel_count + logger.debug(f"Slice {slice_idx}: Mean={slice_mean}, StdDev={slice_stdev}, SNR={slice_snr}") + + return slice_weighted_snr_mean / total_voxel_count - def get_mean_slice_intensitites(self): + def get_mean_slice_intensitites(self, key): - data = self.generate_mask() + data = self.generate_mask(key) mask = np.ma.getmask(data) - dim_x, dim_y, dim_z, dim_t = data.shape + dim_x, dim_y, dim_z, _ = data.shape + + dim_t = key - 4 + + if key > 55: + start = time.time() + differing_slices = self.find_mask_differences(key) + logger.info(f'finding mask differences took {time.time() - start}') - slice_intensity_means = np.zeros( (dim_z,dim_t) ) + + #slice_intensity_means = np.zeros( (dim_z,dim_t) ) slice_voxel_counts = np.zeros( (dim_z), dtype='uint32' ) slice_size = dim_x * dim_y for slice_idx in range(dim_z): slice_voxel_counts[slice_idx] = slice_size - mask[:,:,slice_idx,0].sum() - for volume_idx in range(dim_t): - for slice_idx in range(dim_z): - slice_data = data[:,:,slice_idx,volume_idx] - slice_intensity_means[slice_idx,volume_idx] = slice_data.mean() - - return slice_intensity_means, slice_voxel_counts, data - - - def generate_mask(self): - - mean_data = np.mean(self._fdata_array, axis=3) - - numpy_3d_mask = np.zeros(mean_data.shape, dtype=bool) - - to_mask = (mean_data <= self._slice_means[5]['mask_threshold']) - - mask_lower_count = int(to_mask.sum()) - - numpy_3d_mask = numpy_3d_mask | to_mask - - numpy_4d_mask = np.zeros(self._fdata_array.shape, dtype=bool) - - numpy_4d_mask[numpy_3d_mask] = 1 - - masked_data = np.ma.masked_array(self._fdata_array,mask=numpy_4d_mask) - - return masked_data - - - - - #pdb.set_trace() - ''' + zero_columns = np.where(np.all(self._slice_intensity_means[:,:dim_t] == 0, axis=0))[0].tolist() - if self._slice_means_array.shape[0] != self._mask_array.shape[2]: - raise ValueError("The number of slices in mean_voxel_intensities must equal the number of slices in binary_mask.") + logger.info(f'volumes being calculated: {zero_columns}') - slice_weighted_mean_mean = 0.0 - slice_weighted_stdev_mean = 0.0 - total_voxel_count = 0 - - for slice_idx in range(self._slice_means_array.shape[0]): - slice_voxel_count = np.sum(self._mask_array[:,:,slice_idx] == False) - slice_data = self._slice_means_array[slice_idx, :] + if len(zero_columns) > 20: + for volume_idx in range(dim_t): + for slice_idx in range(dim_z): + slice_data = data[:,:,slice_idx,volume_idx] + self._slice_intensity_means[slice_idx,volume_idx] = slice_data.mean() - if slice_voxel_count > 0: - slice_mean = np.mean(slice_data) - slice_stdev = np.std(slice_data, ddof=1) - - if np.isfinite(slice_mean) and np.isfinite(slice_stdev): - slice_weighted_mean_mean += (slice_mean * slice_voxel_count) - slice_weighted_stdev_mean += (slice_stdev * slice_voxel_count) - total_voxel_count += slice_voxel_count + else: - if total_voxel_count == 0: - raise ValueError("Total voxel count is zero. Check your binary mask.") + for volume_idx in zero_columns: + for slice_idx in range(dim_z): + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - slice_weighted_mean = slice_weighted_mean_mean / total_voxel_count - slice_weighted_stdev = slice_weighted_stdev_mean / total_voxel_count + logger.info(f'recalculating slice means at the following slices: {differing_slices}') + logger.info(f'total of {len(differing_slices)} new slices being computed') - slice_weighted_snr = slice_weighted_mean / slice_weighted_stdev if slice_weighted_stdev != 0 else 0 + if differing_slices: - return slice_weighted_snr - ''' + if key == self._num_vols: + for volume_idx in range(dim_t): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - ''' - slice_wmean_sum = 0 - slice_wstd_sum = 0 - slice_wsnr_sum = 0 - total_masked_voxels = 0 - total_voxel_count = 0 - slice_summary = [] + elif key % 4 == 0: + for volume_idx in range(0, dim_t, 6): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - slice_voxel_counts = self.get_slice_voxel_count() + elif key % 3 == 0: + for volume_idx in range(3, dim_t, 6): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean + ''' + for volume_idx in range(dim_t): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean + ''' - for slice_idx in range(self._slice_means_array.shape[0]): + return self._slice_intensity_means[:, :dim_t], slice_voxel_counts, data - volume_means = self._slice_means_array[slice_idx, :] - # compute the raw slice mean across time - mean = volume_means.mean() - # compute the raw slice std across time - std = volume_means.std(ddof=1) - if std == 0: - snr = 0 - else: - snr = mean / std - # store raw statistics within the slice summary - #slice_summary.append(SliceRow(i + 1, mean, std, snr)) - # count the number of un-masked voxels for the i'th slice + def generate_mask(self, key): - #num_masked_voxels = (self._mask_array[:,:,slice_idx] == True).sum() + mean_data = np.mean(self._fdata_array[...,:key-4], axis=3) - slice_voxel_count = slice_voxel_counts[slice_idx] + numpy_3d_mask = np.zeros(mean_data.shape, dtype=bool) - slice_wsnr_sum += snr * slice_voxel_count + to_mask = (mean_data <= self._mask_threshold) - total_voxel_count += slice_voxel_count + mask_lower_count = int(to_mask.sum()) - # keep a running sum of the number of un-masked voxels + numpy_3d_mask = numpy_3d_mask | to_mask - #total_masked_voxels += num_masked_voxels + numpy_4d_mask = np.zeros(self._fdata_array[..., :key-4].shape, dtype=bool) - # keep running sums of the weighted mean, std, and snr - #slice_wmean_sum += mean * num_masked_voxels - #slice_wstd_sum += std * num_masked_voxels + numpy_4d_mask[numpy_3d_mask] = 1 - #slice_wsnr_sum += snr * num_masked_voxels + masked_data = np.ma.masked_array(self._fdata_array[..., :key-4], mask=numpy_4d_mask) - # compute the weighted slice-based mean, standard deviation, and snr - #wm = slice_wmean_sum / total_masked_voxels - #ws = slice_wstd_sum / total_masked_voxels - wsnr = slice_wsnr_sum / total_voxel_count + mask = np.ma.getmask(masked_data) - logger.debug(wsnr) - - return wsnr - ''' + self._slice_means[key]['mask'] = mask - def get_slice_voxel_count(self, lower_threshold_to_zero=None): - """ - Compute the voxel counts per slice based on a masking array and optional thresholding. - """ + return masked_data - # Apply mask to the data with given mask array - #masked_data = np.ma.masked_array(self._slice_means_array[:, i], mask=self._mask_array[]) + def find_mask_differences(self, key): + num_old_vols = key - 7 + prev_mask = self._slice_means[key-3]['mask'] + current_mask = self._slice_means[key]['mask'] + differences = prev_mask != current_mask[:,:,:,:num_old_vols] + diff_indices = np.where(differences) + differing_slices = [] + for index in zip(*diff_indices): + if int(index[2]) not in differing_slices: + differing_slices.append(int(index[2])) + #slices_to_update = differing_slices.sort() + #pdb.set_trace() + return differing_slices - # Apply threshold if specified - #if lower_threshold_to_zero is not None: - # masked_data = np.ma.where(masked_data > lower_threshold_to_zero, masked_data, 0) - # Number of slices (z-dimension) - num_slices = self._slice_means_array.shape[0] - - # Initialize the array to hold voxel counts for each slice - slice_voxel_counts = np.zeros(num_slices, dtype=int) + def get_mask_threshold(self, ds): + bits_stored = ds.get('BitsStored', None) + receive_coil = self.find_coil(ds) - # Calculate voxels for each slice in the z-dimension - for slice_idx in range(num_slices): - #pdb.set_trace() - #masked_data = np.ma.masked_array(self._slice_means_array[slice_idx,:], mask=self._mask_array[:,:,slice_idx]) - # Counting unmasked (valid) voxels in each slice - slice_voxel_counts[slice_idx] = np.sum(~self._mask_array[slice_idx]) + if bits_stored == 12: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 150.0') + return 150.0, 10 + if bits_stored == 16: + if receive_coil in ['Head_32']: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 1500.0') + return 1500.0, 100 + if receive_coil in ['Head_64', 'HeadNeck_64']: + logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 3000.0') + return 3000.0, 300 + raise MaskThresholdError(f'unexpected bits stored "{bits_stored}" + receive coil "{receive_coil}"') - return slice_voxel_counts + def find_coil(self, ds): + seq = ds[(0x5200, 0x9229)][0] + seq = seq[(0x0018, 0x9042)][0] + return seq[(0x0018, 0x1250)].value def check_snr(self, key): @@ -332,30 +351,3 @@ def check_snr(self, key): return tasks - ## get the index thats snr_interval (10, for example) to the left of the current - ''' - try: - left_index = current_idx - int(self._snr_interval) - for idx in range(left_index, current_idx): - value = self._instances.values()[idx] - logger.debug(f'adding dicom {value["path"]} to snr calculation task') - tasks.append(value) - except IndexError: - pass - - return tasks - - # extract the latest snr value calculation - - if len(self._instances) == self._snr_interval: - prev_snr = 0.0 - else: - for key in reversed(self._instances.keys()): - if self._instances[key]['snr'] is None: - continue - else: - prev_snr = float(self._instances[key]['snr']) - logger.info(f'found previously calculated snr at index {key}') - return tasks, prev_snr - ''' - diff --git a/scanbuddy/proc/snr.py b/scanbuddy/proc/snr.py index bdf81b1..bfec35d 100644 --- a/scanbuddy/proc/snr.py +++ b/scanbuddy/proc/snr.py @@ -10,7 +10,7 @@ import pydicom import subprocess import numpy as np -import nilearn.image +#import nilearn.image import nibabel as nib from pubsub import pub import collections as c @@ -23,10 +23,10 @@ def __init__(self): pub.subscribe(self.listener, 'snr') - def listener(self, tasks): - logger.info('received tasks for snr calculation') + def listener(self, nii_path, tasks): + logger.info('received tasks for fdata extraction') self.snr_tasks = tasks - #self._prev_snr = prev_snr + self._nii_path = nii_path self.run() @@ -37,217 +37,47 @@ def run(self): dcm = self.read_dicoms(self._num_tasks-1) - logger.info(f'calculating slice means for volume {int(dcm.InstanceNumber)}') + instance_num = int(dcm.InstanceNumber) - snr_nii = self.run_dcm2niix(dcm) + logger.info(f'extracting fdata for volume {instance_num}') - #slice_means, mask = self.calculate_snr(dcm, snr_nii) + data_array = self.get_nii_array() - data_array, mask_threshold = self.calculate_snr(dcm, snr_nii) + self.insert_snr(data_array, self.snr_tasks[0], None) - self.insert_snr(data_array, self.snr_tasks[0], mask_threshold) - - self.clean_dir() + self.clean_dir(instance_num) elapsed = time.time() - start - logger.info(f'snr processing took {elapsed} seconds') - - - def calculate_snr(self, first_dcm, snr_nii): - ''' - dcm 1 has already had an snr calculation run for it, - so that value will be used as part of the running average. - Calculate the snr for dcm2 and then add that to dcm 1 and - take the average - ''' - - mask_threshold = self.get_mask_threshold(first_dcm) - - mask, data_array = self.generate_mask(mask_threshold, snr_nii) - - #current_slice_means = self.slice_based_stats(data_array, mask) + logger.info(f'extracting fdata from volume {instance_num} took {elapsed} seconds') - return data_array, mask_threshold - ''' - if self._prev_snr == 0.0: - running_snr = round(current_snr, 2) + def get_nii_array(self): + return nib.load(self._nii_path).get_fdata() - else: - running_snr = float(round((current_snr + self._prev_snr) / 2, 2)) - logger.info(f'value of calculated snr: {current_snr}') - logger.info(f'value of previously calculated snr: {self._prev_snr}') - logger.info(f'value of running snr: {running_snr}') - - return running_snr - ''' def insert_snr(self, slice_means, task, mask): logger.debug(f'state of tasks when inserting {self.snr_tasks}') x, y, z = slice_means.shape data_array_4d = slice_means.reshape(x, y, z, 1) task['slice_means'] = data_array_4d - task['mask_threshold'] = mask - - - def run_dcm2niix(self, dicom): - - self._outdir = Path(os.path.dirname(self.snr_tasks[0]['path']), 'tmp_snr_dir') - os.makedirs(self._outdir, exist_ok=True) - - for idx, dicom in enumerate(self.snr_tasks): - shutil.copy(self.snr_tasks[idx]['path'], self._outdir) - - dcm2niix_cmd = [ - 'dcm2niix', - '-b', 'y', - '-s', 'y' - '-f', f'snr_calc', - '-o', self._outdir, - self._outdir - ] - - output = subprocess.check_output(dcm2niix_cmd, stderr=subprocess.STDOUT) - - logger.debug(f'dcm2niix output: {output}') - - nii_file = self.find_nii(self._outdir) - - return nii_file - - - def find_nii(self, directory): - for file in os.listdir(directory): - if f'tmp_snr_dir' in file and file.endswith('.nii'): - return Path(directory, file) - - def generate_mask(self, mask_threshold, nii): - nii_array = nib.load(nii).get_fdata() - binary_mask = nii_array <= mask_threshold - return binary_mask, nii_array - ''' - keeping this code for posterity. maybe one day we will - want a pre-masked array - - thresholded_nii = nilearn.image.threshold_img(nii, mask_threshold) - output_file_path = Path(os.path.dirname(nii), 'thresholded_img.nii') - nib.save(thresholded_nii, output_file_path) - masked_data_array = nib.load(output_file_path).get_fdata() - return (masked_data_array) - ''' - - def get_mask_threshold(self, ds): - bits_stored = ds.get('BitsStored', None) - receive_coil = self.find_coil(ds) - - if bits_stored == 12: - logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 150.0') - return 150.0 - if bits_stored == 16: - if receive_coil in ['Head_32']: - logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 1500.0') - return 1500.0 - if receive_coil in ['Head_64', 'HeadNeck_64']: - logger.debug(f'scan has "{bits_stored}" bits and receive coil "{receive_coil}", setting mask threshold to 3000.0') - return 3000.0 - raise MaskThresholdError(f'unexpected bits stored "{bits_stored}" + receive coil "{receive_coil}"') - - def slice_based_stats(self, dat, mask): - ''' - Compute slice-based statistics - - :param dat: Data matrix - :type dat: numpy.array - :param mask: Binary mask - :type mask: numpy.array - :returns: Tuple of ['summary', 'mean', 'std', 'snr'] - :rtype: collections.namedtuple - ''' - # get the length of each dimension into named variables - _,_,z = dat.shape - # create slice means for every volume independently - slice_means =np.zeros(z) - volume = np.ma.array(dat[:,:,:], mask=mask[:,:,:], fill_value=0.0) - for j in range(z): - # Get the j'th slice - slice_ij = volume[:, :, j] - # Compute the mean of the slice, automatically handling mask - x = slice_ij.mean() - # Store the mean; if x is masked, it will return the fill_value by default - slice_means[j] = x if not np.ma.is_masked(x) else 0.0 - - return slice_means - - # we are breaking this function into two parts for now. The rest of this - # this function will be in the proc __init__.py file - ''' - - # compute the raw and weighted means, stds, and snrs - slice_wmean_sum = 0 - slice_wstd_sum = 0 - slice_wsnr_sum = 0 - total_masked_voxels = 0 - slice_summary = [] - - for i,volume_means in iter(slice_means.items()): - # cast volume means to a numpy array to make life easier - volume_means = np.array(volume_means) - # compute the raw slice mean across time - mean = volume_means.mean() - # compute the raw slice std across time - std = volume_means.std(ddof=1) - if std == 0: - snr = 0 - else: - snr = mean / std - - # store raw statistics within the slice summary - #slice_summary.append(SliceRow(i + 1, mean, std, snr)) - # count the number of un-masked voxels for the i'th slice - num_masked_voxels = (mask[:,:,i] == False).sum() - # keep a running sum of the number of un-masked voxels - total_masked_voxels += num_masked_voxels - # keep running sums of the weighted mean, std, and snr - slice_wmean_sum += volume_means.mean() * num_masked_voxels - slice_wstd_sum += volume_means.std(ddof=1) * num_masked_voxels - slice_wsnr_sum += snr * num_masked_voxels - # compute the weighted slice-based mean, standard deviation, and snr - wm = slice_wmean_sum / total_masked_voxels - ws = slice_wstd_sum / total_masked_voxels - wsnr = slice_wsnr_sum / total_masked_voxels - - return wsnr - - ''' def read_dicoms(self, last_idx): logger.debug(f'state of tasks when reading dicom: {self.snr_tasks}') dcm1 = self.snr_tasks[0]['path'] - #dcm2 = self.snr_tasks[last_idx]['path'] ds1 = pydicom.dcmread(dcm1, force=True, stop_before_pixels=True) - #ds2 = pydicom.dcmread(dcm2, force=True, stop_before_pixels=True) return ds1 - def find_coil(self, ds): - seq = ds[(0x5200, 0x9229)][0] - seq = seq[(0x0018, 0x9042)][0] - return seq[(0x0018, 0x1250)].value - - def clean_dir(self): - for file in glob.glob(f'{self._outdir}/*.json'): - os.remove(file) - for file in glob.glob(f'{self._outdir}/*.nii'): + def clean_dir(self, instance_num): + if instance_num == 1: + return + for file in glob.glob(f'{os.path.dirname(self._nii_path)}/*.json'): os.remove(file) - for file in glob.glob(f'{self._outdir}/*.dcm'): + for file in glob.glob(f'{os.path.dirname(self._nii_path)}/*.nii'): os.remove(file) def get_num_tasks(self): self._num_tasks = len(self.snr_tasks) - -class MaskThresholdError(Exception): - pass - diff --git a/scanbuddy/proc/volreg.py b/scanbuddy/proc/volreg.py index 103cc09..94398b6 100644 --- a/scanbuddy/proc/volreg.py +++ b/scanbuddy/proc/volreg.py @@ -1,19 +1,22 @@ import os import glob import json +import time import logging import random import pydicom import subprocess import numpy as np from pubsub import pub -import time +from pathlib import Path logger = logging.getLogger(__name__) class VolReg: def __init__(self, mock=False): self._mock = mock + self._dcm1_instance_num = None + self._dcm2_instance_num = None pub.subscribe(self.listener, 'volreg') def listener(self, tasks): @@ -48,9 +51,9 @@ def run(self): arr = self.run_volreg(nii1, nii2, self.out_dir) - self.clean_dir(nii1, nii2, self.out_dir) + #self.clean_dir(nii1, nii2, self.out_dir) - logger.info(f'volreg array from registering volume {int(pydicom.dcmread(dcm2, force=True, stop_before_pixels=True).InstanceNumber)} to volume {int(pydicom.dcmread(dcm1, force=True, stop_before_pixels=True).InstanceNumber)}: {arr}') + logger.info(f'volreg array from registering volume {self._dcm2_instance_num} to volume {self._dcm1_instance_num}: {arr}') self.insert_array(arr, task_idx) @@ -63,11 +66,20 @@ def get_num_tasks(self): self.num_tasks = len(self.tasks) def create_niis(self, task_idx): + dcm1 = self.tasks[task_idx][1]['path'] nii1 = self.run_dcm2niix(dcm1, 1) + self._dcm1_instance_num = int(pydicom.dcmread(dcm1, force=True, stop_before_pixels=True).InstanceNumber) + nii1 = self.run_dcm2niix(dcm1, self._dcm1_instance_num) + if self.tasks[task_idx][1]['nii_path'] is None: + self.tasks[task_idx][1]['nii_path'] = nii1 dcm2 = self.tasks[task_idx][0]['path'] nii2 = self.run_dcm2niix(dcm2, 2) + self._dcm2_instance_num = int(pydicom.dcmread(dcm2, force=True, stop_before_pixels=True).InstanceNumber) + nii2 = self.run_dcm2niix(dcm2, self._dcm2_instance_num) + if self.tasks[task_idx][0]['nii_path'] is None: + self.tasks[task_idx][0]['nii_path'] = nii2 return nii1, nii2, dcm1, dcm2 From e81c770a5b6124198ab0e6fd56d9b566ffacb590 Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Mon, 2 Dec 2024 13:21:49 -0500 Subject: [PATCH 4/8] clean up comments --- scanbuddy/proc/__init__.py | 27 +-------------------------- 1 file changed, 1 insertion(+), 26 deletions(-) diff --git a/scanbuddy/proc/__init__.py b/scanbuddy/proc/__init__.py index 5cf4a6d..0e8fdc3 100644 --- a/scanbuddy/proc/__init__.py +++ b/scanbuddy/proc/__init__.py @@ -45,7 +45,6 @@ def listener(self, ds, path): tasks = self.check_volreg(key) - ### edits need to be made here snr_tasks = self.check_snr(key) logger.debug('publishing message to volreg topic with the following tasks') logger.debug(json.dumps(tasks, indent=2)) @@ -75,7 +74,6 @@ def listener(self, ds, path): x, y, self._z, _ = self._slice_means[key]['slice_means'].shape self._fdata_array = np.empty((x, y, self._z, self._num_vols)) - #self._slice_intensity_means = np.zeros( (self._z, math.ceil(self._num_vols // 2)) ) self._slice_intensity_means = np.zeros( (self._z, self._num_vols) ) @@ -87,21 +85,11 @@ def listener(self, ds, path): insert_position = key - 5 self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'].squeeze() - #if key == 54: - # logger.info('launching calculate and publish snr thread') - - # snr_thread = threading.Thread(target=self.calculate_and_publish_snr, args=(key,)) - # snr_thread.start() - - if key > 53 and (key % 3 == 0) and key < self._num_vols: logger.info('launching calculate and publish snr thread') snr_thread = threading.Thread(target=self.calculate_and_publish_snr, args=(key,)) snr_thread.start() - #snr_metric = round(self.calc_snr(key), 2) - #logger.info(f'running snr metric: {snr_metric}') - #pub.sendMessage('plot_snr', snr_metric=snr_metric) if key == self._num_vols: time.sleep(2) @@ -128,7 +116,6 @@ def calculate_and_publish_snr(self, key): logger.info(f'snr is a nan, decrementing mask threshold by {self._decrement}') self._mask_threshold = self._mask_threshold - self._decrement logger.info(f'new threshold: {self._mask_threshold}') - #self._numpy_4d_mask = np.zeros(self._fdata_array.shape, dtype=bool) self._slice_intensity_means = np.zeros( (self._z, self._num_vols) ) else: pub.sendMessage('plot_snr', snr_metric=snr_metric) @@ -137,10 +124,8 @@ def check_volreg(self, key): tasks = list() current = self._instances[key] - # get numerical index of key O(log n) i = self._instances.bisect_left(key) - # always register current node to left node try: left_index = max(0, i - 1) left = self._instances.values()[left_index] @@ -149,7 +134,6 @@ def check_volreg(self, key): except IndexError: pass - # if there is a right node, re-register to current node try: right_index = i + 1 right = self._instances.values()[right_index] @@ -212,15 +196,12 @@ def get_mean_slice_intensitites(self, key): logger.info(f'finding mask differences took {time.time() - start}') - #slice_intensity_means = np.zeros( (dim_z,dim_t) ) slice_voxel_counts = np.zeros( (dim_z), dtype='uint32' ) slice_size = dim_x * dim_y for slice_idx in range(dim_z): slice_voxel_counts[slice_idx] = slice_size - mask[:,:,slice_idx,0].sum() - #pdb.set_trace() - zero_columns = np.where(np.all(self._slice_intensity_means[:,:dim_t] == 0, axis=0))[0].tolist() logger.info(f'volumes being calculated: {zero_columns}') @@ -265,13 +246,7 @@ def get_mean_slice_intensitites(self, key): slice_data = data[:,:,slice_idx,volume_idx] slice_vol_mean = slice_data.mean() self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - ''' - for volume_idx in range(dim_t): - for slice_idx in differing_slices: - slice_data = data[:,:,slice_idx,volume_idx] - slice_vol_mean = slice_data.mean() - self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - ''' + return self._slice_intensity_means[:, :dim_t], slice_voxel_counts, data From 99e07c84585403c6aae43f4ea3205fe4dedfaab8 Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Mon, 2 Dec 2024 13:24:48 -0500 Subject: [PATCH 5/8] cleanup --- scanbuddy/proc/snr.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scanbuddy/proc/snr.py b/scanbuddy/proc/snr.py index bfec35d..898ab41 100644 --- a/scanbuddy/proc/snr.py +++ b/scanbuddy/proc/snr.py @@ -10,7 +10,6 @@ import pydicom import subprocess import numpy as np -#import nilearn.image import nibabel as nib from pubsub import pub import collections as c From fb221650de6b54cac3a9441de92bd8c35085f219 Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Mon, 2 Dec 2024 15:52:21 -0500 Subject: [PATCH 6/8] modified dockerfile --- Dockerfile | 56 ++++++++++++++++---------- scanbuddy/proc/__init__.py | 80 +++++++++++++++++++++++--------------- 2 files changed, 84 insertions(+), 52 deletions(-) diff --git a/Dockerfile b/Dockerfile index fb0ce0b..3f510ea 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,21 +1,37 @@ FROM rockylinux:8 -# install some base necessities -RUN dnf install -y git vim python39 python39-devel && \ - alternatives --set python /usr/bin/python3 +# Install necessary packages +RUN dnf install -y git vim tcsh libpng15 motif make curl clang zlib-devel \ + libXt-devel libXext-devel expat-devel motif-devel f2c unzip cmake gcc-c++ \ + epel-release && \ + dnf --enablerepo=powertools install -y libstdc++-static -# create a home directory +# Install Mambaforge +WORKDIR /tmp +RUN curl -L https://github.com/conda-forge/miniforge/releases/download/24.9.2-0/Mambaforge-24.9.2-0-Linux-x86_64.sh -o Mambaforge-24.9.2-0.sh && \ + chmod 755 Mambaforge-24.9.2-0.sh && \ + ./Mambaforge-24.9.2-0.sh -b -p /opt/conda && \ + rm Mambaforge-24.9.2-0.sh +ENV PATH=/opt/conda/bin:$PATH + +# Create a Conda environment with Python 3.13 +RUN conda create -n py313 python=3.13 python-freethreading -c conda-forge/label/python_rc -c conda-forge && \ + echo "conda activate py313" >> ~/.bashrc +ENV CONDA_DEFAULT_ENV=py313 +ENV PATH /opt/conda/envs/py313/bin:$PATH + +# Verify Python version +RUN python --version + +# Create a home directory RUN mkdir -p /home/scanbuddy ENV HOME=/home/scanbuddy -# compile and install 3dvolreg from AFNI (linux/amd64 and linux/arm64/v8) +# Compile and install 3dvolreg from AFNI (linux/amd64 and linux/arm64/v8) ARG AFNI_PREFIX="/sw/apps/afni" ARG AFNI_URI="https://github.com/afni/afni" WORKDIR /tmp -RUN dnf install -y epel-release && \ - dnf install -y curl tcsh python2-devel libpng15 motif && \ - dnf install -y make clang zlib-devel libXt-devel libXext-devel expat-devel motif-devel f2c && \ - git clone "${AFNI_URI}" +RUN git clone "${AFNI_URI}" WORKDIR afni/src RUN cp other_builds/Makefile.linux_ubuntu_22_64 Makefile && \ make libmri.a 3dvolreg && \ @@ -23,14 +39,12 @@ RUN cp other_builds/Makefile.linux_ubuntu_22_64 Makefile && \ mv 3dvolreg "${AFNI_PREFIX}" && \ rm -r /tmp/afni -# compile and install dcm2niix (linux/amd64 and linux/arm64/v8) +# Compile and install dcm2niix (linux/amd64 and linux/arm64/v8) ARG D2N_PREFIX="/sw/apps/dcm2niix" ARG D2N_VERSION="1.0.20220720" ARG D2N_URI="https://github.com/rordenlab/dcm2niix/archive/refs/tags/v${D2N_VERSION}.zip" -WORKDIR "/tmp" -RUN dnf install -y unzip cmake gcc-c++ && \ - dnf --enablerepo=powertools install -y libstdc++-static && \ - curl -sL "${D2N_URI}" -o "dcm2niix_src.zip" && \ +WORKDIR /tmp +RUN curl -sL "${D2N_URI}" -o "dcm2niix_src.zip" && \ unzip "dcm2niix_src.zip" && \ rm "dcm2niix_src.zip" && \ mkdir "dcm2niix-${D2N_VERSION}/build" @@ -41,23 +55,23 @@ RUN cmake .. && \ cp bin/dcm2niix "${D2N_PREFIX}" && \ rm -r "/tmp/dcm2niix-${D2N_VERSION}" -# install scanbuddy +# Install scanbuddy ARG SB_PREFIX="/sw/apps/scanbuddy" ARG SB_VERSION="v0.1.7" -RUN python3 -m venv "${SB_PREFIX}" && \ +RUN python -m venv "${SB_PREFIX}" && \ dnf install -y gcc zlib-devel libjpeg-devel python39-tkinter && \ "${SB_PREFIX}/bin/pip" install "git+https://github.com/harvard-nrg/scanbuddy.git@${SB_VERSION}" -# set up afni environment +# Set up AFNI environment ENV PATH="${AFNI_PREFIX}:${PATH}" -# set up dcm2niix environment +# Set up dcm2niix environment ENV PATH="${D2N_PREFIX}:${PATH}" -# set up scanbuddy +# Set up scanbuddy ENV TERM="xterm-256color" -# expose port +# Expose port EXPOSE 11112 -ENTRYPOINT ["/sw/apps/scanbuddy/bin/start.py"] +ENTRYPOINT ["/sw/apps/scanbuddy/bin/start.py"] \ No newline at end of file diff --git a/scanbuddy/proc/__init__.py b/scanbuddy/proc/__init__.py index 0e8fdc3..3577374 100644 --- a/scanbuddy/proc/__init__.py +++ b/scanbuddy/proc/__init__.py @@ -85,7 +85,7 @@ def listener(self, ds, path): insert_position = key - 5 self._fdata_array[:, :, :, insert_position] = self._slice_means[key]['slice_means'].squeeze() - if key > 53 and (key % 3 == 0) and key < self._num_vols: + if key > 53 and (key % 4 == 0) and key < self._num_vols: logger.info('launching calculate and publish snr thread') snr_thread = threading.Thread(target=self.calculate_and_publish_snr, args=(key,)) @@ -190,10 +190,12 @@ def get_mean_slice_intensitites(self, key): dim_t = key - 4 + ''' if key > 55: start = time.time() differing_slices = self.find_mask_differences(key) logger.info(f'finding mask differences took {time.time() - start}') + ''' slice_voxel_counts = np.zeros( (dim_z), dtype='uint32' ) @@ -221,31 +223,48 @@ def get_mean_slice_intensitites(self, key): slice_vol_mean = slice_data.mean() self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - logger.info(f'recalculating slice means at the following slices: {differing_slices}') - logger.info(f'total of {len(differing_slices)} new slices being computed') - - if differing_slices: - - if key == self._num_vols: - for volume_idx in range(dim_t): - for slice_idx in differing_slices: - slice_data = data[:,:,slice_idx,volume_idx] - slice_vol_mean = slice_data.mean() - self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - - elif key % 4 == 0: - for volume_idx in range(0, dim_t, 6): - for slice_idx in differing_slices: - slice_data = data[:,:,slice_idx,volume_idx] - slice_vol_mean = slice_data.mean() - self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean - - elif key % 3 == 0: - for volume_idx in range(3, dim_t, 6): - for slice_idx in differing_slices: - slice_data = data[:,:,slice_idx,volume_idx] - slice_vol_mean = slice_data.mean() - self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean + #logger.info(f'recalculating slice means at the following slices: {differing_slices}') + #logger.info(f'total of {len(differing_slices)} new slices being computed') + + #if differing_slices: + + if key == self._num_vols: + start = time.time() + differing_slices = self.find_mask_differences(key) + logger.info(f'finding mask differences took {time.time() - start}') + logger.info(f'recalculating slice means at the following slices: {differing_slices}') + logger.info(f'total of {len(differing_slices)} new slices being computed') + for volume_idx in range(dim_t): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean + + elif key % 6 == 0: + logger.info(f'inside the even calculation') + start = time.time() + differing_slices = self.find_mask_differences(key) + logger.info(f'finding mask differences took {time.time() - start}') + logger.info(f'recalculating slice means at the following slices: {differing_slices}') + logger.info(f'total of {len(differing_slices)} new slices being computed') + for volume_idx in range(0, dim_t, 8): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean + + elif key % 5 == 0: + logger.info(f'inside the odd calculation') + start = time.time() + differing_slices = self.find_mask_differences(key) + logger.info(f'finding mask differences took {time.time() - start}') + logger.info(f'recalculating slice means at the following slices: {differing_slices}') + logger.info(f'total of {len(differing_slices)} new slices being computed') + for volume_idx in range(5, dim_t, 8): + for slice_idx in differing_slices: + slice_data = data[:,:,slice_idx,volume_idx] + slice_vol_mean = slice_data.mean() + self._slice_intensity_means[slice_idx,volume_idx] = slice_vol_mean return self._slice_intensity_means[:, :dim_t], slice_voxel_counts, data @@ -276,17 +295,16 @@ def generate_mask(self, key): return masked_data def find_mask_differences(self, key): - num_old_vols = key - 7 - prev_mask = self._slice_means[key-3]['mask'] + num_old_vols = key - 8 + last_50 = num_old_vols - 50 + prev_mask = self._slice_means[key-4]['mask'] current_mask = self._slice_means[key]['mask'] - differences = prev_mask != current_mask[:,:,:,:num_old_vols] + differences = prev_mask[:,:,:,-50:] != current_mask[:,:,:,last_50:num_old_vols] diff_indices = np.where(differences) differing_slices = [] for index in zip(*diff_indices): if int(index[2]) not in differing_slices: differing_slices.append(int(index[2])) - #slices_to_update = differing_slices.sort() - #pdb.set_trace() return differing_slices From e72aec3cc50dacdfa0f9e04c87a55ad3ab8a7e4a Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Wed, 4 Dec 2024 13:12:28 -0500 Subject: [PATCH 7/8] updated dockerfile with python 3.13 --- Dockerfile | 70 +++++++++++++++++++++++++----------------------------- 1 file changed, 32 insertions(+), 38 deletions(-) diff --git a/Dockerfile b/Dockerfile index 3f510ea..ec41acd 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,37 +1,28 @@ -FROM rockylinux:8 +FROM rockylinux:9 -# Install necessary packages -RUN dnf install -y git vim tcsh libpng15 motif make curl clang zlib-devel \ - libXt-devel libXext-devel expat-devel motif-devel f2c unzip cmake gcc-c++ \ - epel-release && \ - dnf --enablerepo=powertools install -y libstdc++-static +# install some base necessities +RUN dnf install -y git vim +RUN dnf install -y epel-release +RUN dnf install -y conda +RUN conda create -y -n python3.13t --override-channels -c conda-forge python-freethreading +RUN conda env config vars set PYTHON_GIL=0 -n python3.13t +RUN conda run -n python3.13t --no-capture-output python3 -m pip install git+https://github.com/harvard-nrg/scanbuddy.git +RUN conda run -n python3.13t start.py --help # ENTRYPOINT ["conda", "run", "-n", "python3.13t", "start.py"] +#RUN dnf install -y git vim python3.13 && \ +# alternatives --set python /usr/bin/python3 -# Install Mambaforge -WORKDIR /tmp -RUN curl -L https://github.com/conda-forge/miniforge/releases/download/24.9.2-0/Mambaforge-24.9.2-0-Linux-x86_64.sh -o Mambaforge-24.9.2-0.sh && \ - chmod 755 Mambaforge-24.9.2-0.sh && \ - ./Mambaforge-24.9.2-0.sh -b -p /opt/conda && \ - rm Mambaforge-24.9.2-0.sh -ENV PATH=/opt/conda/bin:$PATH - -# Create a Conda environment with Python 3.13 -RUN conda create -n py313 python=3.13 python-freethreading -c conda-forge/label/python_rc -c conda-forge && \ - echo "conda activate py313" >> ~/.bashrc -ENV CONDA_DEFAULT_ENV=py313 -ENV PATH /opt/conda/envs/py313/bin:$PATH - -# Verify Python version -RUN python --version - -# Create a home directory +# create a home directory RUN mkdir -p /home/scanbuddy ENV HOME=/home/scanbuddy -# Compile and install 3dvolreg from AFNI (linux/amd64 and linux/arm64/v8) +# compile and install 3dvolreg from AFNI (linux/amd64 and linux/arm64/v8) ARG AFNI_PREFIX="/sw/apps/afni" ARG AFNI_URI="https://github.com/afni/afni" WORKDIR /tmp -RUN git clone "${AFNI_URI}" +RUN dnf install -y epel-release && \ + dnf install -y --allowerasing curl tcsh libpng15 motif && \ + dnf install -y make clang zlib-devel libXt-devel libXext-devel expat-devel motif-devel f2c && \ + git clone "${AFNI_URI}" WORKDIR afni/src RUN cp other_builds/Makefile.linux_ubuntu_22_64 Makefile && \ make libmri.a 3dvolreg && \ @@ -39,12 +30,15 @@ RUN cp other_builds/Makefile.linux_ubuntu_22_64 Makefile && \ mv 3dvolreg "${AFNI_PREFIX}" && \ rm -r /tmp/afni -# Compile and install dcm2niix (linux/amd64 and linux/arm64/v8) +# compile and install dcm2niix (linux/amd64 and linux/arm64/v8) ARG D2N_PREFIX="/sw/apps/dcm2niix" ARG D2N_VERSION="1.0.20220720" ARG D2N_URI="https://github.com/rordenlab/dcm2niix/archive/refs/tags/v${D2N_VERSION}.zip" -WORKDIR /tmp -RUN curl -sL "${D2N_URI}" -o "dcm2niix_src.zip" && \ +WORKDIR "/tmp" +RUN dnf install -y unzip cmake gcc-c++ && \ + dnf config-manager --set-enabled crb && \ + dnf install -y libstdc++-static && \ + curl -sL "${D2N_URI}" -o "dcm2niix_src.zip" && \ unzip "dcm2niix_src.zip" && \ rm "dcm2niix_src.zip" && \ mkdir "dcm2niix-${D2N_VERSION}/build" @@ -55,23 +49,23 @@ RUN cmake .. && \ cp bin/dcm2niix "${D2N_PREFIX}" && \ rm -r "/tmp/dcm2niix-${D2N_VERSION}" -# Install scanbuddy +# install scanbuddy ARG SB_PREFIX="/sw/apps/scanbuddy" -ARG SB_VERSION="v0.1.7" -RUN python -m venv "${SB_PREFIX}" && \ - dnf install -y gcc zlib-devel libjpeg-devel python39-tkinter && \ +ARG SB_VERSION="v0.1.10" +RUN python3 -m venv "${SB_PREFIX}" && \ + dnf install -y gcc zlib-devel libjpeg-devel python3-tkinter && \ "${SB_PREFIX}/bin/pip" install "git+https://github.com/harvard-nrg/scanbuddy.git@${SB_VERSION}" -# Set up AFNI environment +# set up afni environment ENV PATH="${AFNI_PREFIX}:${PATH}" -# Set up dcm2niix environment +# set up dcm2niix environment f ENV PATH="${D2N_PREFIX}:${PATH}" -# Set up scanbuddy +# set up scanbuddy ENV TERM="xterm-256color" -# Expose port +# expose port EXPOSE 11112 -ENTRYPOINT ["/sw/apps/scanbuddy/bin/start.py"] \ No newline at end of file +ENTRYPOINT ["/sw/apps/scanbuddy/bin/start.py"] From a0d8142728c4fa3c521c815ede51c4125a6e90cf Mon Sep 17 00:00:00 2001 From: Daniel Asay Date: Wed, 4 Dec 2024 13:12:48 -0500 Subject: [PATCH 8/8] updates to dash ui --- scanbuddy/view/dash.py | 75 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 74 insertions(+), 1 deletion(-) diff --git a/scanbuddy/view/dash.py b/scanbuddy/view/dash.py index 6afa1f0..c23d981 100644 --- a/scanbuddy/view/dash.py +++ b/scanbuddy/view/dash.py @@ -146,6 +146,7 @@ def init_page(self): } ) + ''' metrics_card = dbc.Card( [ dbc.CardBody( @@ -310,6 +311,7 @@ def init_page(self): ), ], style={"border": "1px solid black", "padding": "0"} + '' ) ], style={ @@ -321,11 +323,82 @@ def init_page(self): } ) + ''' + LEFT_COLUMN_WIDTH = "85%" + RIGHT_COLUMN_WIDTH = "85%" + BG_COLOR = "#e0f7fa" # Light blue background + PADDING = "10px" + BORDER_STYLE = "1px solid black" + table_header = [ + html.Thead( + html.Tr( + html.Th( + "Motion Metrics", + colSpan=2, + style={ + 'textAlign': 'center', + 'border': BORDER_STYLE + } + ) + ) + ) + ] + # Metric table rows with ids for each value + row1 = html.Tr([ + html.Td("Number of Volumes", style={'width': LEFT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}), + html.Td(id='number-of-vols', children="0", style={'width': RIGHT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}) + ]) + + row2 = html.Tr([ + html.Td("Movements > .5 mm", style={'width': LEFT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}), + html.Td(id='movements-05mm', children="0", style={'width': RIGHT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}) + ]) + + row3 = html.Tr([ + html.Td("Movements > 1 mm", style={'width': LEFT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}), + html.Td(id='movements-1mm', children="0", style={'width': RIGHT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}) + ]) + + row4 = html.Tr([ + html.Td("Max Abs Motion", style={'width': LEFT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}), + html.Td(id='max-abs-motion', children="0", style={'width': RIGHT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}) + ]) + + row5 = html.Tr([ + html.Td("SNR", style={'width': LEFT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE}), + html.Td(id='snr', children="0.0", style={'width': RIGHT_COLUMN_WIDTH, 'padding': PADDING, 'border': BORDER_STYLE, 'textAlign': 'center'}) + ]) + + table_body = [ + html.Tbody([ + row1, + row2, + row3, + row4, + row5 + ]) + ] + + metrics_table = dbc.Table( + table_header + table_body, + bordered=True, + striped=True, + style={ + 'fontSize': '3.0vh', + 'backgroundColor': BG_COLOR + } + ) + self._app.layout = html.Div([ navbar, dbc.Row( [ - dbc.Col(metrics_card, width=2, style={"marginTop": "50px"}), + dbc.Col( + metrics_table, + style={ + 'margin': '10px' + } + ), dbc.Col( [ displacements_graph,