diff --git a/README.md b/README.md index 1b89b075..fd553ec8 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,9 @@ -# go2com (DICOM and NIFTI image parser) +# go2com (DICOM image parser) ** This package is under active development and can be in a very broken state. Please use the latest released version ** ## TODO -- [ ] Improve NIfTI reader for large file size -- [ ] Improve NIfTI writer to export as NIfTI-2 format -- [ ] Improve DICOM frame parser -- [ ] Support for additional medical image formats -- [ ] Support changing value in NIfTI frame +- [ ] Improve DICOM PixelData frame parser ## Example To parse a DICOM file @@ -15,11 +11,6 @@ To parse a DICOM file ``` -To parse a single NIfTI file: -```go - -``` - ## Supported Transfer Syntaxes ```text ImplicitVRLittleEndian "1.2.840.10008.1.2" diff --git a/dataset.go b/dataset.go new file mode 100644 index 00000000..4cdcb3ef --- /dev/null +++ b/dataset.go @@ -0,0 +1,614 @@ +package go2com + +import ( + "bufio" + "bytes" + "encoding/binary" + "errors" + "fmt" + "github.com/okieraised/go2com/internal/system" + "github.com/okieraised/go2com/internal/utils" + "github.com/okieraised/go2com/pkg/dicom/tag" + "github.com/okieraised/go2com/pkg/dicom/vr" + "io" + "reflect" + "strconv" + "strings" +) + +const ( + VLUndefinedLength uint32 = 0xFFFFFFFF +) + +type Value struct { + RawValue interface{} `json:"raw_value"` +} + +// Element defines the struct for each dicom tag element info. Ordered as below to decrease the memory footprint +type Element struct { + Value Value `json:"value"` + TagName string `json:"tag_name"` + ValueRepresentationStr string `json:"value_representation_str"` + ValueLength uint32 `json:"value_length"` + Tag tag.DicomTag `json:"tag"` + ValueRepresentation vr.VRKind `json:"value_representation"` +} + +type Dataset struct { + Elements []*Element `json:"elements"` +} + +type DicomUID struct { + StudyInstanceUID string `json:"study_instance_uid"` + SeriesInstanceUID string `json:"series_instance_uid"` + SOPInstanceUID string `json:"sop_instance_uid"` +} + +func (ds *Dataset) RetrieveFileUID() (*DicomUID, error) { + res := DicomUID{} + for _, elem := range ds.Elements { + if elem.Tag == tag.SOPInstanceUID { + res.SOPInstanceUID = (elem.Value.RawValue).(string) + } + if elem.Tag == tag.SeriesInstanceUID { + res.SeriesInstanceUID = (elem.Value.RawValue).(string) + } + if elem.Tag == tag.StudyInstanceUID { + res.StudyInstanceUID = (elem.Value.RawValue).(string) + } + } + + if res.StudyInstanceUID == "" || res.SeriesInstanceUID == "" || res.SOPInstanceUID == "" { + return nil, errors.New("missing required UID to identify instance") + } + + return &res, nil +} + +// FindElementByTagStr returns the corresponding element of the input tag. +// Tag must be in 'ggggeeee' or '(gggg,eeee)' format +func (ds *Dataset) FindElementByTagStr(tagStr string) (*Element, error) { + tagStr = utils.FormatTag(tagStr) + for _, elem := range ds.Elements { + if elem.Tag.StringWithoutParentheses() == tagStr { + return elem, nil + } + } + return nil, fmt.Errorf("cannot find tag %s", tagStr) +} + +// FindElementByTagName returns the corresponding element of the input tag name. +func (ds *Dataset) FindElementByTagName(tagName string) (*Element, error) { + tagName = utils.FormatTagName(tagName) + for _, elem := range ds.Elements { + if strings.ToLower(elem.TagName) == tagName { + return elem, nil + } + } + return nil, fmt.Errorf("cannot find tag %s", tagName) +} + +// FindElementByTag returns the corresponding element of the input tag name. +func (ds *Dataset) FindElementByTag(tagName tag.DicomTag) (*Element, error) { + for _, elem := range ds.Elements { + if tagName == elem.Tag { + return elem, nil + } + } + return nil, fmt.Errorf("cannot find tag %s", tagName) +} + +// ReadElement reads the DICOM file tag by tag and returns the pointer to the parsed Element +func ReadElement(r *dcmReader, isImplicit bool, binOrder binary.ByteOrder) (*Element, error) { + tagVal, dcmTagInfo, err := readTag(r) + if err != nil { + return nil, err + } + + if *tagVal == tag.ItemDelimitationItem || *tagVal == tag.Item { + _ = r.skip(4) + return nil, nil + } + + if *tagVal == tag.PixelData && r.SkipPixelData() { + _, err = r.discard(int(r.GetFileSize())) + if err != nil { + return nil, err + } + return nil, nil + } + dmcTagName := dcmTagInfo.Name + dcmVR, err := readVR(r, isImplicit, *tagVal) + if err != nil { + return nil, err + } + + dcmVL, err := readVL(r, isImplicit, *tagVal, dcmVR) + if err != nil { + return nil, err + } + + value, err := readValue(r, *tagVal, dcmVR, dcmVL) + if err != nil { + return nil, err + } + if n, ok := value.([]byte); ok { + dcmVL = uint32(len(n)) + } + + elem := Element{ + Tag: *tagVal, + TagName: dmcTagName, + ValueRepresentationStr: dcmVR, + ValueLength: dcmVL, + Value: Value{RawValue: value}, + } + + return &elem, nil +} + +// readTag returns the tag information +func readTag(r *dcmReader) (*tag.DicomTag, *tag.TagInfo, error) { + group, err := r.readUInt16() + if err != nil { + return nil, nil, err + } + element, err := r.readUInt16() + if err != nil { + return nil, nil, err + } + + t := tag.DicomTag{ + Group: group, + Element: element, + } + + // Check if tag is private. If yes, just return here + // Otherwise, find info about the public tag + if int(group)%2 != 0 { + tagInfo := tag.TagInfo{ + VR: "", + Name: PrivateTag, + VM: "", + Status: "", + } + return &t, &tagInfo, nil + } + + tagInfo, err := tag.Find(t) + if err != nil { + return nil, nil, err + } + return &t, &tagInfo, nil +} + +// readVR returns the value representation of the tag +func readVR(r *dcmReader, isImplicit bool, t tag.DicomTag) (string, error) { + if isImplicit { + record, err := tag.Find(t) + if err != nil { + return vr.Unknown, nil + } + return record.VR, nil + } + return r.readString(2) +} + +// readVL returns the value length of the dicom tag +func readVL(r *dcmReader, isImplicit bool, t tag.DicomTag, valueRepresentation string) (uint32, error) { + if isImplicit { + return r.readUInt32() + } + + switch valueRepresentation { + // if the VR is equal to ‘OB’,’OW’,’OF’,’SQ’,’UI’ or ’UN’, + // the VR is having an extra 2 bytes trailing to it. These 2 bytes trailing to VR are empty and are not decoded. + // When VR is having these 2 extra empty bytes the VL will occupy 4 bytes rather than 2 bytes + case vr.OtherByte, vr.OtherWord, vr.OtherFloat, vr.SequenceOfItems, vr.Unknown, vr.OtherByteOrOtherWord, + strings.ToLower(vr.OtherByteOrOtherWord), vr.UnlimitedText, vr.UniversalResourceIdentifier, + vr.UnlimitedCharacters: + r.skip(2) + valueLength, err := r.readUInt32() + if err != nil { + return 0, err + } + if valueLength == VLUndefinedLength && + (valueRepresentation == vr.UnlimitedCharacters || + valueRepresentation == vr.UniversalResourceIdentifier || + valueRepresentation == vr.UnlimitedText) { + return 0, errors.New("UC, UR and UT must have defined length") + } + return valueLength, nil + default: + valueLength, err := r.readUInt16() + if err != nil { + return 0, err + } + vl := uint32(valueLength) + if vl == 0xffff { + vl = VLUndefinedLength + } + return vl, nil + } +} + +// readValue returns the value of the dicom tag +func readValue(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + // Add this here for a special case when the tag is private and the value representation is UN (unknown) but the + // value is of sequence of items. In this case, we will peek the next 4 bytes and check if it matches the item tag + // If yes then handles like SQ + if valueRepresentation == vr.Unknown && t.Group%2 != 0 { + n, err := r.peek(4) + if err != nil { + return nil, err + } + if binary.BigEndian.Uint32(n) == 0xFFFEE000 || binary.BigEndian.Uint32(n) == 0xFEFF00E0 || binary.BigEndian.Uint32(n) == VLUndefinedLength { + r.SetTransferSyntax(r.ByteOrder(), true) + return readSequence(r, t, valueRepresentation, valueLength) + } + } + vrKind := vr.GetVR(t, valueRepresentation) + switch vrKind { + case vr.VRString, vr.VRDate: + return readStringType(r, t, valueRepresentation, valueLength) + case vr.VRInt16, vr.VRInt32, vr.VRUInt16, vr.VRUInt32, vr.VRTagList: + return readIntType(r, t, valueRepresentation, valueLength) + case vr.VRFloat32, vr.VRFloat64: + return readFloatType(r, t, valueRepresentation, valueLength) + case vr.VRBytes: + return readByteType(r, t, valueRepresentation, valueLength) + case vr.VRPixelData: + return readPixelDataType(r, t, valueRepresentation, valueLength) + case vr.VRSequence: + return readSequence(r, t, valueRepresentation, valueLength) + default: + return readStringType(r, t, valueRepresentation, valueLength) + } +} + +// switchStringToNumeric convert the decimal string to its appropriate value type +func switchStringToNumeric(in interface{}, valueRepresentation string) interface{} { + switch valueRepresentation { + case vr.IntegerString: + switch reflect.ValueOf(in).Kind() { + case reflect.Slice: + ValStrArr, ok := (in).([]string) + if !ok { + return in + } + res := make([]int, 0, len(ValStrArr)) + for _, sub := range ValStrArr { + intVar, err := strconv.Atoi(sub) + if err != nil { + return in + } + res = append(res, intVar) + } + return res + case reflect.String: + valStr, ok := (in).(string) + if !ok { + return in + } + intVal, err := strconv.Atoi(valStr) + if err != nil { + return in + } + return intVal + } + + case vr.DecimalString, vr.OtherFloat, vr.OtherDouble: + switch reflect.ValueOf(in).Kind() { + case reflect.Slice: + ValStrArr, ok := (in).([]string) + if !ok { + return in + } + res := make([]float64, 0, len(ValStrArr)) + for _, sub := range ValStrArr { + flVar, err := strconv.ParseFloat(sub, 64) + if err != nil { + return in + } + res = append(res, flVar) + } + return res + case reflect.String: + valStr, ok := (in).(string) + if !ok { + return in + } + flVar, err := strconv.ParseFloat(valStr, 64) + if err != nil { + return in + } + return flVar + } + default: + } + return in +} + +// readStringType reads the value as string and strips any zero padding +func readStringType(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + sep := "\\" + str, err := r.readString(valueLength) + if err != nil { + return str, err + } + str = strings.Trim(str, " \000") // There is a space " \000", not "\000" + if strings.Contains(str, sep) { + strArr := strings.Split(str, sep) + res := switchStringToNumeric(strArr, valueRepresentation) + return res, nil + } + res := switchStringToNumeric(str, valueRepresentation) + return res, nil +} + +// readPixelDataType reads the raw pixel data +func readPixelDataType(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + if valueLength%2 != 0 && valueLength != VLUndefinedLength { + fmt.Printf("Odd value length encountered for tag: %v with length %d", t.String(), valueLength) + } + byteSize := r.GetFileSize() + if valueLength != VLUndefinedLength { + byteSize = int64(valueLength) + } + + bArr := make([]byte, byteSize) + n, err := io.ReadFull(r, bArr) + sbArr := bArr[:n] + if err != nil { + if err == io.ErrUnexpectedEOF { + return sbArr, nil + } + return nil, err + } + return bArr, nil +} + +// readByteType reads the value as byte array +func readByteType(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + switch valueRepresentation { + case vr.OtherByte, vr.Unknown, vr.OtherByteOrOtherWord, strings.ToLower(vr.OtherByteOrOtherWord): + bArr := make([]byte, valueLength) + n, err := io.ReadFull(r, bArr) + sbArr := bArr[:n] + if err != nil { + if err == io.ErrUnexpectedEOF { + return sbArr, nil + } + return nil, err + } + return bArr, nil + case vr.OtherWord: + if valueLength%2 != 0 { + fmt.Printf("Odd value length encountered for tag: %v with length %d", t.String(), valueLength) + } + buf := bytes.NewBuffer(make([]byte, 0, valueLength)) + numWords := int(valueLength / 2) + for i := 0; i < numWords; i++ { + word, err := r.readUInt16() + if err != nil { + // Handle a case when the actual pixel data is less than the value length. Just return what we can + // read here + if err == io.EOF { + err = binary.Write(buf, system.NativeEndian, word) + if err != nil { + return nil, err + } + r = nil + return buf.Bytes(), nil + } + return nil, err + } + err = binary.Write(buf, system.NativeEndian, word) + if err != nil { + return nil, err + } + + } + return buf.Bytes(), nil + default: + _, err := r.discard(int(valueLength)) + if err != nil { + return nil, err + } + } + return nil, nil +} + +// readIntType reads the value as integer and returns either the value or a slice of value +func readIntType(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + var subVal int + retVal := make([]int, 0, valueLength/2) + n, err := r.peek(int(valueLength)) + if err != nil { + return nil, err + } + subReader := bytes.NewReader(n) + subRd := NewDICOMReader(bufio.NewReader(subReader), WithSkipPixelData(r.SkipPixelData())) + byteRead := 0 + for { + if byteRead >= int(valueLength) { + break + } + switch valueRepresentation { + case vr.UnsignedShort, vr.SignedShortOrUnsignedShort, strings.ToLower(vr.SignedShortOrUnsignedShort): + val, err := subRd.readUInt16() + if err != nil { + return nil, err + } + subVal = int(val) + byteRead += 2 + case vr.AttributeTag: + val, err := subRd.readUInt16() + if err != nil { + return nil, err + } + subVal = int(val) + byteRead += 2 + case vr.UnsignedLong: + val, err := subRd.readUInt32() + if err != nil { + return nil, err + } + subVal = int(val) + byteRead += 4 + case vr.SignedLong: + val, err := subRd.readInt32() + if err != nil { + return nil, err + } + subVal = int(val) + byteRead += 4 + case vr.SignedShort: + val, err := subRd.readInt16() + if err != nil { + return nil, err + } + subVal = int(val) + byteRead += 2 + } + retVal = append(retVal, subVal) + } + _, _ = r.discard(int(valueLength)) + if len(retVal) == 1 { + return retVal[0], nil + } + return retVal, nil +} + +// readFloatType reads the value as float +func readFloatType(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + var subVal float64 + retVal := make([]float64, 0, valueLength/2) + n, err := r.peek(int(valueLength)) + if err != nil { + return nil, err + } + subReader := bytes.NewReader(n) + subRd := NewDICOMReader(bufio.NewReader(subReader), WithSkipPixelData(r.SkipPixelData())) + byteRead := 0 + for { + if byteRead >= int(valueLength) { + break + } + switch valueRepresentation { + case vr.FloatingPointSingle, vr.OtherFloat: + val, err := subRd.readFloat32() + if err != nil { + return nil, err + } + subVal = float64(val) + byteRead += 4 + case vr.FloatingPointDouble: + val, err := subRd.readFloat64() + if err != nil { + return nil, err + } + subVal = val + byteRead += 8 + } + retVal = append(retVal, subVal) + } + _, _ = r.discard(int(valueLength)) + if len(retVal) == 1 { + return retVal[0], nil + } + + return retVal, nil +} + +// readSequence reads the value as sequence of items +func readSequence(r *dcmReader, t tag.DicomTag, valueRepresentation string, valueLength uint32) (interface{}, error) { + var sequences []*Element + // Reference: https://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html + if valueLength == VLUndefinedLength { + for { + subElement, err := ReadElement(r, r.IsImplicit(), r.ByteOrder()) + if err != nil { + return nil, err + } + + if subElement == nil { + continue + } + + if subElement.Tag == tag.SequenceDelimitationItem { + break + } + sequences = append(sequences, subElement) + } + if valueRepresentation == vr.Unknown { + r.SetTransferSyntax(r.ByteOrder(), r.isTrackingImplicit()) + } + } else { + n, err := r.peek(int(valueLength)) + if err != nil { + if err == bufio.ErrBufferFull { + bRaw, err := writeToBuf(r, int(valueLength)) + if err != nil { + return nil, err + } + sequences, err = readDefinedLengthSequences(r, bRaw, valueRepresentation) + if err != nil { + return nil, err + } + return sequences, nil + } + return nil, err + } + sequences, err = readDefinedLengthSequences(r, n, valueRepresentation) + if err != nil { + return nil, err + } + _, _ = r.discard(int(valueLength)) + } + return sequences, nil + +} + +func writeToBuf(r *dcmReader, n int) ([]byte, error) { + buf := bytes.NewBuffer(make([]byte, 0, n)) + + for i := 0; i < n; i++ { + word, err := r.readUInt8() + if err != nil { + return nil, err + } + err = binary.Write(buf, system.NativeEndian, word) + if err != nil { + return nil, err + } + } + return buf.Bytes(), nil +} + +func readDefinedLengthSequences(r *dcmReader, b []byte, valueRepresentation string) ([]*Element, error) { + var sequences []*Element + br := bytes.NewReader(b) + subRd := NewDICOMReader(bufio.NewReaderSize(br, len(b)), WithSkipPixelData(r.SkipPixelData())) + _ = subRd.skip(8) + subRd.SetTransferSyntax(r.ByteOrder(), r.IsImplicit()) + for { + subElement, err := ReadElement(subRd, r.IsImplicit(), r.ByteOrder()) + if err != nil { + if err == io.EOF { + break + } else { + return nil, err + } + } + if subElement == nil { + continue + } + sequences = append(sequences, subElement) + } + if valueRepresentation == vr.Unknown { + r.SetTransferSyntax(r.ByteOrder(), r.isTrackingImplicit()) + } + + return sequences, nil +} diff --git a/dicom_reader_test.go b/dicom_reader_test.go deleted file mode 100644 index 80405f5d..00000000 --- a/dicom_reader_test.go +++ /dev/null @@ -1,150 +0,0 @@ -package go2com - -import ( - "bufio" - "fmt" - "github.com/okieraised/go2com/internal/utils" - "github.com/okieraised/go2com/pkg/dicom/dcm_io" - "github.com/stretchr/testify/assert" - _ "image/jpeg" - "os" - "strings" - "testing" -) - -func Test_ReaderProfiling(t *testing.T) { - assert := assert.New(t) - f, err := os.Open("./test_data/01.dcm") - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - fn := func() { - dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) - err = dcmReader.Parse() - assert.NoError(err) - } - - err = utils.CPUProfilingFunc(fn, "./cpu1.pprof") - assert.NoError(err) -} - -func Test_NewParser1(t *testing.T) { - assert := assert.New(t) - filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_data") - assert.NoError(err) - for _, fPath := range filePaths { - fmt.Println("process:", fPath) - - f, err := os.Open(fPath) - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) - - err = dcmReader.Parse() - if err != nil { - fmt.Println(err) - return - } - assert.NoError(err) - } -} - -func Test_NewParser2(t *testing.T) { - assert := assert.New(t) - filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_full") - assert.NoError(err) - for _, fPath := range filePaths { - fmt.Println("process:", fPath) - f, err := os.Open(fPath) - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) - - err = dcmReader.Parse() - if err != nil { - if strings.Contains(err.Error(), "not in valid dicom format") { - continue - } - } - assert.NoError(err) - } -} - -func Test_NewParser3(t *testing.T) { - assert := assert.New(t) - filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/mammo_dicoms") - assert.NoError(err) - for _, fPath := range filePaths { - fmt.Println("process:", fPath) - - f, err := os.Open(fPath) - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) - - err = dcmReader.Parse() - assert.NoError(err) - - _ = dcmReader.ExportDatasetTags(false) - } -} - -func Test_NewParser4(t *testing.T) { - assert := assert.New(t) - fPath := "/home/tripg/workspace/dicom/mammo_dicoms/1.2.840.113619.2.255.10452022879169.3670200508103440.2701.dicom" - f, err := os.Open(fPath) - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(true), dcm_io.WithSkipDataset(true), dcm_io.WithSetFileSize(fInfo.Size())) - err = parser.Parse() - assert.NoError(err) - - for _, elem := range parser.GetDataset().Elements { - fmt.Println(elem) - } -} - -func Test_NewParser5(t *testing.T) { - assert := assert.New(t) - fPath := "/home/tripg/workspace/10142022/ALI_Technologies/UltraPACS/studies/w0019837/view0001" - f, err := os.Open(fPath) - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) - err = parser.Parse() - assert.NoError(err) - -} - -func Test_NewParser6(t *testing.T) { - assert := assert.New(t) - f, err := os.Open("/home/tripg/workspace/10142022/ALI_Technologies/UltraPACS/studies/w0055053/view0013") - //f, err = os.Open("/home/tripg/workspace/10142022/Acuson/Sequoia/EXAMS/EXAM0000/CLIPS/CLIP0031") - f, err = os.Open("/home/tripg/workspace/10142022/Hamamatsu/Dog_15x15_20x.dcm") - //f, err = os.Open("/home/tripg/workspace/dicom2/PrivateGEImplicitVRBigEndianTransferSyntax16Bits.dcm") - assert.NoError(err) - - fInfo, err := f.Stat() - assert.NoError(err) - - parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) - err = parser.Parse() - assert.NoError(err) -} diff --git a/dicom_reader_test.go.bak b/dicom_reader_test.go.bak new file mode 100644 index 00000000..b3462996 --- /dev/null +++ b/dicom_reader_test.go.bak @@ -0,0 +1,489 @@ +package go2com + +import ( + "bufio" + "fmt" + "github.com/okieraised/go2com/internal/utils" + "github.com/okieraised/go2com/pkg/dicom/dcm_io" + "github.com/stretchr/testify/assert" + _ "image/jpeg" + "os" + "strings" + "testing" +) + +func Test_ReaderProfiling(t *testing.T) { + assert := assert.New(t) + f, err := os.Open("./dicom_test/01.dcm") + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + fn := func() { + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) + err = dcmReader.Parse() + assert.NoError(err) + } + + err = utils.CPUProfilingFunc(fn, "./cpu1.pprof") + assert.NoError(err) +} + +func Test_NewParser_Chest(t *testing.T) { + + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/Downloads/1.3.46.670589.30.1.3.1.1625132584.1673448683234.1.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size()), dcm_io.WithSkipPixelData(true)) + + err = dcmReader.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + + for _, elem := range dcmReader.GetDataset().Elements { + if elem.Tag.String() == "(2627,2730)" || elem.Tag.String() == "(200b,0010)" || elem.Tag.String() == "(6000,3000)" { + continue + } + //if elem.Tag.String() == "(2001,9000)" { + // subElem := elem.Value.RawValue.([]*dcm_io.Element) + // for _, el := range subElem { + // fmt.Println(el) + // } + //} + fmt.Println(elem) + } + //fmt.Println(dcmReader.ExportSeriesTags()) + + } +} + +func Test_NewParser_PDF(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/Downloads/IMG-0003-00001 (1).dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + + //for _, elem := range dcmReader.GetDataset().Elements { + // fmt.Println(elem) + //} + fmt.Println(dcmReader.ExportSeriesTags()) + + } +} + +func Test_NewParser1(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_data") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + } +} + +func Test_NewParser1_1(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_data") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + } +} + +func Test_NewParser2_1(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_full") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + if strings.Contains(err.Error(), "not in valid dicom format") { + continue + } + } + assert.NoError(err) + } +} + +func Test_NewParser2(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_full") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + if err != nil { + if strings.Contains(err.Error(), "not in valid dicom format") { + continue + } + } + assert.NoError(err) + } +} + +func Test_NewParser3(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/mammo_dicoms") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + _ = dcmReader.ExportDatasetTags(false) + } +} + +func Test_NewParser4(t *testing.T) { + assert := assert.New(t) + fPath := "/home/tripg/workspace/dicom/mammo_dicoms/1.2.840.113619.2.255.10452022879169.3670200508103440.2701.dicom" + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(true), dcm_io.WithSkipDataset(true), dcm_io.WithSetFileSize(fInfo.Size())) + err = parser.Parse() + assert.NoError(err) + + for _, elem := range parser.GetDataset().Elements { + fmt.Println(elem) + } +} + +func Test_NewParser5(t *testing.T) { + assert := assert.New(t) + fPath := "/home/tripg/workspace/10142022/ALI_Technologies/UltraPACS/studies/w0019837/view0001" + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) + err = parser.Parse() + assert.NoError(err) + +} + +func Test_NewParser6(t *testing.T) { + assert := assert.New(t) + f, err := os.Open("/home/tripg/workspace/10142022/ALI_Technologies/UltraPACS/studies/w0055053/view0013") + //f, err = os.Open("/home/tripg/workspace/10142022/Acuson/Sequoia/EXAMS/EXAM0000/CLIPS/CLIP0031") + f, err = os.Open("/home/tripg/workspace/10142022/Hamamatsu/Dog_15x15_20x.dcm") + //f, err = os.Open("/home/tripg/workspace/dicom2/PrivateGEImplicitVRBigEndianTransferSyntax16Bits.dcm") + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + parser := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSkipPixelData(false), dcm_io.WithSkipDataset(false), dcm_io.WithSetFileSize(fInfo.Size())) + err = parser.Parse() + assert.NoError(err) +} + +func Test_NewParser7(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + _ = dcmReader.ExportDatasetTags(false) + } +} + +func Test_NewParser7_Single_1(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm/ExplVR_BigEndNoMeta.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + for _, elem := range dcmReader.GetDataset().Elements { + fmt.Println(elem) + } + } +} + +func Test_NewParser7_Single_2(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm/ExplVR_LitEndNoMeta.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + for _, elem := range dcmReader.GetDataset().Elements { + fmt.Println(elem) + } + } +} + +func Test_NewParser7_Single_3(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm/SC_rgb_jpeg.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + for _, elem := range dcmReader.GetDataset().Elements { + fmt.Println(elem) + } + } +} + +func Test_NewParser7_Single_4(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm/no_meta_group_length.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + for _, elem := range dcmReader.GetMetadata().Elements { + fmt.Println(elem) + } + + for _, elem := range dcmReader.GetDataset().Elements { + fmt.Println(elem) + } + } +} + +func Test_NewParser8(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom3/IMAGES/") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + //_ = dcmReader.ExportDatasetTags(false) + //for _, elem := range dcmReader.GetDataset().Elements { + // fmt.Println(elem) + //} + } +} + +func Test_NewParser9(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom3/IMAGES/") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + //_ = dcmReader.ExportDatasetTags(false) + //for _, elem := range dcmReader.GetDataset().Elements { + // fmt.Println(elem) + //} + } +} + +func Test_NewParser10(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_data/File 12943.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + //_ = dcmReader.ExportDatasetTags(false) + //for _, elem := range dcmReader.GetDataset().Elements { + // fmt.Println(elem) + //} + } +} + +func Test_NewParser11(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm/waveform_ecg.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + dcmReader := dcm_io.NewDICOMReader(bufio.NewReader(f), dcm_io.WithSetFileSize(fInfo.Size())) + + err = dcmReader.Parse() + assert.NoError(err) + + for _, elem := range dcmReader.GetDataset().Elements { + fmt.Println(elem) + } + } +} diff --git a/internal/matrix/matrix.go b/internal/matrix/matrix.go deleted file mode 100644 index fd490758..00000000 --- a/internal/matrix/matrix.go +++ /dev/null @@ -1,323 +0,0 @@ -package matrix - -import "C" -import "math" - -type FMat44 struct { - M [4][4]float32 -} - -type FMat33 struct { - M [3][3]float32 -} - -type DMat44 struct { - M [4][4]float64 -} - -type DMat33 struct { - M [3][3]float64 -} - -// Mat33Polar finds the closest orthogonal matrix to input A -// (in both the Frobenius and L2 norms). -// Algorithm is that from NJ Higham, SIAM J Sci Stat Comput, 7:1160-1174. -func Mat33Polar(A DMat33) DMat33 { - var X, Y, Z DMat33 - var alp, bet, gam, gmi float64 - var dif float64 = 1.0 - var k int64 = 0 - - X = A - - gam = Mat33Determinant(X) - for { - if gam > 0.0 { - break - } - gam = 0.00001 * (0.001 + Mat33RowNorm(X)) - X.M[0][0] += gam - X.M[1][1] += gam - X.M[2][2] += gam - gam = Mat33Determinant(X) - } - - for { - Y = Mat33Inverse(X) - if dif > 0.3 { /* far from convergence */ - alp = math.Sqrt(Mat33RowNorm(X) * Mat33ColNorm(X)) - bet = math.Sqrt(Mat33RowNorm(Y) * Mat33ColNorm(Y)) - gam = math.Sqrt(bet / alp) - gmi = 1.0 / gam - } else { - gam = 1.0 - gmi = 1.0 - } - Z.M[0][0] = 0.5 * (gam*X.M[0][0] + gmi*Y.M[0][0]) - Z.M[0][1] = 0.5 * (gam*X.M[0][1] + gmi*Y.M[1][0]) - Z.M[0][2] = 0.5 * (gam*X.M[0][2] + gmi*Y.M[2][0]) - Z.M[1][0] = 0.5 * (gam*X.M[1][0] + gmi*Y.M[0][1]) - Z.M[1][1] = 0.5 * (gam*X.M[1][1] + gmi*Y.M[1][1]) - Z.M[1][2] = 0.5 * (gam*X.M[1][2] + gmi*Y.M[2][1]) - Z.M[2][0] = 0.5 * (gam*X.M[2][0] + gmi*Y.M[0][2]) - Z.M[2][1] = 0.5 * (gam*X.M[2][1] + gmi*Y.M[1][2]) - Z.M[2][2] = 0.5 * (gam*X.M[2][2] + gmi*Y.M[2][2]) - - dif = math.Abs(Z.M[0][0]-X.M[0][0]) + math.Abs(Z.M[0][1]-X.M[0][1]) + math.Abs(Z.M[0][2]-X.M[0][2]) + - math.Abs(Z.M[1][0]-X.M[1][0]) + math.Abs(Z.M[1][1]-X.M[1][1]) + math.Abs(Z.M[1][2]-X.M[1][2]) + - math.Abs(Z.M[2][0]-X.M[2][0]) + math.Abs(Z.M[2][1]-X.M[2][1]) + math.Abs(Z.M[2][2]-X.M[2][2]) - - k = k + 1 - if k > 100 || dif < 3.e-6 { // convergence or exhaustion - break - } - X = Z - } - return Z -} - -// Mat33Determinant computes the determinant of a 3x3 matrix -func Mat33Determinant(R DMat33) float64 { - var r11, r12, r13, r21, r22, r23, r31, r32, r33 float64 - - // [ r11 r12 r13 ] - // [ r21 r22 r23 ] - // [ r31 r32 r33 ] - - r11 = R.M[0][0] - r12 = R.M[0][1] - r13 = R.M[0][2] - r21 = R.M[1][0] - r22 = R.M[1][1] - r23 = R.M[1][2] - r31 = R.M[2][0] - r32 = R.M[2][1] - r33 = R.M[2][2] - - return r11*r22*r33 - r11*r32*r23 - r21*r12*r33 + r21*r32*r13 + r31*r12*r23 - r31*r22*r13 -} - -// Mat33RowNorm computes the max row norm of a 3x3 matrix -func Mat33RowNorm(A DMat33) float64 { - var r1, r2, r3 float64 - - r1 = math.Abs(A.M[0][0]) + math.Abs(A.M[0][1]) + math.Abs(A.M[0][2]) - r2 = math.Abs(A.M[1][0]) + math.Abs(A.M[1][1]) + math.Abs(A.M[1][2]) - r3 = math.Abs(A.M[2][0]) + math.Abs(A.M[2][1]) + math.Abs(A.M[2][2]) - if r1 < r2 { - r1 = r2 - } - if r1 < r3 { - r1 = r3 - } - return r1 -} - -// Mat33ColNorm computes the max column norm of a 3x3 matrix -func Mat33ColNorm(A DMat33) float64 { - var r1, r2, r3 float64 - - r1 = math.Abs(A.M[0][0]) + math.Abs(A.M[1][0]) + math.Abs(A.M[2][0]) - r2 = math.Abs(A.M[0][1]) + math.Abs(A.M[1][1]) + math.Abs(A.M[2][1]) - r3 = math.Abs(A.M[0][2]) + math.Abs(A.M[1][2]) + math.Abs(A.M[2][2]) - - if r1 < r2 { - r1 = r2 - } - if r1 < r3 { - r1 = r3 - } - return r1 -} - -// MatMultiply multiples 2 3x3 matrices -func MatMultiply(A, B DMat33) DMat33 { - var Cm DMat33 - var i, j int64 - - for i = 0; i < 3; i++ { - for j = 0; j < 3; j++ { - Cm.M[i][j] = A.M[i][0]*B.M[0][j] + A.M[i][1]*B.M[1][j] + A.M[i][2]*B.M[2][j] - } - } - return Cm -} - -// Mat33Inverse computes the inverse of a bordered 3x3 matrix -func Mat33Inverse(R DMat33) DMat33 { - var r11, r12, r13, r21, r22, r23, r31, r32, r33, deti float64 - var Q DMat33 - - // INPUT MATRIX: - // [ r11 r12 r13 ] - // [ r21 r22 r23 ] - // [ r31 r32 r33 ] - - r11 = R.M[0][0] - r12 = R.M[0][1] - r13 = R.M[0][2] - r21 = R.M[1][0] - r22 = R.M[1][1] - r23 = R.M[1][2] - r31 = R.M[2][0] - r32 = R.M[2][1] - r33 = R.M[2][2] - - deti = r11*r22*r33 - r11*r32*r23 - r21*r12*r33 + r21*r32*r13 + r31*r12*r23 - r31*r22*r13 - - if deti != 0.0 { - deti = 1.0 / deti - } - - Q.M[0][0] = deti * (r22*r33 - r32*r23) - Q.M[0][1] = deti * (-r12*r33 + r32*r13) - Q.M[0][2] = deti * (r12*r23 - r22*r13) - - Q.M[1][0] = deti * (-r21*r33 + r31*r23) - Q.M[1][1] = deti * (r11*r33 - r31*r13) - Q.M[1][2] = deti * (-r11*r23 + r21*r13) - - Q.M[2][0] = deti * (r21*r32 - r31*r22) - Q.M[2][1] = deti * (-r11*r32 + r31*r12) - Q.M[2][2] = deti * (r11*r22 - r21*r12) - - return Q - -} - -// Mat44Inverse computes the inverse of a bordered 4x4 matrix -func Mat44Inverse(R DMat44) DMat44 { - var r11, r12, r13, r21, r22, r23, r31, r32, r33, v1, v2, v3, deti float64 - var Q DMat44 - - // INPUT MATRIX IS - // [ r11 r12 r13 v1 ] - // [ r21 r22 r23 v2 ] - // [ r31 r32 r33 v3 ] - // [ 0 0 0 1 ] - r11 = R.M[0][0] - r12 = R.M[0][1] - r13 = R.M[0][2] - r21 = R.M[1][0] - r22 = R.M[1][1] - r23 = R.M[1][2] - r31 = R.M[2][0] - r32 = R.M[2][1] - r33 = R.M[2][2] - v1 = R.M[0][3] - v2 = R.M[1][3] - v3 = R.M[2][3] - - deti = r11*r22*r33 - r11*r32*r23 - r21*r12*r33 + r21*r32*r13 + r31*r12*r23 - r31*r22*r13 // determinant - - if deti != 0.0 { - deti = 1.0 / deti - } - - Q.M[0][0] = deti * (r22*r33 - r32*r23) - Q.M[0][1] = deti * (-r12*r33 + r32*r13) - Q.M[0][2] = deti * (r12*r23 - r22*r13) - Q.M[0][3] = deti * (-r12*r23*v3 + r12*v2*r33 + r22*r13*v3 - r22*v1*r33 - r32*r13*v2 + r32*v1*r23) - - Q.M[1][0] = deti * (-r21*r33 + r31*r23) - Q.M[1][1] = deti * (r11*r33 - r31*r13) - Q.M[1][2] = deti * (-r11*r23 + r21*r13) - Q.M[1][3] = deti * (r11*r23*v3 - r11*v2*r33 - r21*r13*v3 + r21*v1*r33 + r31*r13*v2 - r31*v1*r23) - - Q.M[2][0] = deti * (r21*r32 - r31*r22) - Q.M[2][1] = deti * (-r11*r32 + r31*r12) - Q.M[2][2] = deti * (r11*r22 - r21*r12) - Q.M[2][3] = deti * (-r11*r22*v3 + r11*r32*v2 + r21*r12*v3 - r21*r32*v1 - r31*r12*v2 + r31*r22*v1) - - Q.M[3][0] = 0.0 - Q.M[3][1] = 0.0 - Q.M[3][2] = 0.0 - if deti == 0.0 { - Q.M[3][3] = 0.0 - } else { - Q.M[3][3] = 1.0 - } - - return Q - -} - -// MakeOrthoMat44 uses 9 input float64 and make an orthogonal mat44 out of them -func MakeOrthoMat44(r11, r12, r13, r21, r22, r23, r31, r32, r33 float64) DMat44 { - var R DMat44 - var Q, P DMat33 - var val float64 - - R.M[3][0] = 0.0 - R.M[3][1] = 0.0 - R.M[3][2] = 0.0 - R.M[3][3] = 1.0 - - Q.M[0][0] = r11 - Q.M[0][1] = r12 - Q.M[0][2] = r13 - Q.M[1][0] = r21 - Q.M[1][1] = r22 - Q.M[1][2] = r23 - Q.M[2][0] = r31 - Q.M[2][1] = r32 - Q.M[2][2] = r33 - - // normalize row 1 - val = Q.M[0][0]*Q.M[0][0] + Q.M[0][1]*Q.M[0][1] + Q.M[0][2]*Q.M[0][2] - if val > 0.0 { - val = 1.0 / math.Sqrt(val) - Q.M[0][0] *= val - Q.M[0][1] *= val - Q.M[0][2] *= val - } else { - Q.M[0][0] = 1.0 - Q.M[0][1] = 0.0 - Q.M[0][2] = 0.0 - } - - // normalize row 2 - val = Q.M[1][0]*Q.M[1][0] + Q.M[1][1]*Q.M[1][1] + Q.M[1][2]*Q.M[1][2] - if val > 0.0 { - val = 1.0 / math.Sqrt(val) - Q.M[1][0] *= val - Q.M[1][1] *= val - Q.M[1][2] *= val - } else { - Q.M[1][0] = 0.0 - Q.M[1][1] = 1.0 - Q.M[1][2] = 0.0 - } - - // normalize row 3 - val = Q.M[2][0]*Q.M[2][0] + Q.M[2][1]*Q.M[2][1] + Q.M[2][2]*Q.M[2][2] - if val > 0.0 { - val = 1.0 / math.Sqrt(val) - Q.M[2][0] *= val - Q.M[2][1] *= val - Q.M[2][2] *= val - } else { - Q.M[2][0] = Q.M[0][1]*Q.M[1][2] - Q.M[0][2]*Q.M[1][1] - Q.M[2][1] = Q.M[0][2]*Q.M[1][0] - Q.M[0][0]*Q.M[1][2] - Q.M[2][2] = Q.M[0][0]*Q.M[1][1] - Q.M[0][1]*Q.M[1][0] - } - - // P is orthogonal matrix closest to Q - P = Mat33Polar(Q) - - R.M[0][0] = P.M[0][0] - R.M[0][1] = P.M[0][1] - R.M[0][2] = P.M[0][2] - R.M[1][0] = P.M[1][0] - R.M[1][1] = P.M[1][1] - R.M[1][2] = P.M[1][2] - R.M[2][0] = P.M[2][0] - R.M[2][1] = P.M[2][1] - R.M[2][2] = P.M[2][2] - - R.M[0][3] = 0.0 - R.M[1][3] = 0.0 - R.M[2][3] = 0.0 - - return R -} diff --git a/io.go b/io.go new file mode 100644 index 00000000..c6db52bf --- /dev/null +++ b/io.go @@ -0,0 +1,240 @@ +package go2com + +import ( + "bufio" + "encoding/binary" + "fmt" + "github.com/okieraised/go2com/internal/utils" + "github.com/okieraised/go2com/pkg/dicom/tag" + "github.com/okieraised/go2com/pkg/dicom/vr" + "github.com/okieraised/go2com/pkg/plugins/orthanc" + "io" + "strings" +) + +const ( + MagicString = "DICM" + PrivateTag = "PrivateTag" +) + +type dcmReader struct { + reader *bufio.Reader + binaryOrder binary.ByteOrder + dataset Dataset + metadata Dataset + allowNonCompliantDcm bool + isImplicit bool + keepTrackImplicit bool + skipPixelData bool + skipDataset bool + fileSize int64 +} + +// NewDICOMReader returns a new reader +func NewDICOMReader(reader *bufio.Reader, options ...func(*dcmReader)) *dcmReader { + parser := &dcmReader{ + reader: reader, + binaryOrder: binary.LittleEndian, + isImplicit: false, + skipPixelData: false, + skipDataset: false, + } + for _, opt := range options { + opt(parser) + } + return parser +} + +// WithAllowNonCompliantDcm provides option to keep trying to parse the file even if it's not DICOM compliant +// e.g.: Missing header, missing FileMetaInformationGroupLength,... +func WithAllowNonCompliantDcm(allowNonCompliantDcm bool) func(*dcmReader) { + return func(s *dcmReader) { + s.allowNonCompliantDcm = allowNonCompliantDcm + } +} + +// WithSkipPixelData provides option to skip reading pixel data (7FE0,0010). +// If true, pixel data is skipped. If false, pixel data will be read +func WithSkipPixelData(skipPixelData bool) func(*dcmReader) { + return func(s *dcmReader) { + s.skipPixelData = skipPixelData + } +} + +// WithSkipDataset provides option to read only the metadata header. +// If true, only the meta header is read, else, the dataset will be read +func WithSkipDataset(skipDataset bool) func(*dcmReader) { + return func(s *dcmReader) { + s.skipDataset = skipDataset + } +} + +// WithSetFileSize provides option to set the file size to the reader +func WithSetFileSize(fileSize int64) func(*dcmReader) { + return func(s *dcmReader) { + s.fileSize = fileSize + } +} + +// GetMetadata returns the file meta header +func (r *dcmReader) GetMetadata() Dataset { + return r.metadata +} + +// GetDataset returns the dataset +func (r *dcmReader) GetDataset() Dataset { + return r.dataset +} + +func (r *dcmReader) RetrieveFileUID() (*DicomUID, error) { + return r.dataset.RetrieveFileUID() +} + +// Parse reads the DICOM file and parses it into array of elements +func (r *dcmReader) Parse() error { + err := r.parse() + if err != nil { + return err + } + return nil +} + +func (r *dcmReader) SkipPixelData() bool { + return r.skipPixelData +} + +func (r *dcmReader) ByteOrder() binary.ByteOrder { + return r.binaryOrder +} + +// IsValidDICOM checks if the dicom file follows the standard by having 128 bytes preamble followed by the magic string 'DICM' +func (r *dcmReader) IsValidDICOM() error { + preamble, err := r.peek(132) + if err != nil { + return fmt.Errorf("cannot read the first 132 bytes: %v", err) + } + if string(preamble[128:]) != MagicString { + return fmt.Errorf("file is not in valid dicom format") + } + return nil +} + +// GetElementByTagString returns the element value of the input tag +// Tag should be in (gggg,eeee) or ggggeeee format +func (r *dcmReader) GetElementByTagString(tagStr string) (interface{}, error) { + tagStr = utils.FormatTag(tagStr) + + if strings.HasPrefix(tagStr, "0002") { + for _, elem := range r.metadata.Elements { + if tagStr == elem.Tag.StringWithoutParentheses() { + return elem.Value, nil + } + } + return nil, fmt.Errorf("cannot find tag %s", tagStr) + } else { + for _, elem := range r.dataset.Elements { + if tagStr == elem.Tag.StringWithoutParentheses() { + return elem.Value, nil + } + } + return nil, fmt.Errorf("cannot find tag %s", tagStr) + } +} + +// ExportDatasetTags returns the mapped tag/(vr,value) dictionary +func (r *dcmReader) ExportDatasetTags(exportMeta bool) MappedTag { + res := make(MappedTag, len(r.metadata.Elements)+len(r.dataset.Elements)) + if exportMeta { + mt := r.metadata + for _, elem := range mt.Elements { + res.mapElement(elem) + } + } + + ds := r.dataset + //colorImage := false + for _, elem := range ds.Elements { + //if elem.Tag == tag.RedPaletteColorLookupTableData || elem.Tag == tag.BluePaletteColorLookupTableData || elem.Tag == tag.GreenPaletteColorLookupTableData { + // colorImage = true + //} + vrStr := elem.ValueRepresentationStr + if vrStr == "OB" || vrStr == "OW" || vrStr == "UN" || strings.ToLower(vrStr) == "ox" { + continue + } + res.mapElement(elem) + } + + //if colorImage { + // bulkURIMap := createOrthancURI(ds) + // for k, v := range bulkURIMap { + // res[k] = v + // } + //} + + return res +} + +func (r *dcmReader) ExportSeriesTags() MappedTag { + res := make(MappedTag, len(r.dataset.Elements)) + var value interface{} + for _, elem := range r.dataset.Elements { + if _, ok := orthanc.OrthancGetSeriesOfStudyTags[elem.Tag]; ok { + tagStr := elem.Tag.StringWithoutParentheses() + if elem.ValueRepresentationStr == vr.PersonName { + value = utils.AppendToSlice(map[string]interface{}{ + "Alphabetic": elem.Value.RawValue, + }) + } else { + value = utils.AppendToSlice(elem.Value.RawValue) + } + res[tagStr] = tag.TagBrowser{ + VR: elem.ValueRepresentationStr, + Value: value, + } + + } + } + + prefix := "http://127.0.0.1:8042" + firstSeriesURI := fmt.Sprintf("%s/dicom-web/studies/%s/series/%s", prefix, + res[tag.StudyInstanceUID.StringWithoutParentheses()].Value.([]interface{})[0].(string), + res[tag.SeriesInstanceUID.StringWithoutParentheses()].Value.([]interface{})[0].(string), + ) + + res[tag.RetrieveURL.StringWithoutParentheses()] = tag.TagBrowser{ + Value: []interface{}{firstSeriesURI}, + VR: "UR", + } + return res +} + +func (r *dcmReader) Read(p []byte) (int, error) { + return r.reader.Read(p) +} + +func (r *dcmReader) IsImplicit() bool { + return r.isImplicit +} + +func (r *dcmReader) SetTransferSyntax(binaryOrder binary.ByteOrder, isImplicit bool) { + r.binaryOrder = binaryOrder + r.isImplicit = isImplicit +} + +func (r *dcmReader) SetFileSize(fileSize int64) error { + r.fileSize = fileSize + return nil +} + +func (r *dcmReader) GetFileSize() int64 { + return r.fileSize +} + +func (r *dcmReader) readString(n uint32) (string, error) { + data := make([]byte, n) + _, err := io.ReadFull(r, data) + if err != nil { + return "", err + } + return string(data), nil +} diff --git a/io_benchmark_test.go b/io_benchmark_test.go new file mode 100644 index 00000000..3c436f74 --- /dev/null +++ b/io_benchmark_test.go @@ -0,0 +1 @@ +package go2com diff --git a/io_test.go b/io_test.go new file mode 100644 index 00000000..fb7984d4 --- /dev/null +++ b/io_test.go @@ -0,0 +1,220 @@ +package go2com + +import ( + "bufio" + "fmt" + "github.com/okieraised/go2com/internal/utils" + "github.com/stretchr/testify/assert" + "os" + "testing" +) + +func TestNewDICOMReader_1(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_data") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + } +} + +func TestNewDICOMReader_2(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/test_full") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + fmt.Println(err) + continue + } + assert.NoError(err) + } +} + +func TestNewDICOMReader_3(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/mammo_dicoms") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + } +} + +func TestNewDICOMReader_4(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/mammo_dicoms") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + if err != nil { + fmt.Println(err) + return + } + assert.NoError(err) + } +} + +func TestNewDICOMReader_5(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/color_dicom") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + assert.NoError(err) + + _ = rd.ExportDatasetTags(false) + } +} + +func TestNewDICOMReader_6(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/pydicom_dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + assert.NoError(err) + + _ = rd.ExportDatasetTags(false) + } +} + +func TestNewDICOMReader_7(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/dicom/dicoms_mr_func") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + assert.NoError(err) + + _ = rd.ExportDatasetTags(false) + } +} + +func TestNewDICOMReader_8(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/philips") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size())) + + err = rd.Parse() + assert.NoError(err) + + _ = rd.ExportDatasetTags(false) + } +} + +func TestNewDICOMReader_9(t *testing.T) { + assert := assert.New(t) + filePaths, err := utils.ReadDirRecursively("/home/tripg/workspace/img2.dcm") + assert.NoError(err) + for _, fPath := range filePaths { + fmt.Println("process:", fPath) + + f, err := os.Open(fPath) + assert.NoError(err) + + fInfo, err := f.Stat() + assert.NoError(err) + + rd := NewDICOMReader(bufio.NewReader(f), WithSetFileSize(fInfo.Size()), WithSkipPixelData(true)) + + err = rd.Parse() + assert.NoError(err) + + for _, elem := range rd.GetDataset().Elements { + fmt.Println(elem) + } + + _ = rd.ExportDatasetTags(false) + } +} diff --git a/main.py b/main.py index 4616a9e6..aff4b82d 100644 --- a/main.py +++ b/main.py @@ -53,6 +53,17 @@ fpath = '/home/tripg/Downloads/1-1.dcm' fpath = '/home/tripg/workspace/dicom3/winter.dcm' fpath = '/home/tripg/Downloads/PrivateGEImplicitVRBigEndianTransferSyntax16Bits.dcm' +fpath = '/home/tripg/Downloads/IMG-0003-00001 (1).dcm' +fpath = '/home/tripg/Downloads/1.3.46.670589.30.1.3.1.1625132584.1673448683234.1.dcm' +# fpath = '/home/tripg/workspace/pydicom_dcm/SC_rgb_jpeg.dcm' +# fpath = '/home/tripg/workspace/pydicom_dcm/image_dfl.dcm' +fpath = '/home/tripg/workspace/pydicom_dcm/SC_rgb_jpeg.dcm' +fpath = '/home/tripg/Downloads/123.241606668321866.1724728615648318_en.dcm' +fpath = '/home/tripg/workspace/dicom/test_data/File 12943.dcm' +fpath = '/home/tripg/workspace/pydicom_dcm/no_meta_group_length.dcm' +# fpath = '/home/tripg/workspace/pydicom_dcm/waveform_ecg.dcm' +# fpath = '/home/tripg/workspace/pydicom_dcm/meta_missing_tsyntax.dcm' +# fpath = '/home/tripg/workspace/pydicom_dcm/rtstruct.dcm' ds = dcmread(fpath, force=True) print(ds) diff --git a/nifti_sample/avg152T1_LR_nifti.nii b/nifti_sample/avg152T1_LR_nifti.nii deleted file mode 100644 index 342d3acc..00000000 Binary files a/nifti_sample/avg152T1_LR_nifti.nii and /dev/null differ diff --git a/nifti_sample/avg152T1_RL_nifti.nii b/nifti_sample/avg152T1_RL_nifti.nii deleted file mode 100644 index 724e3536..00000000 Binary files a/nifti_sample/avg152T1_RL_nifti.nii and /dev/null differ diff --git a/nifti_sample/zstat1.nii b/nifti_sample/zstat1.nii deleted file mode 100644 index 2aceb1ba..00000000 Binary files a/nifti_sample/zstat1.nii and /dev/null differ diff --git a/nifti_test.go b/nifti_test.go deleted file mode 100644 index e572519b..00000000 --- a/nifti_test.go +++ /dev/null @@ -1,177 +0,0 @@ -package go2com - -import ( - "fmt" - _ "image/jpeg" - "testing" - - "github.com/okieraised/go2com/pkg/nifti/nii_io" - "github.com/stretchr/testify/assert" -) - -func TestNii1(t *testing.T) { - assert := assert.New(t) - - filePath := "/home/tripg/Documents/nifti/Arnow^Corie^Shelvey^OM_segmented.nii" - filePath = "/home/tripg/Documents/nifti/RGB16_4D.nii.gz" - filePath = "/home/tripg/Documents/nifti/someones_anatomy.nii.gz" - filePath = "/home/tripg/Documents/nifti/someones_epi.nii.gz" - //filePath = "/home/tripg/Documents/nifti/RGB8_4D.nii.gz" - //filePath = "/home/tripg/Documents/nifti/jaw.nii.gz" - //filePath = "/home/tripg/Documents/nifti/Arnow^Corie^Shelvey^OM_segmented.nii" - //filePath = "/home/tripg/Documents/nifti/knee.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/fmri.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/structural.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/structural_brain.nii.gz" - //filePath = "/home/tripg/Documents/nifti/JHU_MNI_SS_T1.nii.gz" - //filePath = "/home/tripg/Documents/nifti/avg152T1_LR_nifti2.nii.gz" - //filePath = "/Users/TriPham/Documents/nifti/avg152T1_RL_nifti.nii" - - niiReader, err := nii_io.NewNiiReader(filePath, nii_io.WithInMemory(true)) - assert.NoError(err) - err = niiReader.Parse() - assert.NoError(err) - - shape := niiReader.GetImgShape() - fmt.Println("shape", shape) - fmt.Println(niiReader.GetUnitsOfMeasurements()) - fmt.Println("dtype", niiReader.GetDatatype()) - fmt.Println("nbyper", niiReader.GetNiiData().NByPer) - fmt.Println("GetSFromCode()", niiReader.GetSFormCode()) - fmt.Println("GetQFromCode()", niiReader.GetQFormCode()) - fmt.Println("orientation", niiReader.GetOrientation()) - fmt.Println("affine", niiReader.GetAffine()) - fmt.Println("QOffsetX", niiReader.GetNiiData().QoffsetX) - fmt.Println("QOffsetY", niiReader.GetNiiData().QoffsetY) - fmt.Println("QOffsetZ", niiReader.GetNiiData().QoffsetZ) - - fmt.Println("---------------------------------------------------------------------------------------------------") - - writer, err := nii_io.NewNiiWriter("./out.nii.gz", nii_io.WithNIfTIData(niiReader.GetNiiData()), nii_io.WithCompression(true)) - assert.NoError(err) - err = writer.WriteToFile() - assert.NoError(err) - - filePath = "./out.nii.gz" - niiReader, err = nii_io.NewNiiReader(filePath, nii_io.WithInMemory(true)) - assert.NoError(err) - err = niiReader.Parse() - assert.NoError(err) - - shape = niiReader.GetImgShape() - fmt.Println("shape", shape) - fmt.Println(niiReader.GetUnitsOfMeasurements()) - fmt.Println("dtype", niiReader.GetDatatype()) - fmt.Println("nbyper", niiReader.GetNiiData().NByPer) - fmt.Println("GetSFromCode()", niiReader.GetSFormCode()) - fmt.Println("GetQFromCode()", niiReader.GetQFormCode()) - fmt.Println("orientation", niiReader.GetOrientation()) - fmt.Println("affine", niiReader.GetAffine()) - fmt.Println("QOffsetX", niiReader.GetNiiData().QoffsetX) - fmt.Println("QOffsetY", niiReader.GetNiiData().QoffsetY) - fmt.Println("QOffsetZ", niiReader.GetNiiData().QoffsetZ) - - //res, err := niiReader.GetSlice(1, 0) - //assert.NoError(err) - // - //fmt.Println("len res", len(res)) - //fmt.Println("len res0", len(res[0])) - // - ////fmt.Println(res) - //res2, err := niiReader.GetVolume(0) - //assert.NoError(err) - // - //fmt.Println("len res2", len(res2)) - // - //writer := nii_io.NewNiiWriter("./out.nii", nii_io.WithNIfTIData(niiReader.GetNiiData())) - // - //err = writer.WriteToFile() - //assert.NoError(err) - - //for _, elem := range res { - // fmt.Println(elem) - // - //} - - //fmt.Println(niiReader.GetTimeSeries(26, 30, 16)) - - //for x := 0; x < 57; x++ { - // for y := 0; y < 67; y++ { - // for z := 0; z < 56; z++ { - // fmt.Println(niiReader.GetTimeSeries(int64(x), int64(y), int64(z))) - // } - // } - //} - - shape = niiReader.GetImgShape() - fmt.Println(shape) - -} - -func TestNii2(t *testing.T) { - assert := assert.New(t) - - filePath := "/home/tripg/Documents/nifti/Arnow^Corie^Shelvey^OM_segmented.nii" - filePath = "/home/tripg/Documents/nifti/RGB16_4D.nii.gz" - filePath = "/home/tripg/Documents/nifti/someones_anatomy.nii.gz" - filePath = "/home/tripg/Documents/nifti/someones_epi.nii.gz" - //filePath = "/home/tripg/Documents/nifti/RGB8_4D.nii.gz" - //filePath = "/home/tripg/Documents/nifti/jaw.nii.gz" - //filePath = "/home/tripg/Documents/nifti/Arnow^Corie^Shelvey^OM_segmented.nii" - //filePath = "/home/tripg/Documents/nifti/knee.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/fmri.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/structural.nii.gz" - //filePath = "/home/tripg/Documents/nifti/ExBox11/structural_brain.nii.gz" - //filePath = "/home/tripg/Documents/nifti/JHU_MNI_SS_T1.nii.gz" - //filePath = "/home/tripg/Documents/nifti/avg152T1_LR_nifti2.nii.gz" - filePath = "/home/tripg/Documents/nifti/avg152T1_LR_nifti.img.gz" - - niiReader, err := nii_io.NewNiiReader(filePath, nii_io.WithInMemory(true), nii_io.WithHeaderFile("/home/tripg/Documents/nifti/avg152T1_LR_nifti.hdr.gz")) - assert.NoError(err) - err = niiReader.Parse() - assert.NoError(err) - - shape := niiReader.GetImgShape() - fmt.Println("shape", shape) - fmt.Println(niiReader.GetUnitsOfMeasurements()) - fmt.Println("dtype", niiReader.GetDatatype()) - fmt.Println("nbyper", niiReader.GetNiiData().NByPer) - fmt.Println("GetSFromCode()", niiReader.GetSFormCode()) - fmt.Println("GetQFromCode()", niiReader.GetQFormCode()) - fmt.Println("orientation", niiReader.GetOrientation()) - fmt.Println("affine", niiReader.GetAffine()) - fmt.Println("QOffsetX", niiReader.GetNiiData().QoffsetX) - fmt.Println("QOffsetY", niiReader.GetNiiData().QoffsetY) - fmt.Println("QOffsetZ", niiReader.GetNiiData().QoffsetZ) - fmt.Println("VoxOffset", niiReader.GetNiiData().VoxOffset) - - fmt.Println("---------------------------------------------------------------------------------------------------") - - writer, err := nii_io.NewNiiWriter("./out2.nii.gz", nii_io.WithNIfTIData(niiReader.GetNiiData()), nii_io.WithCompression(true)) - assert.NoError(err) - err = writer.WriteToFile() - assert.NoError(err) - - filePath = "./out2.nii.gz" - niiReader, err = nii_io.NewNiiReader(filePath, nii_io.WithInMemory(true)) - assert.NoError(err) - err = niiReader.Parse() - assert.NoError(err) - - shape = niiReader.GetImgShape() - fmt.Println("shape", shape) - fmt.Println(niiReader.GetUnitsOfMeasurements()) - fmt.Println("dtype", niiReader.GetDatatype()) - fmt.Println("nbyper", niiReader.GetNiiData().NByPer) - fmt.Println("GetSFromCode()", niiReader.GetSFormCode()) - fmt.Println("GetQFromCode()", niiReader.GetQFormCode()) - fmt.Println("orientation", niiReader.GetOrientation()) - fmt.Println("affine", niiReader.GetAffine()) - fmt.Println("QOffsetX", niiReader.GetNiiData().QoffsetX) - fmt.Println("QOffsetY", niiReader.GetNiiData().QoffsetY) - fmt.Println("QOffsetZ", niiReader.GetNiiData().QoffsetZ) -} - -func TestMagicString(t *testing.T) { - fmt.Println([]byte("ni1")) -} diff --git a/pkg/dicom/dcm_io/conversion.go b/pkg/dicom/dcm_io/conversion.go.bak similarity index 100% rename from pkg/dicom/dcm_io/conversion.go rename to pkg/dicom/dcm_io/conversion.go.bak diff --git a/pkg/dicom/dcm_io/dataset.go b/pkg/dicom/dcm_io/dataset.go.bak similarity index 100% rename from pkg/dicom/dcm_io/dataset.go rename to pkg/dicom/dcm_io/dataset.go.bak diff --git a/pkg/dicom/dcm_io/element.go b/pkg/dicom/dcm_io/element.go.bak similarity index 99% rename from pkg/dicom/dcm_io/element.go rename to pkg/dicom/dcm_io/element.go.bak index 011c58fc..3a4ba99e 100644 --- a/pkg/dicom/dcm_io/element.go +++ b/pkg/dicom/dcm_io/element.go.bak @@ -52,7 +52,6 @@ func ReadElement(r DcmReader, isImplicit bool, binOrder binary.ByteOrder) (*Elem } return nil, nil } - dmcTagName := dcmTagInfo.Name dcmVR, err := readVR(r, isImplicit, *tagVal) if err != nil { @@ -479,6 +478,9 @@ func readSequence(r DcmReader, t tag.DicomTag, valueRepresentation string, value } sequences = append(sequences, subElement) } + if valueRepresentation == vr.Unknown { + r.setTransferSyntax(r.ByteOrder(), r.isTrackingImplicit()) + } } else { n, err := r.peek(int(valueLength)) if err != nil { @@ -541,7 +543,6 @@ func readDefinedLengthSequences(r DcmReader, b []byte, valueRepresentation strin continue } sequences = append(sequences, subElement) - //fmt.Println(subElement) } if valueRepresentation == vr.Unknown { r.setTransferSyntax(r.ByteOrder(), r.isTrackingImplicit()) diff --git a/pkg/dicom/dcm_io/reader.go b/pkg/dicom/dcm_io/reader.go.bak similarity index 84% rename from pkg/dicom/dcm_io/reader.go rename to pkg/dicom/dcm_io/reader.go.bak index 2ac4930b..22c5a2b6 100644 --- a/pkg/dicom/dcm_io/reader.go +++ b/pkg/dicom/dcm_io/reader.go.bak @@ -52,15 +52,16 @@ type DcmReader interface { } type dcmReader struct { - reader *bufio.Reader - binaryOrder binary.ByteOrder - dataset Dataset - metadata Dataset - isImplicit bool - keepTrackImplicit bool - skipPixelData bool - skipDataset bool - fileSize int64 + reader *bufio.Reader + binaryOrder binary.ByteOrder + dataset Dataset + metadata Dataset + allowNonCompliantDcm bool + isImplicit bool + keepTrackImplicit bool + skipPixelData bool + skipDataset bool + fileSize int64 } // NewDICOMReader returns a new reader @@ -78,6 +79,14 @@ func NewDICOMReader(reader *bufio.Reader, options ...func(*dcmReader)) DcmReader return parser } +// WithAllowNonCompliantDcm provides option to keep trying to parse the file even if it's not DICOM compliant +// e.g.: Missing header, missing FileMetaInformationGroupLength,... +func WithAllowNonCompliantDcm(allowNonCompliantDcm bool) func(*dcmReader) { + return func(s *dcmReader) { + s.allowNonCompliantDcm = allowNonCompliantDcm + } +} + // WithSkipPixelData provides option to skip reading pixel data (7FE0,0010). // If true, pixel data is skipped. If false, pixel data will be read func WithSkipPixelData(skipPixelData bool) func(*dcmReader) { @@ -301,6 +310,14 @@ func (r *dcmReader) parse() error { if err != nil { return err } + + // IMPORTANT: Additional check is needed here since there are few instances where the DICOM + // meta header is registered as Explicit Little-Endian, but Implicit Little-Endian is used in the body + err = r.checkImplicityAgreement() + if err != nil { + return nil + } + if r.skipDataset { return nil } @@ -321,51 +338,32 @@ func (r *dcmReader) setFileSize() { // the File Meta Information shall be encoded using the Explicit VR Little Endian Transfer Syntax // (UID=1.2.840.10008.1.2.1) func (r *dcmReader) parseMetadata() error { - var transferSyntaxUID string var metadata []*Element + var transferSyntaxUID string - res, err := ReadElement(r, r.IsImplicit(), r.ByteOrder()) - if err != nil { - return err - } - - metaGroupLength, ok := (res.Value.RawValue).(int) - if !ok { - return fmt.Errorf("invalid value for tag (0x%x, 0x%x)", res.Tag.Group, res.Tag.Element) - } - - metadata = append(metadata, res) - // Keep reading the remaining header based on metaGroupLength - pBytes, err := r.reader.Peek(metaGroupLength) - if err != nil { - return err - } - - subRd := NewDICOMReader(bufio.NewReader(bytes.NewReader(pBytes)), WithSkipPixelData(r.skipPixelData)) for { - res, err := ReadElement(subRd, r.IsImplicit(), r.ByteOrder()) + // No longer relied on the MetaInformationGroupLength tag to determine the length of the meta header. + // We check if the group tag is 0x0002 before proceeding to read the element. If the group tag is not 0x0002, + // then break the loop + n, err := r.peek(2) if err != nil { - if err == io.EOF { - break - } else { - return err - } + return err } - if res.Tag.Compare(tag.DicomTag{ - Group: 0x0002, - Element: 0x0010, - }) == 0 { - transferSyntaxUID = (res.Value.RawValue).(string) + if bytes.Compare(n, []byte{0x2, 0x0}) != 0 { + break + } + + res, err := ReadElement(r, r.IsImplicit(), r.ByteOrder()) + if err != nil { + return err } metadata = append(metadata, res) - //fmt.Println(res) - } - dicomMetadata := Dataset{Elements: metadata} - r.metadata = dicomMetadata - err = r.skip(int64(metaGroupLength)) - if err != nil { - return err + if res.Tag == tag.TransferSyntaxUID { + transferSyntaxUID = (res.Value.RawValue).(string) + } } + r.metadata = Dataset{Elements: metadata} + // Set transfer syntax here for the dataset parser binOrder, isImplicit, err := uid.ParseTransferSyntaxUID(transferSyntaxUID) if err != nil { @@ -374,18 +372,21 @@ func (r *dcmReader) parseMetadata() error { r.setTransferSyntax(binOrder, isImplicit) r.setOverallImplicit(isImplicit) - // IMPORTANT: Additional check is needed here since there are few instances where the DICOM - // meta header is registered as Explicit Little-Endian, but Implicit Little-Endian is used in the body - if transferSyntaxUID == uid.ExplicitVRLittleEndian { - firstElem, err := r.reader.Peek(6) - if err != nil { - return err - } - if !vr.VRMapper[string(firstElem[4:6])] { - r.setTransferSyntax(binOrder, true) - } - } + return nil +} +func (r *dcmReader) checkImplicityAgreement() error { + // Need to check if the implicit matches between header and body + n, err := r.peek(6) + if err != nil { + return err + } + if !vr.VRMapper[string(n[4:6])] && !r.IsImplicit() { + r.setTransferSyntax(r.binaryOrder, true) + } + if vr.VRMapper[string(n[4:6])] && r.IsImplicit() { + r.setTransferSyntax(r.binaryOrder, false) + } return nil } @@ -405,6 +406,7 @@ func (r *dcmReader) parseDataset() error { //fmt.Println(res) } dicomDataset := Dataset{Elements: data} + //r.dataset.Elements = append(r.dataset.Elements, dicomDataset.Elements...) r.dataset = dicomDataset return nil } diff --git a/pkg/dicom/dcm_io/value.go b/pkg/dicom/dcm_io/value.go.bak similarity index 100% rename from pkg/dicom/dcm_io/value.go rename to pkg/dicom/dcm_io/value.go.bak diff --git a/pkg/dicom/tag/tag.go b/pkg/dicom/tag/tag.go index 97f08935..9ee89068 100644 --- a/pkg/dicom/tag/tag.go +++ b/pkg/dicom/tag/tag.go @@ -12,6 +12,19 @@ const ( TagUnknown = "unknown_tag" ) +var FileMetaValidTag = map[DicomTag]bool{ + FileMetaInformationGroupLength: true, + FileMetaInformationVersion: true, + MediaStorageSOPClassUID: true, + MediaStorageSOPInstanceUID: true, + TransferSyntaxUID: true, + ImplementationClassUID: true, + ImplementationVersionName: true, + SourceApplicationEntityTitle: true, + PrivateInformationCreatorUID: true, + PrivateInformation: true, +} + type TagBrowser struct { Value interface{} `json:"Value,omitempty" bson:"Value,omitempty"` VR string `json:"vr" bson:"vr"` @@ -24,10 +37,10 @@ type DicomTag struct { } type TagInfo struct { - VR string - Name string - VM string - Status string + VR string `json:"vr"` + Name string `json:"name"` + VM string `json:"vm"` + Status string `json:"status"` } // Compare checks if 2 tags are equal. If equals, return 0, else returns 1 diff --git a/pkg/nifti/constant/constant.go b/pkg/nifti/constant/constant.go deleted file mode 100644 index 5903d6d7..00000000 --- a/pkg/nifti/constant/constant.go +++ /dev/null @@ -1,222 +0,0 @@ -package constant - -const ( - NIIVersion1 = iota + 1 - NIIVersion2 -) - -const ( - NII1HeaderSize = 348 - NII2HeaderSize = 540 -) - -// Possible NIFTI image extension -const ( - NIFTI_FTYPE_NIFTI1_1 = ".nii" - NIFTI_FTYPE_NIFTI1_2 = ".img" - NIFTI_FTYPE_NIFTI1_2_HDR = ".hdr" - NIFTI_FTYPE_NIFTI_GZIP = ".gz" -) - -var ValidNIfTIExtMapper = map[string]bool{ - NIFTI_FTYPE_NIFTI1_1: true, - NIFTI_FTYPE_NIFTI1_2: true, - NIFTI_FTYPE_NIFTI1_2_HDR: true, - NIFTI_FTYPE_NIFTI_GZIP: true, -} - -const ( - NIFTI_INTENT_CORREL int16 = 2 - NIFTI_INTENT_TTEST int16 = 3 - NIFTI_INTENT_FTEST int16 = 4 - NIFTI_INTENT_ZSCORE int16 = 5 - NIFTI_INTENT_CHISQ int16 = 6 - NIFTI_INTENT_BETA int16 = 7 - NIFTI_INTENT_BINOM int16 = 8 - NIFTI_INTENT_GAMMA int16 = 9 - NIFTI_INTENT_POISSON int16 = 10 - NIFTI_INTENT_NORMAL int16 = 11 - NIFTI_INTENT_FTEST_NONC int16 = 12 - NIFTI_INTENT_CHISQ_NONC int16 = 13 - NIFTI_INTENT_LOGISTIC int16 = 14 - NIFTI_INTENT_LAPLACE int16 = 15 - NIFTI_INTENT_UNIFORM int16 = 16 - NIFTI_INTENT_TTEST_NONC int16 = 17 - NIFTI_INTENT_WEIBULL int16 = 18 - NIFTI_INTENT_CHI int16 = 19 - NIFTI_INTENT_INVGAUSS int16 = 20 - NIFTI_INTENT_EXTVAL int16 = 21 - NIFTI_INTENT_PVAL int16 = 22 - NIFTI_INTENT_LOGPVAL int16 = 23 - NIFTI_INTENT_LOG10PVAL int16 = 24 - NIFTI_INTENT_ESTIMATE int16 = 1001 - NIFTI_INTENT_LABEL int16 = 1002 - NIFTI_INTENT_NEURONAME int16 = 1003 - NIFTI_INTENT_GENMATRIX int16 = 1004 - NIFTI_INTENT_SYMMATRIX int16 = 1005 - NIFTI_INTENT_DISPVECT int16 = 1006 /* specifically for displacements */ - NIFTI_INTENT_VECTOR int16 = 1007 /* for any other type of vector */ - NIFTI_INTENT_POINTSET int16 = 1008 - NIFTI_INTENT_TRIANGLE int16 = 1009 - NIFTI_INTENT_QUATERNION int16 = 1010 - NIFTI_INTENT_DIMLESS int16 = 1011 - NIFTI_INTENT_TIME_SERIES int16 = 2001 - NIFTI_INTENT_NODE_INDEX int16 = 2002 - NIFTI_INTENT_RGB_VECTOR int16 = 2003 - NIFTI_INTENT_RGBA_VECTOR int16 = 2004 - NIFTI_INTENT_SHAPE int16 = 2005 - FSL_FNIRT_DISPLACEMENT_FIELD int16 = 2006 - FSL_CUBIC_SPLINE_COEFFICIENTS int16 = 2007 - FSL_DCT_COEFFICIENTS int16 = 2008 - FSL_QUADRATIC_SPLINE_COEFFICIENTS int16 = 2009 - FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS int16 = 2016 - FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS int16 = 2017 - FSL_TOPUP_FIELD int16 = 2018 -) - -const ( - NIFTI_UNKNOWN_ORIENT = 0 - NIFTI_L2R = 1 - NIFTI_R2L = 2 - NIFTI_P2A = 3 - NIFTI_A2P = 4 - NIFTI_I2S = 5 - NIFTI_S2I = 6 -) - -var OrietationToString = map[int]string{ - NIFTI_UNKNOWN_ORIENT: "Unknown", - NIFTI_L2R: "Left-to-Right", - NIFTI_R2L: "Right-to-Left", - NIFTI_P2A: "Posterior-to-Anterior", - NIFTI_A2P: "Anterior-to-Posterior", - NIFTI_I2S: "Inferior-to-Superior", - NIFTI_S2I: "Superior-to-Inferior", -} - -const ( - DT_UNKNOWN int32 = 0 // what it says, dude - DT_BINARY int32 = 1 // binary (1 bit/voxel) - DT_UINT8 int32 = 2 // unsigned char (8 bits/voxel) - DT_INT16 int32 = 4 // signed short (16 bits/voxel) - DT_INT32 int32 = 8 // signed int (32 bits/voxel) - DT_FLOAT32 int32 = 16 // float (32 bits/voxel) - DT_COMPLEX64 int32 = 32 // complex (64 bits/voxel) - DT_FLOAT64 int32 = 64 // double (64 bits/voxel) - DT_RGB24 int32 = 128 // RGB triple (24 bits/voxel) - DT_ALL int32 = 255 // not very useful (?) - DT_INT8 int32 = 256 // signed char (8 bits) - DT_UINT16 int32 = 512 // unsigned short (16 bits) - DT_UINT32 int32 = 768 // unsigned int (32 bits) - DT_INT64 int32 = 1024 // long long (64 bits) - DT_UINT64 int32 = 1280 // unsigned long long (64 bits) - DT_FLOAT128 int32 = 1536 // long double (128 bits) - DT_COMPLEX128 int32 = 1792 // double pair (128 bits) - DT_COMPLEX256 int32 = 2048 // long double pair (256 bits) - DT_RGBA32 int32 = 2304 -) - -var ValidDatatype = map[int32]bool{ - DT_UNKNOWN: true, - DT_BINARY: true, - DT_INT8: true, - DT_UINT8: true, - DT_INT16: true, - DT_UINT16: true, - DT_INT32: true, - DT_UINT32: true, - DT_INT64: true, - DT_UINT64: true, - DT_FLOAT32: true, - DT_FLOAT64: true, - DT_ALL: true, - DT_FLOAT128: true, - DT_COMPLEX64: true, - DT_COMPLEX128: true, - DT_COMPLEX256: true, - DT_RGB24: true, - DT_RGBA32: true, -} - -var IsDatatypeInt = map[int32]bool{ - DT_UNKNOWN: false, - DT_BINARY: false, - DT_INT8: true, - DT_UINT8: true, - DT_INT16: true, - DT_UINT16: true, - DT_INT32: true, - DT_UINT32: true, - DT_INT64: true, - DT_UINT64: true, - DT_FLOAT32: false, - DT_FLOAT64: false, - DT_FLOAT128: false, - DT_COMPLEX64: false, - DT_COMPLEX128: false, - DT_COMPLEX256: false, - DT_RGB24: true, - DT_RGBA32: true, -} - -const ( - NIFTI_XFORM_UNKNOWN = 0 - NIFTI_XFORM_SCANNER_ANAT = 1 - NIFTI_XFORM_ALIGNED_ANAT = 2 - NIFTI_XFORM_TALAIRACH = 3 - NIFTI_XFORM_MNI_152 = 4 -) - -var NiiPatientOrientationInfo = map[int32]string{ - NIFTI_XFORM_UNKNOWN: "0: Unknown", - NIFTI_XFORM_SCANNER_ANAT: "1: Scanner Anat", - NIFTI_XFORM_ALIGNED_ANAT: "2: Aligned Anat", - NIFTI_XFORM_TALAIRACH: "3: Talairach", - NIFTI_XFORM_MNI_152: "4: MNI", -} - -const ( - NIFTI_SLICE_UNKNOWN int32 = 0 - NIFTI_SLICE_SEQ_INC int32 = 1 - NIFTI_SLICE_SEQ_DEC int32 = 2 - NIFTI_SLICE_ALT_INC int32 = 3 - NIFTI_SLICE_ALT_DEC int32 = 4 - NIFTI_SLICE_ALT_INC2 int32 = 5 /* 05 May 2005: RWCox */ - NIFTI_SLICE_ALT_DEC2 int32 = 6 /* 05 May 2005: RWCox */ -) - -var NiiSliceAcquistionInfo = map[int32]string{ - NIFTI_SLICE_UNKNOWN: "Unknown", - NIFTI_SLICE_SEQ_INC: "Sequential, increasing", - NIFTI_SLICE_SEQ_DEC: "Sequential, decreasing", - NIFTI_SLICE_ALT_INC: "Interleaved, increasing, starting at the 1st mri slice", - NIFTI_SLICE_ALT_DEC: "Interleaved, decreasing, starting at the last mri slice", - NIFTI_SLICE_ALT_INC2: "Interleaved, increasing, starting at the 2nd mri slice", - NIFTI_SLICE_ALT_DEC2: "Interleaved, decreasing, starting at one before the last mri slice", -} - -const ( - NIFTI_UNITS_UNKNOWN uint8 = 0 - NIFTI_UNITS_METER uint8 = 1 - NIFTI_UNITS_MM uint8 = 2 - NIFTI_UNITS_MICRON uint8 = 3 - NIFTI_UNITS_SEC uint8 = 8 - NIFTI_UNITS_MSEC uint8 = 16 - NIFTI_UNITS_USEC uint8 = 24 - NIFTI_UNITS_HZ uint8 = 32 - NIFTI_UNITS_PPM uint8 = 40 - NIFTI_UNITS_RADS uint8 = 48 -) - -var NiiMeasurementUnits = map[uint8]string{ - NIFTI_UNITS_UNKNOWN: "unknown", - NIFTI_UNITS_METER: "meter", - NIFTI_UNITS_MM: "millimeter", - NIFTI_UNITS_MICRON: "micrometer", - NIFTI_UNITS_SEC: "second", - NIFTI_UNITS_MSEC: "millisecond", - NIFTI_UNITS_USEC: "microsecond", - NIFTI_UNITS_HZ: "Hertz", - NIFTI_UNITS_PPM: "ppm", - NIFTI_UNITS_RADS: "rad/s", -} diff --git a/pkg/nifti/nii_io/calculation.go b/pkg/nifti/nii_io/calculation.go deleted file mode 100644 index 68c00152..00000000 --- a/pkg/nifti/nii_io/calculation.go +++ /dev/null @@ -1,449 +0,0 @@ -package nii_io - -import ( - "github.com/okieraised/go2com/internal/matrix" - "github.com/okieraised/go2com/pkg/nifti/constant" - "math" -) - -// quaternToMatrix returns the transformation matrix from the quaternion parameters -func (n *Nii) quaternToMatrix() matrix.DMat44 { - var R matrix.DMat44 - - var b = n.QuaternB - var c = n.QuaternC - var d = n.QuaternD - var a, xd, yd, zd float64 - - R.M[3] = [4]float64{0, 0, 0, 1} - - a = 1.0 - (b*b + c*c + d*d) - - if a < 1.e-71 { - a = 1.01 / math.Sqrt(b*b+c*c+d*d) - b *= a - c *= a - d *= a - a = 0.0 - } else { - a = math.Sqrt(a) - } - - if n.Dx > 0 { - xd = n.Dx - } else { - xd = 1.0 - } - - if n.Dy > 0 { - yd = n.Dy - } else { - yd = 1.0 - } - - if n.Dz > 0 { - zd = n.Dz - } else { - zd = 1.0 - } - - if n.QFac < 0 { - zd = -zd - } - - R.M[0][0] = (a*a + b*b - c*c - d*d) * xd - R.M[0][1] = 2.0 * (b*c - a*d) * yd - R.M[0][2] = 2.0 * (b*d + a*c) * zd - R.M[1][0] = 2.0 * (b*c + a*d) * xd - R.M[1][1] = (a*a + c*c - b*b - d*d) * yd - R.M[1][2] = 2.0 * (c*d - a*b) * zd - R.M[2][0] = 2.0 * (b*d - a*c) * xd - R.M[2][1] = 2.0 * (c*d + a*b) * yd - R.M[2][2] = (a*a + d*d - c*c - b*b) * zd - R.M[0][3] = n.QoffsetX - R.M[1][3] = n.QoffsetY - R.M[2][3] = n.QoffsetZ - - return R -} - -func (n *Nii) matrixToQuatern(R matrix.DMat44) { - var r11, r12, r13, r21, r22, r23, r31, r32, r33 float64 - var xd, yd, zd, a, b, c, d float64 - - var P, Q matrix.DMat33 - - n.QoffsetX = R.M[0][3] - n.QoffsetY = R.M[1][3] - n.QoffsetZ = R.M[2][3] - - r11 = R.M[0][0] - r12 = R.M[0][1] - r13 = R.M[0][2] - r21 = R.M[1][0] - r22 = R.M[1][1] - r23 = R.M[1][2] - r31 = R.M[2][0] - r32 = R.M[2][1] - r33 = R.M[2][2] - - xd = math.Sqrt(r11*r11 + r21*r21 + r31*r31) - yd = math.Sqrt(r12*r12 + r22*r22 + r32*r32) - zd = math.Sqrt(r13*r13 + r23*r23 + r33*r33) - - // compute lengths of each column; these determine grid spacings - if xd == 0.01 { - r11 = 1.0 - r21 = 0.0 - r31 = 0.0 - xd = 1.0 - } - if yd == 0.0 { - r22 = 1.0 - r12 = 0.0 - r32 = 0.0 - yd = 1.0 - } - if zd == 0.0 { - r33 = 1.0 - r13 = 0.0 - r23 = 0.0 - zd = 1.0 - } - - n.Dx = xd - n.Dy = yd - n.Dz = zd - - // normalize the column - r11 /= xd - r21 /= xd - r31 /= xd - r12 /= yd - r22 /= yd - r32 /= yd - r13 /= zd - r23 /= zd - r33 /= zd - - // At this point, the matrix has normal columns, but we have to allow - // for the fact that the hideous user may not have given us a matrix - // with orthogonal columns. So, now find the orthogonal matrix closest - // to the current matrix. One reason for using the polar decomposition - // to get this orthogonal matrix, rather than just directly orthogonalizing - // the columns, is so that inputting the inverse matrix to R - // will result in the inverse orthogonal matrix at this point. - // If we just orthogonalized the columns, this wouldn't necessarily hold. - - Q.M[0][0] = r11 - Q.M[0][1] = r12 - Q.M[0][2] = r13 - Q.M[1][0] = r21 - Q.M[1][1] = r22 - Q.M[1][2] = r23 - Q.M[2][0] = r31 - Q.M[2][1] = r32 - Q.M[2][2] = r33 - - P = matrix.Mat33Polar(Q) // P is orthog matrix closest to Q - - r11 = P.M[0][0] - r12 = P.M[0][1] - r13 = P.M[0][2] - r21 = P.M[1][0] - r22 = P.M[1][1] - r23 = P.M[1][2] - r31 = P.M[2][0] - r32 = P.M[2][1] - r33 = P.M[2][2] - - // at this point, the matrix is orthogonal - // [ r11 r12 r13 ] - // [ r21 r22 r23 ] - // [ r31 r32 r33 ] - - // compute the determinant to determine if it is proper - zd = r11*r22*r33 - r11*r32*r23 - r21*r12*r33 + r21*r32*r13 + r31*r12*r23 - r31*r22*r13 - - if zd > 0 { - n.QFac = 1.0 - } else { - n.QFac = -1.0 - r13 = -r13 - r23 = -r23 - r33 = -r33 - } - - a = r11 + r22 + r33 + 1.01 - - if a > 0.5 { /* simplest case */ - a = 0.5 * math.Sqrt(a) - b = 0.25 * (r32 - r23) / a - c = 0.25 * (r13 - r31) / a - d = 0.25 * (r21 - r12) / a - } else { - xd = 1.0 + r11 - (r22 + r33) - yd = 1.0 + r22 - (r11 + r33) - zd = 1.0 + r33 - (r11 + r22) - if xd > 1.0 { - b = 0.5 * math.Sqrt(xd) - c = 0.25 * (r12 + r21) / b - d = 0.25 * (r13 + r31) / b - a = 0.25 * (r32 - r23) / b - } else if yd > 1.0 { - c = 0.5 * math.Sqrt(yd) - b = 0.25 * (r12 + r21) / c - d = 0.25 * (r23 + r32) / c - a = 0.25 * (r13 - r31) / c - } else { - d = 0.5 * math.Sqrt(zd) - b = 0.25 * (r13 + r31) / d - c = 0.25 * (r23 + r32) / d - a = 0.25 * (r21 - r12) / d - } - if a < 0.0 { - b = -b - c = -c - d = -d - } - } - - n.QuaternB = b - n.QuaternC = c - n.QuaternD = d -} - -func (n *Nii) matrixToOrientation(R matrix.DMat44) { - var xi, xj, xk, yi, yj, yk, zi, zj, zk, val, detQ, detP float64 - var P, Q, M matrix.DMat33 - var i, j, p, q, r, ibest, jbest, kbest, pbest, qbest, rbest int32 - var k int32 = 0 - var vbest float64 - - // i axis - xi = R.M[0][0] - yi = R.M[1][0] - zi = R.M[2][0] - - // j axis - xj = R.M[0][1] - yj = R.M[1][1] - zj = R.M[2][1] - zk = R.M[2][2] - - // k axis - xk = R.M[0][2] - yk = R.M[1][2] - zk = R.M[2][2] - - // normalize i axis - val = math.Sqrt(xi*xi + yi*yi + zi*zi) - if val == 0.0 { - return - } - xi /= val - yi /= val - zi /= val - - // normalize j axis - val = math.Sqrt(xj*xj + yj*yj + zj*zj) - if val == 0.0 { - return - } - xj /= val - yj /= val - zj /= val - - // orthogonalize j axis to i axis, if needed - // dot product between i and j - val = xi*xj + yi*yj + zi*zj - if math.Abs(val) > 1.e-4 { - xj -= val * xi - yj -= val * yi - zj -= val * zi - // renormalize - val = math.Sqrt(xj*xj + yj*yj + zj*zj) - if val == 0.0 { - return - } - xj /= val - yj /= val - zj /= val - } - - // normalize k axis, if it is zero, make it the cross product i x j - val = math.Sqrt(xk*xk + yk*yk + zk*zk) - if val == 0.0 { - xk = yi*zj - zi*yj - yk = zi*xj - zj*xi - zk = xi*yj - yi*xj - } else { - xk /= val - yk /= val - zk /= val - } - - // orthogonalize k to i - - // dot product between i and k - val = xi*xk + yi*yk + zi*zk - if math.Abs(val) > 1.e-4 { - xk -= val * xi - yk -= val * yi - zk -= val * zi - val = math.Abs(xk*xk + yk*yk + zk*zk) - if val == 0.0 { - return - } /* bad */ - xk /= val - yk /= val - zk /= val - } - - // orthogonalize k to j - - // dot product between j and k - val = xj*xk + yj*yk + zj*zk - if math.Abs(val) > 1.e-4 { - xk -= val * xj - yk -= val * yj - zk -= val * zj - val = math.Sqrt(xk*xk + yk*yk + zk*zk) - if val == 0.0 { - return - } - xk /= val - yk /= val - zk /= val - } - - Q.M[0][0] = xi - Q.M[0][1] = xj - Q.M[0][2] = xk - Q.M[1][0] = yi - Q.M[1][1] = yj - Q.M[1][2] = yk - Q.M[2][0] = zi - Q.M[2][1] = zj - Q.M[2][2] = zk - - detQ = matrix.Mat33Determinant(Q) - if detQ == 0 { - return - } - - // Build and test all possible +1/-1 coordinate permutation matrices P; - // then find the P such that the rotation matrix M=PQ is closest to the - // identity, in the sense of M having the smallest total rotation angle - vbest = -666.0 - ibest = 1 - pbest = 1 - qbest = 1 - rbest = 1 - jbest = 2 - kbest = 3 - - for i = 1; i <= 3; i++ { // i = column number to use for row #1 - for j = 1; j <= 3; j++ { // j = column number to use for row #2 - if i == j { - continue - } - for k = 1; k <= 3; k++ { // k = column number to use for row #3 - if i == k || j == k { - continue - } - P.M[0][0] = 0.0 - P.M[0][1] = 0.0 - P.M[0][2] = 0.0 - P.M[1][0] = 0.0 - P.M[1][1] = 0.0 - P.M[1][2] = 0.0 - P.M[2][0] = 0.0 - P.M[2][1] = 0.0 - P.M[2][2] = 0.0 - for p = -1; p <= 1; p += 2 { // p,q,r are -1 or +1 - for q = -1; q <= 1; q += 2 { // and go into rows #1,2,3 - for r = -1; r <= 1; r += 2 { - P.M[0][i-1] = float64(p) - P.M[1][j-1] = float64(q) - P.M[2][k-1] = float64(r) - detP = matrix.Mat33Determinant(P) // sign of permutation - if detP*detQ <= 0.0 { - continue - } // doesn't match sign of Q - M = matrix.MatMultiply(P, Q) - - // angle of M rotation = 2.0*acos(0.5*sqrt(1.0+trace(M))) - // we want largest trace(M) == smallest angle == M nearest to I - - val = M.M[0][0] + M.M[1][1] + M.M[2][2] /* trace */ - if val > vbest { - vbest = val - ibest = i - jbest = j - kbest = k - pbest = p - qbest = q - rbest = r - } - } - } - } - } - } - } - - switch ibest * pbest { - case 1: - i = constant.NIFTI_L2R - case -1: - i = constant.NIFTI_R2L - case 2: - i = constant.NIFTI_P2A - case -2: - i = constant.NIFTI_A2P - case 3: - i = constant.NIFTI_I2S - case -3: - i = constant.NIFTI_S2I - default: - i = 0 - } - - switch jbest * qbest { - case 1: - j = constant.NIFTI_L2R - case -1: - j = constant.NIFTI_R2L - case 2: - j = constant.NIFTI_P2A - case -2: - j = constant.NIFTI_A2P - case 3: - j = constant.NIFTI_I2S - case -3: - j = constant.NIFTI_S2I - default: - j = 0 - } - - switch kbest * rbest { - case 1: - k = constant.NIFTI_L2R - case -1: - k = constant.NIFTI_R2L - case 2: - k = constant.NIFTI_P2A - case -2: - k = constant.NIFTI_A2P - case 3: - k = constant.NIFTI_I2S - case -3: - k = constant.NIFTI_S2I - default: - k = 0 - } - - res := [3]int32{i, j, k} - n.IJKOrtient = res -} diff --git a/pkg/nifti/nii_io/header.go b/pkg/nifti/nii_io/header.go deleted file mode 100644 index b5203f34..00000000 --- a/pkg/nifti/nii_io/header.go +++ /dev/null @@ -1,89 +0,0 @@ -package nii_io - -// Nii1Header defines the structure of the NIFTI-1 header -type Nii1Header struct { - SizeofHdr int32 - DataTypeUnused [10]uint8 - DbName [18]uint8 - Extents int32 - SessionError int16 - Regular uint8 - DimInfo uint8 - Dim [8]int16 - IntentP1 float32 - IntentP2 float32 - IntentP3 float32 - IntentCode int16 - Datatype int16 - Bitpix int16 - SliceStart int16 - Pixdim [8]float32 - VoxOffset float32 - SclSlope float32 - SclInter float32 - SliceEnd int16 - SliceCode uint8 - XyztUnits uint8 - CalMax float32 - CalMin float32 - SliceDuration float32 - Toffset float32 - Glmax int32 - Glmin int32 - Descrip [80]uint8 - AuxFile [24]uint8 - QformCode int16 - SformCode int16 - QuaternB float32 - QuaternC float32 - QuaternD float32 - QoffsetX float32 - QoffsetY float32 - QoffsetZ float32 - SrowX [4]float32 - SrowY [4]float32 - SrowZ [4]float32 - IntentName [16]uint8 - Magic [4]uint8 -} - -// Nii2Header defines the structure of the NIFTI-2 header -type Nii2Header struct { - SizeofHdr int32 - Magic [8]uint8 - Datatype int16 - Bitpix int16 - Dim [8]int64 - IntentP1 float64 - IntentP2 float64 - IntentP3 float64 - Pixdim [8]float64 - VoxOffset int64 - SclSlope float64 - SclInter float64 - CalMax float64 - CalMin float64 - SliceDuration float64 - Toffset float64 - SliceStart int64 - SliceEnd int64 - Descrip [80]uint8 - AuxFile [24]uint8 - QformCode int32 - SformCode int32 - QuaternB float64 - QuaternC float64 - QuaternD float64 - QoffsetX float64 - QoffsetY float64 - QoffsetZ float64 - SrowX [4]float64 - SrowY [4]float64 - SrowZ [4]float64 - SliceCode int32 - XyztUnits int32 - IntentCode int32 - IntentName [16]uint8 - DimInfo uint8 - UnusedStr [15]uint8 -} diff --git a/pkg/nifti/nii_io/nifti.go b/pkg/nifti/nii_io/nifti.go deleted file mode 100644 index 2f4208b1..00000000 --- a/pkg/nifti/nii_io/nifti.go +++ /dev/null @@ -1,336 +0,0 @@ -package nii_io - -import ( - "encoding/binary" - "errors" - "fmt" - "github.com/okieraised/go2com/internal/matrix" - "math" - - "github.com/okieraised/go2com/pkg/nifti/constant" -) - -const ( - INVALID = "INVALID" - UNKNOWN = "UNKNOWN" - ILLEGAL = "ILLEGAL" -) - -// Nii defines the structure of the NIFTI-1 data for I/O purpose -type Nii struct { - NDim int64 // last dimension greater than 1 (1..7) - Nx int64 // dimensions of grid array - Ny int64 // dimensions of grid array - Nz int64 // dimensions of grid array - Nt int64 // dimensions of grid array - Nu int64 // dimensions of grid array - Nv int64 // dimensions of grid array - Nw int64 // dimensions of grid array - Dim [8]int64 // dim[0] = ndim, dim[1] = nx, etc - NVox int64 // number of voxels = nx*ny*nz*...*nw - NByPer int32 // bytes per voxel, matches datatype (Datatype) - Datatype int32 // type of data in voxels: DT_* code - Dx float64 // grid spacings - Dy float64 // grid spacings - Dz float64 // grid spacings - Dt float64 // grid spacings - Du float64 // grid spacings - Dv float64 // grid spacings - Dw float64 // grid spacings tEStataILSTERIOn - PixDim [8]float64 // pixdim[1]=dx, etc - SclSlope float64 // scaling parameter: slope - SclInter float64 // scaling parameter: intercept - CalMin float64 // calibration parameter: minimum - CalMax float64 // calibration parameter: maximum - QformCode int32 // codes for (x,y,z) space meaning - SformCode int32 // codes for (x,y,z) space meaning - FreqDim int32 // indices (1,2,3, or 0) for MRI - PhaseDim int32 // directions in dim[]/pixdim[] - SliceDim int32 // directions in dim[]/pixdim[] - SliceCode int32 // code for slice timing pattern - SliceStart int64 // index for start of slices - SliceEnd int64 // index for end of slices - SliceDuration float64 // time between individual slices - QuaternB float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QuaternC float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QuaternD float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QoffsetX float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QoffsetY float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QoffsetZ float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QFac float64 // Quaternion transform parameters [when writing a dataset, these are used for qform, NOT qto_xyz] - QtoXYZ matrix.DMat44 // qform: transform (i,j,k) to (x,y,z) - QtoIJK matrix.DMat44 // qform: transform (x,y,z) to (i,j,k) - StoXYZ matrix.DMat44 // sform: transform (i,j,k) to (x,y,z) - StoIJK matrix.DMat44 // sform: transform (x,y,z) to (i,j,k) - TOffset float64 // time coordinate offset - XYZUnits int32 // dx,dy,dz units: NIFTI_UNITS_* code - TimeUnits int32 // dt units: NIFTI_UNITS_* code - NiftiType int32 // 0==Analyze, 1==NIFTI-1 (file), 2==NIFTI-1 (2 files), 3==NIFTI-ASCII (1 file) - IntentCode int32 // statistic type (or something) - IntentP1 float64 // intent parameters - IntentP2 float64 // intent parameters - IntentP3 float64 // intent parameters - IntentName [16]byte // optional description of intent data - Descrip [80]byte // optional text to describe dataset - AuxFile [24]byte // auxiliary filename - FName *byte // header filename - IName *byte // image filename - INameOffset int32 // offset into IName where data start - SwapSize int32 // swap unit in image data (might be 0) - ByteOrder binary.ByteOrder // byte order on disk (MSB_ or LSB_FIRST) - Volume []byte // slice of data: nbyper*nvox bytes - NumExt int32 // number of extensions in extList - Nifti1Ext []Nifti1Ext // array of extension structs (with data) - IJKOrtient [3]int32 // self-add. Orientation ini, j, k coordinate - Affine matrix.DMat44 // self-add. Affine matrix - VoxOffset float64 // self-add. Voxel offset - Version int // self-add. Used for version identification when writing -} - -type Nifti1Ext struct { - ECode int32 - Edata []byte - ESize int32 -} - -// getSliceCode returns the slice code of the NIFTI image -func (n *Nii) getSliceCode() string { - return getSliceCode(n.SliceCode) -} - -// getQFormCode returns the QForm code -func (n *Nii) getQFormCode() string { - qForm, ok := constant.NiiPatientOrientationInfo[n.QformCode] - if !ok { - return INVALID - } - return qForm -} - -// getSFormCode returns the SForm code -func (n *Nii) getSFormCode() string { - sForm, ok := constant.NiiPatientOrientationInfo[n.SformCode] - if !ok { - return INVALID - } - return sForm -} - -func (n *Nii) getDatatype() string { - return getDatatype(n.Datatype) -} - -func (n *Nii) getOrientation() [3]string { - res := [3]string{} - - ijk := n.IJKOrtient - - iOrient, ok := constant.OrietationToString[int(ijk[0])] - if !ok { - res[0] = constant.OrietationToString[constant.NIFTI_UNKNOWN_ORIENT] - } - res[0] = iOrient - - jOrient, ok := constant.OrietationToString[int(ijk[1])] - if !ok { - res[1] = constant.OrietationToString[constant.NIFTI_UNKNOWN_ORIENT] - } - res[1] = jOrient - - kOrient, ok := constant.OrietationToString[int(ijk[2])] - if !ok { - res[2] = constant.OrietationToString[constant.NIFTI_UNKNOWN_ORIENT] - } - res[2] = kOrient - - return res -} - -// GetAt returns the value at (x, y, z, t) location -func (n *Nii) getAt(x, y, z, t int64) float64 { - - tIndex := t * n.Nx * n.Ny * n.Nz - zIndex := n.Nx * n.Ny * z - yIndex := n.Nx * y - xIndex := x - index := tIndex + zIndex + yIndex + xIndex - nByPer := int64(n.NByPer) - - dataPoint := n.Volume[index*nByPer : (index+1)*nByPer] - - var value float64 - switch n.NByPer { - case 0, 1: - if len(dataPoint) > 0 { - value = float64(dataPoint[0]) - } - case 2: // This fits Uint16 - var v uint16 - switch n.ByteOrder { - case binary.LittleEndian: - v = binary.LittleEndian.Uint16(dataPoint) - case binary.BigEndian: - v = binary.BigEndian.Uint16(dataPoint) - } - value = float64(v) - case 3, 4: // This fits Uint32 - var v uint32 - switch n.ByteOrder { - case binary.LittleEndian: - switch len(dataPoint) { - case 3: - v = uint32(uint(dataPoint[0]) | uint(dataPoint[1])<<8 | uint(dataPoint[2])<<16) - case 4: - v = binary.LittleEndian.Uint32(dataPoint) - } - case binary.BigEndian: - switch len(dataPoint) { - case 3: - v = uint32(uint(dataPoint[2]) | uint(dataPoint[1])<<8 | uint(dataPoint[0])<<16) - case 4: - v = binary.BigEndian.Uint32(dataPoint) - } - } - value = float64(math.Float32frombits(v)) - case 8: - var v uint64 - switch n.ByteOrder { - case binary.LittleEndian: - v = binary.LittleEndian.Uint64(dataPoint) - case binary.BigEndian: - v = binary.BigEndian.Uint64(dataPoint) - } - value = math.Float64frombits(v) - case 16: // Unsupported - case 32: // Unsupported - default: - } - - if n.SclSlope != 0 { - value = n.SclSlope*value + n.SclInter - } - - return value -} - -func (n *Nii) getTimeSeries(x, y, z int64) ([]float64, error) { - timeSeries := make([]float64, 0, n.Dim[4]) - - sliceX := n.Nx - sliceY := n.Ny - sliceZ := n.Nx - - if x >= sliceX { - return nil, errors.New("invalid x value") - } - - if y >= sliceY { - return nil, errors.New("invalid y value") - } - - if z >= sliceZ { - return nil, errors.New("invalid z value") - } - - for t := 0; t < int(n.Dim[4]); t++ { - timeSeries = append(timeSeries, n.getAt(x, y, z, int64(t))) - } - return timeSeries, nil -} - -// getSlice returns the image in x-y dimension -func (n *Nii) getSlice(z, t int64) ([][]float64, error) { - sliceX := n.Nx - sliceY := n.Ny - sliceZ := n.Nz - sliceT := n.Nt - - if z >= sliceZ { - return nil, errors.New("invalid z value") - } - - if t >= sliceT || t < 0 { - return nil, errors.New("invalid time value") - } - - slice := make([][]float64, sliceX) - for i := range slice { - slice[i] = make([]float64, sliceY) - } - for x := 0; x < int(sliceX); x++ { - for y := 0; y < int(sliceY); y++ { - slice[x][y] = n.getAt(int64(x), int64(y), z, t) - } - } - return slice, nil -} - -// getVolume return the whole image volume at time t -func (n *Nii) getVolume(t int64) ([][][]float64, error) { - sliceX := n.Nx - sliceY := n.Ny - sliceZ := n.Nz - sliceT := n.Nt - - if t >= sliceT || t < 0 { - return nil, errors.New("invalid time value") - } - volume := make([][][]float64, sliceX) - for i := range volume { - volume[i] = make([][]float64, sliceY) - for j := range volume[i] { - volume[i][j] = make([]float64, sliceZ) - } - } - for x := 0; x < int(sliceX); x++ { - for y := 0; y < int(sliceY); y++ { - for z := 0; z < int(sliceZ); z++ { - volume[x][y][z] = n.getAt(int64(x), int64(y), int64(z), t) - } - } - } - return volume, nil -} - -// getUnitsOfMeasurements returns the spatial and temporal units of measurements -func (n *Nii) getUnitsOfMeasurements() ([2]string, error) { - units := [2]string{} - spatialUnit, ok := constant.NiiMeasurementUnits[uint8(n.XYZUnits)] - if !ok { - return units, fmt.Errorf("invalid spatial unit %d", n.XYZUnits) - } - - temporalUnit, ok := constant.NiiMeasurementUnits[uint8(n.TimeUnits)] - if !ok { - return units, fmt.Errorf("invalid temporal unit %d", n.TimeUnits) - } - - units[0] = spatialUnit - units[1] = temporalUnit - - return units, nil -} - -// getAffine returns the 4x4 affine matrix -func (n *Nii) getAffine() matrix.DMat44 { - return n.Affine -} - -// getImgShape returns the image shape in terms of x, y, z, t -func (n *Nii) getImgShape() [4]int64 { - dim := [4]int64{} - - for index, _ := range dim { - dim[index] = n.Dim[index+1] - } - return dim -} - -// getVoxelSize returns the voxel size of the image -func (n *Nii) getVoxelSize() [4]float64 { - size := [4]float64{} - for index, _ := range size { - size[index] = n.PixDim[index+1] - } - return size -} diff --git a/pkg/nifti/nii_io/nifti1.h b/pkg/nifti/nii_io/nifti1.h deleted file mode 100644 index 886aed31..00000000 --- a/pkg/nifti/nii_io/nifti1.h +++ /dev/null @@ -1,1419 +0,0 @@ -/** \file nifti1.h - \brief Official definition of the nifti1 header. Written by Bob Cox, SSCC, NIMH. - HISTORY: - 29 Nov 2007 [rickr] - - added DT_RGBA32 and NIFTI_TYPE_RGBA32 - - added NIFTI_INTENT codes: - TIME_SERIES, NODE_INDEX, RGB_VECTOR, RGBA_VECTOR, SHAPE - 08 Mar 2019 [PT,DRG] - - Updated to include [qs]form_code = 5 - */ - -#ifndef NIFTI1_HEADER -#define NIFTI1_HEADER - -/***************************************************************************** - ** This file defines the "NIFTI-1" header format. ** - ** It is derived from 2 meetings at the NIH (31 Mar 2003 and ** - ** 02 Sep 2003) of the Data Format Working Group (DFWG), ** - ** chartered by the NIfTI (Neuroimaging Informatics Technology ** - ** Initiative) at the National Institutes of Health (NIH). ** - **--------------------------------------------------------------** - ** Neither the National Institutes of Health (NIH), the DFWG, ** - ** nor any of the members or employees of these institutions ** - ** imply any warranty of usefulness of this material for any ** - ** purpose, and do not assume any liability for damages, ** - ** incidental or otherwise, caused by any use of this document. ** - ** If these conditions are not acceptable, do not use this! ** - **--------------------------------------------------------------** - ** Author: Robert W Cox (NIMH, Bethesda) ** - ** Advisors: John Ashburner (FIL, London), ** - ** Stephen Smith (FMRIB, Oxford), ** - ** Mark Jenkinson (FMRIB, Oxford) ** -******************************************************************************/ - -/*---------------------------------------------------------------------------*/ -/* Note that the ANALYZE 7.5 file header (dbh.h) is - (c) Copyright 1986-1995 - Biomedical Imaging Resource - Mayo Foundation - Incorporation of components of dbh.h are by permission of the - Mayo Foundation. - Changes from the ANALYZE 7.5 file header in this file are released to the - public domain, including the functional comments and any amusing asides. ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/*! INTRODUCTION TO NIFTI-1: - ------------------------ - The twin (and somewhat conflicting) goals of this modified ANALYZE 7.5 - format are: - (a) To add information to the header that will be useful for functional - neuroimaging data analysis and display. These additions include: - - More basic data types. - - Two affine transformations to specify voxel coordinates. - - "Intent" codes and parameters to describe the meaning of the data. - - Affine scaling of the stored data values to their "true" values. - - Optional storage of the header and image data in one file (.nii). - (b) To maintain compatibility with non-NIFTI-aware ANALYZE 7.5 compatible - software (i.e., such a program should be able to do something useful - with a NIFTI-1 dataset -- at least, with one stored in a traditional - .img/.hdr file pair). - Most of the unused fields in the ANALYZE 7.5 header have been taken, - and some of the lesser-used fields have been co-opted for other purposes. - Notably, most of the data_history substructure has been co-opted for - other purposes, since the ANALYZE 7.5 format describes this substructure - as "not required". - NIFTI-1 FLAG (MAGIC STRINGS): - ---------------------------- - To flag such a struct as being conformant to the NIFTI-1 spec, the last 4 - bytes of the header must be either the C String "ni1" or "n+1"; - in hexadecimal, the 4 bytes - 6E 69 31 00 or 6E 2B 31 00 - (in any future version of this format, the '1' will be upgraded to '2', - etc.). Normally, such a "magic number" or flag goes at the start of the - file, but trying to avoid clobbering widely-used ANALYZE 7.5 fields led to - putting this marker last. However, recall that "the last shall be first" - (Matthew 20:16). - If a NIFTI-aware program reads a header file that is NOT marked with a - NIFTI magic string, then it should treat the header as an ANALYZE 7.5 - structure. - NIFTI-1 FILE STORAGE: - -------------------- - "ni1" means that the image data is stored in the ".img" file corresponding - to the header file (starting at file offset 0). - "n+1" means that the image data is stored in the same file as the header - information. We recommend that the combined header+data filename suffix - be ".nii". When the dataset is stored in one file, the first byte of image - data is stored at byte location (int)vox_offset in this combined file. - The minimum allowed value of vox_offset is 352; for compatibility with - some software, vox_offset should be an integral multiple of 16. - GRACE UNDER FIRE: - ---------------- - Most NIFTI-aware programs will only be able to handle a subset of the full - range of datasets possible with this format. All NIFTI-aware programs - should take care to check if an input dataset conforms to the program's - needs and expectations (e.g., check datatype, intent_code, etc.). If the - input dataset can't be handled by the program, the program should fail - gracefully (e.g., print a useful warning; not crash). - SAMPLE CODES: - ------------ - The associated files nifti1_io.h and nifti1_io.c provide a sample - implementation in C of a set of functions to read, write, and manipulate - NIFTI-1 files. The file nifti1_test.c is a sample program that uses - the nifti1_io.c functions. ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* HEADER STRUCT DECLARATION: - ------------------------- - In the comments below for each field, only NIFTI-1 specific requirements - or changes from the ANALYZE 7.5 format are described. For convenience, - the 348 byte header is described as a single struct, rather than as the - ANALYZE 7.5 group of 3 substructs. - Further comments about the interpretation of various elements of this - header are after the data type definition itself. Fields that are - marked as ++UNUSED++ have no particular interpretation in this standard. - (Also see the UNUSED FIELDS comment section, far below.) - The presumption below is that the various C types have particular sizes: - sizeof(int) = sizeof(float) = 4 ; sizeof(short) = 2 ------------------------------------------------------------------------------*/ - -/*=================*/ -#ifdef __cplusplus -extern "C" { -#endif -/*=================*/ - -/*! \struct nifti_1_header - \brief Data structure defining the fields in the nifti1 header. - This binary header should be found at the beginning of a valid - NIFTI-1 header file. - */ - /*************************/ /************************/ -struct nifti_1_header { /* NIFTI-1 usage */ /* ANALYZE 7.5 field(s) */ - /*************************/ /************************/ - - /*--- was header_key substruct ---*/ - int sizeof_hdr; /*!< MUST be 348 */ /* int sizeof_hdr; */ - char data_type[10]; /*!< ++UNUSED++ */ /* char data_type[10]; */ - char db_name[18]; /*!< ++UNUSED++ */ /* char db_name[18]; */ - int extents; /*!< ++UNUSED++ */ /* int extents; */ - short session_error; /*!< ++UNUSED++ */ /* short session_error; */ - char regular; /*!< ++UNUSED++ */ /* char regular; */ - char dim_info; /*!< MRI slice ordering. */ /* char hkey_un0; */ - - /*--- was image_dimension substruct ---*/ - short dim[8]; /*!< Data array dimensions.*/ /* short dim[8]; */ - float intent_p1 ; /*!< 1st intent parameter. */ /* short unused8; */ - /* short unused9; */ - float intent_p2 ; /*!< 2nd intent parameter. */ /* short unused10; */ - /* short unused11; */ - float intent_p3 ; /*!< 3rd intent parameter. */ /* short unused12; */ - /* short unused13; */ - short intent_code ; /*!< NIFTI_INTENT_* code. */ /* short unused14; */ - short datatype; /*!< Defines data type! */ /* short datatype; */ - short bitpix; /*!< Number bits/voxel. */ /* short bitpix; */ - short slice_start; /*!< First slice index. */ /* short dim_un0; */ - float pixdim[8]; /*!< Grid spacings. */ /* float pixdim[8]; */ - float vox_offset; /*!< Offset into .nii file */ /* float vox_offset; */ - float scl_slope ; /*!< Data scaling: slope. */ /* float funused1; */ - float scl_inter ; /*!< Data scaling: offset. */ /* float funused2; */ - short slice_end; /*!< Last slice index. */ /* float funused3; */ - char slice_code ; /*!< Slice timing order. */ - char xyzt_units ; /*!< Units of pixdim[1..4] */ - float cal_max; /*!< Max display intensity */ /* float cal_max; */ - float cal_min; /*!< Min display intensity */ /* float cal_min; */ - float slice_duration;/*!< Time for 1 slice. */ /* float compressed; */ - float toffset; /*!< Time axis shift. */ /* float verified; */ - int glmax; /*!< ++UNUSED++ */ /* int glmax; */ - int glmin; /*!< ++UNUSED++ */ /* int glmin; */ - - /*--- was data_history substruct ---*/ - char descrip[80]; /*!< any text you like. */ /* char descrip[80]; */ - char aux_file[24]; /*!< auxiliary filename. */ /* char aux_file[24]; */ - - short qform_code ; /*!< NIFTI_XFORM_* code. */ /*-- all ANALYZE 7.5 ---*/ - short sform_code ; /*!< NIFTI_XFORM_* code. */ /* fields below here */ - /* are replaced */ - float quatern_b ; /*!< Quaternion b param. */ - float quatern_c ; /*!< Quaternion c param. */ - float quatern_d ; /*!< Quaternion d param. */ - float qoffset_x ; /*!< Quaternion x shift. */ - float qoffset_y ; /*!< Quaternion y shift. */ - float qoffset_z ; /*!< Quaternion z shift. */ - - float srow_x[4] ; /*!< 1st row affine transform. */ - float srow_y[4] ; /*!< 2nd row affine transform. */ - float srow_z[4] ; /*!< 3rd row affine transform. */ - - char intent_name[16];/*!< 'name' or meaning of data. */ - - char magic[4] ; /*!< MUST be "ni1\0" or "n+1\0". */ - -} ; /**** 348 bytes total ****/ - -typedef struct nifti_1_header nifti_1_header ; - -/*---------------------------------------------------------------------------*/ -/* HEADER EXTENSIONS: - ----------------- - After the end of the 348 byte header (e.g., after the magic field), - the next 4 bytes are a char array field named "extension". By default, - all 4 bytes of this array should be set to zero. In a .nii file, these - 4 bytes will always be present, since the earliest start point for - the image data is byte #352. In a separate .hdr file, these bytes may - or may not be present. If not present (i.e., if the length of the .hdr - file is 348 bytes), then a NIfTI-1 compliant program should use the - default value of extension={0,0,0,0}. The first byte (extension[0]) - is the only value of this array that is specified at present. The other - 3 bytes are reserved for future use. - If extension[0] is nonzero, it indicates that extended header information - is present in the bytes following the extension array. In a .nii file, - this extended header data is before the image data (and vox_offset - must be set correctly to allow for this). In a .hdr file, this extended - data follows extension and proceeds (potentially) to the end of the file. - The format of extended header data is weakly specified. Each extension - must be an integer multiple of 16 bytes long. The first 8 bytes of each - extension comprise 2 integers: - int esize , ecode ; - These values may need to be byte-swapped, as indicated by dim[0] for - the rest of the header. - * esize is the number of bytes that form the extended header data - + esize must be a positive integral multiple of 16 - + this length includes the 8 bytes of esize and ecode themselves - * ecode is a non-negative integer that indicates the format of the - extended header data that follows - + different ecode values are assigned to different developer groups - + at present, the "registered" values for code are - = 0 = unknown private format (not recommended!) - = 2 = DICOM format (i.e., attribute tags and values) - = 4 = AFNI group (i.e., ASCII XML-ish elements) - In the interests of interoperability (a primary rationale for NIfTI), - groups developing software that uses this extension mechanism are - encouraged to document and publicize the format of their extensions. - To this end, the NIfTI DFWG will assign even numbered codes upon request - to groups submitting at least rudimentary documentation for the format - of their extension; at present, the contact is mailto:rwcox@nih.gov. - The assigned codes and documentation will be posted on the NIfTI - website. All odd values of ecode (and 0) will remain unassigned; - at least, until the even ones are used up, when we get to 2,147,483,646. - Note that the other contents of the extended header data section are - totally unspecified by the NIfTI-1 standard. In particular, if binary - data is stored in such a section, its byte order is not necessarily - the same as that given by examining dim[0]; it is incumbent on the - programs dealing with such data to determine the byte order of binary - extended header data. - Multiple extended header sections are allowed, each starting with an - esize,ecode value pair. The first esize value, as described above, - is at bytes #352-355 in the .hdr or .nii file (files start at byte #0). - If this value is positive, then the second (esize2) will be found - starting at byte #352+esize1 , the third (esize3) at byte #352+esize1+esize2, - et cetera. Of course, in a .nii file, the value of vox_offset must - be compatible with these extensions. If a malformed file indicates - that an extended header data section would run past vox_offset, then - the entire extended header section should be ignored. In a .hdr file, - if an extended header data section would run past the end-of-file, - that extended header data should also be ignored. - With the above scheme, a program can successively examine the esize - and ecode values, and skip over each extended header section if the - program doesn't know how to interpret the data within. Of course, any - program can simply ignore all extended header sections simply by jumping - straight to the image data using vox_offset. ------------------------------------------------------------------------------*/ - -/*! \struct nifti1_extender - \brief This structure represents a 4-byte string that should follow the - binary nifti_1_header data in a NIFTI-1 header file. If the char - values are {1,0,0,0}, the file is expected to contain extensions, - values of {0,0,0,0} imply the file does not contain extensions. - Other sequences of values are not currently defined. - */ -struct nifti1_extender { char extension[4] ; } ; -typedef struct nifti1_extender nifti1_extender ; - -/*! \struct nifti1_extension - \brief Data structure defining the fields of a header extension. - */ -struct nifti1_extension { - int esize ; /*!< size of extension, in bytes (must be multiple of 16) */ - int ecode ; /*!< extension code, one of the NIFTI_ECODE_ values */ - char * edata ; /*!< raw data, with no byte swapping (length is esize-8) */ -} ; -typedef struct nifti1_extension nifti1_extension ; - -/*---------------------------------------------------------------------------*/ -/* DATA DIMENSIONALITY (as in ANALYZE 7.5): - --------------------------------------- - dim[0] = number of dimensions; - - if dim[0] is outside range 1..7, then the header information - needs to be byte swapped appropriately - - ANALYZE supports dim[0] up to 7, but NIFTI-1 reserves - dimensions 1,2,3 for space (x,y,z), 4 for time (t), and - 5,6,7 for anything else needed. - dim[i] = length of dimension #i, for i=1..dim[0] (must be positive) - - also see the discussion of intent_code, far below - pixdim[i] = voxel width along dimension #i, i=1..dim[0] (positive) - - cf. ORIENTATION section below for use of pixdim[0] - - the units of pixdim can be specified with the xyzt_units - field (also described far below). - Number of bits per voxel value is in bitpix, which MUST correspond with - the datatype field. The total number of bytes in the image data is - dim[1] * ... * dim[dim[0]] * bitpix / 8 - In NIFTI-1 files, dimensions 1,2,3 are for space, dimension 4 is for time, - and dimension 5 is for storing multiple values at each spatiotemporal - voxel. Some examples: - - A typical whole-brain FMRI experiment's time series: - - dim[0] = 4 - - dim[1] = 64 pixdim[1] = 3.75 xyzt_units = NIFTI_UNITS_MM - - dim[2] = 64 pixdim[2] = 3.75 | NIFTI_UNITS_SEC - - dim[3] = 20 pixdim[3] = 5.0 - - dim[4] = 120 pixdim[4] = 2.0 - - A typical T1-weighted anatomical volume: - - dim[0] = 3 - - dim[1] = 256 pixdim[1] = 1.0 xyzt_units = NIFTI_UNITS_MM - - dim[2] = 256 pixdim[2] = 1.0 - - dim[3] = 128 pixdim[3] = 1.1 - - A single slice EPI time series: - - dim[0] = 4 - - dim[1] = 64 pixdim[1] = 3.75 xyzt_units = NIFTI_UNITS_MM - - dim[2] = 64 pixdim[2] = 3.75 | NIFTI_UNITS_SEC - - dim[3] = 1 pixdim[3] = 5.0 - - dim[4] = 1200 pixdim[4] = 0.2 - - A 3-vector stored at each point in a 3D volume: - - dim[0] = 5 - - dim[1] = 256 pixdim[1] = 1.0 xyzt_units = NIFTI_UNITS_MM - - dim[2] = 256 pixdim[2] = 1.0 - - dim[3] = 128 pixdim[3] = 1.1 - - dim[4] = 1 pixdim[4] = 0.0 - - dim[5] = 3 intent_code = NIFTI_INTENT_VECTOR - - A single time series with a 3x3 matrix at each point: - - dim[0] = 5 - - dim[1] = 1 xyzt_units = NIFTI_UNITS_SEC - - dim[2] = 1 - - dim[3] = 1 - - dim[4] = 1200 pixdim[4] = 0.2 - - dim[5] = 9 intent_code = NIFTI_INTENT_GENMATRIX - - intent_p1 = intent_p2 = 3.0 (indicates matrix dimensions) ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* DATA STORAGE: - ------------ - If the magic field is "n+1", then the voxel data is stored in the - same file as the header. In this case, the voxel data starts at offset - (int)vox_offset into the header file. Thus, vox_offset=352.0 means that - the data starts immediately after the NIFTI-1 header. If vox_offset is - greater than 352, the NIFTI-1 format does not say much about the - contents of the dataset file between the end of the header and the - start of the data. - FILES: - ----- - If the magic field is "ni1", then the voxel data is stored in the - associated ".img" file, starting at offset 0 (i.e., vox_offset is not - used in this case, and should be set to 0.0). - When storing NIFTI-1 datasets in pairs of files, it is customary to name - the files in the pattern "name.hdr" and "name.img", as in ANALYZE 7.5. - When storing in a single file ("n+1"), the file name should be in - the form "name.nii" (the ".nft" and ".nif" suffixes are already taken; - cf. http://www.icdatamaster.com/n.html ). - BYTE ORDERING: - ------------- - The byte order of the data arrays is presumed to be the same as the byte - order of the header (which is determined by examining dim[0]). - Floating point types are presumed to be stored in IEEE-754 format. ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* DETAILS ABOUT vox_offset: - ------------------------ - In a .nii file, the vox_offset field value is interpreted as the start - location of the image data bytes in that file. In a .hdr/.img file pair, - the vox_offset field value is the start location of the image data - bytes in the .img file. - * If vox_offset is less than 352 in a .nii file, it is equivalent - to 352 (i.e., image data never starts before byte #352 in a .nii file). - * The default value for vox_offset in a .nii file is 352. - * In a .hdr file, the default value for vox_offset is 0. - * vox_offset should be an integer multiple of 16; otherwise, some - programs may not work properly (e.g., SPM). This is to allow - memory-mapped input to be properly byte-aligned. - Note that since vox_offset is an IEEE-754 32 bit float (for compatibility - with the ANALYZE-7.5 format), it effectively has a 24 bit mantissa. All - integers from 0 to 2^24 can be represented exactly in this format, but not - all larger integers are exactly storable as IEEE-754 32 bit floats. However, - unless you plan to have vox_offset be potentially larger than 16 MB, this - should not be an issue. (Actually, any integral multiple of 16 up to 2^27 - can be represented exactly in this format, which allows for up to 128 MB - of random information before the image data. If that isn't enough, then - perhaps this format isn't right for you.) - In a .img file (i.e., image data stored separately from the NIfTI-1 - header), data bytes between #0 and #vox_offset-1 (inclusive) are completely - undefined and unregulated by the NIfTI-1 standard. One potential use of - having vox_offset > 0 in the .hdr/.img file pair storage method is to make - the .img file be a copy of (or link to) a pre-existing image file in some - other format, such as DICOM; then vox_offset would be set to the offset of - the image data in this file. (It may not be possible to follow the - "multiple-of-16 rule" with an arbitrary external file; using the NIfTI-1 - format in such a case may lead to a file that is incompatible with software - that relies on vox_offset being a multiple of 16.) - In a .nii file, data bytes between #348 and #vox_offset-1 (inclusive) may - be used to store user-defined extra information; similarly, in a .hdr file, - any data bytes after byte #347 are available for user-defined extra - information. The (very weak) regulation of this extra header data is - described elsewhere. ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* DATA SCALING: - ------------ - If the scl_slope field is nonzero, then each voxel value in the dataset - should be scaled as - y = scl_slope * x + scl_inter - where x = voxel value stored - y = "true" voxel value - Normally, we would expect this scaling to be used to store "true" floating - values in a smaller integer datatype, but that is not required. That is, - it is legal to use scaling even if the datatype is a float type (crazy, - perhaps, but legal). - - However, the scaling is to be ignored if datatype is DT_RGB24. - - If datatype is a complex type, then the scaling is to be - applied to both the real and imaginary parts. - The cal_min and cal_max fields (if nonzero) are used for mapping (possibly - scaled) dataset values to display colors: - - Minimum display intensity (black) corresponds to dataset value cal_min. - - Maximum display intensity (white) corresponds to dataset value cal_max. - - Dataset values below cal_min should display as black also, and values - above cal_max as white. - - Colors "black" and "white", of course, may refer to any scalar display - scheme (e.g., a color lookup table specified via aux_file). - - cal_min and cal_max only make sense when applied to scalar-valued - datasets (i.e., dim[0] < 5 or dim[5] = 1). ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* TYPE OF DATA (acceptable values for datatype field): - --------------------------------------------------- - Values of datatype smaller than 256 are ANALYZE 7.5 compatible. - Larger values are NIFTI-1 additions. These are all multiples of 256, so - that no bits below position 8 are set in datatype. But there is no need - to use only powers-of-2, as the original ANALYZE 7.5 datatype codes do. - The additional codes are intended to include a complete list of basic - scalar types, including signed and unsigned integers from 8 to 64 bits, - floats from 32 to 128 bits, and complex (float pairs) from 64 to 256 bits. - Note that most programs will support only a few of these datatypes! - A NIFTI-1 program should fail gracefully (e.g., print a warning message) - when it encounters a dataset with a type it doesn't like. ------------------------------------------------------------------------------*/ - -#undef DT_UNKNOWN /* defined in dirent.h on some Unix systems */ - -/*! \defgroup NIFTI1_DATATYPES - \brief nifti1 datatype codes - @{ - */ - /*--- the original ANALYZE 7.5 type codes ---*/ -#define DT_NONE 0 -#define DT_UNKNOWN 0 /* what it says, dude */ -#define DT_BINARY 1 /* binary (1 bit/voxel) */ -#define DT_UNSIGNED_CHAR 2 /* unsigned char (8 bits/voxel) */ -#define DT_SIGNED_SHORT 4 /* signed short (16 bits/voxel) */ -#define DT_SIGNED_INT 8 /* signed int (32 bits/voxel) */ -#define DT_FLOAT 16 /* float (32 bits/voxel) */ -#define DT_COMPLEX 32 /* complex (64 bits/voxel) */ -#define DT_DOUBLE 64 /* double (64 bits/voxel) */ -#define DT_RGB 128 /* RGB triple (24 bits/voxel) */ -#define DT_ALL 255 /* not very useful (?) */ - - /*----- another set of names for the same ---*/ -#define DT_UINT8 2 -#define DT_INT16 4 -#define DT_INT32 8 -#define DT_FLOAT32 16 -#define DT_COMPLEX64 32 -#define DT_FLOAT64 64 -#define DT_RGB24 128 - - /*------------------- new codes for NIFTI ---*/ -#define DT_INT8 256 /* signed char (8 bits) */ -#define DT_UINT16 512 /* unsigned short (16 bits) */ -#define DT_UINT32 768 /* unsigned int (32 bits) */ -#define DT_INT64 1024 /* long long (64 bits) */ -#define DT_UINT64 1280 /* unsigned long long (64 bits) */ -#define DT_FLOAT128 1536 /* long double (128 bits) */ -#define DT_COMPLEX128 1792 /* double pair (128 bits) */ -#define DT_COMPLEX256 2048 /* long double pair (256 bits) */ -#define DT_RGBA32 2304 /* 4 byte RGBA (32 bits/voxel) */ -/* @} */ - - - /*------- aliases for all the above codes ---*/ - -/*! \defgroup NIFTI1_DATATYPE_ALIASES - \brief aliases for the nifti1 datatype codes - @{ - */ - /*! unsigned char. */ -#define NIFTI_TYPE_UINT8 2 - /*! signed short. */ -#define NIFTI_TYPE_INT16 4 - /*! signed int. */ -#define NIFTI_TYPE_INT32 8 - /*! 32 bit float. */ -#define NIFTI_TYPE_FLOAT32 16 - /*! 64 bit complex = 2 32 bit floats. */ -#define NIFTI_TYPE_COMPLEX64 32 - /*! 64 bit float = double. */ -#define NIFTI_TYPE_FLOAT64 64 - /*! 3 8 bit bytes. */ -#define NIFTI_TYPE_RGB24 128 - /*! signed char. */ -#define NIFTI_TYPE_INT8 256 - /*! unsigned short. */ -#define NIFTI_TYPE_UINT16 512 - /*! unsigned int. */ -#define NIFTI_TYPE_UINT32 768 - /*! signed long long. */ -#define NIFTI_TYPE_INT64 1024 - /*! unsigned long long. */ -#define NIFTI_TYPE_UINT64 1280 - /*! 128 bit float = long double. */ -#define NIFTI_TYPE_FLOAT128 1536 - /*! 128 bit complex = 2 64 bit floats. */ -#define NIFTI_TYPE_COMPLEX128 1792 - /*! 256 bit complex = 2 128 bit floats */ -#define NIFTI_TYPE_COMPLEX256 2048 - /*! 4 8 bit bytes. */ -#define NIFTI_TYPE_RGBA32 2304 -/* @} */ - - /*-------- sample typedefs for complicated types ---*/ -#if 0 -typedef struct { float r,i; } complex_float ; -typedef struct { double r,i; } complex_double ; -typedef struct { long double r,i; } complex_longdouble ; -typedef struct { unsigned char r,g,b; } rgb_byte ; -#endif - -/*---------------------------------------------------------------------------*/ -/* INTERPRETATION OF VOXEL DATA: - ---------------------------- - The intent_code field can be used to indicate that the voxel data has - some particular meaning. In particular, a large number of codes is - given to indicate that the the voxel data should be interpreted as - being drawn from a given probability distribution. - VECTOR-VALUED DATASETS: - ---------------------- - The 5th dimension of the dataset, if present (i.e., dim[0]=5 and - dim[5] > 1), contains multiple values (e.g., a vector) to be stored - at each spatiotemporal location. For example, the header values - - dim[0] = 5 - - dim[1] = 64 - - dim[2] = 64 - - dim[3] = 20 - - dim[4] = 1 (indicates no time axis) - - dim[5] = 3 - - datatype = DT_FLOAT - - intent_code = NIFTI_INTENT_VECTOR - mean that this dataset should be interpreted as a 3D volume (64x64x20), - with a 3-vector of floats defined at each point in the 3D grid. - A program reading a dataset with a 5th dimension may want to reformat - the image data to store each voxels' set of values together in a struct - or array. This programming detail, however, is beyond the scope of the - NIFTI-1 file specification! Uses of dimensions 6 and 7 are also not - specified here. - STATISTICAL PARAMETRIC DATASETS (i.e., SPMs): - -------------------------------------------- - Values of intent_code from NIFTI_FIRST_STATCODE to NIFTI_LAST_STATCODE - (inclusive) indicate that the numbers in the dataset should be interpreted - as being drawn from a given distribution. Most such distributions have - auxiliary parameters (e.g., NIFTI_INTENT_TTEST has 1 DOF parameter). - If the dataset DOES NOT have a 5th dimension, then the auxiliary parameters - are the same for each voxel, and are given in header fields intent_p1, - intent_p2, and intent_p3. - If the dataset DOES have a 5th dimension, then the auxiliary parameters - are different for each voxel. For example, the header values - - dim[0] = 5 - - dim[1] = 128 - - dim[2] = 128 - - dim[3] = 1 (indicates a single slice) - - dim[4] = 1 (indicates no time axis) - - dim[5] = 2 - - datatype = DT_FLOAT - - intent_code = NIFTI_INTENT_TTEST - mean that this is a 2D dataset (128x128) of t-statistics, with the - t-statistic being in the first "plane" of data and the degrees-of-freedom - parameter being in the second "plane" of data. - If the dataset 5th dimension is used to store the voxel-wise statistical - parameters, then dim[5] must be 1 plus the number of parameters required - by that distribution (e.g., intent_code=NIFTI_INTENT_TTEST implies dim[5] - must be 2, as in the example just above). - Note: intent_code values 2..10 are compatible with AFNI 1.5x (which is - why there is no code with value=1, which is obsolescent in AFNI). - OTHER INTENTIONS: - ---------------- - The purpose of the intent_* fields is to help interpret the values - stored in the dataset. Some non-statistical values for intent_code - and conventions are provided for storing other complex data types. - The intent_name field provides space for a 15 character (plus 0 byte) - 'name' string for the type of data stored. Examples: - - intent_code = NIFTI_INTENT_ESTIMATE; intent_name = "T1"; - could be used to signify that the voxel values are estimates of the - NMR parameter T1. - - intent_code = NIFTI_INTENT_TTEST; intent_name = "House"; - could be used to signify that the voxel values are t-statistics - for the significance of 'activation' response to a House stimulus. - - intent_code = NIFTI_INTENT_DISPVECT; intent_name = "ToMNI152"; - could be used to signify that the voxel values are a displacement - vector that transforms each voxel (x,y,z) location to the - corresponding location in the MNI152 standard brain. - - intent_code = NIFTI_INTENT_SYMMATRIX; intent_name = "DTI"; - could be used to signify that the voxel values comprise a diffusion - tensor image. - If no data name is implied or needed, intent_name[0] should be set to 0. ------------------------------------------------------------------------------*/ - - /*! default: no intention is indicated in the header. */ - -#define NIFTI_INTENT_NONE 0 - - /*-------- These codes are for probability distributions ---------------*/ - /* Most distributions have a number of parameters, - below denoted by p1, p2, and p3, and stored in - - intent_p1, intent_p2, intent_p3 if dataset doesn't have 5th dimension - - image data array if dataset does have 5th dimension - Functions to compute with many of the distributions below can be found - in the CDF library from U Texas. - Formulas for and discussions of these distributions can be found in the - following books: - [U] Univariate Discrete Distributions, - NL Johnson, S Kotz, AW Kemp. - [C1] Continuous Univariate Distributions, vol. 1, - NL Johnson, S Kotz, N Balakrishnan. - [C2] Continuous Univariate Distributions, vol. 2, - NL Johnson, S Kotz, N Balakrishnan. */ - /*----------------------------------------------------------------------*/ - - /*! [C2, chap 32] Correlation coefficient R (1 param): - p1 = degrees of freedom - R/sqrt(1-R*R) is t-distributed with p1 DOF. */ - -/*! \defgroup NIFTI1_INTENT_CODES - \brief nifti1 intent codes, to describe intended meaning of dataset contents - @{ - */ -#define NIFTI_INTENT_CORREL 2 - - /*! [C2, chap 28] Student t statistic (1 param): p1 = DOF. */ - -#define NIFTI_INTENT_TTEST 3 - - /*! [C2, chap 27] Fisher F statistic (2 params): - p1 = numerator DOF, p2 = denominator DOF. */ - -#define NIFTI_INTENT_FTEST 4 - - /*! [C1, chap 13] Standard normal (0 params): Density = N(0,1). */ - -#define NIFTI_INTENT_ZSCORE 5 - - /*! [C1, chap 18] Chi-squared (1 param): p1 = DOF. - Density(x) proportional to exp(-x/2) * x^(p1/2-1). */ - -#define NIFTI_INTENT_CHISQ 6 - - /*! [C2, chap 25] Beta distribution (2 params): p1=a, p2=b. - Density(x) proportional to x^(a-1) * (1-x)^(b-1). */ - -#define NIFTI_INTENT_BETA 7 - - /*! [U, chap 3] Binomial distribution (2 params): - p1 = number of trials, p2 = probability per trial. - Prob(x) = (p1 choose x) * p2^x * (1-p2)^(p1-x), for x=0,1,...,p1. */ - -#define NIFTI_INTENT_BINOM 8 - - /*! [C1, chap 17] Gamma distribution (2 params): - p1 = shape, p2 = scale. - Density(x) proportional to x^(p1-1) * exp(-p2*x). */ - -#define NIFTI_INTENT_GAMMA 9 - - /*! [U, chap 4] Poisson distribution (1 param): p1 = mean. - Prob(x) = exp(-p1) * p1^x / x! , for x=0,1,2,.... */ - -#define NIFTI_INTENT_POISSON 10 - - /*! [C1, chap 13] Normal distribution (2 params): - p1 = mean, p2 = standard deviation. */ - -#define NIFTI_INTENT_NORMAL 11 - - /*! [C2, chap 30] Noncentral F statistic (3 params): - p1 = numerator DOF, p2 = denominator DOF, - p3 = numerator noncentrality parameter. */ - -#define NIFTI_INTENT_FTEST_NONC 12 - - /*! [C2, chap 29] Noncentral chi-squared statistic (2 params): - p1 = DOF, p2 = noncentrality parameter. */ - -#define NIFTI_INTENT_CHISQ_NONC 13 - - /*! [C2, chap 23] Logistic distribution (2 params): - p1 = location, p2 = scale. - Density(x) proportional to sech^2((x-p1)/(2*p2)). */ - -#define NIFTI_INTENT_LOGISTIC 14 - - /*! [C2, chap 24] Laplace distribution (2 params): - p1 = location, p2 = scale. - Density(x) proportional to exp(-abs(x-p1)/p2). */ - -#define NIFTI_INTENT_LAPLACE 15 - - /*! [C2, chap 26] Uniform distribution: p1 = lower end, p2 = upper end. */ - -#define NIFTI_INTENT_UNIFORM 16 - - /*! [C2, chap 31] Noncentral t statistic (2 params): - p1 = DOF, p2 = noncentrality parameter. */ - -#define NIFTI_INTENT_TTEST_NONC 17 - - /*! [C1, chap 21] Weibull distribution (3 params): - p1 = location, p2 = scale, p3 = power. - Density(x) proportional to - ((x-p1)/p2)^(p3-1) * exp(-((x-p1)/p2)^p3) for x > p1. */ - -#define NIFTI_INTENT_WEIBULL 18 - - /*! [C1, chap 18] Chi distribution (1 param): p1 = DOF. - Density(x) proportional to x^(p1-1) * exp(-x^2/2) for x > 0. - p1 = 1 = 'half normal' distribution - p1 = 2 = Rayleigh distribution - p1 = 3 = Maxwell-Boltzmann distribution. */ - -#define NIFTI_INTENT_CHI 19 - - /*! [C1, chap 15] Inverse Gaussian (2 params): - p1 = mu, p2 = lambda - Density(x) proportional to - exp(-p2*(x-p1)^2/(2*p1^2*x)) / x^3 for x > 0. */ - -#define NIFTI_INTENT_INVGAUSS 20 - - /*! [C2, chap 22] Extreme value type I (2 params): - p1 = location, p2 = scale - cdf(x) = exp(-exp(-(x-p1)/p2)). */ - -#define NIFTI_INTENT_EXTVAL 21 - - /*! Data is a 'p-value' (no params). */ - -#define NIFTI_INTENT_PVAL 22 - - /*! Data is ln(p-value) (no params). - To be safe, a program should compute p = exp(-abs(this_value)). - The nifti_stats.c library returns this_value - as positive, so that this_value = -log(p). */ - - -#define NIFTI_INTENT_LOGPVAL 23 - - /*! Data is log10(p-value) (no params). - To be safe, a program should compute p = pow(10.,-abs(this_value)). - The nifti_stats.c library returns this_value - as positive, so that this_value = -log10(p). */ - -#define NIFTI_INTENT_LOG10PVAL 24 - - /*! Smallest intent_code that indicates a statistic. */ - -#define NIFTI_FIRST_STATCODE 2 - - /*! Largest intent_code that indicates a statistic. */ - -#define NIFTI_LAST_STATCODE 24 - - /*---------- these values for intent_code aren't for statistics ----------*/ - - /*! To signify that the value at each voxel is an estimate - of some parameter, set intent_code = NIFTI_INTENT_ESTIMATE. - The name of the parameter may be stored in intent_name. */ - -#define NIFTI_INTENT_ESTIMATE 1001 - - /*! To signify that the value at each voxel is an index into - some set of labels, set intent_code = NIFTI_INTENT_LABEL. - The filename with the labels may stored in aux_file. */ - -#define NIFTI_INTENT_LABEL 1002 - - /*! To signify that the value at each voxel is an index into the - NeuroNames labels set, set intent_code = NIFTI_INTENT_NEURONAME. */ - -#define NIFTI_INTENT_NEURONAME 1003 - - /*! To store an M x N matrix at each voxel: - - dataset must have a 5th dimension (dim[0]=5 and dim[5]>1) - - intent_code must be NIFTI_INTENT_GENMATRIX - - dim[5] must be M*N - - intent_p1 must be M (in float format) - - intent_p2 must be N (ditto) - - the matrix values A[i][[j] are stored in row-order: - - A[0][0] A[0][1] ... A[0][N-1] - - A[1][0] A[1][1] ... A[1][N-1] - - etc., until - - A[M-1][0] A[M-1][1] ... A[M-1][N-1] */ - -#define NIFTI_INTENT_GENMATRIX 1004 - - /*! To store an NxN symmetric matrix at each voxel: - - dataset must have a 5th dimension - - intent_code must be NIFTI_INTENT_SYMMATRIX - - dim[5] must be N*(N+1)/2 - - intent_p1 must be N (in float format) - - the matrix values A[i][[j] are stored in row-order: - - A[0][0] - - A[1][0] A[1][1] - - A[2][0] A[2][1] A[2][2] - - etc.: row-by-row */ - -#define NIFTI_INTENT_SYMMATRIX 1005 - - /*! To signify that the vector value at each voxel is to be taken - as a displacement field or vector: - - dataset must have a 5th dimension - - intent_code must be NIFTI_INTENT_DISPVECT - - dim[5] must be the dimensionality of the displacement - vector (e.g., 3 for spatial displacement, 2 for in-plane) */ - -#define NIFTI_INTENT_DISPVECT 1006 /* specifically for displacements */ -#define NIFTI_INTENT_VECTOR 1007 /* for any other type of vector */ - - /*! To signify that the vector value at each voxel is really a - spatial coordinate (e.g., the vertices or nodes of a surface mesh): - - dataset must have a 5th dimension - - intent_code must be NIFTI_INTENT_POINTSET - - dim[0] = 5 - - dim[1] = number of points - - dim[2] = dim[3] = dim[4] = 1 - - dim[5] must be the dimensionality of space (e.g., 3 => 3D space). - - intent_name may describe the object these points come from - (e.g., "pial", "gray/white" , "EEG", "MEG"). */ - -#define NIFTI_INTENT_POINTSET 1008 - - /*! To signify that the vector value at each voxel is really a triple - of indexes (e.g., forming a triangle) from a pointset dataset: - - dataset must have a 5th dimension - - intent_code must be NIFTI_INTENT_TRIANGLE - - dim[0] = 5 - - dim[1] = number of triangles - - dim[2] = dim[3] = dim[4] = 1 - - dim[5] = 3 - - datatype should be an integer type (preferably DT_INT32) - - the data values are indexes (0,1,...) into a pointset dataset. */ - -#define NIFTI_INTENT_TRIANGLE 1009 - - /*! To signify that the vector value at each voxel is a quaternion: - - dataset must have a 5th dimension - - intent_code must be NIFTI_INTENT_QUATERNION - - dim[0] = 5 - - dim[5] = 4 - - datatype should be a floating point type */ - -#define NIFTI_INTENT_QUATERNION 1010 - - /*! Dimensionless value - no params - although, as in _ESTIMATE - the name of the parameter may be stored in intent_name. */ - -#define NIFTI_INTENT_DIMLESS 1011 - - /*---------- these values apply to GIFTI datasets ----------*/ - - /*! To signify that the value at each location is from a time series. */ - -#define NIFTI_INTENT_TIME_SERIES 2001 - - /*! To signify that the value at each location is a node index, from - a complete surface dataset. */ - -#define NIFTI_INTENT_NODE_INDEX 2002 - - /*! To signify that the vector value at each location is an RGB triplet, - of whatever type. - - dataset must have a 5th dimension - - dim[0] = 5 - - dim[1] = number of nodes - - dim[2] = dim[3] = dim[4] = 1 - - dim[5] = 3 - */ - -#define NIFTI_INTENT_RGB_VECTOR 2003 - - /*! To signify that the vector value at each location is a 4 valued RGBA - vector, of whatever type. - - dataset must have a 5th dimension - - dim[0] = 5 - - dim[1] = number of nodes - - dim[2] = dim[3] = dim[4] = 1 - - dim[5] = 4 - */ - -#define NIFTI_INTENT_RGBA_VECTOR 2004 - - /*! To signify that the value at each location is a shape value, such - as the curvature. */ - -#define NIFTI_INTENT_SHAPE 2005 - - /*! The following intent codes have been used by FSL FNIRT for - displacement/coefficient files. - These codes are included to prevent clashes in community-created - extensions to NIfTI. Encoding and decoding behavior for these - intents is not specified by the standard, and support is OPTIONAL - for conforming implementations. - */ - -#define NIFTI_INTENT_FSL_FNIRT_DISPLACEMENT_FIELD 2006 -#define NIFTI_INTENT_FSL_CUBIC_SPLINE_COEFFICIENTS 2007 -#define NIFTI_INTENT_FSL_DCT_COEFFICIENTS 2008 -#define NIFTI_INTENT_FSL_QUADRATIC_SPLINE_COEFFICIENTS 2009 - - /*! The following intent codes have been used by FSL TOPUP for - displacement/coefficient files. - These codes are included to prevent clashes in community-created - extensions to NIfTI. Encoding and decoding behavior for these - intents is not specified by the standard, and support is OPTIONAL - for conforming implementations. - */ - -#define NIFTI_INTENT_FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS 2016 -#define NIFTI_INTENT_FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS 2017 -#define NIFTI_INTENT_FSL_TOPUP_FIELD 2018 - -/* @} */ - -/*---------------------------------------------------------------------------*/ -/* 3D IMAGE (VOLUME) ORIENTATION AND LOCATION IN SPACE: - --------------------------------------------------- - There are 3 different methods by which continuous coordinates can - attached to voxels. The discussion below emphasizes 3D volumes, and - the continuous coordinates are referred to as (x,y,z). The voxel - index coordinates (i.e., the array indexes) are referred to as (i,j,k), - with valid ranges: - i = 0 .. dim[1]-1 - j = 0 .. dim[2]-1 (if dim[0] >= 2) - k = 0 .. dim[3]-1 (if dim[0] >= 3) - The (x,y,z) coordinates refer to the CENTER of a voxel. In methods - 2 and 3, the (x,y,z) axes refer to a subject-based coordinate system, - with - +x = Right +y = Anterior +z = Superior. - This is a right-handed coordinate system. However, the exact direction - these axes point with respect to the subject depends on qform_code - (Method 2) and sform_code (Method 3). - N.B.: The i index varies most rapidly, j index next, k index slowest. - Thus, voxel (i,j,k) is stored starting at location - (i + j*dim[1] + k*dim[1]*dim[2]) * (bitpix/8) - into the dataset array. - N.B.: The ANALYZE 7.5 coordinate system is - +x = Left +y = Anterior +z = Superior - which is a left-handed coordinate system. This backwardness is - too difficult to tolerate, so this NIFTI-1 standard specifies the - coordinate order which is most common in functional neuroimaging. - N.B.: The 3 methods below all give the locations of the voxel centers - in the (x,y,z) coordinate system. In many cases, programs will wish - to display image data on some other grid. In such a case, the program - will need to convert its desired (x,y,z) values into (i,j,k) values - in order to extract (or interpolate) the image data. This operation - would be done with the inverse transformation to those described below. - N.B.: Method 2 uses a factor 'qfac' which is either -1 or 1; qfac is - stored in the otherwise unused pixdim[0]. If pixdim[0]=0.0 (which - should not occur), we take qfac=1. Of course, pixdim[0] is only used - when reading a NIFTI-1 header, not when reading an ANALYZE 7.5 header. - N.B.: The units of (x,y,z) can be specified using the xyzt_units field. - METHOD 1 (the "old" way, used only when qform_code = 0): - ------------------------------------------------------- - The coordinate mapping from (i,j,k) to (x,y,z) is the ANALYZE - 7.5 way. This is a simple scaling relationship: - x = pixdim[1] * i - y = pixdim[2] * j - z = pixdim[3] * k - No particular spatial orientation is attached to these (x,y,z) - coordinates. (NIFTI-1 does not have the ANALYZE 7.5 orient field, - which is not general and is often not set properly.) This method - is not recommended, and is present mainly for compatibility with - ANALYZE 7.5 files. - METHOD 2 (used when qform_code > 0, which should be the "normal" case): - --------------------------------------------------------------------- - The (x,y,z) coordinates are given by the pixdim[] scales, a rotation - matrix, and a shift. This method is intended to represent - "scanner-anatomical" coordinates, which are often embedded in the - image header (e.g., DICOM fields (0020,0032), (0020,0037), (0028,0030), - and (0018,0050)), and represent the nominal orientation and location of - the data. This method can also be used to represent "aligned" - coordinates, which would typically result from some post-acquisition - alignment of the volume to a standard orientation (e.g., the same - subject on another day, or a rigid rotation to true anatomical - orientation from the tilted position of the subject in the scanner). - The formula for (x,y,z) in terms of header parameters and (i,j,k) is: - [ x ] [ R11 R12 R13 ] [ pixdim[1] * i ] [ qoffset_x ] - [ y ] = [ R21 R22 R23 ] [ pixdim[2] * j ] + [ qoffset_y ] - [ z ] [ R31 R32 R33 ] [ qfac * pixdim[3] * k ] [ qoffset_z ] - The qoffset_* shifts are in the NIFTI-1 header. Note that the center - of the (i,j,k)=(0,0,0) voxel (first value in the dataset array) is - just (x,y,z)=(qoffset_x,qoffset_y,qoffset_z). - The rotation matrix R is calculated from the quatern_* parameters. - This calculation is described below. - The scaling factor qfac is either 1 or -1. The rotation matrix R - defined by the quaternion parameters is "proper" (has determinant 1). - This may not fit the needs of the data; for example, if the image - grid is - i increases from Left-to-Right - j increases from Anterior-to-Posterior - k increases from Inferior-to-Superior - Then (i,j,k) is a left-handed triple. In this example, if qfac=1, - the R matrix would have to be - [ 1 0 0 ] - [ 0 -1 0 ] which is "improper" (determinant = -1). - [ 0 0 1 ] - If we set qfac=-1, then the R matrix would be - [ 1 0 0 ] - [ 0 -1 0 ] which is proper. - [ 0 0 -1 ] - This R matrix is represented by quaternion [a,b,c,d] = [0,1,0,0] - (which encodes a 180 degree rotation about the x-axis). - METHOD 3 (used when sform_code > 0): - ----------------------------------- - The (x,y,z) coordinates are given by a general affine transformation - of the (i,j,k) indexes: - x = srow_x[0] * i + srow_x[1] * j + srow_x[2] * k + srow_x[3] - y = srow_y[0] * i + srow_y[1] * j + srow_y[2] * k + srow_y[3] - z = srow_z[0] * i + srow_z[1] * j + srow_z[2] * k + srow_z[3] - The srow_* vectors are in the NIFTI_1 header. Note that no use is - made of pixdim[] in this method. - WHY 3 METHODS? - -------------- - Method 1 is provided only for backwards compatibility. The intention - is that Method 2 (qform_code > 0) represents the nominal voxel locations - as reported by the scanner, or as rotated to some fiducial orientation and - location. Method 3, if present (sform_code > 0), is to be used to give - the location of the voxels in some standard space. The sform_code - indicates which standard space is present. Both methods 2 and 3 can be - present, and be useful in different contexts (method 2 for displaying the - data on its original grid; method 3 for displaying it on a standard grid). - In this scheme, a dataset would originally be set up so that the - Method 2 coordinates represent what the scanner reported. Later, - a registration to some standard space can be computed and inserted - in the header. Image display software can use either transform, - depending on its purposes and needs. - In Method 2, the origin of coordinates would generally be whatever - the scanner origin is; for example, in MRI, (0,0,0) is the center - of the gradient coil. - In Method 3, the origin of coordinates would depend on the value - of sform_code; for example, for the Talairach coordinate system, - (0,0,0) corresponds to the Anterior Commissure. - QUATERNION REPRESENTATION OF ROTATION MATRIX (METHOD 2) - ------------------------------------------------------- - The orientation of the (x,y,z) axes relative to the (i,j,k) axes - in 3D space is specified using a unit quaternion [a,b,c,d], where - a*a+b*b+c*c+d*d=1. The (b,c,d) values are all that is needed, since - we require that a = sqrt(1.0-(b*b+c*c+d*d)) be nonnegative. The (b,c,d) - values are stored in the (quatern_b,quatern_c,quatern_d) fields. - The quaternion representation is chosen for its compactness in - representing rotations. The (proper) 3x3 rotation matrix that - corresponds to [a,b,c,d] is - [ a*a+b*b-c*c-d*d 2*b*c-2*a*d 2*b*d+2*a*c ] - R = [ 2*b*c+2*a*d a*a+c*c-b*b-d*d 2*c*d-2*a*b ] - [ 2*b*d-2*a*c 2*c*d+2*a*b a*a+d*d-c*c-b*b ] - [ R11 R12 R13 ] - = [ R21 R22 R23 ] - [ R31 R32 R33 ] - If (p,q,r) is a unit 3-vector, then rotation of angle h about that - direction is represented by the quaternion - [a,b,c,d] = [cos(h/2), p*sin(h/2), q*sin(h/2), r*sin(h/2)]. - Requiring a >= 0 is equivalent to requiring -Pi <= h <= Pi. (Note that - [-a,-b,-c,-d] represents the same rotation as [a,b,c,d]; there are 2 - quaternions that can be used to represent a given rotation matrix R.) - To rotate a 3-vector (x,y,z) using quaternions, we compute the - quaternion product - [0,x',y',z'] = [a,b,c,d] * [0,x,y,z] * [a,-b,-c,-d] - which is equivalent to the matrix-vector multiply - [ x' ] [ x ] - [ y' ] = R [ y ] (equivalence depends on a*a+b*b+c*c+d*d=1) - [ z' ] [ z ] - Multiplication of 2 quaternions is defined by the following: - [a,b,c,d] = a*1 + b*I + c*J + d*K - where - I*I = J*J = K*K = -1 (I,J,K are square roots of -1) - I*J = K J*K = I K*I = J - J*I = -K K*J = -I I*K = -J (not commutative!) - For example - [a,b,0,0] * [0,0,0,1] = [0,0,-b,a] - since this expands to - (a+b*I)*(K) = (a*K+b*I*K) = (a*K-b*J). - The above formula shows how to go from quaternion (b,c,d) to - rotation matrix and direction cosines. Conversely, given R, - we can compute the fields for the NIFTI-1 header by - a = 0.5 * sqrt(1+R11+R22+R33) (not stored) - b = 0.25 * (R32-R23) / a => quatern_b - c = 0.25 * (R13-R31) / a => quatern_c - d = 0.25 * (R21-R12) / a => quatern_d - If a=0 (a 180 degree rotation), alternative formulas are needed. - See the nifti1_io.c function mat44_to_quatern() for an implementation - of the various cases in converting R to [a,b,c,d]. - Note that R-transpose (= R-inverse) would lead to the quaternion - [a,-b,-c,-d]. - The choice to specify the qoffset_x (etc.) values in the final - coordinate system is partly to make it easy to convert DICOM images to - this format. The DICOM attribute "Image Position (Patient)" (0020,0032) - stores the (Xd,Yd,Zd) coordinates of the center of the first voxel. - Here, (Xd,Yd,Zd) refer to DICOM coordinates, and Xd=-x, Yd=-y, Zd=z, - where (x,y,z) refers to the NIFTI coordinate system discussed above. - (i.e., DICOM +Xd is Left, +Yd is Posterior, +Zd is Superior, - whereas +x is Right, +y is Anterior , +z is Superior. ) - Thus, if the (0020,0032) DICOM attribute is extracted into (px,py,pz), then - qoffset_x = -px qoffset_y = -py qoffset_z = pz - is a reasonable setting when qform_code=NIFTI_XFORM_SCANNER_ANAT. - That is, DICOM's coordinate system is 180 degrees rotated about the z-axis - from the neuroscience/NIFTI coordinate system. To transform between DICOM - and NIFTI, you just have to negate the x- and y-coordinates. - The DICOM attribute (0020,0037) "Image Orientation (Patient)" gives the - orientation of the x- and y-axes of the image data in terms of 2 3-vectors. - The first vector is a unit vector along the x-axis, and the second is - along the y-axis. If the (0020,0037) attribute is extracted into the - value (xa,xb,xc,ya,yb,yc), then the first two columns of the R matrix - would be - [ -xa -ya ] - [ -xb -yb ] - [ xc yc ] - The negations are because DICOM's x- and y-axes are reversed relative - to NIFTI's. The third column of the R matrix gives the direction of - displacement (relative to the subject) along the slice-wise direction. - This orientation is not encoded in the DICOM standard in a simple way; - DICOM is mostly concerned with 2D images. The third column of R will be - either the cross-product of the first 2 columns or its negative. It is - possible to infer the sign of the 3rd column by examining the coordinates - in DICOM attribute (0020,0032) "Image Position (Patient)" for successive - slices. However, this method occasionally fails for reasons that I - (RW Cox) do not understand. ------------------------------------------------------------------------------*/ - - /* [qs]form_code value: */ /* x,y,z coordinate system refers to: */ - /*-----------------------*/ /*---------------------------------------*/ - -/*! \defgroup NIFTI1_XFORM_CODES - \brief nifti1 xform codes to describe the "standard" coordinate system - @{ - */ - /*! Arbitrary coordinates (Method 1). */ - -#define NIFTI_XFORM_UNKNOWN 0 - - /*! Scanner-based anatomical coordinates */ - -#define NIFTI_XFORM_SCANNER_ANAT 1 - - /*! Coordinates aligned to another file's, - or to anatomical "truth". */ - -#define NIFTI_XFORM_ALIGNED_ANAT 2 - - /*! Coordinates aligned to Talairach- - Tournoux Atlas; (0,0,0)=AC, etc. */ - -#define NIFTI_XFORM_TALAIRACH 3 - - /*! MNI 152 normalized coordinates. */ - -#define NIFTI_XFORM_MNI_152 4 - - /*! Normalized coordinates (for - any general standard template - space). Added March 8, 2019. */ - -#define NIFTI_XFORM_TEMPLATE_OTHER 5 - -/* @} */ - -/*---------------------------------------------------------------------------*/ -/* UNITS OF SPATIAL AND TEMPORAL DIMENSIONS: - ---------------------------------------- - The codes below can be used in xyzt_units to indicate the units of pixdim. - As noted earlier, dimensions 1,2,3 are for x,y,z; dimension 4 is for - time (t). - - If dim[4]=1 or dim[0] < 4, there is no time axis. - - A single time series (no space) would be specified with - - dim[0] = 4 (for scalar data) or dim[0] = 5 (for vector data) - - dim[1] = dim[2] = dim[3] = 1 - - dim[4] = number of time points - - pixdim[4] = time step - - xyzt_units indicates units of pixdim[4] - - dim[5] = number of values stored at each time point - Bits 0..2 of xyzt_units specify the units of pixdim[1..3] - (e.g., spatial units are values 1..7). - Bits 3..5 of xyzt_units specify the units of pixdim[4] - (e.g., temporal units are multiples of 8). - This compression of 2 distinct concepts into 1 byte is due to the - limited space available in the 348 byte ANALYZE 7.5 header. The - macros XYZT_TO_SPACE and XYZT_TO_TIME can be used to mask off the - undesired bits from the xyzt_units fields, leaving "pure" space - and time codes. Inversely, the macro SPACE_TIME_TO_XYZT can be - used to assemble a space code (0,1,2,...,7) with a time code - (0,8,16,32,...,56) into the combined value for xyzt_units. - Note that codes are provided to indicate the "time" axis units are - actually frequency in Hertz (_HZ), in part-per-million (_PPM) - or in radians-per-second (_RADS). - The toffset field can be used to indicate a nonzero start point for - the time axis. That is, time point #m is at t=toffset+m*pixdim[4] - for m=0..dim[4]-1. ------------------------------------------------------------------------------*/ - -/*! \defgroup NIFTI1_UNITS - \brief nifti1 units codes to describe the unit of measurement for - each dimension of the dataset - @{ - */ - /*! NIFTI code for unspecified units. */ -#define NIFTI_UNITS_UNKNOWN 0 - - /** Space codes are multiples of 1. **/ - /*! NIFTI code for meters. */ -#define NIFTI_UNITS_METER 1 - /*! NIFTI code for millimeters. */ -#define NIFTI_UNITS_MM 2 - /*! NIFTI code for micrometers. */ -#define NIFTI_UNITS_MICRON 3 - - /** Time codes are multiples of 8. **/ - /*! NIFTI code for seconds. */ -#define NIFTI_UNITS_SEC 8 - /*! NIFTI code for milliseconds. */ -#define NIFTI_UNITS_MSEC 16 - /*! NIFTI code for microseconds. */ -#define NIFTI_UNITS_USEC 24 - - /*** These units are for spectral data: ***/ - /*! NIFTI code for Hertz. */ -#define NIFTI_UNITS_HZ 32 - /*! NIFTI code for ppm. */ -#define NIFTI_UNITS_PPM 40 - /*! NIFTI code for radians per second. */ -#define NIFTI_UNITS_RADS 48 -/* @} */ - -#undef XYZT_TO_SPACE -#undef XYZT_TO_TIME -#define XYZT_TO_SPACE(xyzt) ( (xyzt) & 0x07 ) -#define XYZT_TO_TIME(xyzt) ( (xyzt) & 0x38 ) - -#undef SPACE_TIME_TO_XYZT -#define SPACE_TIME_TO_XYZT(ss,tt) ( (((char)(ss)) & 0x07) \ - | (((char)(tt)) & 0x38) ) - -/*---------------------------------------------------------------------------*/ -/* MRI-SPECIFIC SPATIAL AND TEMPORAL INFORMATION: - --------------------------------------------- - A few fields are provided to store some extra information - that is sometimes important when storing the image data - from an FMRI time series experiment. (After processing such - data into statistical images, these fields are not likely - to be useful.) - { freq_dim } = These fields encode which spatial dimension (1,2, or 3) - { phase_dim } = corresponds to which acquisition dimension for MRI data. - { slice_dim } = - Examples: - Rectangular scan multi-slice EPI: - freq_dim = 1 phase_dim = 2 slice_dim = 3 (or some permutation) - Spiral scan multi-slice EPI: - freq_dim = phase_dim = 0 slice_dim = 3 - since the concepts of frequency- and phase-encoding directions - don't apply to spiral scan - slice_duration = If this is positive, AND if slice_dim is nonzero, - indicates the amount of time used to acquire 1 slice. - slice_duration*dim[slice_dim] can be less than pixdim[4] - with a clustered acquisition method, for example. - slice_code = If this is nonzero, AND if slice_dim is nonzero, AND - if slice_duration is positive, indicates the timing - pattern of the slice acquisition. The following codes - are defined: - NIFTI_SLICE_SEQ_INC == sequential increasing - NIFTI_SLICE_SEQ_DEC == sequential decreasing - NIFTI_SLICE_ALT_INC == alternating increasing - NIFTI_SLICE_ALT_DEC == alternating decreasing - NIFTI_SLICE_ALT_INC2 == alternating increasing #2 - NIFTI_SLICE_ALT_DEC2 == alternating decreasing #2 - { slice_start } = Indicates the start and end of the slice acquisition - { slice_end } = pattern, when slice_code is nonzero. These values - are present to allow for the possible addition of - "padded" slices at either end of the volume, which - don't fit into the slice timing pattern. If there - are no padding slices, then slice_start=0 and - slice_end=dim[slice_dim]-1 are the correct values. - For these values to be meaningful, slice_start must - be non-negative and slice_end must be greater than - slice_start. Otherwise, they should be ignored. - The following table indicates the slice timing pattern, relative to - time=0 for the first slice acquired, for some sample cases. Here, - dim[slice_dim]=7 (there are 7 slices, labeled 0..6), slice_duration=0.1, - and slice_start=1, slice_end=5 (1 padded slice on each end). - slice - index SEQ_INC SEQ_DEC ALT_INC ALT_DEC ALT_INC2 ALT_DEC2 - 6 : n/a n/a n/a n/a n/a n/a n/a = not applicable - 5 : 0.4 0.0 0.2 0.0 0.4 0.2 (slice time offset - 4 : 0.3 0.1 0.4 0.3 0.1 0.0 doesn't apply to - 3 : 0.2 0.2 0.1 0.1 0.3 0.3 slices outside - 2 : 0.1 0.3 0.3 0.4 0.0 0.1 the range - 1 : 0.0 0.4 0.0 0.2 0.2 0.4 slice_start .. - 0 : n/a n/a n/a n/a n/a n/a slice_end) - The SEQ slice_codes are sequential ordering (uncommon but not unknown), - either increasing in slice number or decreasing (INC or DEC), as - illustrated above. - The ALT slice codes are alternating ordering. The 'standard' way for - these to operate (without the '2' on the end) is for the slice timing - to start at the edge of the slice_start .. slice_end group (at slice_start - for INC and at slice_end for DEC). For the 'ALT_*2' slice_codes, the - slice timing instead starts at the first slice in from the edge (at - slice_start+1 for INC2 and at slice_end-1 for DEC2). This latter - acquisition scheme is found on some Siemens scanners. - The fields freq_dim, phase_dim, slice_dim are all squished into the single - byte field dim_info (2 bits each, since the values for each field are - limited to the range 0..3). This unpleasantness is due to lack of space - in the 348 byte allowance. - The macros DIM_INFO_TO_FREQ_DIM, DIM_INFO_TO_PHASE_DIM, and - DIM_INFO_TO_SLICE_DIM can be used to extract these values from the - dim_info byte. - The macro FPS_INTO_DIM_INFO can be used to put these 3 values - into the dim_info byte. ------------------------------------------------------------------------------*/ - -#undef DIM_INFO_TO_FREQ_DIM -#undef DIM_INFO_TO_PHASE_DIM -#undef DIM_INFO_TO_SLICE_DIM - -#define DIM_INFO_TO_FREQ_DIM(di) ( ((di) ) & 0x03 ) -#define DIM_INFO_TO_PHASE_DIM(di) ( ((di) >> 2) & 0x03 ) -#define DIM_INFO_TO_SLICE_DIM(di) ( ((di) >> 4) & 0x03 ) - -#undef FPS_INTO_DIM_INFO -#define FPS_INTO_DIM_INFO(fd,pd,sd) ( ( ( ((char)(fd)) & 0x03) ) | \ - ( ( ((char)(pd)) & 0x03) << 2 ) | \ - ( ( ((char)(sd)) & 0x03) << 4 ) ) - -/*! \defgroup NIFTI1_SLICE_ORDER - \brief nifti1 slice order codes, describing the acquisition order - of the slices - @{ - */ -#define NIFTI_SLICE_UNKNOWN 0 -#define NIFTI_SLICE_SEQ_INC 1 -#define NIFTI_SLICE_SEQ_DEC 2 -#define NIFTI_SLICE_ALT_INC 3 -#define NIFTI_SLICE_ALT_DEC 4 -#define NIFTI_SLICE_ALT_INC2 5 /* 05 May 2005: RWCox */ -#define NIFTI_SLICE_ALT_DEC2 6 /* 05 May 2005: RWCox */ -/* @} */ - -/*---------------------------------------------------------------------------*/ -/* UNUSED FIELDS: - ------------- - Some of the ANALYZE 7.5 fields marked as ++UNUSED++ may need to be set - to particular values for compatibility with other programs. The issue - of interoperability of ANALYZE 7.5 files is a murky one -- not all - programs require exactly the same set of fields. (Unobscuring this - murkiness is a principal motivation behind NIFTI-1.) - Some of the fields that may need to be set for other (non-NIFTI aware) - software to be happy are: - extents dbh.h says this should be 16384 - regular dbh.h says this should be the character 'r' - glmin, } dbh.h says these values should be the min and max voxel - glmax } values for the entire dataset - It is best to initialize ALL fields in the NIFTI-1 header to 0 - (e.g., with calloc()), then fill in what is needed. ------------------------------------------------------------------------------*/ - -/*---------------------------------------------------------------------------*/ -/* MISCELLANEOUS C MACROS ------------------------------------------------------------------------------*/ - -/*.................*/ -/*! Given a nifti_1_header struct, check if it has a good magic number. - Returns NIFTI version number (1..9) if magic is good, 0 if it is not. */ - -#define NIFTI_VERSION(h) \ - ( ( (h).magic[0]=='n' && (h).magic[3]=='\0' && \ - ( (h).magic[1]=='i' || (h).magic[1]=='+' ) && \ - ( (h).magic[2]>='1' && (h).magic[2]<='9' ) ) \ - ? (h).magic[2]-'0' : 0 ) - -/*.................*/ -/*! Check if a nifti_1_header struct says if the data is stored in the - same file or in a separate file. Returns 1 if the data is in the same - file as the header, 0 if it is not. */ - -#define NIFTI_ONEFILE(h) ( (h).magic[1] == '+' ) - -/*.................*/ -/*! Check if a nifti_1_header struct needs to be byte swapped. - Returns 1 if it needs to be swapped, 0 if it does not. */ - -#define NIFTI_NEEDS_SWAP(h) ( (h).dim[0] < 0 || (h).dim[0] > 7 ) - -/*.................*/ -/*! Check if a nifti_1_header struct contains a 5th (vector) dimension. - Returns size of 5th dimension if > 1, returns 0 otherwise. */ - -#define NIFTI_5TH_DIM(h) ( ((h).dim[0]>4 && (h).dim[5]>1) ? (h).dim[5] : 0 ) - -/*****************************************************************************/ - -/*=================*/ -#ifdef __cplusplus -} -#endif -/*=================*/ - -#endif /* NIFTI1_HEADER */ \ No newline at end of file diff --git a/pkg/nifti/nii_io/nifti2.h b/pkg/nifti/nii_io/nifti2.h deleted file mode 100644 index 341ff4cb..00000000 --- a/pkg/nifti/nii_io/nifti2.h +++ /dev/null @@ -1,113 +0,0 @@ -/** \file nifti2.h - \brief Header structure for NIFTI-2 format. - */ - -#ifndef NIFTI2_HEADER -#define NIFTI2_HEADER - -/*---------------------------------------------------------------------------*/ -/* Changes to the header from NIFTI-1 to NIFTI-2 are intended to allow for - larger and more accurate fields. The changes are as follows: - - short dim[8] -> int64_t dim[8] - - float intent_p1,2,3 -> double intent_p1,2,3 (3 fields) - - float pixdim[8] -> double pixdim[8] - - float vox_offset -> int64_t vox_offset - - float scl_slope -> double scl_slope - - float scl_inter -> double scl_inter - - float cal_max -> double cal_max - - float cal_min -> double cal_min - - float slice_duration -> double slice_duration - - float toffset -> double toffset - - short slice_start -> int64_t slice_start - - short slice_end -> int64_t slice_end - - char slice_code -> int32_t slice_code - - char xyzt_units -> int32_t xyzt_units - - short intent_code -> int32_t intent_code - - short qform_code -> int32_t qform_code - - short sform_code -> int32_t sform_code - - float quatern_b,c,d -> double quatern_b,c,d (3 fields) - - float srow_x,y,z[4] -> double srow_x,y,z[4] (3 fields) - - char magic[4] -> char magic[8] - - char unused_str[15] -> padding added at the end of the header - - previously unused fields have been removed: - data_type, db_name, extents, session_error, regular, glmax, glmin - - the field order has been changed, notably with magic after sizeof_hdr - 2 Jan, 2014 [rickr] ------------------------------------------------------------------------------*/ - -#include - -/*=================*/ -#ifdef __cplusplus -extern "C" { -#endif -/*=================*/ - -/*! \struct nifti_2_header - \brief Data structure defining the fields in the nifti2 header. - This binary header should be found at the beginning of a valid - NIFTI-2 header file. - */ - -/* hopefully cross-platform solution to byte padding added by some compilers */ -#pragma pack(push) -#pragma pack(1) - - /*****************************/ /***********************/ /************/ -struct nifti_2_header { /* NIFTI-2 usage */ /* NIFTI-1 usage */ /* offset */ - /*****************************/ /***********************/ /************/ - int32_t sizeof_hdr; /*!< MUST be 540 */ /* MUST be 348 */ /* 0 */ - char magic[8]; /*!< MUST be valid signature */ /* char magic[4] */ /* 4 */ - int16_t datatype; /*!< Defines data type! */ /* short datatype */ /* 12 */ - int16_t bitpix; /*!< Number bits/voxel */ /* short bitpix */ /* 14 */ - int64_t dim[8]; /*!< Data array dimensions */ /* short dim[8] */ /* 16 */ - double intent_p1; /*!< 1st intent parameter */ /* float intent_p1 */ /* 80 */ - double intent_p2; /*!< 2nd intent parameter */ /* float intent_p2 */ /* 88 */ - double intent_p3; /*!< 3rd intent parameter */ /* float intent_p3 */ /* 96 */ - double pixdim[8]; /*!< Grid spacings */ /* float pixdim[8] */ /* 104 */ - int64_t vox_offset; /*!< Offset into .nii file */ /* float vox_offset */ /* 168 */ - double scl_slope; /*!< Data scaling: slope */ /* float scl_slope */ /* 176 */ - double scl_inter; /*!< Data scaling: offset */ /* float scl_inter */ /* 184 */ - double cal_max; /*!< Max display intensity */ /* float cal_max */ /* 192 */ - double cal_min; /*!< Min display intensity */ /* float cal_min */ /* 200 */ - double slice_duration; /*!< Time for 1 slice */ /* float slice_duration*/ /* 208 */ - double toffset; /*!< Time axis shift */ /* float toffset */ /* 216 */ - int64_t slice_start; /*!< First slice index */ /* short slice_start */ /* 224 */ - int64_t slice_end; /*!< Last slice index */ /* short slice_end */ /* 232 */ - char descrip[80]; /*!< any text you like */ /* char descrip[80] */ /* 240 */ - char aux_file[24]; /*!< auxiliary filename */ /* char aux_file[24] */ /* 320 */ - int32_t qform_code; /*!< NIFTI_XFORM_* code */ /* short qform_code */ /* 344 */ - int32_t sform_code; /*!< NIFTI_XFORM_* code */ /* short sform_code */ /* 348 */ - double quatern_b; /*!< Quaternion b param */ /* float quatern_b */ /* 352 */ - double quatern_c; /*!< Quaternion c param */ /* float quatern_c */ /* 360 */ - double quatern_d; /*!< Quaternion d param */ /* float quatern_d */ /* 368 */ - double qoffset_x; /*!< Quaternion x shift */ /* float qoffset_x */ /* 376 */ - double qoffset_y; /*!< Quaternion y shift */ /* float qoffset_y */ /* 384 */ - double qoffset_z; /*!< Quaternion z shift */ /* float qoffset_z */ /* 392 */ - double srow_x[4]; /*!< 1st row affine transform*/ /* float srow_x[4] */ /* 400 */ - double srow_y[4]; /*!< 2nd row affine transform*/ /* float srow_y[4] */ /* 432 */ - double srow_z[4]; /*!< 3rd row affine transform*/ /* float srow_z[4] */ /* 464 */ - int32_t slice_code; /*!< Slice timing order */ /* char slice_code */ /* 496 */ - int32_t xyzt_units; /*!< Units of pixdim[1..4] */ /* char xyzt_units */ /* 500 */ - int32_t intent_code; /*!< NIFTI_INTENT_* code */ /* short intent_code */ /* 504 */ - char intent_name[16];/*!< name or meaning of data */ /* char intent_name[16]*/ /* 508 */ - char dim_info; /*!< MRI slice ordering */ /* char dim_info */ /* 524 */ - char unused_str[15]; /*!< unused, filled with \0 */ /* 525 */ -}; /****** total bytes: 540 */ -typedef struct nifti_2_header nifti_2_header; - -/* restore packing behavior */ -#pragma pack(pop) - -/* base swap test on the suggested version check, rather than dim[0] - swap4(348)==1543569408, swap4(540)==469893120 */ -#define NIFTI2_NEEDS_SWAP(h) \ - ((h).sizeof_hdr == 1543569408 || (h).sizeof_hdr == 469893120) - -/*=================*/ -#ifdef __cplusplus -} -#endif -/*=================*/ - -#endif /* NIFTI2_HEADER */ \ No newline at end of file diff --git a/private.go b/private.go new file mode 100644 index 00000000..76cd0f61 --- /dev/null +++ b/private.go @@ -0,0 +1,200 @@ +package go2com + +import ( + "bytes" + "encoding/binary" + "github.com/okieraised/go2com/pkg/dicom/tag" + "github.com/okieraised/go2com/pkg/dicom/uid" + "github.com/okieraised/go2com/pkg/dicom/vr" + "io" +) + +func (r *dcmReader) isTrackingImplicit() bool { + return r.keepTrackImplicit +} + +func (r *dcmReader) setOverallImplicit(isImplicit bool) { + r.keepTrackImplicit = isImplicit +} + +func (r *dcmReader) readUInt8() (uint8, error) { + var res uint8 + + err := binary.Read(r, r.binaryOrder, &res) + + return res, err +} + +func (r *dcmReader) readUInt16() (uint16, error) { + var res uint16 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readUInt32() (uint32, error) { + var res uint32 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readUInt64() (uint64, error) { + var res uint64 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readInt8() (int8, error) { + var res int8 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readInt16() (int16, error) { + var res int16 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readInt32() (int32, error) { + var res int32 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readInt64() (int64, error) { + var res int64 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readFloat32() (float32, error) { + var res float32 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) readFloat64() (float64, error) { + var res float64 + err := binary.Read(r, r.binaryOrder, &res) + return res, err +} + +func (r *dcmReader) peek(n int) ([]byte, error) { + return r.reader.Peek(n) +} + +func (r *dcmReader) discard(n int) (int, error) { + return r.reader.Discard(n) +} + +func (r *dcmReader) skip(n int64) error { + _, err := io.CopyN(io.Discard, r, n) + return err +} + +func (r *dcmReader) parse() error { + _ = r.SetFileSize(r.fileSize) + err := r.IsValidDICOM() + if err != nil { + return err + } + _ = r.skip(132) + err = r.parseMetadata() + if err != nil { + return err + } + + // IMPORTANT: Additional check is needed here since there are few instances where the DICOM + // meta header is registered as Explicit Little-Endian, but Implicit Little-Endian is used in the body + err = r.verifyImplicity() + if err != nil { + return nil + } + + if r.skipDataset { + return nil + } + err = r.parseDataset() + if err != nil { + return err + } + return nil +} + +// parseMetadata parses the file meta information according to +// https://dicom.nema.org/dicom/2013/output/chtml/part10/chapter_7.html +// the File Meta Information shall be encoded using the Explicit VR Little Endian Transfer Syntax +// (UID=1.2.840.10008.1.2.1) +func (r *dcmReader) parseMetadata() error { + var metadata []*Element + var transferSyntaxUID string + + for { + // No longer relied on the MetaInformationGroupLength tag to determine the length of the meta header. + // We check if the group tag is 0x0002 before proceeding to read the element. If the group tag is not 0x0002, + // then break the loop + n, err := r.peek(2) + if err != nil { + return err + } + if bytes.Compare(n, []byte{0x2, 0x0}) != 0 { + break + } + + res, err := ReadElement(r, r.IsImplicit(), r.ByteOrder()) + if err != nil { + return err + } + metadata = append(metadata, res) + if res.Tag == tag.TransferSyntaxUID { + transferSyntaxUID = (res.Value.RawValue).(string) + } + } + r.metadata = Dataset{Elements: metadata} + + // Set transfer syntax here for the dataset parser + binOrder, isImplicit, err := uid.ParseTransferSyntaxUID(transferSyntaxUID) + if err != nil { + return err + } + r.SetTransferSyntax(binOrder, isImplicit) + r.setOverallImplicit(isImplicit) + + return nil +} + +func (r *dcmReader) verifyImplicity() error { + // Need to check if the implicit matches between header and body + n, err := r.peek(6) + if err != nil { + return err + } + if !vr.VRMapper[string(n[4:6])] && !r.IsImplicit() { + r.SetTransferSyntax(r.binaryOrder, true) + } + if vr.VRMapper[string(n[4:6])] && r.IsImplicit() { + r.SetTransferSyntax(r.binaryOrder, false) + } + return nil +} + +// parseDataset parses the file dataset after the file meta header +func (r *dcmReader) parseDataset() error { + var data []*Element + for { + res, err := ReadElement(r, r.IsImplicit(), r.ByteOrder()) + if err != nil { + if err == io.EOF { + break + } else { + return err + } + } + data = append(data, res) + //fmt.Println(res) + } + dicomDataset := Dataset{Elements: data} + //r.dataset.Elements = append(r.dataset.Elements, dicomDataset.Elements...) + r.dataset = dicomDataset + return nil +} diff --git a/tag_mapper.go b/tag_mapper.go new file mode 100644 index 00000000..b3ed493e --- /dev/null +++ b/tag_mapper.go @@ -0,0 +1,94 @@ +package go2com + +import ( + "fmt" + "github.com/okieraised/go2com/internal/utils" + "github.com/okieraised/go2com/pkg/dicom/tag" + "github.com/okieraised/go2com/pkg/dicom/vr" + "strings" +) + +type MappedTag map[string]tag.TagBrowser + +// GetElementByTagString returns the element value of the input tag +// Tag should be in (gggg,eeee) or ggggeeee format +func (m MappedTag) GetElementByTagString(tagStr string) (interface{}, error) { + tagStr = utils.FormatTag(tagStr) + + result, ok := m[tagStr] + if !ok { + return nil, fmt.Errorf("tag not found: %s", tagStr) + } + return result.Value, nil +} + +// mapElement returns a map[string]interface{} with key as tag and value as the tag values +func (m MappedTag) mapElement(elem *Element) { + tagStr := elem.Tag.StringWithoutParentheses() + vrStr := elem.ValueRepresentationStr + var value interface{} + + // If VR is SQ then we do type assertion to []*element.Element. If the length of sequence is 0, then do nothing. + // Else, loop through each element in the sequence and extract the info + if vrStr == "SQ" { + subVL := make([]interface{}, 0) + vlArr, ok := (elem.Value.RawValue).([]*Element) + if ok { + if len(vlArr) == 0 { + return + } + groupTag := vlArr[0].Tag.StringWithoutParentheses() + subElemGrp := make(MappedTag) + for index, subVl := range vlArr { + subVRStr := subVl.ValueRepresentationStr + if subVRStr == "OB" || subVRStr == "OW" || subVRStr == "UN" || strings.ToLower(subVRStr) == "ox" { + continue + } + subTag := subVl.Tag.StringWithoutParentheses() + if subTag == groupTag && index > 0 { + subVL = append(subVL, subElemGrp) + subElemGrp = MappedTag{} + } + subElemGrp.mapElement(subVl) + if index == len(vlArr)-1 { + subVL = append(subVL, subElemGrp) + } + } + } + value = subVL + } else { + if elem.ValueRepresentationStr == vr.PersonName { + value = utils.AppendToSlice(map[string]interface{}{ + "Alphabetic": elem.Value.RawValue, + }) + } else { + value = utils.AppendToSlice(elem.Value.RawValue) + } + + } + m[tagStr] = tag.TagBrowser{ + VR: vrStr, + Value: value, + } +} + +func createOrthancURI(ds Dataset) MappedTag { + res := make(MappedTag, 3) + prefix := "http://127.0.0.1:8042" + uids, err := ds.RetrieveFileUID() + if err != nil { + return res + } + for _, elem := range ds.Elements { + var bulkURI string + tagStr := elem.Tag.StringWithoutParentheses() + if elem.Tag == tag.RedPaletteColorLookupTableData || elem.Tag == tag.BluePaletteColorLookupTableData || elem.Tag == tag.GreenPaletteColorLookupTableData { + bulkURI = fmt.Sprintf("%s/dicom-web/studies/%s/series/%s/instances/%s/bulk/%s", prefix, uids.StudyInstanceUID, uids.SeriesInstanceUID, uids.SOPInstanceUID, tagStr) + res[tagStr] = tag.TagBrowser{ + VR: "US", + BulkDataURI: bulkURI, + } + } + } + return res +}