diff --git a/README.md b/README.md index 69f2a56..a47b785 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ Essential quality control analysis of routine phantom stability MRI with single session and trend reporting -### Latest Version : 2021.7.7 +### Latest Version : 2021.7.8 ### External Dependencies FSL 6.0 or later (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FslInstallation) diff --git a/cbicqc/timeseries.py b/cbicqc/timeseries.py index 044627a..dcb352e 100644 --- a/cbicqc/timeseries.py +++ b/cbicqc/timeseries.py @@ -79,8 +79,6 @@ def detrend_timeseries(s_mean_t): :return: """ - np.savetxt("/Users/jmt/sandbox/s_mean_t.csv", s_mean_t, delimiter=",") - # Mean label ROI signal (n_labels x n_timepoints) nl, nt = s_mean_t.shape diff --git a/old-bash/cbicqa.bash b/old-bash/cbicqa.bash deleted file mode 100755 index 5ab61f5..0000000 --- a/old-bash/cbicqa.bash +++ /dev/null @@ -1,72 +0,0 @@ -#!/bin/bash -# -# Single QA analysis script -# -# AUTHOR : Mike Tyszka, Ph.D. -# PLACE : Caltech Brain Imaging Center -# DATES : 10/10/2011 JMT From scratch -# -# Copyright 2011 California Institute of Technology -# All rights reserved. - -# Localization -user_home=/Users/jmt -bin_dir=${user_home}/bin/CBICQA - -if [ $# -lt 3 ]; then - echo "-----------------" - echo "SYNTAX : cbicqa.bash " - echo " has format: YYYYMMDD-HHMM" - echo "This script is designed to be called by cbicqa_all.bash, but can be run manually" - exit -fi - -# Construct QA analysis directory name -qa_data=$1 -qa_date=$2 -qa_overwrite=$3 -qa_dir=${qa_data}/${qa_date} - -# Full paths to commands (for XGrid if used) -cmd_getdicom=${bin_dir}/cbicqa_getdicom.bash -cmd_convert=${bin_dir}/cbicqa_dicom2nifti.bash -cmd_moco=${bin_dir}/cbicqa_moco.bash -cmd_stats=${bin_dir}/cbicqa_stats.bash -cmd_report=${bin_dir}/cbicqa_report.bash - -# Check if directory already exists - if not, get data and analyze -if [ ! -d ${qa_dir} -o ${qa_overwrite} == "Y" ] -then - - echo "" - echo "----------------------------" - echo "NEW QA ANALYSIS : ${qa_date}" - echo "----------------------------" - - if [ -s ${qa_dir}/qa.nii.gz ] - then - - echo " DICOM data has already been retrieved and converted - skipping" - - else - - # Retrieve QA DICOM stack from local OsiriX - $cmd_getdicom $qa_dir $qa_date - - # Convert QA DICOM stack to 4D Nifti-1 volume - $cmd_convert $qa_dir - - fi - - # Motion correct QA series - $cmd_moco $qa_dir - - # Generate descriptive QA stats - $cmd_stats $qa_dir - - # Generate HTML report and summary files - $cmd_report $qa_dir - - echo " Complete" - -fi \ No newline at end of file diff --git a/old-bash/cbicqa_all.bash b/old-bash/cbicqa_all.bash deleted file mode 100755 index 836da92..0000000 --- a/old-bash/cbicqa_all.bash +++ /dev/null @@ -1,71 +0,0 @@ -#!/bin/bash -# -# Master QA analysis script -# -# USAGE: cbicqa_all.bash -# -# Query and retrieve new QA studies from the local OsiriX database -# Use dcmtk command line functions to generate list of new QA data and -# copy DICOM files into QA directory tree ready for analysis. -# -# AUTHOR : Mike Tyszka, Ph.D. -# PLACE : Caltech Brain Imaging Center -# DATES : 10/10/2011 JMT From scratch -# 10/16/2012 JMT Add local OsiriX Q&R to get QA study data -# 10/19/2012 JMT Add optional overwrite argument -# -# Copyright 2011-2012 California Institute of Technology -# All rights reserved. - -# QA data/website directory -qa_data="/Library/Server/Web/Data/Sites/Default/QA" - -if [ $# -lt 1 ]; then - overwrite="N" -else - overwrite=$1 -fi - -echo "CBIC QA Analysis" -echo "----------------" -echo "QA Data Directory : $qa_data" - -# Query all studies on local OsiriX with PatientID = "qa" -echo "Querying local OsiriX database for QA studies" -qa_list=( `findscu -aet QA localhost 11112 -S -k 0008,0052="STUDY" -k PatientID="qa" -k StudyDate= 2>&1 | awk '{ if ( $3 == "DA" ) { print $4 } }'` ) - -# Report length of QA study list -n_qas=${#qa_list[@]} -echo "${n_qas} QA studies found" - -# If no studies are found, likely that OsiriX isn't running or communication is failing -if [ ${n_qas} -lt 1 ] -then - echo "No QA studies found - check that OsiriX is running" - exit -fi - -# -# Loop over all QA studies returned by OsiriX -# - -for q in "${qa_list[@]}": -do - - # Strip any square braces and colons from name - qa_date=${q/[/} - qa_date=${qa_date/]/} - qa_date=${qa_date/:/} - - # Call single study QA analysis - cbicqa.bash ${qa_data} ${qa_date} ${overwrite} - -done - -# Compile summary report from all individual QA reports -# Final report HTML file is in $qa_data/qa_report.html - -cbicqa_report_all.bash $qa_data - -echo "-----------" -echo "Done" diff --git a/old-bash/cbicqa_checklims.bash b/old-bash/cbicqa_checklims.bash deleted file mode 100755 index 058a3ef..0000000 --- a/old-bash/cbicqa_checklims.bash +++ /dev/null @@ -1,20 +0,0 @@ -#!/bin/bash -# -# Return an HTML table row for a given value and acceptable limits with colored check string -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 10/16/2012 JMT From scratch -# -# Copyright 2012 California Institute of Technology -# All rights reserved - -var_name=$1 -var_value=$2 -var_lower=$3 -var_upper=$4 - -# Use awk to generate correct check string -check_str=`echo $2 $3 $4 | awk '{ if ($1 < $2 || $1 > $3) { print "

Fail

" } else { print "

Pass

" } }'` - -# Construct HTML code string to return -echo "${var_name}${var_value}${var_lower}${var_upper}${check_str}" diff --git a/old-bash/cbicqa_detrend.py b/old-bash/cbicqa_detrend.py deleted file mode 100755 index f5a7a09..0000000 --- a/old-bash/cbicqa_detrend.py +++ /dev/null @@ -1,131 +0,0 @@ -#!/opt/local/bin/python -# Detrend signal in a CBIC QA timeseries text file -# -# AUTHOR : Mike Tyszka -# PLACE : Caltech -# DATES : 04/01/2013 JMT From scratch - -import sys -import os -import numpy as np -from pylab import * -from scipy.optimize import curve_fit - -# Main function -def main(): - - # Report python version - # print(sys.version_info) - # print(sp.__version__) - - # QA root directory - qa_root = '/Library/Server/Web/Data/Sites/Default/QA' - - # Get QA date from command line args - if len(sys.argv) >= 2: - qa_dir = sys.argv[1] - else: - qa_dir = os.getcwd() - - print('Detrending ' + qa_dir) - - # Construct timeseries filename - qa_ts_file = os.path.join(qa_root,qa_dir,'qa_timeseries.txt') - if not os.path.isfile(qa_ts_file): - print(qa_ts_file + ' does not exist - exiting') - sys.exit(0) - - # Load timeseries from QA directory - print('Loading timeseries') - x = np.loadtxt(qa_ts_file) - - # Parse timeseries into column vectors - phantom_t = x[:,0] - nyquist_t = x[:,1] - noise_t = x[:,2] - - # Volume number vector - nv = len(phantom_t) - v = np.linspace(0, nv-1, nv) - - # Fit exp + const model to each timeseries - print('Fitting exponential models') - phantom_opt, phantom_cov = curve_fit(expconst, v, phantom_t) - nyquist_opt, nyquist_cov = curve_fit(expconst, v, nyquist_t) - noise_opt, noise_cov = curve_fit(expconst, v, noise_t) - - # Generate fitted curves - phantom_fit = expconst(v, phantom_opt[0], phantom_opt[1], phantom_opt[2], phantom_opt[3]) - nyquist_fit = expconst(v, nyquist_opt[0], nyquist_opt[1], nyquist_opt[2], nyquist_opt[3]) - noise_fit = expconst(v, noise_opt[0], noise_opt[1], noise_opt[2], noise_opt[3]) - - # Fit residuals - phantom_res = phantom_t - phantom_fit - nyquist_res = nyquist_t - nyquist_fit - noise_res = noise_t - noise_fit - - # Residual stats - # Assumes Gaussian white noise + sparse outlier spikes - # Use robust estimate of residual sd (MAD * 1.4826) - phantom_res_sd = np.median(np.abs(phantom_res)) * 1.4826 - nyquist_res_sd = np.median(np.abs(nyquist_res)) * 1.4826 - noise_res_sd = np.median(np.abs(noise_res)) * 1.4826 - - # Count spikes, defined as residuals more than 5 SD from zero - phantom_spikes = (np.abs(phantom_res) > 5 * phantom_res_sd).sum() - nyquist_spikes = (np.abs(nyquist_res) > 5 * nyquist_res_sd).sum() - noise_spikes = (np.abs(noise_res) > 5 * noise_res_sd).sum() - - # Append sd and spike count to fit parameter list - phantom_opt = np.append(phantom_opt, [phantom_res_sd, phantom_spikes]) - nyquist_opt = np.append(nyquist_opt,[nyquist_res_sd, nyquist_spikes]) - noise_opt = np.append(noise_opt,[noise_res_sd, noise_spikes]) - - # Generate detrended timeseries - phantom_0 = phantom_res + phantom_opt[3] - nyquist_0 = nyquist_res + nyquist_opt[3] - noise_0 = noise_res + noise_opt[3] - - # Create array with detrended timeseries in columns (for output by np.savetxt - ts_detrend = np.array([phantom_0, nyquist_0, noise_0]) - ts_detrend = ts_detrend.transpose() - - # Save detrended timeseries - print('Writing detrended timeseries') - ts_detrend_file = os.path.join(qa_root,qa_dir,'qa_timeseries_detrend.txt') - np.savetxt(ts_detrend_file, ts_detrend, delimiter=' ',fmt='%0.6f') - - # Create array with timeseries detrend parameters in columns for each model - # Row order : Exp amp, Exp tau, linear, const, residual sd, spike count - detrend_pars = np.array([phantom_opt, nyquist_opt, noise_opt]) - detrend_pars = detrend_pars.transpose() - - # Save fit parameters - print('Writing detrending parameters') - qa_detrend_pars = os.path.join(qa_root,qa_dir,'qa_detrend.pars') - np.savetxt(qa_detrend_pars, detrend_pars, delimiter=' ',fmt='%0.6f') - - # Plot figures - - fig = figure(figsize = (10,8)) - - subplot(311, title='Phantom') - plot(v, phantom_t, 'ro', v, phantom_fit, '-b', v, phantom_0, '-g') - subplot(312, title='Nyquist Ghost') - plot(v, nyquist_t, 'ro', v, nyquist_fit, '-b', v, nyquist_0, '-g') - subplot(313, title='Noise') - plot(v, noise_t, 'ro', v, noise_fit, '-b', v, noise_0, '-g') - - # Adjust vertical spacing to allow for titles - subplots_adjust(hspace = 0.5) - - # Save figure in QA directory - savefig(os.path.join(qa_root,qa_dir,'qa_detrend.png')) - -# Exponential + constant model -def expconst(t, a, tau, b, c): - return a * np.exp(-t/tau) + b * t + c - -# This is the standard boilerplate that calls the main() function. -if __name__ == '__main__': - main() diff --git a/old-bash/cbicqa_dicom2nifti.bash b/old-bash/cbicqa_dicom2nifti.bash deleted file mode 100755 index 5564572..0000000 --- a/old-bash/cbicqa_dicom2nifti.bash +++ /dev/null @@ -1,80 +0,0 @@ -#!/bin/bash -# -# Convert DICOM QA series to Nifti-1 4D file -# -# Expects ${qa_dir}/DICOM to exist and contain the QA DICOM files -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 09/16/2011 JMT From scratch -# PLACE : Caltech Brain Imaging Center - -if [ $# -lt 1 ]; then - echo "-----------------" - echo "Please provide a QA directory name" - echo "SYNTAX : cbicqa_dicom2nifti " - exit -fi - -qa_dir=$1 - -# Check for QA dir and DICOM subdirectory -if [ ! -d ${qa_dir}/DICOM ] -then - echo "*** ${qa_dir}/DICOM does not exist - exiting" - exit -fi - -echo " Converting DICOM to Nifti" - -# Get the first filename from the DICOM import subdirectory -dc_file=`ls -1 ${qa_dir}/DICOM/* | head -n 1` - -# Generate a series info text file -dicom_info=${qa_dir}/qa_info.txt - -if [ -s ${dicom_info} ]; then - - echo " DICOM info file present - skipping" - -else - - echo " Generating DICOM info file from ${dc_file}" - - # Scanner Serial Number - mri_probedicom --i ${dc_file} --t 0018 1000 >> ${dicom_info} - - # Scan date - mri_probedicom --i ${dc_file} --t 0008 0022 >> ${dicom_info} - - # Scanner Frequency - mri_probedicom --i ${dc_file} --t 0018 0084 >> ${dicom_info} - - # Repetition Time (ms) - mri_probedicom --i ${dc_file} --t 0018 0080 >> ${dicom_info} - -fi - -# Output file name -nifti_file=${qa_dir}/qa.nii.gz - -# Skip if Nifti file already present -if [ -s ${nifti_file} ]; then - - echo " QA Nifti file already exists - skipping" - -else - - # Convert DICOM stack to Nifti 4D - echo " Converting DICOM to Nifti-1 volume" - mri_convert ${dc_file} ${nifti_file} &> ${qa_dir}/qa_mri_convert.log - - # Delete the DICOM subdirectory if conversion successful (ie qa.nii.gz file exists) - if [ -s ${nifti_file} ]; then - echo " Removing temporary DICOM files" - rm -rf ${qa_dir}/DICOM - fi - - # Add number of volumes to info file - fslinfo ${nifti_file} | awk '{ if ($1 == "dim4") { print $2 } }' >> ${dicom_info} - -fi \ No newline at end of file diff --git a/old-bash/cbicqa_getdicom.bash b/old-bash/cbicqa_getdicom.bash deleted file mode 100755 index 387f98e..0000000 --- a/old-bash/cbicqa_getdicom.bash +++ /dev/null @@ -1,50 +0,0 @@ -#!/bin/bash -# -# Single QA analysis script -# -# AUTHOR : Mike Tyszka, Ph.D. -# PLACE : Caltech Brain Imaging Center -# DATES : 10/10/2011 JMT From scratch -# -# Copyright 2011 California Institute of Technology -# All rights reserved. - -# Local OsiriX AE Title -osirix_aet=evendim - -# Local Osirix host name -#osirix_hostname=evendim.caltech.edu -osirix_hostname=127.0.0.1 - -# Local Osirix port number -osirix_port=11112 - -# Local movescu AE Title -movescu_aet=QA - -# Local movescu port number -movescu_port=11113 - -# Full path to study directory for this date -qa_dir=$1 - -# Date string YYYYMMDD of current study -qa_date=$2 - -# Full debug (-d) or quiet movescu -debug=$3 - -# Temporary DICOM import directory -qa_import=${qa_dir}/DICOM -echo " Import directory : ${qa_import}" - -# Create DICOM import directory (and containing directory) -if [ ! -d ${qa_import} ] -then - mkdir -p ${qa_import} -fi - -# Get QA DICOM stack from OsiriX database -echo " Retrieving first QA study on ${qa_date} from OsiriX database" -#movescu -aet ${movescu_aet} -aem ${movescu_aet} --port ${movescu_port} -aec ${osirix_aet} -S ${debug} -k 0008,0052="STUDY" -k PatientID="qa" -k StudyDate=${qa_date} -od ${qa_import} ${osirix_hostname} ${osirix_port} -movescu --port ${movescu_port} -S ${debug} -k 0008,0052="STUDY" -k PatientID="qa" -k StudyDate=${qa_date} -od ${qa_import} ${osirix_hostname} ${osirix_port} diff --git a/old-bash/cbicqa_meansd.bash b/old-bash/cbicqa_meansd.bash deleted file mode 100755 index 69d05cf..0000000 --- a/old-bash/cbicqa_meansd.bash +++ /dev/null @@ -1,19 +0,0 @@ -#!/bin/bash -# -# Calculate column mean and sd in a text file of floats -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 10/25/2012 JMT From scratch -# -# Copyright 2012 California Institute of Technology -# All rights reserved - -# Arguments -filename=$1 -column=$2 - -awk -v c=$column ' - BEGIN { sum = 0; sum2 = 0 } - { sum += $c; sum2 += $c*$c } - END { printf "%0.5f %0.5f", sum/NR,sqrt( (sum2 - sum*sum/NR) / NR ) } -' ${filename} \ No newline at end of file diff --git a/old-bash/cbicqa_moco.bash b/old-bash/cbicqa_moco.bash deleted file mode 100755 index 04e6908..0000000 --- a/old-bash/cbicqa_moco.bash +++ /dev/null @@ -1,29 +0,0 @@ -#!/bin/bash -# -# Register all volumes over time (use mcflirt) -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 09/16/2011 JMT From scratch - -if [ $# -lt 1 ]; then - echo "Please provide a QA directory name" - echo "SYNTAX : cbicqa_moco.bash " - exit -fi - -# Key directories and files -qa_dir=$1 -qa_nifti=${qa_dir}/qa -qa_mcf=${qa_dir}/qa_mcf -qa_mcf_par=${qa_dir}/qa_mcf.par - -# Check whether MOCO has already been performed -if [ -s ${qa_mcf}.nii.gz ] && [ -s ${qa_mcf_par} ]; then - echo " Motion correction has already been done - skipping" -else - echo " Motion correcting" - # Run mcflirt on QA Nifti file - # Output file defaults to qa_mcf.nii.gz, parameters in qa_mcf.par - mcflirt -in ${qa_nifti} -out ${qa_mcf} -refvol 0 -plots - -fi \ No newline at end of file diff --git a/old-bash/cbicqa_plot.bash b/old-bash/cbicqa_plot.bash deleted file mode 100755 index 835320c..0000000 --- a/old-bash/cbicqa_plot.bash +++ /dev/null @@ -1,41 +0,0 @@ -#!/bin/bash -# -# Create plot image of column in qa.dat -# -# AUTHOR : Mike Tyszka, Ph.D. -# PLACE : Caltech Brain Imaging Center -# DATES : 10/10/2011 JMT From scratch -# -# Copyright 2011 California Institute of Technology -# All rights reserved. - -if [ $# -lt 3 ]; then - echo "-----------------" - echo "SYNTAX : cbicqa_pot.bash " - exit -fi - -data=$1 -col=$2 -out=$3 -ymin=$4 -ymax=$5 - -# Plot graph as PBM image -gnuplot << ENDPLOT -set output '${out}.pbm' -set term pbm size 800,200 color -set xdata time -set yrange [${ymin}:${ymax}] -set timefmt '%Y-%m-%d' # YYYY-MM-DD formate dates -set format y '%10.3f' -set datafile separator ',' # CSV File -# Plot values from CSV file, starting at second line (skip header row) -plot '${data}' using 2:${col} every ::2 notitle with points lt 3 pt 6 ps 3 -ENDPLOT - -# Convert to PNG -convert ${out}.pbm ${out}.png - -# Remove PBM image -rm -rf ${out}.pbm diff --git a/old-bash/cbicqa_report.bash b/old-bash/cbicqa_report.bash deleted file mode 100755 index 5a8f4bc..0000000 --- a/old-bash/cbicqa_report.bash +++ /dev/null @@ -1,198 +0,0 @@ -#!/bin/bash -# -# Generate HTML report from analysis results in a QA directory -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 09/16/2011 JMT From scratch -# -# Copyright 2012 California Institute of Technology -# All rights reserved - -if [ $# -lt 1 ]; then - echo "Please provide a QA directory name" - echo "SYNTAX : cbicqa_report.bash " - exit -fi - -# Directories -qa_dir=$1 - -echo " Generating report" - -# HTML output file -qa_report=${qa_dir}/qa_report.html - -# Info summary file -qa_summary=${qa_dir}/qa_summary.txt - -# Text files -qa_timeseries=${qa_dir}/qa_timeseries.txt - -qa_signal_volume=${qa_dir}/qa_signal_volume.txt -qa_signal_mean=${qa_dir}/qa_signal_mean.txt -qa_noise_mean=${qa_dir}/qa_noise_mean.txt -qa_noise_sd=${qa_dir}/qa_noise_sd.txt -qa_snr=${qa_dir}/qa_snr.txt - -qa_mcf_par=${qa_dir}/qa_mcf.par - -# Read scan info from file -info=$(<${qa_dir}/qa_info.txt) -set -- $info -scanner_serno=$1 -acq_date=$2 -scanner_freq=$3 -TR_ms=$4 -num_volumes=$5 - -# Analysis date and time -analysis_date=`date "+%Y%m%d-%H%M%S"` - -# -# Summary statistics -# - -# Timeseries file has three columns for mean signal within: -# 1. Phantom -# 2. Nyquist ghost -# 3. Noise (non-phantom, non-ghost) - -# Temporal mean and sd of basic timeseries -tmeansd_phantom=`cbicqa_meansd.bash ${qa_timeseries} 1` -tmeansd_nyquist=`cbicqa_meansd.bash ${qa_timeseries} 2` -tmeansd_noise=`cbicqa_meansd.bash ${qa_timeseries} 3` - -# Extract temporal means -tmean_phantom=`echo $tmeansd_phantom | awk '{ printf "%0.2f", $1 }'` -tmean_nyquist=`echo $tmeansd_nyquist | awk '{ printf "%0.2f", $1 }'` -tmean_noise=`echo $tmeansd_noise | awk '{ printf "%0.2f", $1 }'` - -# Extract temporal SDs -tsd_phantom=`echo $tmeansd_phantom | awk '{ printf "%0.2f", $2 }'` -tsd_nyquist=`echo $tmeansd_nyquist | awk '{ printf "%0.2f", $2 }'` -tsd_noise=`echo $tmeansd_noise | awk '{ printf "%0.2f", $2 }'` - -# Derived measures from basic timeseries: SNR and drift -tmean_snr=`awk '{sum+=$1/$3} END {printf "%0.1f",sum/NR}' ${qa_timeseries}` -sig_drift_perc=`awk 'BEGIN {max=0; min=9999} $1>max {max=$1} $1amax {amax=sqrt($4*$4)} END { printf "%.1f", amax * 1000 }' ${qa_mcf_par}` -max_ady=`awk 'BEGIN {amax=0} sqrt($5*$5)>amax {amax=sqrt($5*$5)} END { printf "%.1f", amax * 1000 }' ${qa_mcf_par}` -max_adz=`awk 'BEGIN {amax=0} sqrt($6*$6)>amax {amax=sqrt($6*$6)} END { printf "%.1f", amax * 1000 }' ${qa_mcf_par}` - -# Create and clear report HTML file -echo "" > $qa_report - -# -# Construct single QA report page -# - -# Write HTML header to report file -echo " - - - -" >> $qa_report - -# Add QA acquisition information to HTML page -echo "

CBIC Quality Assurance Report

- - - - - - - -" >> $qa_report - -echo `cbicqa_checklims.bash "Mean Signal" ${tmean_phantom} 550 1000` >> $qa_report -echo `cbicqa_checklims.bash "Mean Nyquist" ${tmean_nyquist} 0 50` >> $qa_report -echo `cbicqa_checklims.bash "Mean Noise" ${tmean_noise} 0 30` >> $qa_report -echo `cbicqa_checklims.bash "SD Signal" ${tsd_phantom} 0 1` >> $qa_report -echo `cbicqa_checklims.bash "SD Nyquist" ${tsd_nyquist} 0 1` >> $qa_report -echo `cbicqa_checklims.bash "SD Noise" ${tsd_noise} 0 0.1` >> $qa_report -echo `cbicqa_checklims.bash "Mean SNR" ${tmean_snr} 25 100` >> $qa_report -echo `cbicqa_checklims.bash "Signal Drift (%)" ${sig_drift_perc} 0 10` >> $qa_report -echo `cbicqa_checklims.bash "CofM x (mm)" ${com_mean_x} -20 20` >> $qa_report -echo `cbicqa_checklims.bash "CofM y (mm)" ${com_mean_y} -20 20` >> $qa_report -echo `cbicqa_checklims.bash "CofM z (mm)" ${com_mean_z} -20 20` >> $qa_report -echo `cbicqa_checklims.bash "Max |dx| (um)" ${max_adx} 0 1000` >> $qa_report -echo `cbicqa_checklims.bash "Max |dy| (um)" ${max_ady} 0 1000` >> $qa_report -echo `cbicqa_checklims.bash "Max |dz| (um)" ${max_adz} 0 1000` >> $qa_report - -echo "
Scanner Serial No $scanner_serno
Acquisition Date $acq_date
Analysis Date $analysis_date
Frequency ${scanner_freq} MHz
Volumes ${num_volumes}
TR ${TR_ms} ms
" >> $qa_report - -# Create orthogonal slice views of mean and sd images -scale_factor=2 -slicer ${qa_dir}/qa_mean -s ${scale_factor} -a ${qa_dir}/qa_mean_ortho.png -slicer ${qa_dir}/qa_sd -s ${scale_factor} -a ${qa_dir}/qa_sd_ortho.png -slicer ${qa_dir}/qa_mask -s ${scale_factor} -a ${qa_dir}/qa_mask_ortho.png - -# Insert mean and sd orthoslice images into HTML page -echo "

Temporal Mean Image


" >> $qa_report -echo "

Temporal SD Image


" >> $qa_report -echo "

Region Mask


" >> $qa_report - -# Create graph images of mean regional signal timeseries -# TODO -fsl_tsplot -i ${qa_dir}/qa_timeseries.txt --start=1 --finish=1 -a Phantom -o ${qa_dir}/qa_phantom_ts.png -fsl_tsplot -i ${qa_dir}/qa_timeseries.txt --start=2 --finish=2 -a Nyquist -o ${qa_dir}/qa_nyquist_ts.png -fsl_tsplot -i ${qa_dir}/qa_timeseries.txt --start=3 --finish=3 -a Noise -o ${qa_dir}/qa_noise_ts.png - -# Add graphs to page -echo "

Signal and Noise

-
-
-
" >> $qa_report - -# Create translation parameter graph -fsl_tsplot -i ${qa_dir}/qa_mcf.par -t "Translation (mm)" -a x,y,z --start=4 --finish=6 -o ${qa_dir}/qa_trans.png - -# Add motion parameter graphs to page -echo "

Apparent rigid body motion

-
" >> $qa_report - -# Close out HTML report -echo "" >> $qa_report - -# Write QA summary to text file for gobal reporting -echo "ScannerSerNo $scanner_serno -acq_date $acq_date -analysis_date $analysis_date -scanner_freq $scanner_freq -tmean_phantom $tmean_phantom -tmean_nyquist $tmean_nyquist -tmean_noise $tmean_noise -tsd_phantom $tsd_phantom -tsd_nyquist $tsd_nyquist -tsd_noise $tsd_noise -tmean_snr $tmean_snr -sig_drift_perc $sig_drift_perc -com_mean_x $com_mean_x -com_mean_y $com_mean_y -com_mean_z $com_mean_z -max_adx $max_adx -max_ady $max_ady -max_adz $max_adz" > $qa_summary diff --git a/old-bash/cbicqa_report.py b/old-bash/cbicqa_report.py deleted file mode 100755 index f5a7a09..0000000 --- a/old-bash/cbicqa_report.py +++ /dev/null @@ -1,131 +0,0 @@ -#!/opt/local/bin/python -# Detrend signal in a CBIC QA timeseries text file -# -# AUTHOR : Mike Tyszka -# PLACE : Caltech -# DATES : 04/01/2013 JMT From scratch - -import sys -import os -import numpy as np -from pylab import * -from scipy.optimize import curve_fit - -# Main function -def main(): - - # Report python version - # print(sys.version_info) - # print(sp.__version__) - - # QA root directory - qa_root = '/Library/Server/Web/Data/Sites/Default/QA' - - # Get QA date from command line args - if len(sys.argv) >= 2: - qa_dir = sys.argv[1] - else: - qa_dir = os.getcwd() - - print('Detrending ' + qa_dir) - - # Construct timeseries filename - qa_ts_file = os.path.join(qa_root,qa_dir,'qa_timeseries.txt') - if not os.path.isfile(qa_ts_file): - print(qa_ts_file + ' does not exist - exiting') - sys.exit(0) - - # Load timeseries from QA directory - print('Loading timeseries') - x = np.loadtxt(qa_ts_file) - - # Parse timeseries into column vectors - phantom_t = x[:,0] - nyquist_t = x[:,1] - noise_t = x[:,2] - - # Volume number vector - nv = len(phantom_t) - v = np.linspace(0, nv-1, nv) - - # Fit exp + const model to each timeseries - print('Fitting exponential models') - phantom_opt, phantom_cov = curve_fit(expconst, v, phantom_t) - nyquist_opt, nyquist_cov = curve_fit(expconst, v, nyquist_t) - noise_opt, noise_cov = curve_fit(expconst, v, noise_t) - - # Generate fitted curves - phantom_fit = expconst(v, phantom_opt[0], phantom_opt[1], phantom_opt[2], phantom_opt[3]) - nyquist_fit = expconst(v, nyquist_opt[0], nyquist_opt[1], nyquist_opt[2], nyquist_opt[3]) - noise_fit = expconst(v, noise_opt[0], noise_opt[1], noise_opt[2], noise_opt[3]) - - # Fit residuals - phantom_res = phantom_t - phantom_fit - nyquist_res = nyquist_t - nyquist_fit - noise_res = noise_t - noise_fit - - # Residual stats - # Assumes Gaussian white noise + sparse outlier spikes - # Use robust estimate of residual sd (MAD * 1.4826) - phantom_res_sd = np.median(np.abs(phantom_res)) * 1.4826 - nyquist_res_sd = np.median(np.abs(nyquist_res)) * 1.4826 - noise_res_sd = np.median(np.abs(noise_res)) * 1.4826 - - # Count spikes, defined as residuals more than 5 SD from zero - phantom_spikes = (np.abs(phantom_res) > 5 * phantom_res_sd).sum() - nyquist_spikes = (np.abs(nyquist_res) > 5 * nyquist_res_sd).sum() - noise_spikes = (np.abs(noise_res) > 5 * noise_res_sd).sum() - - # Append sd and spike count to fit parameter list - phantom_opt = np.append(phantom_opt, [phantom_res_sd, phantom_spikes]) - nyquist_opt = np.append(nyquist_opt,[nyquist_res_sd, nyquist_spikes]) - noise_opt = np.append(noise_opt,[noise_res_sd, noise_spikes]) - - # Generate detrended timeseries - phantom_0 = phantom_res + phantom_opt[3] - nyquist_0 = nyquist_res + nyquist_opt[3] - noise_0 = noise_res + noise_opt[3] - - # Create array with detrended timeseries in columns (for output by np.savetxt - ts_detrend = np.array([phantom_0, nyquist_0, noise_0]) - ts_detrend = ts_detrend.transpose() - - # Save detrended timeseries - print('Writing detrended timeseries') - ts_detrend_file = os.path.join(qa_root,qa_dir,'qa_timeseries_detrend.txt') - np.savetxt(ts_detrend_file, ts_detrend, delimiter=' ',fmt='%0.6f') - - # Create array with timeseries detrend parameters in columns for each model - # Row order : Exp amp, Exp tau, linear, const, residual sd, spike count - detrend_pars = np.array([phantom_opt, nyquist_opt, noise_opt]) - detrend_pars = detrend_pars.transpose() - - # Save fit parameters - print('Writing detrending parameters') - qa_detrend_pars = os.path.join(qa_root,qa_dir,'qa_detrend.pars') - np.savetxt(qa_detrend_pars, detrend_pars, delimiter=' ',fmt='%0.6f') - - # Plot figures - - fig = figure(figsize = (10,8)) - - subplot(311, title='Phantom') - plot(v, phantom_t, 'ro', v, phantom_fit, '-b', v, phantom_0, '-g') - subplot(312, title='Nyquist Ghost') - plot(v, nyquist_t, 'ro', v, nyquist_fit, '-b', v, nyquist_0, '-g') - subplot(313, title='Noise') - plot(v, noise_t, 'ro', v, noise_fit, '-b', v, noise_0, '-g') - - # Adjust vertical spacing to allow for titles - subplots_adjust(hspace = 0.5) - - # Save figure in QA directory - savefig(os.path.join(qa_root,qa_dir,'qa_detrend.png')) - -# Exponential + constant model -def expconst(t, a, tau, b, c): - return a * np.exp(-t/tau) + b * t + c - -# This is the standard boilerplate that calls the main() function. -if __name__ == '__main__': - main() diff --git a/old-bash/cbicqa_report_all.bash b/old-bash/cbicqa_report_all.bash deleted file mode 100755 index 81ce24b..0000000 --- a/old-bash/cbicqa_report_all.bash +++ /dev/null @@ -1,148 +0,0 @@ -#!/bin/bash -# -# Generate HTML report from analysis results in a QA directory -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 09/16/2011 JMT From scratch - -echo "" -echo "-------------------------" -echo "GENERATING SUMMARY REPORT" -echo "-------------------------" - -if [ $# -lt 1 ]; then - echo "Please provide a QA data directory name" - echo "SYNTAX : cbicqa_report_all.bash " - exit -fi - -# QA Data Directory -qa_data=$1 - -# QA Web Site Directory -qa_website=/Library/Server/Web/Data/Sites/Default/QA - -# Create website directory if necessary -if [ ! -d ${qa_website} ] -then - echo " Creating QA Website Directory" - mkdir -p ${qa_website} -fi - -# Output report HTML file -qa_report=${qa_website}/index.html - -# Output aggregated summary file in CSV format -qa_csv=${qa_website}/qa.csv - -# Clear summary QA data file -rm -rf $qa_csv - -# -# Create HTML report page for all QA series -# - -# Clear and create new report HTML file -rm -rf $qa_report - -echo "Creating HTML report in ${qa_report}" - -# Write HTML header to report file -echo " - - - -" >> $qa_report - -echo "

CBIC Quality Assurance

" >> $qa_report - -# Add timestamp -timestamp=`date` -echo "Autogenerated : $timestamp" >> $qa_report - -# Create QA table -echo "
" >> $qa_report - -# Table header row -echo " and and /} - #tmp=${tmp//<\/tr>/} - #tmp=${tmp//
Date SNR Drift % Nyquist tMean Noise tMean Nyquist tSD Noise tSD Report Status" >> $qa_report - -# Write column headers to CSV file -echo "Date,SNR,Drift,Nyquist_tMean,Noise_tMean,Nyquist_tSD,Noise_tSD," >> $qa_csv - -# Loop over all QA series in the data directory -# QA directories have form YYYYMMDD - 8 characters - -# Collect directory names in reverse lexical order (ie most recent to oldest) -qa_dirs=`ls -1dr ${qa_data}/[0-2][0-9]*` - -for qa_dir in ${qa_dirs} -do - - if [ -d $qa_dir ]; then - - # Study (directory) name - study_name=`basename ${qa_dir}` - - # Summary file for all QA studies analyzed - summary=${qa_dir}/qa_summary.txt - - # Parse summary fields - acq_date=`awk '$1 == "acq_date" {print $2}' $summary` - tmean_phantom=`awk '$1 == "tmean_phantom" {print $2}' $summary` - tmean_nyquist=`awk '$1 == "tmean_nyquist" {print $2}' $summary` - tmean_noise=`awk '$1 == "tmean_noise" {print $2}' $summary` - tsd_nyquist=`awk '$1 == "tsd_nyquist" {print $2}' $summary` - tsd_noise=`awk '$1 == "tsd_noise" {print $2}' $summary` - tmean_snr=`awk '$1 == "tmean_snr" {print $2}' $summary` - sig_drift_perc=`awk '$1 == "sig_drift_perc" {print $2}' $summary` - - # Convert acquisition date to YYYY-MM-DD format (compatible with GNUPLOT and Excel) - YYYY=${acq_date:0:4} - MM=${acq_date:4:2} - DD=${acq_date:6:2} - acq_date="${YYYY}-${MM}-${DD}" - - # Add table row - echo "
${acq_date} ${tmean_snr} ${sig_drift_perc} ${tmean_nyquist} ${tmean_noise} ${tsd_nyquist} ${tsd_noise}" >> $qa_report - - # Add link to daily QA report - daily_report_url="${study_name}/qa_report.html" - echo "Report" >> $qa_report - - # Get QA fail rows (if any) from QA report page - fails=`fgrep Fail ${qa_dir}/qa_report.html` - - # Add colored QA status flag - if [ -z "${fails}" ] - then - - echo "

Pass

" >> $qa_report - - else - - # Remove all
- #tmp=${fails//
/} - - # Write fail cell to table - echo "

Check

" >> $qa_report - - fi - - # Add line to CSV file - echo "$acq_date,$tmean_snr,$sig_drift_perc,$tmean_nyquist,$tmean_noise,$tsd_nyquist,$tsd_noise," >> $qa_csv - - fi - -done - -# Close out HTML report -echo "
" >> $qa_report diff --git a/old-bash/cbicqa_stats.bash b/old-bash/cbicqa_stats.bash deleted file mode 100755 index a3b8935..0000000 --- a/old-bash/cbicqa_stats.bash +++ /dev/null @@ -1,124 +0,0 @@ -#!/bin/bash -# -# Calculate important QA stats for phantom QA dataset -# -# AUTHOR : Mike Tyszka, Ph.D. -# DATES : 09/16/2011 JMT From scratch -# 10/03/2012 JMT Replace corner airspace mask with systematic noise estimate in artifact mask -# 10/12/2012 JMT Use signal volume for estimating artifact volume given known phantom volume -# -# Copyright 2011-2012 California Institute of Technology. -# All rights reserved. - -# CBIC FBIRN phantom volume in ml -phantom_volume_ml=2500 - -# Check arguments -if [ $# -lt 1 ]; then - echo "Please provide a QA directory name" - echo "SYNTAX : cbicqa_stats.bash " - exit -fi - -# Directory containing QA study -qa_dir=$1 - -echo " Generating descriptive QA statistics" - -# Image filenames -qa_mcf=${qa_dir}/qa_mcf -qa_mean=${qa_dir}/qa_mean -qa_sd=${qa_dir}/qa_sd -qa_mask=${qa_dir}/qa_mask - -# Timeseries text file -qa_timeseries=${qa_dir}/qa_timeseries.txt - -# Temporal mean of registered images -if [ -s ${qa_mean}.nii.gz ] -then - echo " Mean image exists - skipping" -else - echo " Calculating temporal mean image" - fslmaths $qa_mcf -Tmean $qa_mean -fi - -# Temporal SD of registered images -if [ -s ${qa_sd}.nii.gz ] -then - echo " SD image exists - skipping" -else - echo " Calculating temporal SD image" - fslmaths $qa_mcf -Tstd $qa_sd -fi - -# Create regional mask for phantom, Nyquist ghost and noise -if [ -s ${qa_mask}.nii.gz ] -then - echo " Mask image exists - skipping" -else - - # Temporary Nifti files - tmp_signal=${qa_dir}/tmp_signal - tmp_signal_dil=${qa_dir}/tmp_signal_dil - tmp_phantom=${qa_dir}/tmp_phantom - tmp_nyquist=${qa_dir}/tmp_nyquist - tmp_upper=${qa_dir}/tmp_upper - tmp_lower=${qa_dir}/tmp_lower - tmp_noise=${qa_dir}/tmp_noise - - # Signal threshold for basic segmentation - signal_threshold=`fslstats $qa_mean -p 99 | awk '{ print $1 * 0.1 }'` - - # Signal mask - echo " Creating signal mask (threshold = ${signal_threshold})" - fslmaths ${qa_mean} -thr ${signal_threshold} -bin ${tmp_signal} - - # Erode signal mask (6 mm radius sphere) to create phantom mask - echo " Creating phantom mask" - fslmaths ${tmp_signal} -kernel sphere 6.0 -ero ${tmp_phantom} - - # Dilate phantom mask (6 mm radius sphere) *twice* to create dilated signal mask - echo " Creating dilated signal mask" - fslmaths ${tmp_phantom} -kernel sphere 6.0 -dilF -dilF ${tmp_signal_dil} - - echo " Creating Nyquist mask" - # Extract upper and lower halves of dilated volume mask in Y dimension (PE) - fslroi ${tmp_signal_dil} ${tmp_lower} 0 -1 0 32 0 -1 0 -1 - fslroi ${tmp_signal_dil} ${tmp_upper} 0 -1 32 -1 0 -1 0 -1 - - # Create shifted (Nyqusit) mask from swapped upper and lower masks - fslmerge -y ${tmp_nyquist} ${tmp_upper} ${tmp_lower} - - # Correct y offset in sform matrix - sform=`fslorient -getsform ${tmp_signal}` - fslorient -setsform ${sform} ${tmp_nyquist} - - # XOR Nyquist and dilated signal masks - fslmaths ${tmp_nyquist} -mul ${tmp_signal_dil} -mul -1.0 -add ${tmp_nyquist} ${tmp_nyquist} - - echo " Creating noise mask" - # Create noise mask by subtracting Nyquist mask from NOT dilated signal mask - fslmaths ${tmp_signal_dil} -binv -sub ${tmp_nyquist} ${tmp_noise} - - # Finally merge all three masks into an indexed file - # Phantom = 1 - # Nyquist = 2 - # Noise = 3 - fslmaths ${tmp_nyquist} -mul 2 ${tmp_nyquist} - fslmaths ${tmp_noise} -mul 3 ${tmp_noise} - fslmaths ${tmp_phantom} -add ${tmp_nyquist} -add ${tmp_noise} ${qa_mask} - - # Clean up temporary images - rm -rf ${qa_dir}/tmp*.* - -fi - -# Extract time-series stats within each ROI -if [ -s ${qa_timeseries} ]; then - echo " Signal timecourses exist - skipping" -else - echo " Extracting mean signal timecourses for each region" - # Timecourse of mean signal within each mask label - fslmeants -i ${qa_mcf} -o ${qa_timeseries} --label=${qa_mask} -fi \ No newline at end of file diff --git a/setup.py b/setup.py index c7b7bdf..955c4cc 100644 --- a/setup.py +++ b/setup.py @@ -45,7 +45,7 @@ # For a discussion on single-sourcing the version across setup.py and the # project code, see # https://packaging.python.org/en/latest/single_source_version.html - version='2021.7.7', # Required + version='2021.7.8', # Required # This is a one-line description or tagline of what your project does. This # corresponds to the "Summary" metadata field: