Merge pull request #3 from NPACore/shim_and_station

Shim and station

acq2sqlite.py

@@ -7,6 +7,7 @@
 import re
 import sqlite3
 import sys
+from datetime import datetime, timedelta
 from typing import Optional
 
 from dcmmeta2tsv import NULLVAL, TagValues
@@ -121,7 +122,7 @@ def __init__(self, sql=None):
         #  only add if not already in the DB
         acq_uniq_col = set(self.all_columns) - set(self.CONSTS) - set(["filename"])
         assert acq_uniq_col == set(
-            ["AcqTime", "AcqDate", "SeriesNumber", "SubID", "Operator"]
+            ["AcqTime", "AcqDate", "SeriesNumber", "SubID", "Operator", "Shims", "Station"]
         )
         # TODO: include station?
 
@@ -269,6 +270,24 @@ def get_template(self, pname: str, seqname: str) -> sqlite3.Row:
         logging.debug("found template: %s", res)
         return res
 
+    def find_acquisitions_since(self, since_date: Optional[str] = None):
+        """
+        Retrieve all acquisitions with AcqDate greater than the specified date.
+
+        :param since_date: Date string in 'YYYY-MM-DD' format; defaults to yesterday if None.
+        :return: List of acquisition rows with AcqDate > since_date.
+        """
+
+        # Default to yesterday if since_date is None
+        if since_date is None:
+            since_date = (datetime.now() - timedelta(days=1)).strftime("%Y-%m-%d")
+
+        query = "select * from acq where AcqDate > ?"
+        logging.info("Finding acquisitions since %s", since_date)
+
+        cur = self.sql.execute(query, (since_date,))
+        return cur.fetchall()
+
 
 def have_pipe_data():
     return os.isatty(sys.stdout.fileno())

compliance_check_html.py

@@ -0,0 +1,49 @@
+from template_checker import TemplateChecker, CheckResult
+from jinja2 import Template
+
+def load_template(template_path: str) -> Template:
+    """
+    Load an HTML template from the template.html file
+
+    :param template_path: Path to the HTML template file.
+    :returns: A Jinja2 Template object.
+    """
+    with open(template_path, 'r') as file:
+        template_content = file.read()
+    return Template(template_content)
+
+def generate_html_report(check_result: CheckResult, template_path: str) -> str:
+    """
+    Generate an HTML report of DICOM header comparison, highlighting mismatches.
+
+    :param check_result: Output from the check_header function, containing the comparison results.
+    :returns: HTML string with results formatted using a Jinja2 template.
+    """
+
+    # Headers to check
+    headers_to_check = ["Project", "SequenceName", "TR", "TE", "FA", "iPAT", "Comments"]
+
+    # Initialize the rows list
+    rows = []
+
+    # Add rows for each header parameter
+    for header in headers_to_check:
+        expected_value = check_result["template"].get(header, "N/A")
+        actual_value = check_result["input"].get(header, "N/A")
+        class_name = "mismatch" if header in check_result["errors"] else "match"
+        rows.append({
+            "header": header, 
+            "expected_value": expected_value,
+            "actual_value": actual_value,
+            "class_name": class_name
+        })
+
+    # Load the template from the file
+    template = load_template(template_path)
+
+    # Render the template with the data
+    html = template.render(rows=rows)
+
+    return html
+
+

copyDicom.sh

@@ -0,0 +1,25 @@
+#!/usr/bin/bash
+#Copy the first two DICOM files from rest, task, and dwi folders
+#
+output_dir="$HOME/src/mrrc-hdr-qa/dicoms/"
+
+mkdir -p "$output_dir"
+
+copy_first_two() {
+	local sessions=("$@")
+	for session in "${sessions[@]}"; do
+		# copy the first two dicoms in the session
+		files=($(ls "$session" | head -n2))
+		for file in "${files[@]}"; do
+			cp "$session/$file" "$output_dir/$(basename "$file")"
+		done
+	done
+}
+
+rest_sessions=($(printf "%s\n" /Volumes/Hera/Raw/MRprojects/Habit/2*/1*/Resting-state_ME_4*/ | head -n2))
+task_sessions=($(printf "%s\n" /Volumes/Hera/Raw/MRprojects/Habit/2*/1*/HabitTask_704*/ | head -n2))
+dwi_sessions=($(printf "%s\n" /Volumes/Hera/Raw/MRprojects/Habit/2*/1*/dMRI_b0_AP_1*/ | head -n2))
+
+copy_first_two "${rest_sessions[@]}"
+copy_first_two "${task_sessions[@]}"
+copy_first_two "${dwi_sessions[@]}"

dcm2nii_check.bash

@@ -1,5 +1,9 @@
 # use dcm2niix to make bids sidecar json ('-b o' is json only, no nifti) files from dicom header
+<<<<<<< HEAD
 # want to check acq dir there in addition to dicom tag in case we're looking at the wrong one
+=======
+# want to check acq dir there in addition to dicom tag in case we're looking at the wrong one
+>>>>>>> 48cf3dc (first commit for audit.py)
 echo /Volumes/Hera/Raw/MRprojects/Habit/20*-*/1*_2*/dMRI_*/ |
     head -n 15 |
     xargs -I{} -n1 dcm2niix -o example_jsons/ -b o {}

dcmmeta2tsv.bash

@@ -2,14 +2,16 @@
 # quick pass at building minimal text database of dicom headers
 # 20240907 WF - init
 #
-export TAG_ARGS=$(cut -f2 taglist.txt|grep -v csa\ header | sed '1d;/#/d;s/^/-tag /;'|paste -sd' ')
+export TAG_ARGS=$(cut -f2 taglist.txt|grep -Pv 'csa\ header|ASCCOV' | sed '1d;/#/d;s/^/-tag /;'|paste -sd' ')
 dcmmeta2tsv(){
  declare -g TAG_ARGS
  #echo "# $1" >&2
  gdcmdump -dC "$1" |
    perl -ne 'BEGIN{%a=(Phase=>"NA", ucPAT=>"NA")}
    $a{substr($1,0,5)} = $2 if m/(PhaseEncodingDirectionPositive.*Data..|ucPATMode\s+=\s+)(\d+)/;
    END {print join("\t", @a{qw/Phase ucPAT/}), "\t"}'
+ # WARNING  as of 20241119 dcmmeta2tsv.py diverged:
+ # WARNING  ASCCOV shim values not extracted here
  dicom_hinfo -sepstr $'\t' -last -full_entry $TAG_ARGS "$@"
 }
 

dcmmeta2tsv.py

@@ -7,20 +7,25 @@
 import re
 import sys
 import warnings
-from typing import TypedDict
+from typing import Optional, TypedDict
 
 import pydicom
 
 with warnings.catch_warnings():
     warnings.simplefilter("ignore")
     # UserWarning: The DICOM readers are highly experimental...
-    import nibabel.nicom.csareader as csareader
+    from nibabel.nicom import csareader
 
 logging.basicConfig(level=os.environ.get("LOGLEVEL", logging.INFO))
 
-#: object that has obj.value for when a dicom tag does not exist
-#: using 'null' to match AFNI's dicom_hinfo
-NULLVAL = type("", (object,), {"value": "null"})()
+
+class NULLVAL:
+    """Container to imitate ``pydicom.dcmread``.
+    object that has ``obj.value`` for when a dicom tag does not exist.
+    Using "null" to match AFNI's dicom_hinfo missing text"""
+
+    value: str = "null"
+
 
 TagTuple = tuple[int, int]
 
@@ -67,6 +72,8 @@ def read_known_tags(tagfile="taglist.txt") -> TagDicts:
         if re.search("^[0-9]{4},", tags[i]["tag"]):
             tags[i]["tag"] = tagpair_to_hex(tags[i]["tag"])
             tags[i]["loc"] = "header"
+        elif tags[i]["name"] == "shims":
+            tags[i]["loc"] = "asccov"
         else:
             tags[i]["loc"] = "csa"
 
@@ -93,11 +100,54 @@ def csa_fetch(csa_tr: dict, itemname: str) -> str:
     return val
 
 
-def read_csa(csa) -> list[str]:
+def read_shims(csa_s: Optional[dict]) -> list:
+    """
+    read current shim parameters from ASCCONV
+    from 0x0029,0x1020 CSA Series Header Info
+
+    csa_s = dcmmeta2tsv.read_csa(dcm.get(())
+
+    CHM maltab code concats
+      sAdjData.uiAdjShimMode
+      sGRADSPEC.asGPAData[0].lOffset{X,Y,Z}
+      sGRADSPEC.alShimCurrent[0:4]
+      sTXSPEC.asNucleusInfo[0].lFrequency
+    >>> csa_s = pydicom.dcmread('example_dicoms/RewardedAnti_good.dcm').get((0x0029, 0x1020))
+    >>> ",".join(read_shims(read_csa(csa_s)))
+    '1174,-2475,4575,531,-20,59,54,-8,123160323,4'
+    >>> read_shims(None)
+    [None]*10
+    """
+
+    if csa_s is None:
+        csa_s = {}
+    try:
+        asccov = csa_s["tags"]["MrPhoenixProtocol"]["items"][0]
+    except KeyError:
+        return [NULLVAL.value] * 10
+
+    key = "|".join(
+        [
+            "sAdjData.uiAdjShimMode",
+            "sGRADSPEC.asGPAData\\[0\\].lOffset[XYZ]",
+            "sGRADSPEC.alShimCurrent\\[[0-4]\\]",
+            "sTXSPEC.asNucleusInfo\\[0\\].lFrequency",
+        ]
+    )
+
+    # keys are like
+    #   sGRADSPEC.asGPAData[0].lOffsetX\t = \t1174
+    reg = re.compile(f"({key})\\s*=\\s*([^\\s]+)")
+    res = reg.findall(asccov)
+    # could be more rigerous about order by moving tuple results into dict
+    return [x[1] for x in res]
+
+
+def read_csa(csa) -> Optional[dict]:
     """
     extract parameters from siemens CSA
     :param csa: content of siemens private tag (0x0029, 0x1010)
-    :return: [pepd, ipat] is phase encode positive direction and GRAPA iPATModeText
+    :return: nibabel's csareader dictionary or None if cannot read
 
     >>> read_csa(None) is None
     True
@@ -142,7 +192,12 @@ def read_tags(dcm_path: os.PathLike, tags: TagDicts) -> TagValues:
     csa = read_csa(dcm.get((0x0029, 0x1010)))
     for tag in tags:
         k = tag["name"]
-        if tag["loc"] == "csa":
+        if k == "Shims":
+            # 20241118: add shims
+            csa_s = read_csa(dcm.get((0x0029, 0x1020)))
+            shims = read_shims(csa_s)
+            out[k] = ",".join(shims)
+        elif tag["loc"] == "csa":
             out[k] = csa_fetch(csa, tag["tag"]) if csa is not None else NULLVAL.value
         else:
             out[k] = dcm.get(tag["tag"], NULLVAL).value

generated_report.html

@@ -0,0 +1,74 @@
+
+<html>
+<head>
+    <style>
+        table {
+            border-collapse: collapse;
+            width: 100%;
+        }
+        th, td {
+            border: 1px solid black;
+            padding: 8px;
+            text-align: center;
+        }
+        th {
+            background-color: #f2f2f2;
+        }
+        .match { background-color: #d4edda; } /* green for matches */
+        .mismatch { background-color: #f8d7da; } /* red for mismatches */
+    </style>
+</head>
+<body>
+    <h2>DICOM Header Compliance Report</h2>
+    <table>
+        <tr>
+            <th>Header</th>
+            <th>Expected Value</th>
+            <th>Actual Value</th>
+        </tr>
+
+            <tr class="match">
+            <td>Project</td>
+            <td>Brain^wpc-8620</td>
+            <td>Brain^wpc-8620</td>
+        </tr>
+
+            <tr class="match">
+            <td>SequenceName</td>
+            <td>HabitTask</td>
+            <td>HabitTask</td>
+        </tr>
+
+            <tr class="mismatch">
+            <td>TR</td>
+            <td>1300</td>
+            <td>1400</td>
+        </tr>
+
+            <tr class="match">
+            <td>TE</td>
+            <td>30</td>
+            <td>30</td>
+        </tr>
+
+            <tr class="mismatch">
+            <td>FA</td>
+            <td>60</td>
+            <td>70</td>
+        </tr>
+
+            <tr class="match">
+            <td>iPAT</td>
+            <td>GRAPPA</td>
+            <td>GRAPPA</td>
+        </tr>
+
+            <tr class="match">
+            <td>Comments</td>
+            <td>Unaliased MB3/PE4/LB SENSE1</td>
+            <td>Unaliased MB3/PE4/LB SENS1</td>
+        </tr>
+
+    </table>
+</body>
+</html>