shims are signed (#608)

Canon/README.md

@@ -6,18 +6,26 @@ dcm2niix can convert Canon (né Toshiba) DICOM format images to NIfTI. This page
 
 Users of Canon MRI equipment are strongly advised to export data from their scanners as enhanced DICOM (with all images from the series stored as a single file) rather than classic DICOM (each 2D slice stored as a separate file). Limitations of the Canon classic DICOMs are described [here](https://github.com/rordenlab/dcm2niix/issues/495) and [here](https://github.com/neurolabusc/dcm_qa_canon).
 
+## Avoid Enhanced DICOM
+
+Users of Canon MRI equipment are strongly advised to export data from their scanners as classic DICOM (with each slice from the series stored as a separate files) rather than enhanced DICOM. The enhanced 4D sequences such as fMRI incorrectly sets the public tag TemporalPositionIndex (0020,9128) as 1 for all volumes: both (0020,9157) and (0020,9128). This is a violation of the DICOM standard and interferes with attempts to sort data into the correct temporal order. Any software or user who assumes these tags are truthful will fail.
+
+## The Enhanced versus Classic DICOM dilemma
+
+The prior two sections provide conflicting advice, due to limitations in both the classic and enhanced datasets generated by Canon instruments. Users of Canon equipment should lobby the manufacturer to honor their DICOM conformance statement. DICOM data from Canon instruments may not be handled correctly by any software tool including dcm2niix. This reflects a limitation in the source DICOM data, not dcm2niix.
+
 ## Diffusion Weighted Imaging Notes
 
 In contrast to several other vendors, Toshiba used public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. These coordinates are with respect to the scanner bore, not image space.
 
-Since the acquisition by Canon, these public tags are no longer populated for images saved in classic 2D DICOM format. The diffusion gradient directions are now stored in the ASCII Image Comments tag. Like GE (but unlike [Siemens, GE and Toshiba](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI)), these directions are with respect to the image space, not the scanner bore. Further, gradient direction is not adjusted for phase encoding polarity, and it is impossible to determine phase encoding polarity. For detailed discussion and a validation dataset that exhibits these attributes please see [dcm_qa_canon](https://github.com/neurolabusc/dcm_qa_canon). A Canon classic DICOM DWI image may report:
+Be aware that Canon software V6.0 stored diffusion directions in the ASCII Image Comments tag when exporting to classic DICOM. For detailed discussion and a validation dataset that exhibits these attributes please see [dcm_qa_canon](https://github.com/neurolabusc/dcm_qa_canon). A Canon classic DICOM DWI image may report:
 
 ```
  (0018,9087) FD 1500                                     #   8, 1 DiffusionBValue
  (0020,4000) LT [b=1500(0.445,0.000,0.895)]              #  26, 1 ImageComments
 ```
 
-In contrast, when exporting images as enhanced (4D) DICOM, information is stored in public tags and does appear to compensate for phase encode polarity. These coordinates are with respect to the scanner bore, not image space. A Canon classic DICOM DWI image may report:
+In contrast, Canon software [V6.1](https://github.com/neurolabusc/dcm_qa_canon_61) uses pulbic tags for both classic and enhanced DICOMs. A Canon V6.1 DICOM DWI image may report:
 
 ```
 (0018,9087) FD 1500                                     #   8, 1 DiffusionBValue

Philips/README.md

@@ -148,7 +148,7 @@ The BIDS tag [PartialFourier](https://bids-specification.readthedocs.io/en/stabl
 
 Philips DICOMs do not allow one to determine the temporal order of volumes for diffusion (DWI, DTI) and arterial spin labelling (ASL) sequences. This can hinder tools that attempt to model motion motion or T1 effects. For ASL data, dcm2niix will attempt to store volumes in the order phase < control/label < repeat. [For ASL examples and a table describing the ordering hierarchy, see this dataset](https://github.com/neurolabusc/dcm_qa_philips_asl). For diffusion images, dcm2niix may order Philips volumes differently depending on if the data is stored as classic or enhanced DICOM. For enhanced DICOM, dcm2niix follows the hierarchy specified by DimensionIndexValues (0020,9157). On the other hand, for classic DICOMs, volumes with identical b-value index (2005,1412) will be stored sequentially, with ties sorted based on gradient direction number (2005,1413). [For diffusion examples, see series 22 and 23 from this dataset](https://github.com/neurolabusc/dcm_qa_philips_enh). These two different sorting methods may not necessarily sort data identically.
 
-[Slice timing correction](https://www.mccauslandcenter.sc.edu/crnl/tools/stc) can account for some variability in fMRI datasets. Unfortunately, Philips DICOM data [does not explicitly encode slice timing information](https://neurostars.org/t/heudiconv-no-extraction-of-slice-timing-data-based-on-philips-dicoms/2201/4). Therefore, dcm2niix is unable to populate the "SliceTiming" BIDS field. However, the [philips_order.py script](https://neurostars.org/t/how-dcm2niix-handles-different-imaging-types/22697/7) may be able to determine slice timing using an undocumented and unverified technique. Additionally, one can typically infer slice timing by recording the [mode and number of packages](https://en.wikibooks.org/w/index.php?title=SPM/Slice_Timing&stable=0#Philips_scanners) reported for the sequence on the scanner console or the [sequence PDF](http://adni.loni.usc.edu/wp-content/uploads/2017/09/ADNI-3-Basic-Philips-R5.pdf). For precise timing, it is also worth knowing if equidistant "temporal slice spacing" is set and whether "prospect. motion corr." is on or off (if on, a short delay occurs between volumes).
+[Slice timing correction](https://www.mccauslandcenter.sc.edu/crnl/tools/stc) can account for some variability in fMRI datasets. Unfortunately, Philips DICOM data [does not explicitly encode slice timing information](https://neurostars.org/t/heudiconv-no-extraction-of-slice-timing-data-based-on-philips-dicoms/2201/4). Therefore, dcm2niix is unable to populate the "SliceTiming" BIDS field. Note, the [philips_order.py script](https://neurostars.org/t/how-dcm2niix-handles-different-imaging-types/22697/7) may be able to determine slice timing using an undocumented and unverified technique that works in some instances but certainly fails with others. Additionally, one can typically infer slice timing by recording the [mode and number of packages](https://en.wikibooks.org/w/index.php?title=SPM/Slice_Timing&stable=0#Philips_scanners) reported for the sequence on the scanner console or the [sequence PDF](http://adni.loni.usc.edu/wp-content/uploads/2017/09/ADNI-3-Basic-Philips-R5.pdf). For precise timing, it is also worth knowing if equidistant "temporal slice spacing" is set and whether "prospect. motion corr." is on or off (if on, a short delay occurs between volumes).
 
 Likewise, the BIDS tag "PhaseEncodingDirection" allows tools like [eddy](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy) and [TOPUP](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup) to undistort images. While the Philips DICOM header distinguishes the phase encoding axis (e.g. anterior-posterior vs left-right) it does not encode the polarity (A->P vs P->A).
 

console/nii_dicom.cpp

@@ -1138,6 +1138,14 @@ double dcmFloatDouble(const size_t lByteLength, const unsigned char lBuffer[], c
 } //dcmFloatDouble()
 #endif
 
+// SS/IS datatype
+int16_t dcmIntSS(int lByteLength, unsigned char lBuffer[], bool littleEndian) { //read binary 16 or 32 bit integer
+	if (littleEndian)
+		return (uint16_t)lBuffer[0] | ((uint16_t)lBuffer[1] << 8); //shortint vs word?
+	return (uint16_t)lBuffer[1] | ((uint16_t)lBuffer[1] << 0); //shortint vs word?
+} //dcmInt()
+
+//UL/US unsigned integer
 int dcmInt(int lByteLength, unsigned char lBuffer[], bool littleEndian) { //read binary 16 or 32 bit integer
 	if (littleEndian) {
 		if (lByteLength <= 3)
@@ -4166,7 +4174,7 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
 #define kDirectoryRecordSequence 0x0004 + (0x1220 << 16)
 //#define kSpecificCharacterSet 0x0008+(0x0005 << 16 ) //someday we should handle foreign characters...
 #define kImageTypeTag 0x0008 + (0x0008 << 16)
-//#define kSOPInstanceUID 0x0008+(0x0018 << 16 ) //Philips inserts time as last item, e.g. ?.?.?.YYYYMMDDHHmmSS.SSSS
+#define kSOPInstanceUID 0x0008+(0x0018 << 16 ) //Philips inserts time as last item, e.g. ?.?.?.YYYYMMDDHHmmSS.SSSS
 // not reliable https://neurostars.org/t/heudiconv-no-extraction-of-slice-timing-data-based-on-philips-dicoms/2201/21
 #define kStudyDate 0x0008 + (0x0020 << 16)
 #define kAcquisitionDate 0x0008 + (0x0022 << 16)
@@ -4557,6 +4565,7 @@ const uint32_t kEffectiveTE = 0x0018 + (0x9082 << 16);
 	float vRLPhilips = 0.0;
 	float vAPPhilips = 0.0;
 	float vFHPhilips = 0.0;
+	double acquisitionTimePhilips = -1.0;
 	bool isPhase = false;
 	bool isReal = false;
 	bool isImaginary = false;
@@ -4667,6 +4676,7 @@ const uint32_t kEffectiveTE = 0x0018 + (0x9082 << 16);
 				//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
 				if ((nDimIndxVal > 1) && (volumeNumber > 0) && (inStackPositionNumber > 0) && ((d.aslFlags == kASL_FLAG_PHILIPS_LABEL) || (d.aslFlags == kASL_FLAG_PHILIPS_CONTROL))) {
+
 					isKludgeIssue533 = true;
 					for (int i = 0; i < nDimIndxVal; i++)
 						d.dimensionIndexValues[i] = 0;
@@ -4679,6 +4689,11 @@ const uint32_t kEffectiveTE = 0x0018 + (0x9082 << 16);
 					nDimIndxVal = 4; //slice < phase < control/label < volume
 					//printf("slice %d phase %d control/label %d repeat %d\n", inStackPositionNumber, d.phaseNumber, d.aslFlags == kASL_FLAG_PHILIPS_LABEL, volumeNumber);
 				}
+				if ((volumeNumber == 1) && (acquisitionTimePhilips >= 0.0) && (inStackPositionNumber > 0)) {
+					d.CSA.sliceTiming[inStackPositionNumber - 1] = acquisitionTimePhilips;
+					printf("%d\t%f\n", inStackPositionNumber, acquisitionTimePhilips);
+					acquisitionTimePhilips = - 1.0;
+				}
 				int ndim = nDimIndxVal;
 				if (inStackPositionNumber > 0) {
 					//for images without SliceNumberMrPhilips (2001,100A)
@@ -5222,6 +5237,32 @@ const uint32_t kEffectiveTE = 0x0018 + (0x9082 << 16);
 			dcmStr(lLength, &buffer[lPos], acquisitionDateTimeTxt);
 			//printMessage("%s\n",acquisitionDateTimeTxt);
 			break;
+/* Failed attempt to infer slice timing for Philips MRI
+https://neurostars.org/t/how-dcm2niix-handles-different-imaging-types/22697/6
+
+		case kSOPInstanceUID: {
+			if (d.manufacturer != kMANUFACTURER_PHILIPS) break;
+			char uid[kDICOMStrLarge];
+			dcmStr(lLength, &buffer[lPos], uid, true);
+			char *timeStr = strrchr(uid, '.');
+			//nb Manufactuer (0008,0070) comes AFTER (0008,0018) SOPInstanceUID.
+			//format of (0008,0018) UI 
+			//[1.23.4.2019051416101221842
+			//       .YYYYMMDDHHmmssxxxxx
+			timeStr++; //skip "."
+			if ((strlen(timeStr) != 19) || (strlen(d.studyDate) < 8)) break;
+			bool sameDay = true;
+			for (int z = 0; z < 8; z++)
+				if (timeStr[z] != d.studyDate[z]) sameDay = false;
+			if (!sameDay)
+				printf("SOPInstanceUID does not match StudyDate: assuming study cross midnight\n");
+			char *hourStr = timeStr + 8; //Skip 8 charactersYear,Month,Day YYYYMMDD
+			acquisitionTimePhilips = (double) atof(hourStr) * (double) 0.00001;
+			//printf("   %s  %s  %f\n", timeStr, hourStr, acquisitionTimePhilips);
+
+			break;
+		}
+*/
 		case kStudyDate:
 			dcmStr(lLength, &buffer[lPos], d.studyDate);
 			if (((int)strlen(d.studyDate) > 7) && (strstr(d.studyDate, "19000101") != NULL))
@@ -6650,17 +6691,17 @@ const uint32_t kEffectiveTE = 0x0018 + (0x9082 << 16);
 		case kShimGradientX: 
 			if (d.manufacturer != kMANUFACTURER_GE)
 				break;
-			d.shimGradientX = dcmInt(lLength, &buffer[lPos], d.isLittleEndian);
+			d.shimGradientX = dcmIntSS(lLength, &buffer[lPos], d.isLittleEndian);
 			break;
 		case kShimGradientY: 
 			if (d.manufacturer != kMANUFACTURER_GE)
 				break;
-			d.shimGradientY = dcmInt(lLength, &buffer[lPos], d.isLittleEndian);
+			d.shimGradientY = dcmIntSS(lLength, &buffer[lPos], d.isLittleEndian);
 			break;
 		case kShimGradientZ: 
 			if (d.manufacturer != kMANUFACTURER_GE)
 				break;
-			d.shimGradientZ = dcmInt(lLength, &buffer[lPos], d.isLittleEndian);
+			d.shimGradientZ = dcmIntSS(lLength, &buffer[lPos], d.isLittleEndian);
 			break;
 		case kPrescanReuseString: //LO
 			if (d.manufacturer != kMANUFACTURER_GE)

console/nii_dicom_batch.cpp

@@ -5705,6 +5705,8 @@ void checkSliceTiming(struct TDICOMdata *d, struct TDICOMdata *d1, int verbose,
 	//modified 20190704: this function now ensures all slice times are in msec
 	if ((d->TR < 0.0) || (d->CSA.sliceTiming[0] < 0.0))
 		return; //no slice timing
+	if (d->manufacturer == kMANUFACTURER_PHILIPS)
+		return; //Philips does not provide slice timing details in DICOM
 	if (d->manufacturer == kMANUFACTURER_GE)
 		return; //compute directly from Protocol Block
 	if (d->modality == kMODALITY_PT)