diff --git a/.appveyor.yml b/.appveyor.yml index a1aa35b6..e0a05a46 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -78,8 +78,8 @@ for: - mkdir -p build/bin - lipo -create -output build/bin/dcm2niix intel/bin/dcm2niix apple/bin/dcm2niix - strip -Sx build/bin/* - - 7z a dcm2niix_macos.zip ./build/bin/* &>/dev/null - - appveyor PushArtifact dcm2niix_macos.zip + - 7z a dcm2niix_mac.zip ./build/bin/* &>/dev/null + - appveyor PushArtifact dcm2niix_mac.zip deploy: - provider: GitHub diff --git a/BIDS/README.md b/BIDS/README.md index 37200c58..c48d5a43 100644 --- a/BIDS/README.md +++ b/BIDS/README.md @@ -263,7 +263,7 @@ Note that [many of the fields listed by BIDS](https://bids-specification.readthe |---------------------------------|--------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|------------| | ArterialSpinLabelingType | | `PASL` or `PCASL` | B | | LabelOffset | | 2D pCASL [ep2d_pcasl](http://www.loft-lab.org) | D | -| PostLabelDelay | s | 2D pCASL [ep2d_pcasl](http://www.loft-lab.org), [FAIREST](https://pubmed.ncbi.nlm.nih.gov/11746944/) (`ep2d_fairest`) | D | +| PostLabelingDelay | s | 2D pCASL [ep2d_pcasl](http://www.loft-lab.org), [FAIREST](https://pubmed.ncbi.nlm.nih.gov/11746944/) (`ep2d_fairest`) | B | | PostInversionDelay | s | 2D [FAIREST](https://pubmed.ncbi.nlm.nih.gov/11746944/) (`ep2d_fairest`) | D | | NumRFBlocks | | 2D pCASL [ep2d_pcasl](http://www.loft-lab.org) | D | | RFGap | | 2,3D pCASL [ep2d_pcasl, tgse_pcasl](http://www.loft-lab.org) | D | @@ -334,7 +334,7 @@ Fields specific to [Siemens XA-series](https://github.com/rordenlab/dcm2niix/tre | ReceiveCoilActiveElements | | DICOM tag 0021,114F | B | | BandwidthPerPixelPhaseEncode | Hz | DICOM tag 0021,1153 | D | | ScanningSequence | | DICOM tag 0021,105a | D | -| PostLabelDelay | s | DICOM tag 0018,9258 | D | +| PostLabelingDelay | s | DICOM tag 0018,9258 | B | | NonlinearGradientCorrection | b | 0008,0008 or 0021,1175 | B | | PhaseEncodingDirection | | polarity from 0021,111c | B | | SpoilingState | | 0021,105B | B | @@ -349,3 +349,7 @@ Fields specific to United Imaging Healthcare systems (e.g. uMR 770). |--------------------------------|------|--------------------------|------------| | BandwidthPerPixelPhaseEncode | Hz | DICOM tag 0019,1028 | B | | ParallelReductionFactorInPlane | | DICOM tag 0065,100D | B | + +### Notes + + - Prior to v1.0.20230411, dcm2niix created the internal tag `PostLabelDelay`, subsequent versions use the newly introduced BIDS tag `PostLabelingDelay` \ No newline at end of file diff --git a/BidsGuess/BidsGuess.png b/BidsGuess/BidsGuess.png new file mode 100644 index 00000000..dbc7102d Binary files /dev/null and b/BidsGuess/BidsGuess.png differ diff --git a/BidsGuess/README.md b/BidsGuess/README.md new file mode 100644 index 00000000..07ac8a52 --- /dev/null +++ b/BidsGuess/README.md @@ -0,0 +1,63 @@ +## About + +**The BidsGuess feature is intended for wrapper developers only and should not be used in isolation of a wrapper** + +dcm2niix version v1.0.20230731 and later will insert the field `BidsGuess` into the [BIDS](https://bids-specification.readthedocs.io/en/stable/) JSON sidecar files. This experimental feature is designed to aid wrappers that use dcm2niix to create BIDS compatible datasets. BIDS aids reproducible, reusable and automatic analysis of neuroimaging data. This feature was inspired by the automatic BIDS conversion wrappers [niix2bids](https://github.com/benoitberanger/niix2bids) and [ezBIDS](https://brainlife.io/ezbids/). The `BidsGuess` field can be leveraged by wrappers to fully automate conversion (though this will likely require a ezBIDS style user validation) or to flag improbable naming from user configuration files. + +## Compiling the development branch + +The BidsGuess feature is currently only available in the development branch, so you will need to compile and run this version (v1.0.20230731 and later). + +``` +git clone --branch development https://github.com/rordenlab/dcm2niix.git +cd dcm2niix/console +make +./dcm2niix +``` + +## A concrete example + +The BidsGuess feature is designed to be leveraged by dcm2niix wrappers, and not used to directly create BIDS format files. Specifically, bidsGuess converts each DICOM series in isolation, and has no information about the user intention. Therefore, it is unable to resolve fieldmap [IntendedFor](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#using-intendedfor-metadata), unable to distinguish fMRI tasks from resting state, BIDS subject ID, BIDS session number, or create meaningful [dataset_description](https://bids-specification.readthedocs.io/en/stable/glossary.html#dataset_description-files) or [readme](https://bids-specification.readthedocs.io/en/stable/glossary.html#readme-files) files. + +For the developers of wrappers, you can use the hazardous file naming argument (`-f $h`) to create a minimal BIDS structure for the [bids-validator](https://github.com/bids-standard/bids-validator). Note that this mode will always claim that the data is from `sub-1` and that data is only from a single session. Here is a simple example: + +``` +cd ~ +mkdir bids +git clone git@github.com:neurolabusc/dcm_qa_pdt2.git +dcm2niix -f %h -w 1 -i y -o ~/bids ~/dcm_qa_pdt2 +bids-validator ~/bids +``` +You can see that the bids-validator is happy with the results and that the data appears organized correctly: + +![BidsGuess](BidsGuess.png) + +Inspecting the JSON files, we can see that dcm2niix has suggested a likely <[datatype](https://bids-specification.readthedocs.io/en/stable/schema/index.html#bids-filenames)> (`anat`) and <[entities](https://bids-specification.readthedocs.io/en/stable/schema/index.html#bids-filenames)>. + +``` + "BidsGuess": ["anat","_acq-tse2_run-3_PDw"], + "BidsGuess": ["anat","_acq-tse2_run-3_T2w"], +``` + +Developers can use the `hazardous` file naming to validate and extend the modality detection. However, production quality wrappers should use the `BidsGuess` in the JSON file combined by a file naming scheme that avoids name clashes between different participants and sessions (e.g. one can segment data by datetime, series and protocol name with `-f %t/%s_%p`). + +## The acq entity + +The `BidsGuess` creates a meaningful [`_acq-` entity](https://bids-specification.readthedocs.io/en/stable/appendices/entities.html#acq) that can provide consistency across wrappers, minimizes the risk of naming clashes and allows users to quickly detect sequence details. The first part of this reveals the manufacturer's name for the sequence. For example, a Siemens 2D turbo spin-echo will report `tse2`, while a 2D echo-planar spin-echo will report `epse2` and a 3D turbo flash will report `tfl3`. If acceleration was used, this will be reported next. For example, a 2D echo-planar gradient-echo with x3 in-plane and x4 mult-iband acceleration would be reported as `epfid2p3m4`. + +## The run entity + +The `BidsGuess` will typically include a [`_run-` entity](https://bids-specification.readthedocs.io/en/stable/appendices/entities.html#run) that reports the DICOM series number of an image. This is often useful for determining the temporal order of images (e.g. if the first attempt to acquire the data was due to head motion). This entity also avoids naming clashes. Note that wrapper developers may want to remove this entity from the BIDS guess - it is redundant with the `SeriesNumber` stored in the sidecar JSON. Finally, dcm2niix will not append a `_run` entity for fieldmaps: the BIDS validation tool expects that the different images associated with a fieldmap have identical file names except the suffix (e.g. `magnitude1` and `phasediff`). + +## The dir entity + +The `BidsGuess` will typically include a [`_dir-` entity](https://bids-specification.readthedocs.io/en/stable/appendices/entities.html#dir) for modalities when required by the BIDS validator. This field is redundant with the JSON `PhaseEncodingDirection` field. Note that the JSON uses the values `j`, `j-`, `k` and `k-` to specify row versus column and polarity. In contrast, the BidsGuess will use the `AP`, `PA`, `LR`, and `RL` tags though note these tags will only be correct for axial acquisitions. Note that it is impossible to infer phase encoding polarity for Philips data, so this entity will not be populated leading to errors from the bids-validator. + +## Limitations + +This feature is very experimental, and is currently provided to get feedback from wrapper developers and to get community support to enhance the guessing accuracy. + + - This feature currently only supports images from GE, Philips and Siemens MR scanners. + - Multi-Echo MP-RAGE where individual echos (rather than mean) are saved will include the [_echo](https://bids-specification.readthedocs.io/en/stable/appendices/entities.html#echo) entity to distinguish them, e.g. `_acq-tflme3p2_run-5_echo-2_T1w`, `_acq-tflme3p2_run-5_echo-1_T1w`. This will [cause issues](https://github.com/bids-standard/bids-specification/issues/654) with the current bids-validator. + - ASL datasets will generate errors with the bids-validator. The ASL BEP introduced many required tags for converted data without explaining how these are determined or even if they exist in the core DICOM images. The validation dataset [only provides the desired BIDS translation and not the source DICOMs](https://osf.io/yru2q/). + - Philips DICOMs are underspecified for BIDS conversion. Beyond the previously noted issue with phaseEncodingDirection, the `SliceTiming` is also unknown. This is a limitation of Philips DICOM, not dcm2niix. diff --git a/CONTRIBUTE.md b/CONTRIBUTE.md index dcb0e593..6d84a0a6 100644 --- a/CONTRIBUTE.md +++ b/CONTRIBUTE.md @@ -7,7 +7,7 @@ Like the [Brain Imaging Data Structure](https://bids.neuroimaging.io/get_involve The easiest way to contribute to dcm2niix is to ask questions you have by [generating Github issues](https://github.com/rordenlab/dcm2niix/issues) or [asking a question on the NITRC forum](https://www.nitrc.org/forum/?group_id=880). The code is open source, and you can share your improvements by [creating a pull request](https://github.com/rordenlab/dcm2niix/pulls). -dcm2niix is a community project that has benefitted from many [contrbutors](https://github.com/rordenlab/dcm2niix/graphs/contributors). +dcm2niix is a community project that has benefitted from many [contributors](https://github.com/rordenlab/dcm2niix/graphs/contributors). The INCF suggests indicating who is responsible for maintaining software for [stability and support](https://incf.org/incf-standards-review-criteria-v20). Therefore, below we indicate several active contributors and their primary domain of expertise. However, this list is not comprehensive, and it is noted that the project has been supported by contributions from many users. This list does not reflect magnitude of prior contributions, rather it is a non-exhaustive list of members who are actively maintaining the project. diff --git a/FILENAMING.md b/FILENAMING.md index 890796a0..436114a1 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -39,12 +39,13 @@ In general dcm2niix creates images with 3D dimensions, or 4 dimensions when the - _cNx.._cNz where C* refers to the coil name (typically only seen for uncombined data, where a separate image is generated for each antenna) - _e1..eN echo number for multi-echo sequences - _Eq is commonly seen in [CT scans](https://github.com/neurolabusc/dcm_qa_ct). For example, CT scans of the brain often have many slices closely packed near the brain stem and only a few slices spread far apart near the top of the head. Variable between-slice spacing is rarer in MRI, and if you see this from a MRI sequence you should ensure that [all of the acquired slices have been provided to dcm2niix](https://neurostars.org/t/field-mapping-siemens-scanners-dcm2niix-output-2-bids/2075/7). NIfTI assumes all 2D slices that form a 3D stack are equidistant. Therefore, dcm2niix reslices the input data to generate an equidistant volume. - - _ph phase map - - _iN appended image number for non-parallel slices + - _i1..N appended image number for non-parallel slices - _imaginary imaginary component of complex image - _MoCo is appended to the ProtocolName if Image Type (0008,0008) includes the term 'MOCO'. This helps disambiguate Siemens fMRI runs where both motion corrected and raw data is stored for a single session. - - _real real component of complex image + - _ph phase map - _phMag rare case where phase and magnitude are saved as the 4th dimension + - _r1..rN repetition number for sequences with multiple repetition times + - _real real component of complex image - _t If the trigger delay time (0020,9153) or trigger time (0018,1060) is non-zero, it will be recorded in the file name. For example, the files "T1_t178.nii" and "T1_t511" suggests that the T1 scan was acquired with two cardiac trigger delays (178 and 511ms after the last R-peak). - _Tilt is specific to [CT scans](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Computed_Tomography_.28CT.2C_CAT.29). These scans can be acquired with a gantry tilt that causes a skew that can not be stored in a NIfTI qForm. Therefore, the slices are resampled to remove the effect of tilt. @@ -69,7 +70,7 @@ dcm2niix will attempt to write your image using the naming scheme you specify wi ## Special Characters -[Some characters are not permitted](https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names) in file names. The following characters will be replaced with underscorces (`_`). Note that the forbidden characters vary between operating systems (Linux only forbids the forward slash, MacOS forbids forward slash and colon, while Windows forbids any of the characters listed below). To ensure that files can be easily copied between file systems, [dcm2niix restricts file names to characters allowed by Windows](https://github.com/rordenlab/dcm2niix/issues/237). While technically legal in all filesystems, the semicolon can wreak havoc in [Windows](https://stackoverflow.com/questions/3869594/semi-colons-in-windows-filenames) and [Linux](https://forums.plex.tv/t/linux-hates-semicolons-in-file-names/49098/2). +[Some characters are not permitted](https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names) in file names. The following characters will be replaced with underscorces (`_`). Note that the forbidden characters vary between operating systems (Linux only forbids the forward slash, MacOS forbids forward slash and colon, while Windows forbids any of the characters listed below). To ensure that files can be easily copied between file systems, [dcm2niix restricts file names to characters allowed by Windows](https://github.com/rordenlab/dcm2niix/issues/237). While technically legal in all filesystems, the semicolon can wreak havoc in [Windows](https://stackoverflow.com/questions/3869594/semi-colons-in-windows-filenames) and [Linux](https://forums.plex.tv/t/linux-hates-semicolons-in-file-names/49098/2). Likewise, the [dollar sign can cause issues](https://github.com/rordenlab/dcm2niix/issues/749). ### List of Forbidden Characters ``` @@ -82,7 +83,8 @@ dcm2niix will attempt to write your image using the naming scheme you specify wi | (vertical bar or pipe) ? (question mark) * (asterisk) -; (semicolon) +; (semicolon) +$ (dollar sign) ``` [Control characters](https://en.wikipedia.org/wiki/ASCII#Control_characters) like backspace and tab are also forbidden. diff --git a/GE/README.md b/GE/README.md index 7efbdc92..5bfbfceb 100644 --- a/GE/README.md +++ b/GE/README.md @@ -50,15 +50,25 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6 ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Specifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first echo to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor` 2 (Full Fourier) or 4 (Partial Fourier)). +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Specifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first echo to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). Note that the Effective Echo Spacing (0043,102C) in GE DICOMS is defined as the time between two consecutives acquired phase encoding lines divided by the number of shots (usually, this is equal to 1 for fMRI and Diffusion). +While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor` 2 (Full Fourier) or 4 (Partial Fourier)). -The formula for FSL's definition of TotalReadoutTime (in seconds) is: +Let `NotPhysicalNumberOfAcquiredPELinesGE` be the number of acquired phase encoding lines if there was no partial Fourier and `NotPhysicalTotalReadOutTimeGE` be the physical total read-out time if there was no partial Fourier. Please, note that these two intermediate variables do not take partial Fourier into account. These two variables can be computed as ``` -TotalReadoutTime = ( ( ceil ((1/Round_factor) * PE_AcquisitionMatrix / Asset_R_factor ) * Round_factor) - 1 ] * EchoSpacing * 0.000001 -EffectiveEchoSpacing = TotalReadoutTime/ (reconMatrixPE - 1) +NotPhysicalNumberOfAcquiredPELinesGE = (ceil((1/Round_factor) * PE_AcquisitionMatrix / Asset_R_factor) * Round_factor) +NotPhysicalTotalReadOutTimeGE = (NotPhysicalNumberOfAcquiredPELinesGE - 1) * EchoSpacing * 0.000001 ``` +with `EchoSpacing` referring to the GE private DICOM field "(0043,102C) Effective Echo Spacing". +Then, the formula for FSL's definition of `EffectiveEchoSpacing` and `TotalReadoutTime` (in seconds) are: + +``` +EffectiveEchoSpacing = NotPhysicalTotalReadOutTimeGE / (PE_AcquisitionMatrix - 1) +TotalReadoutTime = EffectiveEchoSpacing * (ReconMatrixPE - 1) +``` +When there is no image interpolation (i.e. `ReconMatrixPE = PE_AcquisitionMatrix`) the `TotalReadoutTime` has the same value as `NotPhysicalTotalReadOutTimeGE`. In other words, `NotPhysicalTotalReadOutTimeGE` is the `TotalReadoutTime` without any image interpolation. + Consider an example: ``` @@ -73,6 +83,7 @@ From this we can derive: ``` ASSET= 1.5 PE_AcquisitionMatrix= 128 EchoSpacing= 636 Round_Factor= 4 TotalReadoutTime= 0.055332 ``` +For debugging purposes, the intermediate variables `NotPhysicalTotalReadOutTimeGE`, `NotPhysicalNumberOfAcquiredPELinesGE` and `EchoSpacing` (renamed `EchoSpacingMicroSecondsGE`) are written to the BIDS-sidecar JSON file when dcm2niix is compiled with the flag `MY_DEBUG`. ## Image Acceleration diff --git a/README.md b/README.md index 1a9f26ff..9a15240b 100644 --- a/README.md +++ b/README.md @@ -141,7 +141,6 @@ The following tools exploit dcm2niix - [BraTS-Preprocessor](https://neuronflow.github.io/BraTS-Preprocessor/) uses dcm2niix to import files for [Brain Tumor Segmentation](https://www.frontiersin.org/articles/10.3389/fnins.2020.00125/full). - [CardioNIfTI](https://github.com/UK-Digital-Heart-Project/CardioNIfTI) processes cardiac MR DICOM datasets and converts them to NIfTI. - [clinica](https://github.com/aramis-lab/clinica) is a software platform for clinical neuroimaging studies that uses dcm2niix to convert DICOM images. - - [clinical_dicom2bids_smk](https://github.com/greydongilmore/clinical_dicom2bids_smk) Snakemake workflow to convert a clinical dicom directory into BIDS structure. - [clpipe](https://github.com/cohenlabUNC/clpipe) uses dcm2bids for DICOM import. - [conversion](https://github.com/pnlbwh/conversion) is a Python library that can convert dcm2niix created NIfTI files to the popular NRRD format (including DWI gradient tables). Note, recent versions of dcm2niix can directly convert DICOM images to NRRD. - [convert_source](https://github.com/AdebayoBraimah/convert_source) to convert DICOM to BIDS directory layout. @@ -162,7 +161,7 @@ The following tools exploit dcm2niix - [dicom2nifti_batch](https://github.com/scanUCLA/dicom2nifti_batch) is a Matlab script for automating dcm2niix. - [dicomConversionToNifti](https://github.com/bsmarine/dicomConversionToNifti) converts, de-identifies and assigns standardized naming convention to medical imaging. - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. - - [DPABI Data Processing & Analysis for Brain Imaging](http://rfmri.org/dpabi) includes dcm2niix. + - [DPABI](https://github.com/Chaogan-Yan/DPABI) [Data Processing & Analysis for Brain Imaging](https://rfmri.org/DPABI) includes dcm2niix. - [ExploreASL](https://sites.google.com/view/exploreasl/exploreasl) uses dcm2niix to import images. - [ExploreASL-GUI](https://github.com/MauricePasternak/ExploreASL-GUI) uses dcm2niix for image conversion. - [ezBIDS](https://github.com/brainlife/ezbids) is a [web service](https://brainlife.io/ezbids/) for converting directory full of DICOM images into BIDS without users having to learn python nor custom configuration file. @@ -196,6 +195,7 @@ The following tools exploit dcm2niix - [pydra-dcm2bids](https://github.com/aramis-lab/pydra-dcm2bids) supports Pydra tasks for dcm2bids. - [pydra-dcm2niix](https://github.com/nipype/pydra-dcm2niix) is a contains Pydra task interface for dcm2niix. - [qsm](https://github.com/CAIsr/qsm) Quantitative Susceptibility Mapping software. + - [QSMxT](https://github.com/QSMxT/QSMxT) is an end-to-end software toolbox for Quantitative Susceptibility Mapping. - [reproin](https://github.com/ReproNim/reproin) is a setup for automatic generation of shareable, version-controlled BIDS datasets from MR scanners. - [Retina_OCT_dcm2nii](https://github.com/Choupan/Retina_OCT_dcm2nii) converts optical coherence tomography (OCT) data to NIfTI. - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) Flywheel Gear (docker). diff --git a/RENAMING.md b/RENAMING.md index 5d115d9a..d2c9f429 100644 --- a/RENAMING.md +++ b/RENAMING.md @@ -4,7 +4,7 @@ dcm2niix is primarily designed to convert DICOM images into NIfTI images. Howeve ## Limitation -dcm2niix renames and copies the DICOM images, but the current version does not copy or create a new [DICOMDIR](http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_F.2.2.2.html) file. Most users ignore these files. However, you should not use this featire if you wish to preserve your DICOMDIR files. +dcm2niix renames and copies the DICOM images, but the current version does not copy or create a new [DICOMDIR](http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_F.2.2.2.html) file. Most users ignore these files. However, you should not use this feature if you wish to preserve your DICOMDIR files. Note that this feature only copies your DICOM images with a new filename. It does not modify the contents of the DICOM header. This means it will not compress or anonymize your files. Free tools for these functions include [dcmcjpeg](https://dicom.offis.de/dcmtk.php.en), [gdcmanon](http://gdcm.sourceforge.net/html/gdcmanon.html) and [gdcmconv](http://gdcm.sourceforge.net/html/gdcmconv.html). diff --git a/SuperBuild/External-CLOUDFLARE-ZLIB.cmake b/SuperBuild/External-CLOUDFLARE-ZLIB.cmake index 4c6a6acb..74168530 100644 --- a/SuperBuild/External-CLOUDFLARE-ZLIB.cmake +++ b/SuperBuild/External-CLOUDFLARE-ZLIB.cmake @@ -11,7 +11,6 @@ ExternalProject_Add(zlib ${OSX_ARCHITECTURES} # Compiler settings -DCMAKE_C_COMPILER:FILEPATH=${CMAKE_C_COMPILER} - -DCMAKE_CXX_COMPILER:FILEPATH=${CMAKE_CXX_COMPILER} # Install directories -DCMAKE_INSTALL_PREFIX:PATH=${DEP_INSTALL_DIR} ) diff --git a/SuperBuild/External-OPENJPEG.cmake b/SuperBuild/External-OPENJPEG.cmake index 1d99ac28..a311fa0b 100644 --- a/SuperBuild/External-OPENJPEG.cmake +++ b/SuperBuild/External-OPENJPEG.cmake @@ -12,7 +12,6 @@ ExternalProject_Add(openjpeg ${OSX_ARCHITECTURES} # Compiler settings -DCMAKE_C_COMPILER:FILEPATH=${CMAKE_C_COMPILER} - # Not used -DCMAKE_CXX_COMPILER:FILEPATH=${CMAKE_CXX_COMPILER} # Install directories -DCMAKE_INSTALL_PREFIX:PATH=${DEP_INSTALL_DIR} ) diff --git a/UIH/README.md b/UIH/README.md index 2a2d9d64..02c471d1 100644 --- a/UIH/README.md +++ b/UIH/README.md @@ -17,7 +17,7 @@ UIH supports two ways of archiving the DWI/DTI and fMRI data. One way is one DIC |0065,1002 | ImageOrientationDisplayed | SH | 1 | Image Orientation Displayed | Sag or Sag>Cor| |0065,1003 | ReceiveCoil | LO | 1 | Receive Coil Information | H 8| |0065,1004 | Interpolation | SH | 1 | Interpolation | I| -|0065,1005 | PE Direction Displayed | SH | 1 | Phase encoding diretion displayed | A->P or H->F| +|0065,1005 | PE Direction Displayed | SH | 1 | Phase encoding direction displayed | A->P or H->F| |0065,1006 | Slice Group ID | IS | 1 | Slice Group ID | 1| |0065,1007 | Uprotocol | OB | 1 | Uprotocol value |  | |0065,1009 | BActualValue | FD | 1 | Actual B-Value from sequence | 1000.0| diff --git a/console/CMakeLists.txt b/console/CMakeLists.txt index aa52bd7c..79fcd38b 100644 --- a/console/CMakeLists.txt +++ b/console/CMakeLists.txt @@ -91,9 +91,9 @@ set(DCM2NIIX_SRCS nii_dicom_batch.cpp) -option(USE_JNIfTI "Build with JNIfTI support" ON) +option(USE_JNIFTI "Build with JNIfTI support" ON) if(USE_JNIFTI) - add_definitions(-DmyEnableJNIfTI) + add_definitions(-DmyEnableJNIFTI) set(DCM2NIIX_SRCS ${DCM2NIIX_SRCS} cJSON.cpp base64.cpp) endif() diff --git a/console/charls/README.md b/console/charls/README.md index 4ff0549d..4407d62e 100644 --- a/console/charls/README.md +++ b/console/charls/README.md @@ -64,7 +64,7 @@ According to preliminary test results published on http://imagecompression.info/ * No support for (optional) JPEG restart markers (RST). These markers are rarely used in practice. * No support for the SPIFF file header. * No support for oversize image dimension. Maximum supported image dimensions are [1, 65535] by [1, 65535]. -* After releasing the original baseline standrd 14495-1:1999, ISO released an extension to the JPEG-LS standard called ISO/IEC 14495-2:2003: "Lossless and near-lossless compression of continuous-tone still images: Extensions". CharLS doesn't support these extensions. +* After releasing the original baseline standard 14495-1:1999, ISO released an extension to the JPEG-LS standard called ISO/IEC 14495-2:2003: "Lossless and near-lossless compression of continuous-tone still images: Extensions". CharLS doesn't support these extensions. ## Supported platforms diff --git a/console/dcm2niix_fswrapper.cpp b/console/dcm2niix_fswrapper.cpp index e8eedfc6..1286fa37 100644 --- a/console/dcm2niix_fswrapper.cpp +++ b/console/dcm2niix_fswrapper.cpp @@ -1,4 +1,5 @@ #include +#include #include "nii_dicom.h" #include "dcm2niix_fswrapper.h" @@ -52,7 +53,7 @@ numSeries = 0 */ // set TDCMopts defaults, overwrite settings to output in mgz orientation -void dcm2niix_fswrapper::setOpts(const char* dcmindir, const char* niioutdir, bool createBIDS) +void dcm2niix_fswrapper::setOpts(const char* dcmindir, const char* niioutdir, bool createBIDS, int ForceStackSameSeries) { memset(&tdcmOpts, 0, sizeof(tdcmOpts)); setDefaultOpts(&tdcmOpts, NULL); @@ -62,7 +63,9 @@ void dcm2niix_fswrapper::setOpts(const char* dcmindir, const char* niioutdir, bo if (niioutdir != NULL) strcpy(tdcmOpts.outdir, niioutdir); - strcpy(tdcmOpts.filename, "%4s.%p"); + // dcmunpack actually uses seriesDescription, set FName = `printf %04d.$descr $series` + // change it from "%4s.%p" to "%4s.%d" + strcpy(tdcmOpts.filename, "%4s.%d"); // set the options for freesurfer mgz orientation tdcmOpts.isRotate3DAcq = false; @@ -70,7 +73,7 @@ void dcm2niix_fswrapper::setOpts(const char* dcmindir, const char* niioutdir, bo tdcmOpts.isIgnoreSeriesInstanceUID = true; tdcmOpts.isCreateBIDS = createBIDS; tdcmOpts.isGz = false; - tdcmOpts.isForceStackSameSeries = 1; // merge 2D slice '-m y' + tdcmOpts.isForceStackSameSeries = ForceStackSameSeries; // merge 2D slice '-m y', tdcmOpts.isForceStackSameSeries = 1 tdcmOpts.isForceStackDCE = false; //tdcmOpts.isForceOnsetTimes = false; } @@ -107,6 +110,12 @@ MRIFSSTRUCT* dcm2niix_fswrapper::getMrifsStruct(void) return nii_getMrifsStruct(); } +// interface to nii_getMrifsStructVector() +std::vector* dcm2niix_fswrapper::getMrifsStructVector(void) +{ + return nii_getMrifsStructVector(); +} + // return nifti header saved in MRIFSSTRUCT nifti_1_header* dcm2niix_fswrapper::getNiiHeader(void) { @@ -121,11 +130,226 @@ const unsigned char* dcm2niix_fswrapper::getMRIimg(void) return mrifsStruct->imgM; } -void dcm2niix_fswrapper::dicomDump(const char* dicomdir) +void dcm2niix_fswrapper::dicomDump(const char* dicomdir, const char *series_info, bool max, bool extrainfo) { strcpy(tdcmOpts.indir, dicomdir); tdcmOpts.isDumpNotConvert = true; nii_loadDirCore(tdcmOpts.indir, &tdcmOpts); + char fnamecopy[2048] = {'\0'}; + memcpy(fnamecopy, series_info, strlen(series_info)); + char *logdir = dirname(fnamecopy); + + FILE *fpout = fopen(series_info, "w"); + + std::vector *mrifsStruct_vector = dcm2niix_fswrapper::getMrifsStructVector(); + int nitems = (*mrifsStruct_vector).size(); + for (int n = 0; n < nitems; n++) + { + struct TDICOMdata *tdicomData = &(*mrifsStruct_vector)[n].tdicomData; + + // output the dicom list for the series into seriesNum-dicomflst.txt + char dicomflst[2048] = {'\0'}; + sprintf(dicomflst, "%s/%ld-dicomflst.txt", logdir, tdicomData->seriesNum); + FILE *fp_dcmLst = fopen(dicomflst, "w"); + for (int nDcm = 0; nDcm < (*mrifsStruct_vector)[n].nDcm; nDcm++) + fprintf(fp_dcmLst, "%s\n", (*mrifsStruct_vector)[n].dicomlst[nDcm]); + fclose(fp_dcmLst); + + // output series_info + fprintf(fpout, "%ld %s %f %f %f %f\\%f %c %f %s %s", + tdicomData->seriesNum, tdicomData->seriesDescription, + tdicomData->TE, tdicomData->TR, tdicomData->flipAngle, tdicomData->xyzMM[1], tdicomData->xyzMM[2], + tdicomData->phaseEncodingRC, tdicomData->pixelBandwidth, (*mrifsStruct_vector)[n].dicomfile, tdicomData->imageType); +#if 0 + if (max) + { + //$max kImageStart 0x7FE0 + (0x0010 << 16) + fprintf(fpout, " max-value"); + } +#endif + if (extrainfo) + fprintf(fpout, " %s %s %s %s %f %s", tdicomData->patientName, tdicomData->studyDate, mfrCode2Str(tdicomData->manufacturer), tdicomData->manufacturersModelName, tdicomData->fieldStrength, tdicomData->deviceSerialNumber); + + fprintf(fpout, "\n"); + } + fclose(fpout); + return; } + +const char *dcm2niix_fswrapper::mfrCode2Str(int code) +{ + if (code == kMANUFACTURER_SIEMENS) + return "SIEMENS"; + else if (code == kMANUFACTURER_GE) + return "GE"; + else if (code == kMANUFACTURER_PHILIPS) + return "PHILIPS"; + else if (code == kMANUFACTURER_TOSHIBA) + return "TOSHIBA"; + else if (code == kMANUFACTURER_UIH) + return "UIH"; + else if (code == kMANUFACTURER_BRUKER) + return "BRUKER"; + else if (code == kMANUFACTURER_HITACHI) + return "HITACHI"; + else if (code == kMANUFACTURER_CANON) + return "CANON"; + else if (code == kMANUFACTURER_MEDISO) + return "MEDISO"; + else + return "UNKNOWN"; +} + + +void dcm2niix_fswrapper::seriesInfoDump(FILE *fpdump, const MRIFSSTRUCT *pmrifsStruct) +{ + // print dicom-info for the series converted using (*mrifsStruct_vector)[0] + // see output from mri_probedicom --i $f --no-siemens-ascii + + fprintf(fpdump, "###### %s ######\n", pmrifsStruct->dicomfile); + + const struct TDICOMdata *tdicomData = &(pmrifsStruct->tdicomData); + + // kManufacturer 0x0008 + (0x0070 << 16) + fprintf(fpdump, "Manufacturer %s\n", dcm2niix_fswrapper::mfrCode2Str(tdicomData->manufacturer)); + + // kManufacturersModelName 0x0008 + (0x1090 << 16) + fprintf(fpdump, "ScannerModel %s\n", tdicomData->manufacturersModelName); + + // kSoftwareVersions 0x0018 + (0x1020 << 16) //LO + fprintf(fpdump, "SoftwareVersion %s\n", tdicomData->softwareVersions); + + // kDeviceSerialNumber 0x0018 + (0x1000 << 16) //LO + fprintf(fpdump, "ScannerSerialNo %s\n", tdicomData->deviceSerialNumber); + + // kInstitutionName 0x0008 + (0x0080 << 16) + fprintf(fpdump, "Institution %s\n", tdicomData->institutionName); + + // kSeriesDescription 0x0008 + (0x103E << 16) // '0008' '103E' 'LO' 'SeriesDescription' + fprintf(fpdump, "SeriesDescription %s\n", tdicomData->seriesDescription); + + // kStudyInstanceUID 0x0020 + (0x000D << 16) + fprintf(fpdump, "StudyUID %s\n", tdicomData->studyInstanceUID); + + // kStudyDate 0x0008 + (0x0020 << 16) + fprintf(fpdump, "StudyDate %s\n", tdicomData->studyDate); + + // kStudyTime 0x0008 + (0x0030 << 16) + fprintf(fpdump, "StudyTime %s\n", tdicomData->studyTime); + + // kPatientName 0x0010 + (0x0010 << 16) + fprintf(fpdump, "PatientName %s\n", tdicomData->patientName); + + // kSeriesNum 0x0020 + (0x0011 << 16) + fprintf(fpdump, "SeriesNo %ld\n", tdicomData->seriesNum); + + // kImageNum 0x0020 + (0x0013 << 16) + fprintf(fpdump, "ImageNo %d\n", tdicomData->imageNum); + + // kMRAcquisitionType 0x0018 + (0x0023 << 16) + fprintf(fpdump, "AcquisitionType %s\n", (tdicomData->is3DAcq) ? "3D" : ((tdicomData->is3DAcq) ? "2D" : "unknown")); //dcm2niix has two values: d.is2DAcq, d.is3DAcq + + // kImageTypeTag 0x0008 + (0x0008 << 16) + fprintf(fpdump, "ImageType %s\n", tdicomData->imageType); + + // kImagingFrequency 0x0018 + (0x0084 << 16) //DS + fprintf(fpdump, "ImagingFrequency %f\n", tdicomData->imagingFrequency); + + // kPixelBandwidth 0x0018 + (0x0095 << 16) //'DS' 'PixelBandwidth' + fprintf(fpdump, "PixelFrequency %f\n", tdicomData->pixelBandwidth); + + // dcm2niix doesn't seem to retrieve this 0x18, 0x85 + //fprintf(fpdump, "ImagedNucleus %s\n",e->d.string); + + // kEchoNum 0x0018 + (0x0086 << 16) //IS + fprintf(fpdump, "EchoNumber %d\n", tdicomData->echoNum); + + // kMagneticFieldStrength 0x0018 + (0x0087 << 16) //DS + fprintf(fpdump, "FieldStrength %f\n", tdicomData->fieldStrength); + + // kSequenceName 0x0018 + (0x0024 << 16) + fprintf(fpdump, "PulseSequence %s\n", tdicomData->sequenceName); + + // kProtocolName 0x0018 + (0x1030 << 16) + fprintf(fpdump, "ProtocolName %s\n", tdicomData->protocolName); + + // kScanningSequence 0x0018 + (0x0020 << 16) + fprintf(fpdump, "ScanningSequence %s\n", tdicomData->scanningSequence); + + // dcm2niix doesn't seem to retrieve this + // kTransmitCoilName 0x0018 + (0x1251 << 16) // SH issue527 + //fprintf(fpdump, "TransmittingCoil %s\n",e->d.string); + + // kPatientOrient 0x0018 + (0x5100 << 16) //0018,5100. patient orientation - 'HFS' + fprintf(fpdump, "PatientPosition %s\n", tdicomData->patientOrient); + + // kFlipAngle 0x0018 + (0x1314 << 16) + fprintf(fpdump, "FlipAngle %f\n", tdicomData->flipAngle); + + // kTE 0x0018 + (0x0081 << 16) + fprintf(fpdump, "EchoTime %f\n", tdicomData->TE); + + // kTR 0x0018 + (0x0080 << 16) + fprintf(fpdump, "RepetitionTime %f\n", tdicomData->TR); + + // kTI 0x0018 + (0x0082 << 16) // Inversion time + fprintf(fpdump, "InversionTime %f\n", tdicomData->TI); + + // kPhaseEncodingSteps 0x0018 + (0x0089 << 16) //'IS' + fprintf(fpdump, "NPhaseEnc %d\n", tdicomData->phaseEncodingSteps); + + // kInPlanePhaseEncodingDirection 0x0018 + (0x1312 << 16) //CS + fprintf(fpdump, "PhaseEncDir %c\n", tdicomData->phaseEncodingRC); + + // kZSpacing 0x0018 + (0x0088 << 16) //'DS' 'SpacingBetweenSlices' + fprintf(fpdump, "SliceDistance %f\n", tdicomData->zSpacing); + + // kZThick 0x0018 + (0x0050 << 16) + fprintf(fpdump, "SliceThickness %f\n", tdicomData->zThick); // d.zThick=d.xyzMM[3] + + // kXYSpacing 0x0028 + (0x0030 << 16) //DS 'PixelSpacing' + fprintf(fpdump, "PixelSpacing %f\\%f\n", tdicomData->xyzMM[1], tdicomData->xyzMM[2]); + + // kDim2 0x0028 + (0x0010 << 16) + fprintf(fpdump, "NRows %d\n", tdicomData->xyzDim[2]); + + // kDim1 0x0028 + (0x0011 << 16) + fprintf(fpdump, "NCols %d\n", tdicomData->xyzDim[1]); + + // kBitsAllocated 0x0028 + (0x0100 << 16) + fprintf(fpdump, "BitsPerPixel %d\n", tdicomData->bitsAllocated); + + // dcm2niix doesn't seem to retrieve this + //fprintf(fpdump, "HighBit %d\n",*(e->d.us)); + + // dcm2niix doesn't seem to retrieve this + //fprintf(fpdump, "SmallestValue %d\n",*(e->d.us)); + + // dcm2niix doesn't seem to retrieve this + //fprintf(fpdump, "LargestValue %d\n",*(e->d.us)); + + // kOrientation 0x0020 + (0x0037 << 16) + //fprintf(fpdump, "ImageOrientation %f\\%f\\%f\\%f\\%f\\%f\n", + fprintf(fpdump, "ImageOrientation %g\\%g\\%g\\%g\\%g\\%g\n", + tdicomData->orient[1], tdicomData->orient[2], tdicomData->orient[3], + tdicomData->orient[4], tdicomData->orient[5], tdicomData->orient[6]); // orient[7] ??? + + // kImagePositionPatient 0x0020 + (0x0032 << 16) // Actually ! patientPosition[4] ??? + fprintf(fpdump, "ImagePosition %f\\%f\\%f\n", + tdicomData->patientPosition[1], tdicomData->patientPosition[2], tdicomData->patientPosition[3]); //two values: d.patientPosition, d.patientPositionLast + + // dcm2niix doesn't seem to retrieve this + // kSliceLocation 0x0020+(0x1041 << 16 ) //DS would be useful if not optional type 3 + //case kSliceLocation : //optional so useless, infer from image position patient (0020,0032) and image orientation (0020,0037) + // sliceLocation = dcmStrFloat(lLength, &buffer[lPos]); + // break; + //fprintf(fpdump, "SliceLocation %s\n",e->d.string); + + // kTransferSyntax 0x0002 + (0x0010 << 16) + fprintf(fpdump, "TransferSyntax %s\n", tdicomData->transferSyntax); + + // dcm2niix doesn't seem to retrieve this 0x51, 0x1016 + //fprintf(fpdump, "SiemensCrit %s\n",e->d.string); +} diff --git a/console/dcm2niix_fswrapper.h b/console/dcm2niix_fswrapper.h index 8cb0189f..0a71a8da 100644 --- a/console/dcm2niix_fswrapper.h +++ b/console/dcm2niix_fswrapper.h @@ -17,7 +17,7 @@ class dcm2niix_fswrapper { public: // set TDCMopts defaults, overwrite settings to output in mgz orientation. - static void setOpts(const char* dcmindir, const char* niioutdir=NULL, bool createBIDS=false); + static void setOpts(const char* dcmindir, const char* niioutdir=NULL, bool createBIDS=false, int ForceStackSameSeries=1); // interface to isDICOMfile() in nii_dicom.cpp static bool isDICOM(const char* file); @@ -32,10 +32,18 @@ class dcm2niix_fswrapper // return nifti header saved in MRIFSSTRUCT static nifti_1_header* getNiiHeader(void); + // interface to nii_getMrifsStructVector() + static std::vector* getMrifsStructVector(void); + // return image data saved in MRIFSSTRUCT static const unsigned char* getMRIimg(void); - static void dicomDump(const char* dicomdir); + static void dicomDump(const char* dicomdir, const char *series_info, bool max=false, bool extrainfo=false); + + // + static void seriesInfoDump(FILE *fpdump, const MRIFSSTRUCT *pmrifsStruct); + + static const char *mfrCode2Str(int code); private: static struct TDCMopts tdcmOpts; diff --git a/console/main_console.cpp b/console/main_console.cpp index 539fd276..8865955b 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -80,7 +80,7 @@ void showHelp(const char *argv[], struct TDCMopts opts) { printf(" -ba : anonymize BIDS (y/n, default %c)\n", bool2Char(opts.isAnonymizeBIDS)); printf(" -c : comment stored in NIfTI aux_file (up to 24 characters e.g. '-c VIP', empty to anonymize e.g. 0020,4000 e.g. '-c \"\"')\n"); printf(" -d : directory search depth. Convert DICOMs in sub-folders of in_folder? (0..9, default %d)\n", opts.dirSearchDepth); -#ifdef myEnableJNIfTI +#ifdef myEnableJNIFTI printf(" -e : export as NRRD (y) or MGH (o) or JSON/JNIfTI (j) or BJNIfTI (b) instead of NIfTI (y/n/o/j/b, default n)\n"); #else printf(" -e : export as NRRD (y) or MGH (o) instead of NIfTI (y/n/o/j/b, default n)\n"); @@ -347,6 +347,12 @@ int main(int argc, const char *argv[]) { opts.isAnonymizeBIDS = false; else opts.isAnonymizeBIDS = true; + } else if (argv[i][2] == 'i') { //"-bi M2022" provide BIDS subject ID + i++; + snprintf(opts.bidsSubject, kOptsStr-1, "%s", argv[i]); + } else if (argv[i][2] == 'v') { //"-bv 1222" provide BIDS subject visit + i++; + snprintf(opts.bidsSession, kOptsStr-1, "%s", argv[i]); } else printf("Error: Unknown command line argument: '%s'\n", argv[i]); } else if ((argv[i][1] == 'c') && ((i + 1) < argc)) { @@ -370,7 +376,7 @@ int main(int argc, const char *argv[]) { opts.saveFormat = kSaveFormatJNII; if ((argv[i][0] == 'b') || (argv[i][0] == 'B') || (argv[i][0] == '4')) opts.saveFormat = kSaveFormatBNII; - #ifndef myEnableJNIfTI + #ifndef myEnableJNIFTI if ((opts.saveFormat == kSaveFormatJNII) || (opts.saveFormat == kSaveFormatBNII)) { printf("Recompile for JNIfTI support.\n"); return EXIT_SUCCESS; diff --git a/console/makefile b/console/makefile index 1f35e23c..33131075 100644 --- a/console/makefile +++ b/console/makefile @@ -25,7 +25,7 @@ endif #run "JNIfTI=0 make" to disable JNIFTI build JSFLAGS= ifneq "$(JNIfTI)" "0" - JSFLAGS=-DmyEnableJNIfTI base64.cpp cJSON.cpp + JSFLAGS=-DmyEnableJNIFTI base64.cpp cJSON.cpp endif ifneq ($(OS),Windows_NT) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index a14ea39c..21b496e1 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -37,7 +37,9 @@ #include #include // discriminate files from folders #include +#ifndef USING_DCM2NIIXFSWRAPPER #include +#endif #ifdef USING_R #undef isnan @@ -893,6 +895,7 @@ struct TDICOMdata clear_dicom_data() { d.isValid = false; d.isXRay = false; d.isMultiEcho = false; + d.numberOfTR = 1; d.isSigned = false; //default is unsigned! d.isFloat = false; //default is for integers, not single or double precision d.isResampled = false; //assume data not resliced to remove gantry tilt problems @@ -951,6 +954,9 @@ struct TDICOMdata clear_dicom_data() { d.CSA.SeriesHeader_offset = 0; d.CSA.SeriesHeader_length = 0; d.CSA.coilNumber = -1; + strcpy(d.CSA.bidsDataType, ""); + strcpy(d.CSA.bidsEntitySuffix, ""); + strcpy(d.CSA.bidsTask, ""); return d; } //clear_dicom_data() @@ -1462,7 +1468,7 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i CSA->coilNumber = csaICEdims(&buff[lPos]); else if (strcmp(tagCSA.name, "NumberOfImagesInMosaic") == 0) CSA->mosaicSlices = (int)round(csaMultiFloat(&buff[lPos], 1, lFloats, &itemsOK)); - else if (strcmp(tagCSA.name, "B_value") == 0) { + else if (strcmp(tagCSA.name, "B_value") == 0) { CSA->dtiV[0] = csaMultiFloat(&buff[lPos], 1, lFloats, &itemsOK); if (CSA->dtiV[0] < 0.0) { printWarning("(Corrupt) CSA reports negative b-value! %g\n", CSA->dtiV[0]); @@ -1798,7 +1804,7 @@ void cleanStr(char *lOut) { int isSameFloatGE(float a, float b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! - //return (a == b); //niave approach does not have any tolerance for rounding errors + //return (a == b); //naive approach does not have any tolerance for rounding errors return (fabs(a - b) <= 0.0001); } @@ -2042,6 +2048,9 @@ struct TDICOMdata nii_readParRec(char *parname, int isVerbose, struct TDTI4D *dt strcat(d.seriesDescription, Comment[6]); strcat(d.seriesDescription, Comment[7]); cleanStr(d.seriesDescription); +#ifdef USING_DCM2NIIXFSWRAPPER + remove_specialchars(d.seriesDescription); +#endif } if ((strcmp(Comment[0], "Examination") == 0) && (strcmp(Comment[1], "date/time") == 0)) { if ((strlen(Comment[3]) >= 10) && (strlen(Comment[5]) >= 8)) { @@ -2362,7 +2371,12 @@ struct TDICOMdata nii_readParRec(char *parname, int isVerbose, struct TDTI4D *dt //offset images by type: mag+0,real+1, imag+2,phase+3 //if (cols[kImageType] != 0) //yikes - phase maps! // slice = slice + numExpected; - maxSlice2D = std::max(slice, maxSlice2D); +#ifdef USING_DCM2NIIXFSWRAPPER + // fix ubuntu18 compiler error + maxSlice2D = (slice > maxSlice2D) ? slice : maxSlice2D; +#else + maxSlice2D = std::max(slice, maxSlice2D); +#endif if ((slice >= 0) && (slice < kMaxSlice2D) && (numSlice2D < kMaxSlice2D) && (numSlice2D >= 0)) { dti4D->sliceOrder[slice - 1] = numSlice2D; //printMessage("%d\t%d\t%d\n", numSlice2D, slice, (int)cols[kSlice],(int)vol); @@ -2428,7 +2442,12 @@ struct TDICOMdata nii_readParRec(char *parname, int isVerbose, struct TDTI4D *dt if (d.isHasImaginary) nType ++; if (d.isHasPhase) nType ++; if (d.isHasReal) nType ++; - nType = std::max(nType, 1); +#ifdef USING_DCM2NIIXFSWRAPPER + // fix ubuntu18 compiler error + nType = (nType > 1) ? nType : 1; +#else + nType = std::max(nType, 1); +#endif if (slice != numSlice2D) { printError("Catastrophic error: found %d but expected %d slices. %s\n", slice, numSlice2D, parname); printMessage(" slices*grad*bval*cardiac*echo*dynamic*mix*labels*types = %d*%d*%d*%d*%d*%d*%d*%d*%d\n", @@ -2590,7 +2609,7 @@ struct TDICOMdata nii_readParRec(char *parname, int isVerbose, struct TDTI4D *dt y.v[2] = (float)d.xyzDim[3] - 1.0f; y.v[3] = 1.0f; y = nifti_vect44mat44_mul(y, R44); - int iOri = 2; //for axial, slices are 3rd dimenson (indexed from 0) (k) + int iOri = 2; //for axial, slices are 3rd dimension (indexed from 0) (k) if (d.sliceOrient == kSliceOrientSag) iOri = 0; //for sagittal, slices are 1st dimension (i) if (d.sliceOrient == kSliceOrientCor) @@ -3133,12 +3152,25 @@ unsigned char *nii_byteswap(unsigned char *img, struct nifti_1_header *hdr) { #ifdef myEnableJasper unsigned char *nii_loadImgCoreJasper(char *imgname, struct nifti_1_header hdr, struct TDICOMdata dcm, int compressFlag) { +#if defined(JAS_VERSION_MAJOR) && JAS_VERSION_MAJOR >= 3 + jas_conf_clear(); + jas_conf_set_debug_level(0); + jas_conf_set_max_mem_usage(1 << 30); // 1 GiB + if (jas_init_library()) { + return NULL; + } + if (jas_init_thread()) { + jas_cleanup_library(); + return NULL; + } +#else if (jas_init()) { return NULL; } + jas_setdbglevel(0); +#endif jas_stream_t *in; jas_image_t *image; - jas_setdbglevel(0); if (!(in = jas_stream_fopen(imgname, "rb"))) { printError("Cannot open input image file %s\n", imgname); return NULL; @@ -3265,6 +3297,10 @@ unsigned char *nii_loadImgCoreJasper(char *imgname, struct nifti_1_header hdr, s jas_matrix_destroy(data); jas_image_destroy(image); jas_image_clearfmts(); +#if defined(JAS_VERSION_MAJOR) && JAS_VERSION_MAJOR >= 3 + jas_cleanup_thread(); + jas_cleanup_library(); +#endif return img; } //nii_loadImgCoreJasper() #endif @@ -3801,14 +3837,19 @@ unsigned char *nii_loadImgXL(char *imgname, struct nifti_1_header *hdr, struct T return img; } //nii_loadImgXL() -int isSQ(uint32_t groupElement) { //Detect sequence VR ("SQ") for implicit tags +int isSQ(uint32_t groupElement, bool isPhilips) { //Detect sequence VR ("SQ") for implicit tags static const int array_size = 35; - uint32_t array[array_size] = {0x2005 + (uint32_t(0x140F) << 16), 0x0008 + (uint32_t(0x1111) << 16), 0x0008 + (uint32_t(0x1115) << 16), 0x0008 + (uint32_t(0x1140) << 16), 0x0008 + (uint32_t(0x1199) << 16), 0x0008 + (uint32_t(0x2218) << 16), 0x0008 + (uint32_t(0x9092) << 16), 0x0018 + (uint32_t(0x9006) << 16), 0x0018 + (uint32_t(0x9042) << 16), 0x0018 + (uint32_t(0x9045) << 16), 0x0018 + (uint32_t(0x9049) << 16), 0x0018 + (uint32_t(0x9112) << 16), 0x0018 + (uint32_t(0x9114) << 16), 0x0018 + (uint32_t(0x9115) << 16), 0x0018 + (uint32_t(0x9117) << 16), 0x0018 + (uint32_t(0x9119) << 16), 0x0018 + (uint32_t(0x9125) << 16), 0x0018 + (uint32_t(0x9152) << 16), 0x0018 + (uint32_t(0x9176) << 16), 0x0018 + (uint32_t(0x9226) << 16), 0x0018 + (uint32_t(0x9239) << 16), 0x0020 + (uint32_t(0x9071) << 16), 0x0020 + (uint32_t(0x9111) << 16), 0x0020 + (uint32_t(0x9113) << 16), 0x0020 + (uint32_t(0x9116) << 16), 0x0020 + (uint32_t(0x9221) << 16), 0x0020 + (uint32_t(0x9222) << 16), 0x0028 + (uint32_t(0x9110) << 16), 0x0028 + (uint32_t(0x9132) << 16), 0x0028 + (uint32_t(0x9145) << 16), 0x0040 + (uint32_t(0x0260) << 16), 0x0040 + (uint32_t(0x0555) << 16), 0x0040 + (uint32_t(0xa170) << 16), 0x5200 + (uint32_t(0x9229) << 16), 0x5200 + (uint32_t(0x9230) << 16)}; + uint32_t array[array_size] = {0x0008 + (uint32_t(0x1111) << 16), 0x0008 + (uint32_t(0x1115) << 16), 0x0008 + (uint32_t(0x1140) << 16), 0x0008 + (uint32_t(0x1199) << 16), 0x0008 + (uint32_t(0x2218) << 16), 0x0008 + (uint32_t(0x9092) << 16), 0x0018 + (uint32_t(0x9006) << 16), 0x0018 + (uint32_t(0x9042) << 16), 0x0018 + (uint32_t(0x9045) << 16), 0x0018 + (uint32_t(0x9049) << 16), 0x0018 + (uint32_t(0x9112) << 16), 0x0018 + (uint32_t(0x9114) << 16), 0x0018 + (uint32_t(0x9115) << 16), 0x0018 + (uint32_t(0x9117) << 16), 0x0018 + (uint32_t(0x9119) << 16), 0x0018 + (uint32_t(0x9125) << 16), 0x0018 + (uint32_t(0x9152) << 16), 0x0018 + (uint32_t(0x9176) << 16), 0x0018 + (uint32_t(0x9226) << 16), 0x0018 + (uint32_t(0x9239) << 16), 0x0020 + (uint32_t(0x9071) << 16), 0x0020 + (uint32_t(0x9111) << 16), 0x0020 + (uint32_t(0x9113) << 16), 0x0020 + (uint32_t(0x9116) << 16), 0x0020 + (uint32_t(0x9221) << 16), 0x0020 + (uint32_t(0x9222) << 16), 0x0028 + (uint32_t(0x9110) << 16), 0x0028 + (uint32_t(0x9132) << 16), 0x0028 + (uint32_t(0x9145) << 16), 0x0040 + (uint32_t(0x0260) << 16), 0x0040 + (uint32_t(0x0555) << 16), 0x0040 + (uint32_t(0xa170) << 16), 0x5200 + (uint32_t(0x9229) << 16), 0x5200 + (uint32_t(0x9230) << 16)}; + //uint32_t array[array_size] = {0x2005 + (uint32_t(0x140F) << 16), 0x0008 + (uint32_t(0x1111) << 16), 0x0008 + (uint32_t(0x1115) << 16), 0x0008 + (uint32_t(0x1140) << 16), 0x0008 + (uint32_t(0x1199) << 16), 0x0008 + (uint32_t(0x2218) << 16), 0x0008 + (uint32_t(0x9092) << 16), 0x0018 + (uint32_t(0x9006) << 16), 0x0018 + (uint32_t(0x9042) << 16), 0x0018 + (uint32_t(0x9045) << 16), 0x0018 + (uint32_t(0x9049) << 16), 0x0018 + (uint32_t(0x9112) << 16), 0x0018 + (uint32_t(0x9114) << 16), 0x0018 + (uint32_t(0x9115) << 16), 0x0018 + (uint32_t(0x9117) << 16), 0x0018 + (uint32_t(0x9119) << 16), 0x0018 + (uint32_t(0x9125) << 16), 0x0018 + (uint32_t(0x9152) << 16), 0x0018 + (uint32_t(0x9176) << 16), 0x0018 + (uint32_t(0x9226) << 16), 0x0018 + (uint32_t(0x9239) << 16), 0x0020 + (uint32_t(0x9071) << 16), 0x0020 + (uint32_t(0x9111) << 16), 0x0020 + (uint32_t(0x9113) << 16), 0x0020 + (uint32_t(0x9116) << 16), 0x0020 + (uint32_t(0x9221) << 16), 0x0020 + (uint32_t(0x9222) << 16), 0x0028 + (uint32_t(0x9110) << 16), 0x0028 + (uint32_t(0x9132) << 16), 0x0028 + (uint32_t(0x9145) << 16), 0x0040 + (uint32_t(0x0260) << 16), 0x0040 + (uint32_t(0x0555) << 16), 0x0040 + (uint32_t(0xa170) << 16), 0x5200 + (uint32_t(0x9229) << 16), 0x5200 + (uint32_t(0x9230) << 16)}; for (int i = 0; i < array_size; i++) { //if (array[i] == groupElement) printMessage(" implicitSQ %04x,%04x\n", groupElement & 65535,groupElement>>16); if (array[i] == groupElement) return 1; } + //issue 775: Canon/Toshiba use 2005,140f but not as a SQ? + uint32_t kPhilipsSQ = 0x2005 + (uint32_t(0x140F) << 16); + if ((isPhilips) && (kPhilipsSQ == groupElement)) + return 1; return 0; } //isSQ() @@ -3963,6 +4004,13 @@ void dcmMultiFloatDouble(size_t lByteLength, unsigned char lBuffer[], size_t lnF lFloats[i] = dcmFloatDouble((int)floatlen, lBuffer + i * floatlen, isLittleEndian); } //dcmMultiFloatDouble() +void dcmMultiFloatSingle(size_t lByteLength, unsigned char lBuffer[], size_t lnFloats, float *lFloats, bool isLittleEndian) { +// dcmFloat kTRArray + size_t floatlen = lByteLength / lnFloats; + for (size_t i = 0; i < lnFloats; ++i) + lFloats[i] = dcmFloat((int)floatlen, lBuffer + i * floatlen, isLittleEndian); +} //dcmMultiFloatSingle() + void set_orientation0018_9089(struct TVolumeDiffusion *ptvd, int lLength, unsigned char *inbuf, bool isLittleEndian) { if (ptvd->_isPhilipsNonDirectional) { for (int i = 1; i < 4; ++i) // Deliberately ignore inbuf; it might be nonsense. @@ -4014,7 +4062,13 @@ void _update_tvd(struct TVolumeDiffusion *ptvd) { return; //no B=0 if (isReady) { for (int i = 1; i < 4; ++i) { - if (ptvd->_dtiV[i] > 1) { + // Check that _dtiV is not still at its default of [-1, 2, 2, 2] from + // clear_volume(struct TVolumeDiffusion *ptvd). This is mostly avoided + // because of the ptvd->_dtiV[0] >= 0 check above, but was supposedly + // needed at some point. + // issue 769: some bvec components may be slightly (~5%) > 1. AFAIK, + // the relevant value to guard against would be +2. + if (ptvd->_dtiV[i] > 1.5) { isReady = false; break; } @@ -4100,7 +4154,7 @@ void _update_tvd(struct TVolumeDiffusion *ptvd) { struct TDCMdim { //DimensionIndexValues uint32_t dimIdx[MAX_NUMBER_OF_DIMENSIONS]; uint32_t diskPos; - float triggerDelayTime, TE, intenScale, intenIntercept, intenScalePhilips, RWVScale, RWVIntercept; + float triggerDelayTime, TE, TR, intenScale, intenIntercept, intenScalePhilips, RWVScale, RWVIntercept; float V[4]; bool isPhase; bool isReal; @@ -4302,6 +4356,7 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D #define kSeriesTime 0x0008 + (0x0031 << 16) #define kAcquisitionTime 0x0008 + (0x0032 << 16) //TM //#define kContentTime 0x0008+(0x0033 << 16 ) //TM +#define kAccessionNumber 0x0008 + (0x0050 << 16) #define kModality 0x0008 + (0x0060 << 16) //CS #define kManufacturer 0x0008 + (0x0070 << 16) #define kInstitutionName 0x0008 + (0x0080 << 16) @@ -4315,9 +4370,9 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D #define kReferencedImageEvidenceSQ (uint32_t)0x0008 + (0x9092 << 16) #define kComplexImageComponent (uint32_t)0x0008 + (0x9208 << 16) //'0008' '9208' 'CS' 'ComplexImageComponent' #define kAcquisitionContrast (uint32_t)0x0008 + (0x9209 << 16) //'0008' '9209' 'CS' 'AcquisitionContrast' +#define kIconSQ 0x0009 + (0x1110 << 16) #define kPatientName 0x0010 + (0x0010 << 16) #define kPatientID 0x0010 + (0x0020 << 16) -#define kAccessionNumber 0x0008 + (0x0050 << 16) #define kPatientBirthDate 0x0010 + (0x0030 << 16) #define kPatientSex 0x0010 + (0x0040 << 16) #define kPatientAge 0x0010 + (0x1010 << 16) @@ -4338,7 +4393,6 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D #define kTI 0x0018 + (0x0082 << 16) // Inversion time #define kNumberOfAverages 0x0018 + (0x0083 << 16) //DS #define kImagingFrequency 0x0018 + (0x0084 << 16) //DS -//#define kEffectiveTE 0x0018+(0x9082 << 16 ) #define kEchoNum 0x0018 + (0x0086 << 16) //IS #define kMagneticFieldStrength 0x0018 + (0x0087 << 16) //DS #define kZSpacing 0x0018 + (0x0088 << 16) //'DS' 'SpacingBetweenSlices' @@ -4386,6 +4440,7 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D #define kPartialFourier 0x0018 + uint32_t(0x9081 << 16) //CS const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //#define kDiffusionBFactorSiemens 0x0019+(0x100C<< 16 ) // 0019;000C;SIEMENS MR HEADER;B_value +#define kDiffusionGradientDirectionSQ 0x0018 + uint32_t(0x9076 << 16) // SQ #define kDiffusion_bValue 0x0018 + uint32_t(0x9087 << 16) // FD #define kDiffusionOrientation 0x0018 + uint32_t(0x9089 << 16) // FD, seen in enhanced DICOM from Philips 5.* and Siemens XA10. #define kImagingFrequencyFD 0x0018 + uint32_t(0x9098 << 16) //FD @@ -4454,6 +4509,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD #define kScanOptionsSiemens 0x0021 + (0x105C << 16) //CS Siemens ONLY #define kPATModeText 0x0021 + (0x1009 << 16) //LO, see kImaPATModeText #define kCSASeriesHeaderInfoXA 0x0021 + (0x1019 << 16) +#define kCSASeriesHeaderInfoXA2 0x0021 + (0x11FE << 16) #define kTimeAfterStart 0x0021 + (0x1104 << 16) //DS #define kICE_dims 0x0021 + (0x1106 << 16) //LO [X_4_1_1_1_1_160_1_1_1_1_1_277] #define kPhaseEncodingDirectionPositiveSiemens 0x0021 + (0x111C << 16) //IS @@ -4461,6 +4517,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //#define kPATModeText2 0x0021 + (0x1156 << 16) //LO, always same as 0021,1009 #define kBandwidthPerPixelPhaseEncode21 0x0021 + (0x1153 << 16) //FD #define kCoilElements 0x0021 + (0x114F << 16) //LO +#define kICE_dimsLong 0x0021 + (0x118e << 16) #define kAcquisitionMatrixText21 0x0021 + (0x1158 << 16) //SH #define kImageTypeText 0x0021 + (0x1175 << 16) //CS #define kDeepLearningText 0x0021 + (0x1176 << 16) //LO @@ -4529,6 +4586,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD // #define kSeriesType 0x0054+(0x1000 << 16 ) #define kDoseCalibrationFactor 0x0054 + (0x1322 << 16) #define kPETImageIndex 0x0054 + (0x1330 << 16) +#define kReferencedSegmentNumber 0x0062 + (0x000B << 16) //US #define kPEDirectionDisplayedUIH 0x0065 + (0x1005 << 16) //SH #define kDiffusion_bValueUIH 0x0065 + (0x1009 << 16) //FD #define kParallelInformationUIH 0x0065 + (0x100D << 16) //SH @@ -4558,6 +4616,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //#define kStackSliceNumber 0x2001+(0x1035 << 16 )//? Potential way to determine slice order for Philips? #define kNumberOfDynamicScans 0x2001 + (0x1081 << 16) //'2001' '1081' 'IS' 'NumberOfDynamicScans' //#define kTRPhilips 0x2005 + (0x1030 << 16) //(2005,1030) FL 30\150 +#define kTRArray 0x2005 + (0x1030 << 16) //FL Array, e.g. "15\75", see issue 743 #define kMRfMRIStatusIndicationPhilips 0x2005 + (0x1063 << 16) #define kMRAcquisitionTypePhilips 0x2005 + (0x106F << 16) //SS #define kAngulationAP 0x2005 + (0x1071 << 16) //'2005' '1071' 'FL' 'MRStackAngulationAP' @@ -4645,6 +4704,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD int temporalPositionIndex = 0; int minTemporalPositionIndex = 65535; int maxTemporalPositionIndex = 0; + int minReferencedSegmentNumber = 65535; + int maxReferencedSegmentNumber = 0; //int temporalPositionIdentifier = 0; int locationsInAcquisitionPhilips = 0; int imagesInAcquisition = 0; @@ -4654,6 +4715,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD int volumeNumber = -1; int gradientOrientationNumberPhilips = -1; int numberOfFrames = 0; + int numberOfFramesICEdims = 0; //int MRImageGradientOrientationNumber = 0; //int minGradNum = kMaxDTI4D + 1; //int maxGradNum = -1; @@ -4689,6 +4751,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD char imageType1st[kDICOMStr] = ""; bool isEncapsulatedData = false; int diffusionDirectionTypeGE = 0; //issue690 + int seriesdiffusionDirectionTypeGE = 0; //issue690, 777 int multiBandFactor = 0; int frequencyRows = 0; int numberOfImagesInMosaic = 0; @@ -4819,7 +4882,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD if (sqDepth < 0) sqDepth = 0; //should not happen, but protect for faulty anonymization //if we leave the folder MREchoSequence 0018,9114 - if ((nDimIndxVal > 0) && ((d.manufacturer == kMANUFACTURER_CANON) || (d.manufacturer == kMANUFACTURER_BRUKER) || (d.manufacturer == kMANUFACTURER_PHILIPS)) && (sqDepth00189114 >= sqDepth)) { + if ((nDimIndxVal > 0) && ((d.manufacturer == kMANUFACTURER_UNKNOWN) || (d.manufacturer == kMANUFACTURER_CANON) || (d.manufacturer == kMANUFACTURER_BRUKER) || (d.manufacturer == kMANUFACTURER_PHILIPS)) && (sqDepth00189114 >= sqDepth)) { sqDepth00189114 = -1; //triggered //printf("slice %d---> 0020,9157 = %d %d %d\n", inStackPositionNumber, d.dimensionIndexValues[0], d.dimensionIndexValues[1], d.dimensionIndexValues[2]); // d.aslFlags = kASL_FLAG_PHILIPS_LABEL; kASL_FLAG_PHILIPS_LABEL @@ -4839,7 +4902,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD } if ((volumeNumber == 1) && (acquisitionTimePhilips >= 0.0) && (inStackPositionNumber > 0)) { d.CSA.sliceTiming[inStackPositionNumber - 1] = acquisitionTimePhilips; - printf("%d\t%f\n", inStackPositionNumber, acquisitionTimePhilips); + //printf("%d\t%f\n", inStackPositionNumber, acquisitionTimePhilips); acquisitionTimePhilips = - 1.0; } int ndim = nDimIndxVal; @@ -4884,7 +4947,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD for (int i = 0; i < ndim; i++) dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[dimensionIndexOrder[i]]; } - dcmDim[numDimensionIndexValues].TE = TE; + dcmDim[numDimensionIndexValues].TE = d.TE; + dcmDim[numDimensionIndexValues].TR = d.TR; dcmDim[numDimensionIndexValues].intenScale = d.intenScale; dcmDim[numDimensionIndexValues].intenIntercept = d.intenIntercept; dcmDim[numDimensionIndexValues].isPhase = isPhase; @@ -4972,14 +5036,15 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD nDimIndxVal = -1; //we need DimensionIndexValues } //record dimensionIndexValues slice information } //groupElement == kItemDelimitationTag : delimit item exits folder + uint32_t itemTagLength = 0; if (groupElement == kItemTag) { - uint32_t slen = dcmInt(4, &buffer[lPos], d.isLittleEndian); + itemTagLength = dcmInt(4, &buffer[lPos], d.isLittleEndian); uint32_t kUndefinedLen = 0xFFFFFFFF; - if (slen != kUndefinedLen) { + if (itemTagLength != kUndefinedLen) { nNestPos++; if (nNestPos >= kMaxNestPost) nNestPos = kMaxNestPost - 1; - nestPos[nNestPos] = slen + lFileOffset + lPos; + nestPos[nNestPos] = itemTagLength + lFileOffset + lPos; } lLength = 4; sqDepth++; @@ -5023,17 +5088,25 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD else lLength = buffer[lPos + 3] | (buffer[lPos + 2] << 8) | (buffer[lPos + 1] << 16) | (buffer[lPos] << 24); lPos += 4; //we have loaded the 32-bit length - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isSQ(groupElement))) { //https://github.com/rordenlab/dcm2niix/issues/144 + if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isSQ(groupElement, true))) { //https://github.com/rordenlab/dcm2niix/issues/144 vr[0] = 'S'; vr[1] = 'Q'; lLength = 0; //Do not skip kItemTag - required to determine nesting of Philips Enhanced } - if ((d.manufacturer != kMANUFACTURER_PHILIPS) && (isSQ(groupElement))) { //https://github.com/rordenlab/dcm2niix/issues/144 + if ((d.manufacturer != kMANUFACTURER_PHILIPS) && (isSQ(groupElement, false))) { //https://github.com/rordenlab/dcm2niix/issues/144 + vr[0] = 'S'; + vr[1] = 'Q'; + lLength = 0; //Do not skip kItemTag - required to determine nesting of Philips Enhanced + } + if (groupElement == kDiffusionGradientDirectionSQ) { //https://github.com/rordenlab/dcm2niix/issues/144 vr[0] = 'S'; vr[1] = 'Q'; lLength = 0; //Do not skip kItemTag - required to determine nesting of Philips Enhanced } } //if explicit else implicit VR + if ((lLength == 0xFFFFFFFF) && (vr[0] == 'O') && (vr[1] == 'B') && (isIconImageSequence)) { + lLength = 0; //encapuslated data of unspecified length + } //issue713 if (lLength == 0xFFFFFFFF) { lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length //09032018 - do not count these as SQs: Horos does not count even groups @@ -5044,6 +5117,9 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD vr[1] = 'Q'; //} } + if ((groupElement == kItemTag) && (vr[0] == 'O') && (vr[1] == 'B') && (isIconImageSequence)) + lLength += itemTagLength; //issue713 + //printf("issue713::%c%c %04x,%04x %d@%lu\n", vr[0], vr[1], groupElement & 65535, groupElement >> 16, lLength, lPos); if ((groupElement == kItemTag) && (isEncapsulatedData)) { //use this to find image fragment for compressed datasets, e.g. JPEG transfer syntax d.imageBytes = dcmInt(4, &buffer[lPos], d.isLittleEndian); lPos = lPos + 4; @@ -5213,6 +5289,9 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD char transferSyntax[kDICOMStr]; strcpy(transferSyntax, ""); dcmStr(lLength, &buffer[lPos], transferSyntax); +#ifdef USING_DCM2NIIXFSWRAPPER + strcpy(d.transferSyntax, transferSyntax); +#endif if (strcmp(transferSyntax, "1.2.840.10008.1.2.1") == 0) ; //default isExplicitVR=true; //d.isLittleEndian=true else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.50") == 0) { @@ -5510,6 +5589,13 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD if (strlen(d.studyTime) < 2) dcmStr(lLength, &buffer[lPos], d.studyTime); break; + case kIconSQ: {//issue760 GEIIS icon strikes again + if ((vr[0] != 'S') || (vr[1] != 'Q') || (lLength != 8)) break; //only risk kludge for explicit VR with length + isIconImageSequence = true; + if (sqDepthIcon < 0) + sqDepthIcon = sqDepth; + break; + } case kPatientName: dcmStr(lLength, &buffer[lPos], d.patientName); break; @@ -5562,6 +5648,9 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD break; case kSeriesDescription: dcmStr(lLength, &buffer[lPos], d.seriesDescription); +#ifdef USING_DCM2NIIXFSWRAPPER + remove_specialchars(d.seriesDescription); +#endif break; case kInstitutionalDepartmentName: dcmStr(lLength, &buffer[lPos], d.institutionalDepartmentName); @@ -5590,6 +5679,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD d.isXA10A = true; if ((slen > 4) && (strstr(d.softwareVersions, "XA30") != NULL)) d.isXA10A = true; + if ((slen > 4) && (strstr(d.softwareVersions, "XA31") != NULL)) + d.isXA10A = true; if ((slen < 5) || (strstr(d.softwareVersions, "XA10") == NULL)) break; d.isXA10A = true; @@ -5730,9 +5821,13 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD case kMRAcquisitionPhaseEncodingStepsOutOfPlane: d.phaseEncodingStepsOutOfPlane = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); break; - case kGradientEchoTrainLength: - d.echoTrainLength = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); + case kGradientEchoTrainLength: { + int etl = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); + if (etl < 2) //issue 717 + break; + d.echoTrainLength = etl; break; + } case kNumberOfImagesInMosaic: if (d.manufacturer == kMANUFACTURER_SIEMENS) numberOfImagesInMosaic = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); @@ -5857,6 +5952,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD break; dcmStr(lLength, &buffer[lPos], d.pulseSequenceName); if (strstr( d.pulseSequenceName, "epi_pepolar") != NULL) { + //if ((strstr( d.pulseSequenceName, "epi_pepolar") != NULL) || (strstr( d.pulseSequenceName, "epi2_pepolar") != NULL)){ d.epiVersionGE = kGE_EPI_PEPOLAR_FWD; //n.b. combine with 0019,10B3 } else if (strstr( d.pulseSequenceName, "epi2") != NULL) { d.epiVersionGE = kGE_EPI_EPI2; //-1 = not epi, 0 = epi, 1 = epiRT, 2 = epi2 @@ -5865,7 +5961,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD } else if (strstr( d.pulseSequenceName, "epi") != NULL) { d.epiVersionGE = kGE_EPI_EPI; //-1 = not epi, 0 = epi, 1 = epiRT } - if (strcmp( d.pulseSequenceName, "3db0map") == 0) { + if (strstr( d.pulseSequenceName, "b0map")) { isGEfieldMap = true; //issue501 } break; @@ -5875,6 +5971,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD break; char epiStr[kDICOMStr]; dcmStr(lLength, &buffer[lPos], epiStr); + strcat(d.phaseEncodingDirectionDisplayedUIH, epiStr); //overload if ((d.epiVersionGE < kGE_EPI_PEPOLAR_FWD) && (strcmp(epiStr, "EPI") == 0)) { d.internalepiVersionGE = 1; //-1 = not EPI, 1 = EPI, 2 = EPI2 if (d.epiVersionGE != 1) { // 1 = epiRT by kEpiRTGroupDelayGE or kPulseSequenceNameGE @@ -5896,7 +5993,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD dcmStr(lLength, &buffer[lPos], d.studyInstanceUID); break; case kSeriesInstanceUID: // 0020,000E - if (sqDepth > seriesInstanceUIDsqDepth) break; //issue655 + if (sqDepth > seriesInstanceUIDsqDepth) break; //issue655 seriesInstanceUIDsqDepth = sqDepth; dcmStr(lLength, &buffer[lPos], d.seriesInstanceUID); d.seriesUidCrc = mz_crc32X((unsigned char *)&d.seriesInstanceUID, strlen(d.seriesInstanceUID)); @@ -5989,6 +6086,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD minTemporalPositionIndex = temporalPositionIndex; break; case kDimensionIndexPointer: + if (dimensionIndexPointerCounter >= MAX_NUMBER_OF_DIMENSIONS) break; dimensionIndexPointer[dimensionIndexPointerCounter++] = dcmAttributeTag(&buffer[lPos], d.isLittleEndian); break; case kFrameContentSequence: @@ -6078,6 +6176,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //printMessage("p%gs%d\n", d.accelFactPE, multiBandFactor); break; } +// case kCSASeriesHeaderInfoXA2: +// printf("do something profound\n"); case kCSASeriesHeaderInfoXA: if (d.manufacturer != kMANUFACTURER_SIEMENS) break; @@ -6123,7 +6223,22 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //printMessage("%d:%d:'%s'\n", d.echoNum, echo, iceStr); if (iceStr != end) d.echoNum = echo; - //printMessage("%d:'%s'\n", echo, echoStr); + break; + } + case kICE_dimsLong: { //issue568: LO (0021,1106) [X_4_1_1_1_1_160_1_1_1_1_1_277] + if ((d.manufacturer != kMANUFACTURER_SIEMENS) || (numberOfFramesICEdims > 0)) + break; + char iceStr[kDICOMStr]; + dcmStr(lLength, &buffer[lPos], iceStr); + int idx = 0; + char * pch = strtok (iceStr,"_"); + char *end; + while (pch != NULL) { + if (idx ==20) + numberOfFramesICEdims = (int)strtol(pch, &end, 10); + idx ++; + pch = strtok (NULL, "_"); + } break; } case kPhaseEncodingDirectionPositiveSiemens: { @@ -6189,6 +6304,14 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD case kPETImageIndex: PETImageIndex = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); break; + case kReferencedSegmentNumber: { //US issue706 + int segn = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); + if (segn < minReferencedSegmentNumber) + minReferencedSegmentNumber = segn; + if (segn > maxReferencedSegmentNumber) + maxReferencedSegmentNumber = segn; + break; + } case kPEDirectionDisplayedUIH: if (d.manufacturer != kMANUFACTURER_UIH) break; @@ -6555,8 +6678,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD case kScanningSequenceSiemens: if (d.manufacturer == kMANUFACTURER_SIEMENS) dcmStr(lLength, &buffer[lPos], scanningSequenceSiemens); - if (d.manufacturer == kMANUFACTURER_GE) //issue690 - diffusionDirectionTypeGE = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); + if (d.manufacturer == kMANUFACTURER_GE) //issue690, series-level 16=DFAXDTI + seriesdiffusionDirectionTypeGE = dcmInt(lLength, &buffer[lPos], d.isLittleEndian); break; case kSequenceVariant21: if (d.manufacturer != kMANUFACTURER_SIEMENS) @@ -6575,6 +6698,23 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD dcmStr(lLength, &buffer[lPos], d.sequenceName); break; } + case kTRArray: { + if (d.manufacturer != kMANUFACTURER_PHILIPS) + break; + //support multiple TRs: issue 743 + const int kMax = 4; + float v[kMax] = { 0.0 }; + int nFloat32 = (lLength / 4); + if (nFloat32 < 2) break; + if (nFloat32 > kMax) + nFloat32 = kMax; + dcmMultiFloatSingle(nFloat32 * 4, &buffer[lPos], nFloat32, v, d.isLittleEndian); + d.numberOfTR = 0; + //count number of TRs > 0, see issue 749 + for(int i = 0; i< nFloat32; i++) + if (v[i] > 0.0) + d.numberOfTR ++; + } case kMRfMRIStatusIndicationPhilips: {//fmri volume number if (d.manufacturer != kMANUFACTURER_PHILIPS) break; @@ -6946,7 +7086,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD break; d.shimGradientZ = dcmIntSS(lLength, &buffer[lPos], d.isLittleEndian); break; - case kVasCollapseFlagGE: //SS issue 690 16=DiffusionDtiDicomValue + case kVasCollapseFlagGE: //SS issue 690 image-level 16=DiffusionDtiDicomValue or 14=DiffusionT2DicomValue (initial b0) if (d.manufacturer != kMANUFACTURER_GE) break; diffusionDirectionTypeGE = dcmIntSS(lLength, &buffer[lPos], d.isLittleEndian); @@ -7070,6 +7210,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD if (d.manufacturer == kMANUFACTURER_GE) { d.CSA.dtiV[0] = (float)set_bValGE(&volDiffusion, lLength, &buffer[lPos]); d.CSA.numDti = 1; + d.isDiffusion = true; } break; case kEpiRTGroupDelayGE: //FL @@ -7154,7 +7295,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD isIconImageSequence = true; if (sqDepthIcon < 0) sqDepthIcon = sqDepth; - //geiisBug = true; //compressed thumbnails do not follow transfer syntax! GE should not re-use pulbic tags for these proprietary images http://sonca.kasshin.net/gdcm/Doc/GE_ImageThumbnails + //geiisBug = true; //compressed thumbnails do not follow transfer syntax! GE should not reuse pulbic tags for these proprietary images http://sonca.kasshin.net/gdcm/Doc/GE_ImageThumbnails } break; case kStudyComments: { @@ -7454,9 +7595,10 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD } if (isGEfieldMap) { //issue501 : to do check zip factor //Volume 1) derived phase field map [Hz] and 2) magnitude volume. - d.isDerived = (d.imageNum <= locationsInAcquisitionGE); //first volume - d.isRealIsPhaseMapHz = d.isDerived; - d.isHasReal = d.isDerived; + //issue 777 while a fieldmap is technically derived, do not exclude with -i y + bool isDerived = (d.imageNum <= locationsInAcquisitionGE); //first volume + d.isRealIsPhaseMapHz = isDerived; + d.isHasReal = isDerived; } /* SAH.end */ if (locationsInAcquisitionGE < d.locationsInAcquisition) { @@ -7553,6 +7695,11 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD d.xyzDim[3] = d.xyzDim[3] / numberOfDynamicScans; d.xyzDim[4] = numberOfDynamicScans; } + int nSegment = maxReferencedSegmentNumber - minReferencedSegmentNumber + 1; + if ((nSegment > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % nSegment) == 0)) { //issue706 + d.xyzDim[3] = d.xyzDim[3] / nSegment; + d.xyzDim[4] = nSegment; + } if ((maxInStackPositionNumber > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % maxInStackPositionNumber) == 0)) { d.xyzDim[4] = d.xyzDim[3] / maxInStackPositionNumber; d.xyzDim[3] = maxInStackPositionNumber; @@ -7605,7 +7752,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD qsort(objects, numberOfFrames, sizeof(struct fidx), fcmp); numDimensionIndexValues = numberOfFrames; for (int i = 0; i < numberOfFrames; i++) { - // printf("%d > %g\n", objects[i].index, objects[i].value); + //printf("%d > %g\n", objects[i].index, objects[i].value); dcmDim[objects[i].index].dimIdx[0] = i; } for (int i = 0; i < 4; i++) { @@ -7664,7 +7811,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD int stackPositionItem = 0; if (dimensionIndexPointerCounter > 0) for (size_t i = 0; i < dimensionIndexPointerCounter; i++) - if (dimensionIndexPointer[i] == kInStackPositionNumber) + if ((dimensionIndexPointer[i] == kInStackPositionNumber) || (dimensionIndexPointer[i] == kImagePositionPatient)) //issue706 stackPositionItem = i; if ((d.manufacturer == kMANUFACTURER_CANON) && (nVariableItems == 1) && (d.xyzDim[4] > 1)) { //WARNING: Canon CANON V6.0SP2001* (0018,9005) = "AX fMRI" strangely sets TemporalPositionIndex(0020,9128) as 1 for all volumes: (0020,9157) and (0020,9128) are INCORRECT! @@ -7724,7 +7871,10 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD if (dti4D->intenScalePhilips[i] != dti4D->intenScalePhilips[0]) d.isScaleVariesEnh = true; } - if (!(d.manufacturer == kMANUFACTURER_BRUKER && d.isDiffusion) && (d.xyzDim[4] > 1) && (d.xyzDim[4] < kMaxDTI4D)) { //record variations in TE + //Revert Bruker change as new Bruker data uses DimensionIndexValues 0020,9157 for diffusion vectors + // https://github.com/rordenlab/dcm2niix/commit/8207877de984dab59e0374906c7ad212756f85a6#diff-120bb215d28dcbddeca8ea56dbb97e237501901ac2dfa1fbe4182282a8dd2b75 + //if (!(d.manufacturer == kMANUFACTURER_BRUKER && d.isDiffusion) && (d.xyzDim[4] > 1) && (d.xyzDim[4] < kMaxDTI4D)) { //record variations in TE + if ((d.xyzDim[4] > 1) && (d.xyzDim[4] < kMaxDTI4D)) { //record variations in TE d.isScaleOrTEVaries = false; bool isTEvaries = false; bool isScaleVaries = false; @@ -7734,6 +7884,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD int slice = j + (i * d.xyzDim[3]); //dti4D->gradDynVol[i] = 0; //only PAR/REC dti4D->TE[i] = dcmDim[slice].TE; + dti4D->TR[i] = dcmDim[slice].TR; dti4D->isPhase[i] = dcmDim[slice].isPhase; dti4D->isReal[i] = dcmDim[slice].isReal; dti4D->isImaginary[i] = dcmDim[slice].isImaginary; @@ -7814,6 +7965,13 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //if (contentTime != 0.0) && (numDimensionIndexValues < (MAX_NUMBER_OF_DIMENSIONS - 1)){ // uint_32t timeCRC = mz_crc32X((unsigned char*) &contentTime, sizeof(double)); //} + if (numberOfFramesICEdims < 2) //issue751: icedims[20] frames for EPI only + numberOfFramesICEdims = 0; + if ((numberOfFramesICEdims > 0) && (d.xyzDim[3] != numberOfFramesICEdims)) { + printWarning("Series %ld includes partial volume (issue 742): %d slices acquired but ICE dims (0021,118e) specifies %d \n", d.seriesNum, d.xyzDim[3], numberOfFramesICEdims); + d.seriesNum += 1000; + d.isDerived = true; + } //if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strcmp(d.sequenceName, "fldyn3d1")== 0)) { if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "fldyn3d1") != NULL)) { //combine DCE series https://github.com/rordenlab/dcm2niix/issues/252 @@ -7884,8 +8042,13 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD } d.seriesUidCrc = mz_crc32X((unsigned char *)&d.seriesInstanceUID, strlen(d.seriesInstanceUID)); } - if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "fl2d1") != NULL)) { - d.isLocalizer = true; + if (d.manufacturer == kMANUFACTURER_SIEMENS) { + if (strstr(d.seriesDescription, "AAHScout") != NULL) + d.isLocalizer = true; + if (strstr(d.protocolName, "AAHScout") != NULL) + d.isLocalizer = true; + if (strstr(d.sequenceName, "fl2d1") != NULL) + d.isLocalizer = true; } // detect GE diffusion gradient cycling mode (see issue 635) // GE diffusion epi @@ -7922,6 +8085,9 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD d.accelFactPE = accelFactPE; //printf("Determining accelFactPE from 0018,9069 not 0021,1009 or 0051,1011\n"); } + //avoid false positives: non-EPI GE scans can report b-value = 0 + if ((d.isDiffusion) && (d.manufacturer == kMANUFACTURER_GE)) + d.isDiffusion = ((d.internalepiVersionGE == kGE_EPI_EPI2) || (d.epiVersionGE == kGE_EPI_EPI2)); //detect pepolar https://github.com/nipy/heudiconv/issues/479 if ((d.epiVersionGE == kGE_EPI_PEPOLAR_FWD) && (userData12GE == 1)) d.epiVersionGE = kGE_EPI_PEPOLAR_REV; @@ -7946,7 +8112,7 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD if ((d.manufacturer == kMANUFACTURER_UIH) && (strstr(d.sequenceName, "gre_fsp") != NULL)) d.echoTrainLength = 0; //printf(">>%s\n", d.sequenceName); d.isValid = false; - // Andrey Fedorov has requested keeping GE bvalues, see issue 264 + // Andrey Fedorov has requested keeping GE bvaecues, see issue 264 //if ((d.CSA.numDti > 0) && (d.manufacturer == kMANUFACTURER_GE) && (d.numberOfDiffusionDirectionGE < 1)) // d.CSA.numDti = 0; //https://github.com/rordenlab/dcm2niix/issues/264 if ((!d.isLocalizer) && (isInterpolated) && (d.imageNum <= 1)) @@ -7978,8 +8144,11 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD //in practice 0020,0110 not used //https://github.com/bids-standard/bep001/blob/repetitiontime/Proposal_RepetitionTime.md } - //issue690 - if ((d.manufacturer == kMANUFACTURER_GE) && (diffusionDirectionTypeGE > 0) && (diffusionDirectionTypeGE != 16)) + //issue690, 777 + // detect non-DTI for GE + if ((d.manufacturer == kMANUFACTURER_GE) && (diffusionDirectionTypeGE > 0) && (diffusionDirectionTypeGE != 16) && (diffusionDirectionTypeGE != 14)) + d.numberOfDiffusionDirectionGE = 0; + if ((d.manufacturer == kMANUFACTURER_GE) && (seriesdiffusionDirectionTypeGE > 0) && (seriesdiffusionDirectionTypeGE != 16)) d.numberOfDiffusionDirectionGE = 0; //issue 542 if ((d.manufacturer == kMANUFACTURER_GE) && (isNeologica) && (!isSameFloat(d.CSA.dtiV[0], 0.0f)) && ((isSameFloat(d.CSA.dtiV[1], 0.0f)) && (isSameFloat(d.CSA.dtiV[2], 0.0f)) && (isSameFloat(d.CSA.dtiV[3], 0.0f)) ) ) @@ -7999,6 +8168,13 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD d.rawDataRunNumber = philMRImageDiffVolumeNumber; d.phaseNumber = 0; } + // Phase encoding polarity + //next conditionals updated: make GE match Siemens + // it also rescues Siemens XA20 images without CSA header but with 0021,111C + if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) + d.CSA.phaseEncodingDirectionPositive = 1; + if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_FLIPPED) + d.CSA.phaseEncodingDirectionPositive = 0; //issue 668: several SAR levels for different regions (IEC_HEAD, IEC_LOCAL, etc) d.SAR = fmax(maxSAR, d.SAR); // d.rawDataRunNumber = (d.rawDataRunNumber > d.phaseNumber) ? d.rawDataRunNumber : d.phaseNumber; //will not work: conflict for MultiPhase ASL with multiple averages @@ -8016,6 +8192,8 @@ const uint32_t kEffectiveTE = 0x0018 + uint32_t(0x9082 << 16); //FD d.isValid = false; } //printf("%g\t%g\t%s\n", d.intenIntercept, d.intenScale, fname); + if ((d.isLocalizer) && (strstr(d.seriesDescription, "b1map"))) //issue751 b1map uses same base as scout + d.isLocalizer = false; return d; } // readDICOMx() @@ -8040,3 +8218,32 @@ struct TDICOMdata readDICOM(char *fname) { free(dti4D); return ret; } // readDICOM() + + +#ifdef USING_DCM2NIIXFSWRAPPER +// remove spaces, '[', ']' in given buf +// ??? also remove <, >, & +void remove_specialchars(char *buf) +{ + if (strchr(buf, ' ') == NULL) + return; + + int buflen = strlen(buf)+1; + char *newbuf = new char[buflen]; + memset(newbuf, '\0' , buflen); + + char *ptr_newbuf = &newbuf[0]; + for (int i = 0; i < buflen; i++) + { + // skip spaces, '[', ']' + if (buf[i] == ' ' || buf[i] == '[' || buf[i] == ']') + continue; + + *ptr_newbuf = buf[i]; + ptr_newbuf++; + } + + memcpy(buf, newbuf, ptr_newbuf-newbuf); + free(newbuf); +} +#endif diff --git a/console/nii_dicom.h b/console/nii_dicom.h index aa86ef41..2fbb776e 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20230411" +#define kDCMdate "v1.0.20240202" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic @@ -64,6 +64,7 @@ static const int kMaxDTI4D = kMaxSlice2D; //issue460: maximum number of DTI dire #define kDICOMStr 66 //64 characters plus NULL https://github.com/rordenlab/dcm2niix/issues/268 #define kDICOMStrLarge 256 +#define kDICOMStrExtraLarge 65536 // for Siemens WipMemBlock only #define kMANUFACTURER_UNKNOWN 0 #define kMANUFACTURER_SIEMENS 1 @@ -192,7 +193,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; int sliceOrder[kMaxSlice2D]; // [7,3,2] means the first slice on disk should be moved to 7th position int gradDynVol[kMaxDTI4D]; //used to parse dimensions of Philips data, e.g. file with multiple dynamics, echoes, phase+magnitude //int fragmentOffset[kMaxDTI4D], fragmentLength[kMaxDTI4D]; //for images with multiple compressed fragments - float frameReferenceTime[kMaxDTI4D], frameDuration[kMaxDTI4D], decayFactor[kMaxDTI4D], volumeOnsetTime[kMaxDTI4D], triggerDelayTime[kMaxDTI4D], TE[kMaxDTI4D], RWVScale[kMaxDTI4D], RWVIntercept[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; + float frameReferenceTime[kMaxDTI4D], frameDuration[kMaxDTI4D], decayFactor[kMaxDTI4D], volumeOnsetTime[kMaxDTI4D], triggerDelayTime[kMaxDTI4D], TE[kMaxDTI4D], TR[kMaxDTI4D], RWVScale[kMaxDTI4D], RWVIntercept[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; bool isReal[kMaxDTI4D]; bool isImaginary[kMaxDTI4D]; bool isPhase[kMaxDTI4D]; @@ -230,20 +231,23 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; float sliceTiming[kMaxEPI3D], dtiV[4], sliceNormV[4], bandwidthPerPixelPhaseEncode, sliceMeasurementDuration; int coilNumber, numDti, SeriesHeader_offset, SeriesHeader_length, multiBandFactor, sliceOrder, slice_start, slice_end, mosaicSlices, protocolSliceNumber1, phaseEncodingDirectionPositive; bool isPhaseMap; + char bidsDataType[kDICOMStr]; //anat, func, dwi + char bidsEntitySuffix[kDICOMStrLarge]; //anat, func, dwi + char bidsTask[kDICOMStr]; //rest, naming40 }; struct TDICOMdata { long seriesNum; int xyzDim[5]; uint32_t coilCrc, seriesUidCrc, instanceUidCrc; int overlayStart[kMaxOverlay]; - int postLabelDelay, shimGradientX, shimGradientY, shimGradientZ, phaseNumber, spoiling, mtState, partialFourierDirection, interp3D, aslFlags, durationLabelPulseGE, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfImagesInGridUIH, numberOfDiffusionT2GE, numberOfDiffusionDirectionGE, tensorFileGE, diffCyclingModeGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, frequencyEncodingSteps, phaseEncodingStepsOutOfPlane, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme; + int postLabelDelay, shimGradientX, shimGradientY, shimGradientZ, phaseNumber, spoiling, mtState, partialFourierDirection, interp3D, aslFlags, durationLabelPulseGE, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfTR, numberOfImagesInGridUIH, numberOfDiffusionT2GE, numberOfDiffusionDirectionGE, tensorFileGE, diffCyclingModeGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, frequencyEncodingSteps, phaseEncodingStepsOutOfPlane, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme; float compressedSensingFactor, xRayTubeCurrent, exposureTimeMs, numberOfExcitations, numberOfArms, numberOfPointsPerArm, groupDelay, decayFactor, percentSampling,waterFatShift, numberOfAverages, patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, accelFactOOP, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4], velocityEncodeScaleGE; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float frameReferenceTime, frameDuration, ecat_isotope_halflife, ecat_dosage; float pixelPaddingValue; // used for both FloatPixelPaddingValue (0028, 0122) and PixelPaddingValue (0028, 0120); NaN if not present. double imagingFrequency, acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode; - char parallelAcquisitionTechnique[kDICOMStr], radiopharmaceutical[kDICOMStr], convolutionKernel[kDICOMStr], unitsPT[kDICOMStr], decayCorrection[kDICOMStr], attenuationCorrectionMethod[kDICOMStr],reconstructionMethod[kDICOMStr]; + char parallelAcquisitionTechnique[kDICOMStr], radiopharmaceutical[kDICOMStr], convolutionKernel[kDICOMStr], unitsPT[kDICOMStr], decayCorrection[kDICOMStr], attenuationCorrectionMethod[kDICOMStr],reconstructionMethod[kDICOMStr], transferSyntax[kDICOMStr]; char prescanReuseString[kDICOMStr], imageOrientationText[kDICOMStr], pulseSequenceName[kDICOMStr], coilElements[kDICOMStr], coilName[kDICOMStr], phaseEncodingDirectionDisplayedUIH[kDICOMStr], imageBaseName[kDICOMStr], stationName[kDICOMStr], softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], instanceUID[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageTypeText[kDICOMStr], imageType[kDICOMStr], institutionalDepartmentName[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr], accessionNumber[kDICOMStr], seriesDescription[kDICOMStr], studyID[kDICOMStr], sequenceName[kDICOMStr], protocolName[kDICOMStr],sequenceVariant[kDICOMStr],scanningSequence[kDICOMStr], patientBirthDate[kDICOMStr], patientAge[kDICOMStr], studyDate[kDICOMStr],studyTime[kDICOMStr]; char deepLearningText[kDICOMStrLarge], scanOptions[kDICOMStrLarge], institutionAddress[kDICOMStrLarge], imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; @@ -276,6 +280,9 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int isVerbose); int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeSForm) ; unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose, struct TDTI4D *dti4D); +#ifdef USING_DCM2NIIXFSWRAPPER + void remove_specialchars(char *buf); +#endif #ifdef __cplusplus } #endif diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index f913f2ff..3b2eb676 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -38,7 +38,7 @@ #endif #include "nii_dicom.h" #include "nii_ortho.h" -#ifdef myEnableJNIfTI +#ifdef myEnableJNIFTI #include "base64.h" #include "cJSON.h" #endif @@ -113,6 +113,7 @@ const char kFileSep[2] = "/"; // create the struct to save nifti header, image data, TDICOMdata, & TDTI information. // no .nii, .bval, .bvec are created. MRIFSSTRUCT mrifsStruct; +std::vector mrifsStruct_vector; // retrieve the struct MRIFSSTRUCT* nii_getMrifsStruct() @@ -125,6 +126,36 @@ void nii_clrMrifsStruct() { free(mrifsStruct.imgM); free(mrifsStruct.tdti); + + for (int n = 0; n < mrifsStruct.nDcm; n++) + free(mrifsStruct.dicomlst[n]); + + if (mrifsStruct.dicomlst != NULL) + free(mrifsStruct.dicomlst); + +} + +// retrieve the struct +std::vector* nii_getMrifsStructVector() +{ + return &mrifsStruct_vector; +} + +// free the memory used for the image and dti +void nii_clrMrifsStructVector() +{ + int nitem = mrifsStruct_vector.size(); + for (int n = 0; n < nitem; n++) + { + free(mrifsStruct_vector[n].imgM); + free(mrifsStruct_vector[n].tdti); + + for (int n = 0; n < mrifsStruct.nDcm; n++) + free(mrifsStruct_vector[n].dicomlst[n]); + + if (mrifsStruct_vector[n].dicomlst != NULL) + free(mrifsStruct_vector[n].dicomlst); + } } #endif @@ -266,8 +297,7 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV if ((toupper(d->patientOrient[0]) == 'H') && (toupper(d->patientOrient[1]) == 'F') && (toupper(d->patientOrient[2]) == 'S')) ; //participant was head first supine else { - printMessage("GE DTI directions require head first supine acquisition\n"); - return; + printWarning("Limited validation for non-HFS (Head first supine) GE DTI: confirm gradient vector transformation\n"); } bool col = false; if (d->phaseEncodingRC == 'C') @@ -587,7 +617,7 @@ float readKeyFloat(const char *key, char *buffer, int remLength) { //look for te return atof(str); } //readKeyFloat() -void readKeyStr(const char *key, char *buffer, int remLength, char *outStr) { +void readKeyStrLen(const char *key, char *buffer, int remLength, char *outStr, int outStrLen) { //if key is CoilElementID.tCoilID the string 'CoilElementID.tCoilID = ""Head_32""' returns 'Head32' strcpy(outStr, ""); char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); @@ -599,7 +629,7 @@ void readKeyStr(const char *key, char *buffer, int remLength, char *outStr) { tmpstr[1] = 0; bool isQuote = false; while ((i < remLength) && (keyPos[i] != 0x0A)) { - if ((isQuote) && (keyPos[i] != '"') && (outLen < kDICOMStrLarge)) { + if ((isQuote) && (keyPos[i] != '"') && (outLen < outStrLen)) { tmpstr[0] = keyPos[i]; strcat(outStr, tmpstr); outLen++; @@ -613,6 +643,10 @@ void readKeyStr(const char *key, char *buffer, int remLength, char *outStr) { } } //readKeyStr() +void readKeyStr(const char *key, char *buffer, int remLength, char *outStr) { + readKeyStrLen(key, buffer, remLength, outStr, kDICOMStrLarge); +} //readKeyStr() + int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) { //returns offset to ASCII Phoenix data if (lLength < 36) @@ -650,14 +684,16 @@ int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) { } //phoenixOffsetCSASeriesHeader() #define kMaxWipFree 64 +#define freeDiffusionMaxN 512 typedef struct { float TE0, TE1, delayTimeInTR, phaseOversampling, phaseResolution, txRefAmp, accelFactTotal; - int phaseEncodingLines, existUcImageNumb, ucMode, baseResolution, interp, partialFourier, echoSpacing, - difBipolar, parallelReductionFactorInPlane, refLinesPE, combineMode, patMode, ucMTC, accelFact3D; + int lInvContrasts, lContrasts ,phaseEncodingLines, existUcImageNumb, ucMode, baseResolution, interp, partialFourier, echoSpacing, + difBipolar, parallelReductionFactorInPlane, refLinesPE, combineMode, patMode, ucMTC, accelFact3D, freeDiffusionN; float alFree[kMaxWipFree]; float adFree[kMaxWipFree]; float alTI[kMaxWipFree]; float sPostLabelingDelay, ulLabelingDuration, dAveragesDouble, dThickness, ulShape, sPositionDTra, sNormalDTra; + vec3 freeDiffusionVec[freeDiffusionMaxN]; } TCsaAscii; void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, int csaLength, float *shimSetting, char *coilID, char *consistencyInfo, char *coilElements, char *pulseSequenceDetails, char *fmriExternalInfo, char *protocolName, char *wipMemBlock) { @@ -678,6 +714,8 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i csaAscii->interp = 0; csaAscii->partialFourier = 0; csaAscii->echoSpacing = 0; + csaAscii->lInvContrasts = 0; + csaAscii->lContrasts = 0; csaAscii->difBipolar = 0; //0=not assigned,1=bipolar,2=monopolar csaAscii->parallelReductionFactorInPlane = 0; csaAscii->accelFact3D = 0;//lAccelFact3D @@ -686,6 +724,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i csaAscii->combineMode = 0; csaAscii->patMode = 0; csaAscii->ucMTC = 0; + csaAscii->freeDiffusionN = 0; for (int i = 0; i < 8; i++) shimSetting[i] = 0.0; strcpy(coilID, ""); @@ -727,7 +766,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i char keyStr[] = "### ASCCONV BEGIN"; //skip to start of ASCII often "### ASCCONV BEGIN ###" but also "### ASCCONV BEGIN object=MrProtDataImpl@MrProtocolData" char *keyPos = (char *)memmem(bufferTrim, csaLengthTrim, keyStr, strlen(keyStr)); if (keyPos) { - //We could detect multi-echo MPRAGE here, e.g. "lContrasts = 4"- but ideally we want an earlier detection + //We can detect multi-echo MPRAGE here, e.g. "lContrasts = 4"- but ideally we want an earlier detection csaLengthTrim -= (keyPos - bufferTrim); //FmriExternalInfo listed AFTER AscConvEnd and uses different delimiter || // char keyStrExt[] = "FmriExternalInfo"; @@ -741,7 +780,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i #endif char keyStrLns[] = "sKSpace.lPhaseEncodingLines"; csaAscii->phaseEncodingLines = readKey(keyStrLns, keyPos, csaLengthTrim); - char keyStrUcImg[] = "sSliceArray.ucImageNumb"; + char keyStrUcImg[] = "sSliceArray.ucImageNumb"; //some non-mosaics like ToF include "sSliceArray.ucImageNumbSag" csaAscii->existUcImageNumb = readKey(keyStrUcImg, keyPos, csaLengthTrim); char keyStrUcMode[] = "sSliceArray.ucMode"; csaAscii->ucMode = readKeyN1(keyStrUcMode, keyPos, csaLengthTrim); @@ -753,6 +792,10 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i csaAscii->partialFourier = readKey(keyStrPF, keyPos, csaLengthTrim); char keyStrES[] = "sFastImaging.lEchoSpacing"; csaAscii->echoSpacing = readKey(keyStrES, keyPos, csaLengthTrim); + char keyStrNumInv[] = "lInvContrasts"; + csaAscii->lInvContrasts = readKey(keyStrNumInv, keyPos, csaLengthTrim); + char keyStrNumEcho[] = "lContrasts"; + csaAscii->lContrasts = readKey(keyStrNumEcho, keyPos, csaLengthTrim); char keyStrDS[] = "sDiffusion.dsScheme"; csaAscii->difBipolar = readKey(keyStrDS, keyPos, csaLengthTrim); if (csaAscii->difBipolar == 0) { @@ -796,7 +839,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i char keyStrSeq[] = "tSequenceFileName"; readKeyStr(keyStrSeq, keyPos, csaLengthTrim, pulseSequenceDetails); char keyStrWipMemBlock[] = "sWipMemBlock.tFree"; - readKeyStr(keyStrWipMemBlock, keyPos, csaLengthTrim, wipMemBlock); + readKeyStrLen(keyStrWipMemBlock, keyPos, csaLengthTrim, wipMemBlock, kDICOMStrExtraLarge); char keyStrPn[] = "tProtocolName"; readKeyStr(keyStrPn, keyPos, csaLengthTrim, protocolName); char keyStrTE0[] = "alTE[0]"; @@ -903,6 +946,29 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i shimSetting[6] = readKeyFloat(keyStrSh6, keyPos, csaLengthTrim); char keyStrSh7[] = "sGRADSPEC.alShimCurrent[4]"; shimSetting[7] = readKeyFloat(keyStrSh7, keyPos, csaLengthTrim); + // pull out the directions in the DVI + char keyStrDVIn[] = "sDiffusion.sFreeDiffusionData.lDiffDirections"; + int nDiffDir = readKey(keyStrDVIn, keyPos, csaLengthTrim); + csaAscii->freeDiffusionN = min(nDiffDir, freeDiffusionMaxN); + + //printMessage("Free diffusion: %i\n", csaAscii->freeDiffusionN); + + for (int k = 0; k < csaAscii->freeDiffusionN; k++) { + char txt[128]; + + snprintf(txt, 128, "sDiffusion.sFreeDiffusionData.asDiffDirVector[%i].dSag", k); + float x = readKeyFloat(txt, keyPos, csaLengthTrim); + + snprintf(txt, 128, "sDiffusion.sFreeDiffusionData.asDiffDirVector[%i].dCor", k); + float y = readKeyFloat(txt, keyPos, csaLengthTrim); + + snprintf(txt, 128, "sDiffusion.sFreeDiffusionData.asDiffDirVector[%i].dTra", k); + float z = readKeyFloat(txt, keyPos, csaLengthTrim); + + csaAscii->freeDiffusionVec[k].v[0] = x; + csaAscii->freeDiffusionVec[k].v[1] = y; + csaAscii->freeDiffusionVec[k].v[2] = z; + } } free(buffer); return; @@ -915,7 +981,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i #define myReadGeProtocolBlock #endif #ifdef myReadGeProtocolBlock -int geProtocolBlock(const char *filename, int geOffset, int geLength, int isVerbose, int *sliceOrder, int *viewOrder, int *mbAccel, int *nSlices, float *groupDelay, char ioptGE[]) { +int geProtocolBlock(const char *filename, int geOffset, int geLength, int isVerbose, int *sliceOrder, int *viewOrder, int *mbAccel, int *nSlices, float *groupDelay, char ioptGE[], char seqStr[]) { *sliceOrder = -1; *viewOrder = 0; *mbAccel = 0; @@ -994,7 +1060,7 @@ int geProtocolBlock(const char *filename, int geOffset, int geLength, int isVerb // if ((pUnCmp[0] == '<') && (pUnCmp[1] == '?')) printWarning("New XML-based GE Protocol Block is not yet supported: please report issue on dcm2niix Github page\n"); - char keyStrSO[] = "SLICEORDER"; + char keyStrSO[] = "\nSLICEORDER"; *sliceOrder = readKeyN1(keyStrSO, (char *)pUnCmp, unCmpSz); char keyStrVO[] = "VIEWORDER"; *viewOrder = readKey(keyStrVO, (char *)pUnCmp, unCmpSz); @@ -1007,6 +1073,9 @@ int geProtocolBlock(const char *filename, int geOffset, int geLength, int isVerb readKeyStr(keyStrDELACQ, (char *)pUnCmp, unCmpSz, DELACQ); char keyStrGD[] = "DELACQNOAV"; *groupDelay = readKeyFloat(keyStrGD, (char *)pUnCmp, unCmpSz); + + char keyStrPSEQ[] = "PSEQ"; + readKeyStr(keyStrPSEQ, (char *)pUnCmp, unCmpSz, seqStr); char keyStrIOPT[] = "IOPT"; readKeyStr(keyStrIOPT, (char *)pUnCmp, unCmpSz, ioptGE); char PHASEDELAYS1[10000]; @@ -1129,7 +1198,7 @@ void rescueProtocolName(struct TDICOMdata *d, const char *filename) { return; #ifdef myReadAsciiCsa float shimSetting[8]; - char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrLarge]; + char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrExtraLarge]; TCsaAscii csaAscii; siemensCsaAscii(filename, &csaAscii, d->CSA.SeriesHeader_offset, d->CSA.SeriesHeader_length, shimSetting, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo, protocolName, wipMemBlock); if (strlen(protocolName) >= kDICOMStr) @@ -1152,7 +1221,18 @@ void nii_SaveBIDSX(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts char txtname[2048] = {""}; strcpy(txtname, pathoutname); strcat(txtname, ".json"); +#ifdef USING_R + FILE *fp = NULL; + if (opts.isImageInMemory) + { + ImageList *images = (ImageList *)opts.imageList; + fp = images->jsonHandle(); + } + if (fp == NULL) + fp = fopen(txtname, "w"); +#else FILE *fp = fopen(txtname, "w"); +#endif fprintf(fp, "{\n"); switch (d.modality) { case kMODALITY_CR: @@ -1219,7 +1299,7 @@ tse3d: T2*/ fprintf(fp, "\t\"MagneticFieldStrength\": %g,\n", d.fieldStrength); //Imaging Frequency (0018,0084) can be useful https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth // however, UIH stores 128176031 not 128.176031 https://github.com/rordenlab/dcm2niix/issues/225 - if (d.imagingFrequency < 9000000) + if ((d.imagingFrequency < 9000000.0) && (d.imagingFrequency > 0.0)) fprintf(fp, "\t\"ImagingFrequency\": %.10g,\n", d.imagingFrequency); switch (d.manufacturer) { case kMANUFACTURER_BRUKER: @@ -1294,6 +1374,14 @@ tse3d: T2*/ json_Float(fp, "\t\"PatientWeight\": %g,\n", d.patientWeight); //d.patientBirthDate //convert from DICOM YYYYMMDD to JSON //d.patientAge //4-digit Age String: nnnD, nnnW, nnnM, nnnY; + //issue 763 following BIDS standard, unit for Age is YEARS + int ageLen = strlen(d.patientAge); + if ((ageLen > 1) && (d.patientAge[ageLen-1] == 'Y')) { + char ageStr[kDICOMStr]; + strcpy(ageStr, d.patientAge); + ageStr[ageLen -1] = '\0'; + fprintf(fp, "\t\"PatientAge\": %d,\n", atoi(ageStr)); + } } if (d.isQuadruped) json_Bool(fp, "\t\"Quadruped\": %s,\n", true); // BIDS suggests 0018,9020 but Siemens V-series do not populate this, alternatives are CSA or (0018,0021) CS [SK\MTC\SP] @@ -1406,7 +1494,7 @@ tse3d: T2*/ fprintf(fp, "\t\"UsePhilipsFloatNotDisplayScaling\": %d,\n", opts.isPhilipsFloatNotDisplayScaling); } //https://bids-specification--622.org.readthedocs.build/en/622/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#case-3-direct-field-mapping - if ((d.isRealIsPhaseMapHz) && (d.isHasReal)) + if ((d.isRealIsPhaseMapHz) && ((d.manufacturer == kMANUFACTURER_GE) || (d.isHasReal))) fprintf(fp, "\t\"Units\": \"Hz\",\n"); // //PET ISOTOPE MODULE ATTRIBUTES json_Str(fp, "\t\"Radiopharmaceutical\": \"%s\",\n", d.radiopharmaceutical); @@ -1526,7 +1614,8 @@ tse3d: T2*/ bool interp = false; //2D interpolation float phaseOversampling = 0.0; //n.b. https://neurostars.org/t/getting-missing-ge-information-required-by-bids-for-common-preprocessing/1357/7 - json_Str(fp, "\t\"PhaseEncodingDirectionDisplayed\": \"%s\",\n", d.phaseEncodingDirectionDisplayedUIH); + if (d.manufacturer == kMANUFACTURER_UIH) + json_Str(fp, "\t\"PhaseEncodingDirectionDisplayed\": \"%s\",\n", d.phaseEncodingDirectionDisplayedUIH); if ((d.manufacturer == kMANUFACTURER_GE) && (d.phaseEncodingGE != kGE_PHASE_ENCODING_POLARITY_UNKNOWN)) { //only set for GE if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) fprintf(fp, "\t\"PhaseEncodingPolarityGE\": \"Unflipped\",\n"); @@ -1568,7 +1657,7 @@ tse3d: T2*/ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) { float pf = 1.0f; //partial fourier float shimSetting[8]; - char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrLarge]; + char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrExtraLarge]; TCsaAscii csaAscii; siemensCsaAscii(filename, &csaAscii, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, shimSetting, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo, protocolName, wipMemBlock); if ((d.phaseEncodingLines < 1) && (csaAscii.phaseEncodingLines > 0)) @@ -1584,7 +1673,7 @@ tse3d: T2*/ if (csaAscii.dAveragesDouble > 1.0) //*spcR_44ns fractional and independent of (0018,0083) DS NumberOfAverages, e.g. 0018,0083=2, dAveragesDouble = 1.4? json_Float(fp, "\t\"AveragesDouble\": %g,\n", csaAscii.dAveragesDouble); phaseOversampling = csaAscii.phaseOversampling; - if (csaAscii.existUcImageNumb > 0) { + if ((d.CSA.mosaicSlices > 1) && (csaAscii.existUcImageNumb > 0)) { if (d.CSA.protocolSliceNumber1 < 2) { printWarning("Assuming mosaics saved in reverse order due to 'sSliceArray.ucImageNumb'\n"); //never seen such an image in the wild.... sliceDir may need to be reversed @@ -1610,12 +1699,16 @@ tse3d: T2*/ json_FloatNotNan(fp, "\t\"LabelingDuration\": %g,\n", csaAscii.ulLabelingDuration * (1.0 / 1000000.0)); //usec->sec json_FloatNotNan(fp, "\t\"PostLabelingDelay\": %g,\n", csaAscii.sPostLabelingDelay * (1.0 / 1000000.0)); //usec -> sec } + //ASL specific tags - 2D pCASL Danny J.J. Wang http://www.loft-lab.org if ((strstr(pulseSequenceDetails, "ep2d_pcasl")) || (strstr(pulseSequenceDetails, "ep2d_pcasl_UI_PHC"))) { isPCASL = true; repetitionTimePreparation = d.TR; json_FloatNotNan(fp, "\t\"LabelOffset\": %g,\n", csaAscii.adFree[1]); //mm - json_FloatNotNan(fp, "\t\"PostLabelDelay\": %g,\n", csaAscii.adFree[2] * (1.0 / 1000000.0)); //usec -> sec - json_FloatNotNan(fp, "\t\"NumRFBlocks\": %g,\n", csaAscii.adFree[3]); + json_FloatNotNan(fp, "\t\"PostLabelingDelay\": %g,\n", csaAscii.adFree[2] * (1.0 / 1000000.0)); //usec -> sec + float num_RF_Block = csaAscii.adFree[3]; + json_FloatNotNan(fp, "\t\"NumRFBlocks\": %g,\n", num_RF_Block); + // Sep 5, 2023, at 7:56 PM, Danny JJ Wang the labeling duration is (0.92*20*Num_RF_Block) ms + json_FloatNotNan(fp, "\t\"LabelingDuration\": %g,\n", (0.92 *20.0 *num_RF_Block)/1000.0); //in seconds json_FloatNotNan(fp, "\t\"RFGap\": %g,\n", csaAscii.adFree[4] * (1.0 / 1000000.0)); //usec -> sec json_FloatNotNan(fp, "\t\"MeanGzx10\": %g,\n", csaAscii.adFree[10]); json_FloatNotNan(fp, "\t\"PhiAdjust\": %g,\n", csaAscii.adFree[11]); // percent @@ -1627,6 +1720,10 @@ tse3d: T2*/ json_FloatNotNan(fp, "\t\"RFGap\": %g,\n", csaAscii.adFree[4] * (1.0 / 1000000.0)); //usec -> sec json_FloatNotNan(fp, "\t\"MeanGzx10\": %g,\n", csaAscii.adFree[10]); // mT/m json_FloatNotNan(fp, "\t\"T1\": %g,\n", csaAscii.adFree[12] * (1.0 / 1000000.0)); //usec -> sec + float num_RF_Block = csaAscii.adFree[3]; + json_FloatNotNan(fp, "\t\"NumRFBlocks\": %g,\n", num_RF_Block); + // Sep 5, 2023, at 7:56 PM, Danny JJ Wang the labeling duration is (0.92*20*Num_RF_Block) ms + json_FloatNotNan(fp, "\t\"LabelingDuration\": %g,\n", (0.92 *20.0 *num_RF_Block)/1000.0); //in seconds } //ASL specific tags - 2D PASL Siemens Product if (strstr(pulseSequenceDetails, "ep2d_pasl")) { @@ -1645,7 +1742,7 @@ tse3d: T2*/ if (strstr(pulseSequenceDetails, "ep2d_fairest")) { isPASL = true; json_FloatNotNan(fp, "\t\"PostInversionDelay\": %g,\n", csaAscii.adFree[2] * (1.0 / 1000.0)); //usec->sec - json_FloatNotNan(fp, "\t\"PostLabelDelay\": %g,\n", csaAscii.adFree[4] * (1.0 / 1000.0)); //usec -> sec + json_FloatNotNan(fp, "\t\"PostLabelingDelay\": %g,\n", csaAscii.adFree[4] * (1.0 / 1000.0)); //usec -> sec } //ASL specific tags - Oxford (Thomas OKell) bool isOxfordASL = false; @@ -1747,6 +1844,33 @@ tse3d: T2*/ // https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#common-metadata-fields-applicable-to-both-pcasl-and-pasl if (((isPASL) || (isPCASL)) && (csaAscii.interp <= 0)) fprintf(fp, "\t\"AcquisitionVoxelSize\": [\n\t\t%g,\n\t\t%g,\n\t\t%g\t],\n", d.xyzMM[1], d.xyzMM[2], d.zThick); + // lund free waveform sequence, see https://github.com/filip-szczepankiewicz/fwf_sequence_tools + if ( (strstr(pulseSequenceDetails, "ep2d_diff_fwf") != 0) || (strstr(pulseSequenceDetails, "ep2d_diff_sms_fwf_simple") != 0)) { + for (int i = 0; i < kMaxWipFree; i++) { + if (!isnan(csaAscii.adFree[i])) + fprintf(fp, "\t\"FWF_adFree[%i]\": %g,\n", i, csaAscii.adFree[i]); + } + + for (int i = 0; i < kMaxWipFree; i++) { + if (!isnan(csaAscii.alFree[i])) + fprintf(fp, "\t\"FWF_alFree[%i]\": %g,\n", i, csaAscii.alFree[i]); + } + + for (int d = 0; d < 3; d++) + { + char str [4096]; + strcpy(str, ""); + for (int i = 0; i < csaAscii.freeDiffusionN;i++) { + char tmp[10]; + snprintf(tmp, 10, "%1.4f", csaAscii.freeDiffusionVec[i].v[d]); + strcat(str, tmp); + if ( (i+1) < csaAscii.freeDiffusionN) + strcat(str, ", "); + } + char dchar = 'x' + d; + fprintf(fp, "\t\"FWF_DVS%c\": [ %s ],\n", dchar, str); + } + } //general properties if ((csaAscii.partialFourier > 0) && ((d.modality == kMODALITY_MR))) { //check MR, e.g. do not report for Siemens PET //https://github.com/ismrmrd/siemens_to_ismrmrd/blob/master/parameter_maps/IsmrmrdParameterMap_Siemens_EPI_FLASHREF.xsl @@ -1878,7 +2002,7 @@ tse3d: T2*/ */ //generic public ASL tags if (d.postLabelDelay > 0) - json_Float(fp, "\t\"PostLabelDelay\": %g,\n", float(d.postLabelDelay) / 1000.0); + json_Float(fp, "\t\"PostLabelingDelay\": %g,\n", float(d.postLabelDelay) / 1000.0); //ASL BIDS tags if ((d.aslFlags & kASL_FLAG_GE_CONTINUOUS) || (d.aslFlags & kASL_FLAG_GE_3DCASL)) fprintf(fp, "\t\"ArterialSpinLabelingType\": \"CASL\",\n"); @@ -1888,7 +2012,7 @@ tse3d: T2*/ fprintf(fp, "\t\"ArterialSpinLabelingType\": \"PASL\",\n"); if (d.durationLabelPulseGE > 0) { json_Float(fp, "\t\"LabelingDuration\": %g,\n", d.durationLabelPulseGE / 1000.0); - json_Float(fp, "\t\"PostLabelDelay\": %g,\n", d.TI / 1000.0); //For GE ASL: InversionTime -> Post-label delay + json_Float(fp, "\t\"PostLabelingDelay\": %g,\n", d.TI / 1000.0); //For GE ASL: InversionTime -> Post-label delay json_Float(fp, "\t\"NumberOfPointsPerArm\": %g,\n", d.numberOfPointsPerArm); json_Float(fp, "\t\"NumberOfArms\": %g,\n", d.numberOfArms); } @@ -2001,10 +2125,18 @@ tse3d: T2*/ float roundFactor = 2.0; if (d.isPartialFourier) roundFactor = 4.0; - float totalReadoutTime = ((ceil(1 / roundFactor * d.phaseEncodingLines / d.accelFactPE) * roundFactor) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; - //printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); + float NotPhysicalNumberOfAcquiredPELinesGE = (ceil(1 / roundFactor * d.phaseEncodingLines / d.accelFactPE) * roundFactor); + float NotPhysicalTotalReadOutTimeGE = ( NotPhysicalNumberOfAcquiredPELinesGE - 1.0) * d.effectiveEchoSpacingGE * 0.000001; + // printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTimeGE= %g NumKyLineGE= %d\n", + // d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, NotPhysicalTotalReadOutTimeGE, (int)NotPhysicalNumberOfAcquiredPELinesGE); //json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); - effectiveEchoSpacing = totalReadoutTime / (reconMatrixPE - 1); + effectiveEchoSpacing = NotPhysicalTotalReadOutTimeGE / (d.phaseEncodingLines - 1); + // if this is considered acceptable, meaningful intermediate variables can be written, this might help the end-user. +#ifdef MY_DEBUG + fprintf(fp, "\t\"EchoSpacingMicroSecondsGE\": %d,\n", d.effectiveEchoSpacingGE); + fprintf(fp, "\t\"NotPhysicalNumberOfAcquiredPELinesGE\": %d,\n", (int)(NotPhysicalNumberOfAcquiredPELinesGE)); + json_Float(fp, "\t\"NotPhysicalTotalReadOutTimeGE\": %g,\n", NotPhysicalTotalReadOutTimeGE); +#endif } json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); // Calculate true echo spacing (should match what Siemens reports on the console) @@ -2039,16 +2171,20 @@ tse3d: T2*/ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0)) fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9); // Phase encoding polarity +/* +following logic now occurs earlier to aid bids guess int phPos = d.CSA.phaseEncodingDirectionPositive; //next two conditionals updated: make GE match Siemens if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) phPos = 1; if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_FLIPPED) phPos = 0; +*/ //if ((phPos >= 0) && (d.phaseEncodingRC == 'R') && (d.manufacturer == kMANUFACTURER_UIH)) phPos = 1 - phPos; //issue410 bool isSkipPhaseEncodingAxis = d.is3DAcq; if (d.echoTrainLength > 1) isSkipPhaseEncodingAxis = false; //issue 371: ignore phaseEncoding for 3D MP-RAGE/SPACE, but report for 3D EPI + int phPos = d.CSA.phaseEncodingDirectionPositive; if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && (phPos < 0)) { //when phase encoding axis is known but we do not know phase encoding polarity // https://github.com/rordenlab/dcm2niix/issues/163 @@ -2109,6 +2245,10 @@ tse3d: T2*/ fprintf(fp, "\t\"InPlanePhaseEncodingDirectionDICOM\": \"COL\",\n"); if (d.phaseEncodingRC == 'R') fprintf(fp, "\t\"InPlanePhaseEncodingDirectionDICOM\": \"ROW\",\n"); + if ((opts.isGuessBidsFilename) && (strlen(d.CSA.bidsDataType)) && (strlen(d.CSA.bidsDataType))) + fprintf(fp, "\t\"BidsGuess\": [\"%s\",\"%s\"],\n",d.CSA.bidsDataType, d.CSA.bidsEntitySuffix); + //json_Str(fp, "\t\"StationName\": \"%s\",\n", d.stationName); + // Finish up with info on the conversion tool fprintf(fp, "\t\"ConversionSoftware\": \"dcm2niix\",\n"); fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMdate); @@ -2117,14 +2257,12 @@ tse3d: T2*/ fclose(fp); } // nii_SaveBIDSX() -#ifndef USING_R - void swapEndian(struct nifti_1_header *hdr, unsigned char *im, bool isNative) { //swap endian from big->little or little->big // must be told which is native to detect datatype and number of voxels // one could also auto-detect: hdr->sizeof_hdr==348 if (!isNative) -#ifdef USING_MGH_NIFTI_IO +#if defined(USING_MGH_NIFTI_IO) || defined(USING_R) swap_nifti_header(hdr, 1); #else swap_nifti_header(hdr); @@ -2136,7 +2274,7 @@ void swapEndian(struct nifti_1_header *hdr, unsigned char *im, bool isNative) { int bitpix = hdr->bitpix; int datatype = hdr->datatype; if (isNative) -#ifdef USING_MGH_NIFTI_IO +#if defined(USING_MGH_NIFTI_IO) || defined(USING_R) swap_nifti_header(hdr, 1); #else swap_nifti_header(hdr); @@ -2152,6 +2290,8 @@ void swapEndian(struct nifti_1_header *hdr, unsigned char *im, bool isNative) { nifti_swap_8bytes(nVox, im); } +#ifndef USING_R + void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct nifti_1_header *h, const char *filename) { struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); dti4D->sliceOrder[0] = -1; @@ -2267,8 +2407,13 @@ int *nii_saveDTI(char pathoutname[], int nConvert, struct TDCMsort dcmSort[], st allB0 = false; if (numDti > 1) allB0 = false; - if (nConvert > 1) - allB0 = false; + if (nConvert > 1) { + for (int i = 0; i < nConvert; i++) { //for each image + float b0 = dcmList[dcmSort[i].indx].CSA.dtiV[0]; + if (!isSameFloat(b0, 0.0)) + allB0 = false; + } + } if ((numDti == 1) && (dti4D->S[0].V[0] > 50.0)) allB0 = false; if (allB0) { @@ -2401,7 +2546,7 @@ int *nii_saveDTI(char pathoutname[], int nConvert, struct TDCMsort dcmSort[], st *numADC = 0; bvals = (float *)malloc(numDti * sizeof(float)); int numGEwarn = 0; - bool isGEADC = (dcmList[indx0].numberOfDiffusionDirectionGE == 0); + bool isGEADC = (dcmList[indx0].numberOfDiffusionDirectionGE == 0); // GE non-DTI for (int i = 0; i < numDti; i++) { bvals[i] = vx[i].V[0]; //printMessage("---bxyz %g %g %g %g\n",vx[i].V[0],vx[i].V[1],vx[i].V[2],vx[i].V[3]); @@ -2423,8 +2568,9 @@ int *nii_saveDTI(char pathoutname[], int nConvert, struct TDCMsort dcmSort[], st } bvals[i] = bvals[i] + (0.5 * i / numDti); //add a small bias so ties are kept in sequential order } - if (numGEwarn > 0) - printWarning("Some images had bval>0 but bvec=0 (either Trace or b=0, see issue 245)\n"); + // See issue 777: removed the warning because GE DTI b=0 with bval>0 but bvec=0 (prior to version 29.1) will be handled by geCorrectBvecs() + // if (numGEwarn > 0) + // printWarning("Some images had bval>0 but bvec=0 (either Trace or b=0, see issue 245)\n"); /*if ((*numADC == numDti) || (numGEwarn == numDti)) { //issue 405: we now save bvals file for isotropic series //all isotropic/ADC images - no valid bvecs *numADC = 0; @@ -2768,7 +2914,7 @@ float computeGantryTiltPrecise(struct TDICOMdata d1, struct TDICOMdata d2, int i float dotX = dotProduct(slice_vector90, slice_vector); float cosX = dotX / (len * len90); float degX = acos(cosX) * (180.0 / M_PI); //arccos, radian -> degrees - if (!isSameFloat(cosX, 1.0)) + if (!isSameFloatGE(cosX, 1.0)) ret = degX; if ((isSameFloat(ret, 0.0)) && (isSameFloat(ret, d1.gantryTilt))) return 0.0; @@ -3175,6 +3321,51 @@ void cleanISO8859(char *cString) { cString[i] = 'y'; } } + +void createDummyBidsBoilerplate(char *pth, bool isFunc) { + //https://remi-gau.github.io/bids_cookbook/#starters + char pathSep[2] = {"a"}; + pathSep[0] = kPathSeparator; + char descfnm[PATH_MAX] = {""}; + char taskfnm[PATH_MAX] = {""}; + char fnm[PATH_MAX] = {""}; + strcat(fnm, pth); + strcat(fnm, pathSep); + strcat(taskfnm, fnm); + strcat(descfnm, fnm); + snprintf(fnm + strlen(fnm), PATH_MAX - strlen(fnm), "%s", "README.md"); + if (!is_fileexists(fnm)) { + FILE *fp = fopen(fnm, "w"); + static const char readmePre[] = "Generated using dcm2niix ("; + static const char readmePost[] = ")\n\nDescribe your dataset here. This file was generated by dcm2niix in a single pass. Details like IntendedFor, Subject ID, Session and tasks are not defined."; + + if (fp != NULL) + fprintf(fp, readmePre); + fprintf(fp, kDCMdate); + fprintf(fp, readmePost); + fclose(fp); + } + snprintf(descfnm + strlen(descfnm), PATH_MAX - strlen(descfnm), "%s", "dataset_description.json"); + if (!is_fileexists(descfnm)) { + FILE *fp = fopen(descfnm, "w"); + static const char readme[] = "{\n \"Name\": \"dcm2niix dummy dataset\",\n \"Authors\": [\"Chris Rorden\", \"Alex Teghipco\"],\n \"BIDSVersion\": \"1.6.0\"\n}\n"; + if (fp != NULL) + fprintf(fp, readme); + fclose(fp); + } + if (!isFunc) + return; //only functional data gets a task file + snprintf(taskfnm + strlen(taskfnm), PATH_MAX - strlen(taskfnm), "%s", "task-rest_bold.json"); + if (!is_fileexists(taskfnm)) { + FILE *fp = fopen(taskfnm, "w"); + static const char taskRest[] = "{\n\"TaskName\": \"rest\",\n\"CogAtlasID\": \"https://www.cognitiveatlas.org/task/id/trm_4c8a834779883/\"\n}\n"; + if (fp != NULL) + fprintf(fp, taskRest); + fclose(fp); + } + +} + int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts opts) { char pth[PATH_MAX] = {""}; if (strlen(opts.outdir) > 0) { @@ -3210,6 +3401,8 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts strcpy(inname, "T%t_N%n_S%s"); } const char kTempPathSeparator = '\a'; + char pathSep[2] = {"a"}; + pathSep[0] = kTempPathSeparator; for (size_t pos = 0; pos < strlen(inname); pos++) if ((inname[pos] == '\\') || (inname[pos] == '/')) inname[pos] = kTempPathSeparator; @@ -3252,6 +3445,49 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts strcat(outname, opts.indirParent); if (f == 'G') strcat(outname, dcm.accessionNumber); + if (f == 'H') { + printWarning("hazardous (%%h) bids naming experimental\n"); + char bidsSubject[kOptsStr] = "sub-"; + if (strlen(opts.bidsSubject) <= 0) + strcat(bidsSubject, "1"); + else + strcat(bidsSubject, opts.bidsSubject); + char bidsSession[kOptsStr] = "ses-"; + if (strlen(opts.bidsSession) <= 0) + strcat(bidsSession, "1"); + else + strcat(bidsSession, opts.bidsSession); + createDummyBidsBoilerplate(pth, (strstr(dcm.CSA.bidsDataType, "func") != NULL)); + if (strlen(dcm.CSA.bidsDataType) < 1) { + strcat(outname, "Unknown"); + snprintf(newstr, PATH_MAX, "%c", kTempPathSeparator); + strcat(outname, newstr); + snprintf(newstr, PATH_MAX, "%ld", dcm.seriesNum); + strcat(outname, newstr); + strcat(outname, "_"); + strcat(outname, dcm.protocolName); + + } else { + isAddNamePostFixes = false; + strcat(outname, bidsSubject); + strcat(outname, pathSep); + strcat(outname, bidsSession); + strcat(outname, pathSep); + strcat(outname, dcm.CSA.bidsDataType); + strcat(outname, pathSep); + strcat(outname, bidsSubject); + strcat(outname, "_"); + strcat(outname, bidsSession); + if (strstr(dcm.CSA.bidsDataType, "func") != NULL) { + strcat(outname, "_task-"); + if (strlen(dcm.CSA.bidsTask) > 0) + strcat(outname, dcm.CSA.bidsTask); + else + strcat(outname, "rest"); + } + strcat(outname, dcm.CSA.bidsEntitySuffix); + } + } if (f == 'I') strcat(outname, dcm.patientID); if (f == 'J') @@ -3330,21 +3566,37 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts if (f == 'V') { if (dcm.manufacturer == kMANUFACTURER_BRUKER) strcat(outname, "Bruker"); + else if (dcm.manufacturer == kMANUFACTURER_CANON) + strcat(outname, "Canon"); else if (dcm.manufacturer == kMANUFACTURER_GE) strcat(outname, "GE"); + else if (dcm.manufacturer == kMANUFACTURER_HYPERFINE) + strcat(outname, "Hyperfine"); + else if (dcm.manufacturer == kMANUFACTURER_MEDISO) + strcat(outname, "Mediso"); + else if (dcm.manufacturer == kMANUFACTURER_MRSOLUTIONS) + strcat(outname, "MRsolutions"); else if (dcm.manufacturer == kMANUFACTURER_PHILIPS) strcat(outname, "Philips"); else if (dcm.manufacturer == kMANUFACTURER_SIEMENS) strcat(outname, "Siemens"); else if (dcm.manufacturer == kMANUFACTURER_TOSHIBA) strcat(outname, "Toshiba"); - else if (dcm.manufacturer == kMANUFACTURER_CANON) - strcat(outname, "Canon"); else if (dcm.manufacturer == kMANUFACTURER_UIH) strcat(outname, "UIH"); else strcat(outname, "NA"); } + if (f == 'W') {//Weird includes personal data in filename patientWeight + snprintf(newstr, PATH_MAX, "dob%sg%cwt%d", dcm.patientBirthDate, dcm.patientSex, (int)round(dcm.patientWeight)); + if (strstr(dcm.institutionName, "Richland")) + strcat(newstr, "R"); + strcat(outname, newstr); + } + if ((f == 'Y') && (dcm.rawDataRunNumber >= 0)) { + snprintf(newstr, PATH_MAX, "%d", dcm.rawDataRunNumber); //GE (0019,10A2) else (0020,0100) + strcat(outname, newstr); + } if (f == 'X') strcat(outname, dcm.studyID); if ((f == 'Y') && (dcm.rawDataRunNumber >= 0)) { @@ -3391,6 +3643,13 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts //strcat (outname,newstr); strcat(outname, "_c"); strcat(outname, dcm.coilName); +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s_c%s", mrifsStruct.namePostFixes, dcm.coilName); +#endif + } + if ((isAddNamePostFixes) && (dcm.numberOfTR > 1) && (dcm.TR > 0)) { + snprintf(newstr, PATH_MAX, "_r%g", dcm.TR); + strcat(outname, newstr); } // myMultiEchoFilenameSkipEcho1 https://github.com/rordenlab/dcm2niix/issues/237 #ifdef myMultiEchoFilenameSkipEcho1 @@ -3404,15 +3663,24 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts snprintf(newstr, PATH_MAX, "_e%d", dcm.echoNum); strcat(outname, newstr); isEchoReported = true; +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s%s", mrifsStruct.namePostFixes, newstr); +#endif } if ((isAddNamePostFixes) && (!isSeriesReported) && (!isEchoReported) && (dcm.echoNum > 1)) { //last resort: user provided no method to disambiguate echo number in filename snprintf(newstr, PATH_MAX, "_e%d", dcm.echoNum); strcat(outname, newstr); isEchoReported = true; +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s%s", mrifsStruct.namePostFixes, newstr); +#endif } if ((dcm.isNonParallelSlices) && (!isImageNumReported)) { snprintf(newstr, PATH_MAX, "_i%05d", dcm.imageNum); strcat(outname, newstr); +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s%s", mrifsStruct.namePostFixes, newstr); +#endif } /*if (dcm.maxGradDynVol > 0) { //Philips segmented snprintf(newstr, PATH_MAX, "_v%04d", dcm.gradDynVol+1); //+1 as indexed from zero @@ -3420,21 +3688,36 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts }*/ if ((isAddNamePostFixes) && (dcm.isHasImaginary)) { strcat(outname, "_imaginary"); //has phase map +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s_imaginary", mrifsStruct.namePostFixes); +#endif } if ((isAddNamePostFixes) && (dcm.isHasReal) && (dcm.isRealIsPhaseMapHz)) { strcat(outname, "_fieldmaphz"); //has field map +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s_fieldmaphz", mrifsStruct.namePostFixes); +#endif } if ((isAddNamePostFixes) && (dcm.isHasReal) && (!dcm.isRealIsPhaseMapHz)) { strcat(outname, "_real"); //has phase map +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s_real", mrifsStruct.namePostFixes); +#endif } if ((isAddNamePostFixes) && (dcm.isHasPhase)) { strcat(outname, "_ph"); //has phase map if (dcm.isHasMagnitude) strcat(outname, "Mag"); //Philips enhanced with BOTH phase and Magnitude in single file +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s_ph", mrifsStruct.namePostFixes); +#endif } if ((isAddNamePostFixes) && (dcm.aslFlags == kASL_FLAG_NONE) && (dcm.triggerDelayTime >= 1) && (dcm.manufacturer != kMANUFACTURER_GE)) { //issue 336 GE uses this for slice timing snprintf(newstr, PATH_MAX, "_t%d", (int)roundf(dcm.triggerDelayTime)); strcat(outname, newstr); +#ifdef USING_DCM2NIIXFSWRAPPER + sprintf(mrifsStruct.namePostFixes, "%s%s", mrifsStruct.namePostFixes, newstr); +#endif } //could add (isAddNamePostFixes) to these next two, but consequences could be catastrophic if (dcm.isRawDataStorage) //avoid name clash for Philips XX_ files @@ -3458,12 +3741,18 @@ int nii_createFilename(struct TDICOMdata dcm, char *niiFilename, struct TDCMopts for (size_t pos = 0; pos < strlen(outname); pos++) if ((outname[pos] == '\\') || (outname[pos] == '/') || (outname[pos] == ' ') || (outname[pos] == '<') || (outname[pos] == '>') || (outname[pos] == ':') || (outname[pos] == ';') || (outname[pos] == '"') // || (outname[pos] == '/') || (outname[pos] == '\\') //|| (outname[pos] == '^') issue398 + || (outname[pos] == '$') //issue749 || (outname[pos] == '*') || (outname[pos] == '|') || (outname[pos] == '?')) outname[pos] = '_'; #else for (size_t pos = 0; pos < strlen(outname); pos++) if (outname[pos] == ':') //not allowed by MacOS outname[pos] = '_'; +#endif +#if !defined(_WIN64) || !defined(_WIN32) + for (size_t pos = 0; pos < strlen(outname); pos++) + if (outname[pos] == '`' || outname[pos] == '$') // unix shell expansion characters + outname[pos] = '_'; #endif cleanISO8859(outname); //re-insert explicit path separators: -f %t/%s_%p will have folder for time, but will not segment a protocol named "fMRI\bold" @@ -3583,7 +3872,7 @@ void nii_createDummyFilename(char *niiFilename, struct TDCMopts opts) { strcat(niiFilename, ".nhdr'"); else strcat(niiFilename, ".nrrd'"); - #ifdef myEnableJNIfTI + #ifdef myEnableJNIFTI } else if (opts.saveFormat == kSaveFormatJNII) { strcat(niiFilename, ".jnii'"); } else if (opts.saveFormat == kSaveFormatBNII) { @@ -3710,175 +3999,7 @@ void writeNiiGz(char *baseName, struct nifti_1_header hdr, unsigned char *src_bu } //writeNiiGz() #endif -#ifdef USING_R - -// Version of nii_saveNII() for R/divest: create nifti_image pointer and push onto stack -int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, struct TDCMopts opts, struct TDICOMdata d) { - hdr.vox_offset = 352; - // Extract the basename from the full file path - char *start = niiFilename + strlen(niiFilename); - while (start >= niiFilename && *start != '/' && *start != kPathSeparator) - start--; - std::string name(++start); - nifti_image *image = nifti_convert_nhdr2nim(hdr, niiFilename); - if (image == NULL) - return EXIT_FAILURE; - image->data = (void *)im; - ImageList *images = (ImageList *)opts.imageList; - images->append(image, name); - free(image); - return EXIT_SUCCESS; -} - -void nii_saveAttributes(struct TDICOMdata &data, struct nifti_1_header &header, struct TDCMopts &opts, const char *filename) { - ImageList *images = (ImageList *)opts.imageList; - switch (data.modality) { - case kMODALITY_CR: - images->addAttribute("modality", "CR"); - break; - case kMODALITY_CT: - images->addAttribute("modality", "CT"); - break; - case kMODALITY_MR: - images->addAttribute("modality", "MR"); - break; - case kMODALITY_PT: - images->addAttribute("modality", "PT"); - break; - case kMODALITY_US: - images->addAttribute("modality", "US"); - break; - } - switch (data.manufacturer) { - case kMANUFACTURER_SIEMENS: - images->addAttribute("manufacturer", "Siemens"); - break; - case kMANUFACTURER_GE: - images->addAttribute("manufacturer", "GE"); - break; - case kMANUFACTURER_MEDISO: - images->addAttribute("manufacturer", "Mediso"); - break; - case kMANUFACTURER_PHILIPS: - images->addAttribute("manufacturer", "Philips"); - break; - case kMANUFACTURER_TOSHIBA: - images->addAttribute("manufacturer", "Toshiba"); - break; - case kMANUFACTURER_UIH: - images->addAttribute("manufacturer", "UIH"); - break; - case kMANUFACTURER_BRUKER: - images->addAttribute("manufacturer", "Bruker"); - break; - case kMANUFACTURER_HITACHI: - images->addAttribute("manufacturer", "Hitachi"); - break; - case kMANUFACTURER_CANON: - images->addAttribute("manufacturer", "Canon"); - break; - case kMANUFACTURER_MRSOLUTIONS: - images->addAttribute("manufacturer", "MRSolutions"); - break; - case kMANUFACTURER_HYPERFINE: - images->addAttribute("manufacturer", "Hyperfine"); - break; - } - images->addAttribute("scannerModelName", data.manufacturersModelName); - images->addAttribute("imageType", data.imageType); - if (data.seriesNum > 0) - images->addAttribute("seriesNumber", int(data.seriesNum)); - images->addAttribute("seriesDescription", data.seriesDescription); - images->addAttribute("sequenceName", data.sequenceName); - images->addAttribute("protocolName", data.protocolName); - images->addDateAttribute("studyDate", data.studyDate); - images->addTimeAttribute("studyTime", data.studyTime); - images->addAttribute("fieldStrength", data.fieldStrength); - images->addAttribute("flipAngle", data.flipAngle); - images->addAttribute("echoTime", data.TE); - images->addAttribute("repetitionTime", data.TR); - images->addAttribute("inversionTime", data.TI); - if (!data.isXRay) { - images->addAttribute("sliceThickness", data.zThick); - images->addAttribute("sliceSpacing", data.zSpacing); - } - if (data.CSA.multiBandFactor > 1) - images->addAttribute("multibandFactor", data.CSA.multiBandFactor); - if (data.phaseEncodingSteps > 0) - images->addAttribute("phaseEncodingSteps", data.phaseEncodingSteps); - if (data.phaseEncodingLines > 0) - images->addAttribute("phaseEncodingLines", data.phaseEncodingLines); - // Calculations relating to the reconstruction in the phase encode direction, - // which are needed to derive effective echo spacing and readout time below. - // See the nii_SaveBIDS() function for details - int reconMatrixPE = data.phaseEncodingLines; - if ((header.dim[2] > 0) && (header.dim[1] > 0)) { - if (header.dim[1] == header.dim[2]) //phase encoding does not matter - reconMatrixPE = header.dim[2]; - else if (data.phaseEncodingRC == 'C') - reconMatrixPE = header.dim[2]; - else if (data.phaseEncodingRC == 'R') - reconMatrixPE = header.dim[1]; - } - double bandwidthPerPixelPhaseEncode = data.bandwidthPerPixelPhaseEncode; - if (bandwidthPerPixelPhaseEncode == 0.0) - bandwidthPerPixelPhaseEncode = data.CSA.bandwidthPerPixelPhaseEncode; - double effectiveEchoSpacing = 0.0; - if ((reconMatrixPE > 0) && (bandwidthPerPixelPhaseEncode > 0.0)) - effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * reconMatrixPE); - if (data.effectiveEchoSpacingGE > 0.0) { - double roundFactor = data.isPartialFourier ? 4.0 : 2.0; - double totalReadoutTime = ((ceil(1.0 / roundFactor * data.phaseEncodingLines / data.accelFactPE) * roundFactor) - 1.0) * data.effectiveEchoSpacingGE * 0.000001; - effectiveEchoSpacing = totalReadoutTime / (reconMatrixPE - 1); - } - images->addAttribute("effectiveEchoSpacing", effectiveEchoSpacing); - if (data.manufacturer == kMANUFACTURER_UIH) - images->addAttribute("effectiveReadoutTime", data.acquisitionDuration / 1000.0); - else if ((reconMatrixPE > 0) && (effectiveEchoSpacing > 0.0)) - images->addAttribute("effectiveReadoutTime", effectiveEchoSpacing * (reconMatrixPE - 1.0)); - images->addAttribute("pixelBandwidth", data.pixelBandwidth); - if ((data.manufacturer == kMANUFACTURER_SIEMENS) && (data.dwellTime > 0)) - images->addAttribute("dwellTime", data.dwellTime * 1e-9); - // Phase encoding polarity - // We only save these attributes if both direction and polarity are known - bool isSkipPhaseEncodingAxis = data.is3DAcq; - if (data.echoTrainLength > 1) - isSkipPhaseEncodingAxis = false; //issue 371: ignore phaseEncoding for 3D MP-RAGE/SPACE, but report for 3D EPI - if (((data.phaseEncodingRC == 'R') || (data.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && ((data.CSA.phaseEncodingDirectionPositive == 1) || (data.CSA.phaseEncodingDirectionPositive == 0))) { - if (data.phaseEncodingRC == 'C') { - images->addAttribute("phaseEncodingDirection", "j"); - // Notice the XOR (^): the sense of phaseEncodingDirectionPositive - // is reversed if we are flipping the y-axis - images->addAttribute("phaseEncodingSign", ((data.CSA.phaseEncodingDirectionPositive == 0) ^ opts.isFlipY) ? -1 : 1); - } else if (data.phaseEncodingRC == 'R') { - images->addAttribute("phaseEncodingDirection", "i"); - images->addAttribute("phaseEncodingSign", data.CSA.phaseEncodingDirectionPositive == 0 ? -1 : 1); - } - } - // Slice timing (stored in seconds) - if (data.CSA.sliceTiming[0] >= 0.0 && (data.manufacturer == kMANUFACTURER_UIH || data.manufacturer == kMANUFACTURER_GE || (data.manufacturer == kMANUFACTURER_SIEMENS && !data.isXA10A))) { - std::vector sliceTimes; - for (int i = 0; i < header.dim[3]; i++) { - if (data.CSA.sliceTiming[i] < 0.0) - break; - sliceTimes.push_back(data.CSA.sliceTiming[i] / 1000.0); - } - images->addAttribute("sliceTiming", sliceTimes); - } - images->addAttribute("patientIdentifier", data.patientID); - images->addAttribute("patientName", data.patientName); - images->addDateAttribute("patientBirthDate", data.patientBirthDate); - if (strlen(data.patientAge) > 0 && strcmp(data.patientAge, "000Y") != 0) - images->addAttribute("patientAge", data.patientAge); - if (data.patientSex == 'F') - images->addAttribute("patientSex", "F"); - else if (data.patientSex == 'M') - images->addAttribute("patientSex", "M"); - images->addAttribute("patientWeight", data.patientWeight); - images->addAttribute("comments", data.imageComments); -} - -#else +#ifndef USING_R int pigz_File(char *fname, struct TDCMopts opts, size_t imgsz) { //given "/dir/file.nii" creates "/dir/file.nii.gz" @@ -3893,13 +4014,17 @@ int pigz_File(char *fname, struct TDCMopts opts, size_t imgsz) { if ((opts.gzLevel > 0) && (opts.gzLevel < 12)) { char newstr[256]; snprintf(newstr, 256, "\"%s -n -f -%d \"", blockSize, opts.gzLevel); + //749 snprintf(newstr, 256, "\"%s -n -f -%d '", blockSize, opts.gzLevel); strcat(command, newstr); } else { char newstr[256]; snprintf(newstr, 256, "\"%s -n \"", blockSize); + //749 snprintf(newstr, 256, "\"%s -n '", blockSize); strcat(command, newstr); } strcat(command, fname); + //issue749 use single quote to prevent expansion of $ + //749 strcat(command, "'"); //add quotes in case spaces in filename 'pigz "c:\my dir\img.nii"' strcat(command, "\""); //add quotes in case spaces in filename 'pigz "c:\my dir\img.nii"' #if defined(_WIN64) || defined(_WIN32) //using CreateProcess instead of system to run in background (avoids screen flicker) DWORD exitCode; @@ -3974,7 +4099,7 @@ void writeMghGz(char *baseName, Tmgh hdr, TmghFooter footer, unsigned char *src_ //add image strm.avail_in = (unsigned int)src_len; // size of input strm.next_in = (uint8_t *)src_buffer; // input image -- TPX strm.next_in = (Bytef *)src_buffer; - deflate(&strm, Z_FINISH); + deflate(&strm, Z_FINISH); //add footer strm.avail_in = (unsigned int)sizeof(footer); // size of input strm.next_in = (uint8_t *) &footer.TR; @@ -4107,7 +4232,7 @@ int nii_saveMGH(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, #endif if (isGz) { strcat(fname, ".mgz"); - writeMghGz(fname, mgh, footer, im, imgsz, opts.gzLevel); + writeMghGz(fname, mgh, footer, im, imgsz, opts.gzLevel); } else { strcat(fname, ".mgh"); FILE *fp = fopen(fname, "wb"); @@ -4390,7 +4515,7 @@ int nii_saveNRRD(char *niiFilename, struct nifti_1_header hdr, unsigned char *im enum TZipMethod {zmZlib, zmGzip, zmBase64}; -#ifdef myEnableJNIfTI +#ifdef myEnableJNIFTI #ifdef Z_DEFLATED int zmat_run(const size_t inputsize, unsigned char *inputstr, size_t *outputsize, unsigned char **outputbuf, const int zipid, int *ret, const int iscompress){ @@ -4956,12 +5081,12 @@ int nii_savejnii(char *niiFilename, struct nifti_1_header hdr, unsigned char *im cJSON_Delete(root); return EXIT_SUCCESS; } // nii_savejnii() -#endif //#ifdef myEnableJNIfTI +#endif //#ifdef myEnableJNIFTI int nii_saveForeign(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, struct TDCMopts opts, struct TDICOMdata d, struct TDTI4D *dti4D, int numDTI) { if (opts.saveFormat == kSaveFormatMGH) return nii_saveMGH(niiFilename, hdr, im, opts, d, dti4D, numDTI); - #ifdef myEnableJNIfTI + #ifdef myEnableJNIFTI else if (opts.saveFormat == kSaveFormatJNII || opts.saveFormat == kSaveFormatBNII) return nii_savejnii(niiFilename, hdr, im, opts, d, dti4D, numDTI); #endif @@ -5026,9 +5151,25 @@ void removeSclSlopeInter(struct nifti_1_header *hdr, unsigned char *img) { //printWarning("NRRD unable to record scl_slope/scl_inter %g/%g\n", hdr->scl_slope, hdr->scl_inter); } -#ifndef USING_R - int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, struct TDCMopts opts, struct TDICOMdata d) { +#ifdef USING_R + ImageList *images = (ImageList *)opts.imageList; + if (opts.isImageInMemory) { + hdr.vox_offset = 352; + // Extract the basename from the full file path + char *start = niiFilename + strlen(niiFilename); + while (start >= niiFilename && *start != '/' && *start != kPathSeparator) + start--; + std::string name(++start); + nifti_image *image = nifti_convert_nhdr2nim(hdr, niiFilename); + if (image == NULL) + return EXIT_FAILURE; + image->data = (void *)im; + images->append(image, name); + free(image); + return EXIT_SUCCESS; + } +#else if (opts.isOnlyBIDS) return EXIT_SUCCESS; if (opts.saveFormat != kSaveFormatNIfTI) { @@ -5037,6 +5178,7 @@ int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, free(dti4D); return ret; } +#endif hdr.vox_offset = 352; size_t imgsz = nii_ImgBytes(hdr); if (imgsz < 1) { @@ -5065,6 +5207,9 @@ int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, true); //byte-swap endian (e.g. little->big) writeNiiGz(niiFilename, hdr, im, imgsz, opts.gzLevel, false); +#ifdef USING_R + images->appendPath(std::string(niiFilename) + ".nii.gz"); +#endif if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, false); //unbyte-swap endian (e.g. big->little) return EXIT_SUCCESS; @@ -5088,12 +5233,15 @@ int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, if ((opts.gzLevel > 0) && (opts.gzLevel < 12)) { char newstr[256]; snprintf(newstr, 256, "\" -n -f -%d > \"", opts.gzLevel); + //749 snprintf(newstr, 256, "\" -n -f -%d > '", opts.gzLevel); strcat(command, newstr); } else strcat(command, "\" -n -f > \""); //current versions of pigz (2.3) built on Windows can hang if the filename is included, presumably because it is not finding the path characters ':\' + //749 strcat(command, "\" -n -f > '"); //current versions of pigz (2.3) built on Windows can hang if the filename is included, presumably because it is not finding the path characters ':\' strcat(command, fname); + //issue749 single not double quotes so $ character does not cause issues + //749 strcat(command, ".gz'"); //add quotes in case spaces in filename 'pigz "c:\my dir\img.nii"' strcat(command, ".gz\""); //add quotes in case spaces in filename 'pigz "c:\my dir\img.nii"' - //strcat(command, "x.gz\""); //add quotes in case spaces in filename 'pigz "c:\my dir\img.nii"' if (opts.isVerbose) printMessage("Compress: %s\n", command); FILE *pigzPipe; @@ -5125,10 +5273,14 @@ int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, fwrite(&pad, sizeof(pad), 1, fp); fwrite(&im[0], imgsz, 1, fp); fclose(fp); -#endif if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, false); //unbyte-swap endian (e.g. big->little) +#endif + +#ifdef USING_R + images->appendPath(fname); +#else if ((opts.isGz) && (strlen(opts.pigzname) > 0)) { #ifndef myDisableGzSizeLimits if ((imgsz + hdr.vox_offset) > kMaxPigz) { @@ -5138,11 +5290,10 @@ int nii_saveNII(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, #endif return pigz_File(fname, opts, imgsz); } +#endif return EXIT_SUCCESS; } // nii_saveNII() -#endif - int nii_saveNIIx(char *niiFilename, struct nifti_1_header hdr, unsigned char *im, struct TDCMopts opts) { struct TDICOMdata dcm = clear_dicom_data(); return nii_saveNII(niiFilename, hdr, im, opts, dcm); @@ -6253,18 +6404,666 @@ void sliceTimingGE_Testx0021x105E(struct TDICOMdata *d, struct TDCMopts opts, st printMessage("\t%g\t%g\n", d->CSA.sliceTiming[v], sliceTiming[v]); } -void reportProtocolBlockGE(struct TDICOMdata *d, const char *filename) { +void reportProtocolBlockGE(struct TDICOMdata *d, const char *filename, int isVerbose) { #ifdef myReadGeProtocolBlock - if ((d->manufacturer != kMANUFACTURER_GE) || (d->protocolBlockStartGE < 1) || (d->protocolBlockLengthGE < 19)) + if ((d->manufacturer != kMANUFACTURER_GE) || (d->modality != kMODALITY_MR)) + return; + if ((d->protocolBlockStartGE < 1) || (d->protocolBlockLengthGE < 19)) { + //if (isVerbose) + printWarning("Missing GE protocol data block (0025,101B)\n"); return; + } int viewOrderGE = -1; int sliceOrderGE = -1; int mbAccel = -1; int nSlices = -1; float groupDelay = 0.0; char ioptGE[3000] = ""; - geProtocolBlock(filename, d->protocolBlockStartGE, d->protocolBlockLengthGE, 2, &sliceOrderGE, &viewOrderGE, &mbAccel, &nSlices, &groupDelay, ioptGE); + char seqName[kDICOMStr] = ""; + geProtocolBlock(filename, d->protocolBlockStartGE, d->protocolBlockLengthGE, isVerbose, &sliceOrderGE, &viewOrderGE, &mbAccel, &nSlices, &groupDelay, ioptGE, seqName); + //if (strlen(d->procedureStepDescription) < 2) + // strcpy(d->procedureStepDescription, seqName); + strcat(d->procedureStepDescription, seqName); //issue790 #endif +} //bidsGE + +void setBidsSiemens(struct TDICOMdata *d, int nConvert, int isVerbose, const char *filename) { + char seqDetails[kDICOMStrLarge] = ""; + char acqStr[kDICOMStrLarge] = ""; + char preAcqStr[kDICOMStrLarge] = ""; + float inv1 = NAN; + float inv2 = NAN; + bool isDualTI = false; + int lContrasts = 0; //detect me-mprage + #ifdef myReadAsciiCsa + if ((d->CSA.SeriesHeader_offset > 0) && (d->CSA.SeriesHeader_length > 0)) { + float pf = 1.0f; //partial fourier + float shimSetting[8]; + char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], wipMemBlock[kDICOMStrExtraLarge]; + TCsaAscii csaAscii; + siemensCsaAscii(filename, &csaAscii, d->CSA.SeriesHeader_offset, d->CSA.SeriesHeader_length, shimSetting, coilID, consistencyInfo, coilElements, seqDetails, fmriExternalInfo, protocolName, wipMemBlock); + inv1 = csaAscii.alTI[0] / 1000.0; + inv2 = csaAscii.alTI[1] / 1000.0; + lContrasts = csaAscii.lContrasts; + //If parameter lInvContrasts exists in the protocol, a value of 1 indicates MPRAGE and a value of 2 MP2RAGE. Note that lInvContrasts is different from lContrasts and that only lInvContrasts must be considered. + //If parameter lInvContrasts does not exist, then the presence of alTI[1] indicates that this is an MP2RAGE protocol. An MPRAGE protocol will only contain alTI[0]. + if (csaAscii.lInvContrasts == 1) //explicitly reports one TI + inv2 = NAN; + if ((!isnan(inv1)) && (!isnan(inv2)) && (inv1 > 0.0) && (inv2 > 0.0)) + isDualTI = true; + + } + #endif // myReadAsciiCsa + char *dataTypeBIDS = d->CSA.bidsDataType; + strcpy(dataTypeBIDS, ""); + char *suffixBIDS = d->CSA.bidsEntitySuffix; + strcpy(suffixBIDS, ""); + char modalityBIDS[kDICOMStrLarge] = ""; //T1w, bold, dwi + char recBIDS[kDICOMStrLarge] = ""; //_desc-preproc + if (d->manufacturer != kMANUFACTURER_SIEMENS) + return; + bool isReportEcho = true; //do not report _echo for PD/T2 pair or fieldmap + bool isMultiEcho = false; + bool isDerived = d->isDerived; //report phase encoding direction + bool isDirLabel = false; + bool isPart = false; + if (d->isHasPhase) + isPart = true; + bool isAddSeriesToRun = true; //except phasemap, where phasediff and magnitude have different series numbers but must share acq + char seqName[kDICOMStrLarge]; + strcpy(seqName, d->sequenceName); + if (strlen(d->sequenceName) < 2) //e.g. XA uses 0018,9005 while VE uses 0018,0024 + strcpy(seqName, d->pulseSequenceName); + if (strlen(seqDetails) < 2) + strcpy(seqDetails, seqName); + if (strstr(d->imageType, "DERIVED")) { + isDerived = true; //to do: respond to derived images + } + if (d->modality != kMODALITY_MR) + return; + if (((d->xyzDim[3] < 2) && (nConvert < 1)) || (d->isLocalizer)) { //need nConvert or nifti header + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "localizer"); + } else if (strstr(seqDetails,"b1map")) { + //issue 751 nb both T1 and b1map can use tfl base + //https://bids-specification.readthedocs.io/en/stable/appendices/qmri.html#tb1tfl-and-tb1rfm-specific-notes + strcpy(dataTypeBIDS, "fmap"); + strcpy(modalityBIDS, "TB1TFL"); + if ((strstr(d->imageType, "FLIP ANGLE MAP")) || (strstr(d->imageType, "FLIP ANGLE MAP"))) + strcpy(preAcqStr, "famp"); + else + strcpy(preAcqStr, "anat"); + } else if ((strstr(seqDetails, "tfl") != NULL) || (strstr(seqDetails, "mp2rage") != NULL) || (strstr(seqDetails, "wip925") != NULL)) { //prog_mprage + strcpy(dataTypeBIDS, "anat"); + if (isDualTI) + strcpy(modalityBIDS, "MP2RAGE"); + else { //bizarrely mprage can populate both alTI[0] alTI[1] see dcm_qa_xa30 + strcpy(modalityBIDS, "T1w"); + //bork inv2 = NAN; + } + if (strstr(d->imageType, "T1 MAP") != NULL) { + strcpy(modalityBIDS, "T1map"); + isDualTI = false; //issue 750 derived from two images + } + if (strstr(d->imageType, "_UNI") != NULL) { + isDualTI = false; //issue 750 derived from two images + strcpy(modalityBIDS, "UNIT1"); + if (strstr(d->imageComments, "DENOISED IMAGE") != NULL) + strcat(recBIDS, "denoise"); + } + isDerived = false; //issue750 Siemens E11 considers UNIT1 derived, but BIDS does not + } else if ((d->CSA.numDti > 0) || (strstr(seqDetails, "_diff") != NULL) || (strstr(seqDetails, "resolve") != NULL) || (strstr(seqDetails, "PtkSmsVB13ADwDualSpinEchoEpi") != NULL) || (strstr(seqDetails, "ep2d_stejskal_386") != NULL)) { //prog_diff + strcpy(dataTypeBIDS, "dwi"); + strcpy(modalityBIDS, "dwi"); + if (strstr(d->seriesDescription, "_SBRef") != NULL) + strcpy(modalityBIDS, "sbref"); + //if seriesDesc trace", "fa", "adc" isDerived = true; + isDirLabel = true; + } else if ((strstr(seqDetails, "fairest")) || (strstr(seqDetails, "_asl") != NULL) || (strstr(seqDetails, "_pasl") != NULL) || (strstr(seqDetails, "pcasl") != NULL) || (strstr(seqDetails, "PCASL") != NULL)) { //prog_asl + strcpy(dataTypeBIDS, "perf"); + strcpy(modalityBIDS, "asl"); + if (strstr(d->seriesDescription, "_m0") != NULL) + strcpy(modalityBIDS, "m0scan"); + } else if (strstr(d->pulseSequenceName, "spcR") != NULL) { + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "T2w"); + } else if (strstr(seqDetails, "tse_vfl") != NULL) { //prog_tse_vfl + strcpy(dataTypeBIDS, "anat"); + if ((strstr(seqDetails, "spcir") != NULL) || (strstr(d->sequenceName, "spcir") != NULL)) + strcpy(modalityBIDS, "FLAIR"); + else + strcpy(modalityBIDS, "T2w"); + } else if (strstr(seqDetails, "tse") != NULL) { //prog_tse + isReportEcho = false; //do not report _echo for T2/PDw pair + strcpy(dataTypeBIDS, "anat"); + if (strstr(d->sequenceName, "tir") != NULL) + strcpy(modalityBIDS, "FLAIR"); + if ((strstr(d->sequenceName, "tse2d") != NULL) || (strstr(d->pulseSequenceName, "tse2d") != NULL)) { + if (d->TE < 50) + strcpy(modalityBIDS, "PDw"); + else + strcpy(modalityBIDS, "T2w"); + } + } else if ((strstr(seqDetails, "ep2d_ase") != NULL)) { //prog_ep2d_se + // oxygen extraction fraction(OEF) Asymmetric Spin Echo (ASE) + //printWarning("BIDS does not yet specify `ase` data\n"); + //n.b. we do not specify dataTypeBIDS as not yet defined + strcpy(modalityBIDS, "oef_ase"); + } else if ((strstr(seqDetails, "ep2d_se") != NULL)) { //prog_ep2d_se + strcpy(dataTypeBIDS, "fmap"); + strcpy(modalityBIDS, "epi"); + isDirLabel = true; + } else if ((strstr(seqDetails, "gre_field_mapping") != NULL)) { //prog_fmap + isReportEcho = false; //echo encoded in "_magnitude" + isAddSeriesToRun = false; + isPart = false; + strcpy(dataTypeBIDS, "fmap"); + if (d->isHasPhase) + strcpy(modalityBIDS, "phasediff"); + if (d->isHasMagnitude) + snprintf(modalityBIDS, kDICOMStrLarge - strlen(modalityBIDS), "magnitude%d", d->echoNum); + //printf("magnitude%d\n", d->echoNum); + //} else if (( seqName,"_tfl2d1") == 0) || (strcmp(dcmList[indx].sequenceName, "_fl3d1_ns") == 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1" + } else if ((strstr(seqDetails, "\\trufi") != NULL) || (strstr(seqName, "fl3d1_ns") != NULL) || (strstr(seqName, "fl2d1") != NULL) || (strstr(seqName, "tfl2d1") != NULL)) { //localizers + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "localizer"); + } else if ((strstr(seqName, "fl3d1r") != NULL) || (strstr(seqName, "fl2d1") != NULL) || (strstr(seqName, "tfl2d1") != NULL)) { //localizers + //nb ToF fl3d1r_ts fl3d1r_t70 fl3d1r7t but check for fl3d1r SWI + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "angio"); + } else if (strstr(seqDetails, "ep_seg_fid") != NULL) { + //n.b. large echoTrainLength even for single echo acquisition + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "T2starw"); + isPart = true; + if (strstr(d->seriesDescription, "mIP") != NULL) { //derived minimum intensity + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "mIP"); + } + if (strstr(d->seriesDescription, "SWI_Images") != NULL) { //derived SWI + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "SWI_Images"); + } + } else if (strstr(seqDetails, "gre") != NULL) { //prog_gre + //printf("GRE T2starw or if MEGRE\n"); + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "T2starw"); + if ((d->echoNum > 1) || (lContrasts > 1)) { + strcpy(modalityBIDS, "MEGRE"); + isMultiEcho = true; + } + isPart = true; + } else if ((strstr(seqDetails, "AALScout") != NULL) || (strstr(seqDetails, "haste") != NULL)) { //localizer: unused + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "localizer"); + } else if ((strstr(seqDetails, "_bold")) || (strstr(seqDetails, "pace")) + || (strstr(d->imageType, "FMRI")) || (strstr(seqDetails, "ep2d_fid"))) { //prog_bold + //n.b. "Space" is not "pace" + strcpy(dataTypeBIDS, "func"); + strcpy(modalityBIDS, "bold"); + isDirLabel = true; + char* p = strstr(d->protocolName, "func_"); + if (p == d->protocolName) + //todo issue753 infer task from d->protocolName + if (strstr(d->seriesDescription, "_SBRef") != NULL) + strcpy(modalityBIDS, "sbref"); + } else if (strstr(d->sequenceName, "*epse2d") != NULL) { + //pepolar? + strcpy(dataTypeBIDS, "fmap"); + strcpy(modalityBIDS, "epi"); + isDirLabel = true; + } + strcpy(acqStr, "_acq-"); + strcat(acqStr, preAcqStr); + int len = strlen(seqName); + if (len > 0) { + for (int i = 0; i < len; i++) { + char ch = seqName[i]; + if (ch == '*') continue; + if ((ch >= '0' && ch <= '9') || (ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z')) + strncat(acqStr, &ch, 1); + if (isdigit(ch)) break; + } + } //len > 0 + if (d->accelFactPE > 1) + snprintf(acqStr + strlen(acqStr), kDICOMStrLarge - strlen(acqStr), "p%d", (int)round(d->accelFactPE)); + if (d->CSA.multiBandFactor > 1) + snprintf(acqStr + strlen(acqStr), kDICOMStrLarge - strlen(acqStr), "m%d", (int)round(d->CSA.multiBandFactor)); + strcat(suffixBIDS, acqStr); + //add _rec + // https://bids-specification.readthedocs.io/en/stable/05-derivatives/02-common-data-types.html#common-file-level-metadata-fields + if (strlen(recBIDS) > 0) { + strcat(suffixBIDS, "_rec-"); + strcat(suffixBIDS, recBIDS); + } + //add dir + //nb for func dir comes before echo + int phPos = d->CSA.phaseEncodingDirectionPositive; + if (isDirLabel) { + char dirLabel[kDICOMStrLarge] = "_dir-"; + if (d->phaseEncodingRC == 'C') { + if (phPos) + strcat(dirLabel, "AP"); + else + strcat(dirLabel, "PA"); + } else { + if (phPos) + strcat(dirLabel, "RL"); + else + strcat(dirLabel, "LR"); + } + strcat(suffixBIDS, dirLabel); + } + //add run + if (isAddSeriesToRun) + snprintf(suffixBIDS + strlen(suffixBIDS), kDICOMStrLarge - strlen(suffixBIDS), "_run-%ld", d->seriesNum); + //add echo + if (isReportEcho) + if ((d->echoNum > 1) || (isMultiEcho) || ((d->isMultiEcho) && (d->echoNum > 0))) + snprintf(suffixBIDS + strlen(suffixBIDS), kDICOMStrLarge - strlen(suffixBIDS), "_echo-%d", d->echoNum); + //add inv + // https://bids-specification.readthedocs.io/en/stable/appendices/entities.html#inv + if (isDualTI) { + //n.b. Siemens uses alTI[0]/alTI[1] for inversion times: alTI[2] used by ASL + int invIdx = 1; + if (isSameFloatGE(d->TI, inv2)) + invIdx = 2; + snprintf(suffixBIDS + strlen(suffixBIDS), kDICOMStrLarge - strlen(suffixBIDS), "_inv-%d", invIdx); + } + //add part + if (isPart) { + if (d->isHasPhase) + strcat(suffixBIDS, "_part-phase"); + if (d->isHasMagnitude) + strcat(suffixBIDS, "_part-mag"); + } + if (strlen(modalityBIDS) > 0) { + strcat(suffixBIDS, "_"); + strcat(suffixBIDS, modalityBIDS); + } + if ((isVerbose > 0) || (strlen(dataTypeBIDS) < 1)) + printf("::autoBids:Siemens CSAseqFname:'%s' pulseSeq:'%s' seqName:'%s'\n", + seqDetails, d->pulseSequenceName, d->sequenceName); + if (isDerived) + strcpy(dataTypeBIDS, "derived"); + //bork - ARC data follows + /* + if (strstr(dataTypeBIDS, "dwi")) { + if (strstr(d->protocolName, "12 dirs")) + strcpy(dataTypeBIDS, "discard"); + } + if (strstr(dataTypeBIDS, "fmap")) + strcpy(dataTypeBIDS, "discard"); + if (strstr(dataTypeBIDS, "perf")) + strcpy(dataTypeBIDS, "discard"); + if (strstr(dataTypeBIDS, "func")) { + if (d->TR > 9999) { + if (nConvert < 60) + strcpy(dataTypeBIDS, "discard"); + strcpy(d->CSA.bidsTask, "naming40"); + } else { + if (!strcasestr(d->protocolName, "REST")) + strcpy(dataTypeBIDS, "discard"); + } + if (nConvert < 10) + strcpy(dataTypeBIDS, "discard"); + }*/ +} // setBidsSiemens() + +void setBidsPhilips(struct TDICOMdata *d, int nConvert, int isVerbose) { + if (d->manufacturer != kMANUFACTURER_PHILIPS) + return; + char *dataTypeBIDS = d->CSA.bidsDataType; + strcpy(dataTypeBIDS, ""); + char *suffixBIDS = d->CSA.bidsEntitySuffix; + strcpy(suffixBIDS, ""); + if (d->modality != kMODALITY_MR) + return; + char seqName[kDICOMStr] = ""; + strcpy(seqName, d->sequenceVariant); + char modalityBIDS[kDICOMStrLarge] = ""; //_acq-FL3p2m2 + bool isReportEcho = true; + bool isDirLabel = false; //report phase encoding direction + bool isDerived = d->isDerived; //report phase encoding direction + bool isAddSeriesToRun = true; + bool isPart = false; + //d->pulseSequenceName, d->scanningSequence, d->sequenceVariant) + if (((d->xyzDim[3] < 4) && (nConvert < 4)) || (d->isLocalizer)) { + strcpy(dataTypeBIDS, "discard"); + strcpy(modalityBIDS, "localizer"); + } else if (strstr(seqName, "MP") != NULL) { + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "T1w"); + } else if ((d->isDiffusion) && (strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "SE") != NULL) ) { + strcpy(dataTypeBIDS, "dwi"); + strcpy(modalityBIDS, "dwi"); + } else if (strstr(d->imageType, "PERFUSION") != NULL) { + //scanSeq:'GR' seqVariant:'SK' + strcpy(dataTypeBIDS, "perf"); + strcpy(modalityBIDS, "asl"); + } else if ((strstr(d->pulseSequenceName, "SEEPI") != NULL) && (!d->isDiffusion) && (strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "SE") != NULL) ) { + //pepolar? d->pulseSequenceName + isAddSeriesToRun = false; + strcpy(dataTypeBIDS, "fmap"); + strcpy(modalityBIDS, "epi"); + printWarning("Unable to estimate BIDS `_dir` for fmap epi as Philips DICOMs do not report phase encoding polarity\n"); + isDirLabel = true; + } else if ((!d->isDiffusion) && (strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "SE") != NULL) ) { + strcpy(dataTypeBIDS, "anat"); + if (false)//((strstr(d->scanningSequence, "IR") != NULL)) + strcpy(modalityBIDS, "FLAIR"); + else if (d->TE < 40) + strcpy(modalityBIDS, "PDw"); + else + strcpy(modalityBIDS, "T2w"); + } else if ((strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "IR") != NULL) ) { + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "FLAIR"); + } else if ((strstr(d->imageType, "PERFUSION") != NULL) && (strstr(d->pulseSequenceName, "FEEPI") != NULL) && (strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "GR") != NULL) ) { + strcpy(dataTypeBIDS, "perf"); + strcpy(modalityBIDS, "asl"); + } else if (((d->rawDataRunNumber >= 1) || (strstr(d->pulseSequenceName, "FEEPI") != NULL)) && (strstr(seqName, "SK") != NULL) && (strstr(d->scanningSequence, "GR") != NULL) ) { + //nb older Philips data does not record EPI so use 2005,1063 fMRIStatusIndication + strcpy(dataTypeBIDS, "func"); + strcpy(modalityBIDS, "bold"); + isDirLabel = true; + } else if ((strstr(seqName, "SS") != NULL) && (strstr(d->scanningSequence, "GR") != NULL) ) { + strcpy(dataTypeBIDS, "anat"); + strcpy(modalityBIDS, "T2starw"); + isPart = true; + } else if ((d->isRealIsPhaseMapHz) && (strstr(seqName, "SS") != NULL) && (strstr(d->scanningSequence, "RM") != NULL)) { + //https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#expressing-the-mr-protocol-intent-for-fieldmaps + //isRealIsPhaseMapHz for 'Case 3: Direct field mapping' + //sub-