forked from henghuang/nifti
/
nifti.go
484 lines (429 loc) · 16 KB
/
nifti.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
package nifti
import (
"encoding/binary"
"fmt"
"io"
"io/ioutil"
"math"
"os"
pkgPath "path/filepath"
"strings"
gzip "github.com/klauspost/pgzip"
)
type Nifti1Header struct {
SizeofHdr int32 /*!< MUST be 348 */ /* int32 sizeof_hdr; */
DataType [10]byte /*!< ++UNUSED++ */ /* char data_type[10]; */
DbName [18]byte /*!< ++UNUSED++ */ /* char db_name[18]; */
Extents int32 /*!< ++UNUSED++ */ /* int32 extents; */
SessionError int16 /*!< ++UNUSED++ */ /* short session_error; */
Regular byte /*!< ++UNUSED++ */ /* char regular; */
DimInfo byte /*!< MRI slice ordering. */ /* char hkey_un0; */
/*--- was image_dimension substruct ---*/
Dim [8]int16 /*!< Data array dimensions.*/ /* short dim[8]; */
IntentP1 float32 /*!< 1st intent parameter. */ /* short unused8; */
/* short unused9; */
IntentP2 float32 /*!< 2nd intent parameter. */ /* short unused10; */
/* short unused11; */
IntentP3 float32 /*!< 3rd intent parameter. */ /* short unused12; */
/* short unused13; */
IntentCode int16 /*!< NIFTI_INTENT_* code. */ /* short unused14; */
Datatype int16 /*!< Defines data type! */ /* short datatype; */
Bitpix int16 /*!< Number bits/voxel. */ /* short bitpix; */
SliceStart int16 /*!< First slice index. */ /* short dim_un0; */
Pixdim [8]float32 /*!< Grid spacings. */ /* float32 pixdim[8]; */
VoxOffset float32 /*!< Offset into .nii file */ /* float32 vox_offset; */
SclSlope float32 /*!< Data scaling: slope. */ /* float32 funused1; */
SclInter float32 /*!< Data scaling: offset. */ /* float32 funused2; */
SliceEnd int16 /*!< Last slice index. */ /* float32 funused3; */
SliceCode byte /*!< Slice timing order. */
XyztUnits byte /*!< Units of pixdim[1..4] */
CalMax float32 /*!< Max display intensity */ /* float32 cal_max; */
CalMin float32 /*!< Min display intensity */ /* float32 cal_min; */
SliceDuration float32 /*!< Time for 1 slice. */ /* float32 compressed; */
Toffset float32 /*!< Time axis shift. */ /* float32 verified; */
Glmax int32 /*!< ++UNUSED++ */ /* int32 glmax; */
Glmin int32 /*!< ++UNUSED++ */ /* int32 glmin; */
/*--- was data_history substruct ---*/
Descrip [80]byte /*!< any text you like. */ /* char descrip[80]; */
AuxFile [24]byte /*!< auxiliary filename. */ /* char aux_file[24]; */
QformCode int16 /*!< NIFTI_XFORM_* code. */ /*-- all ANALYZE 7.5 ---*/
SformCode int16 /*!< NIFTI_XFORM_* code. */ /* fields below here */
/* are replaced */
QuaternB float32 /*!< Quaternion b param. */
QuaternC float32 /*!< Quaternion c param. */
QuaternD float32 /*!< Quaternion d param. */
QoffsetX float32 /*!< Quaternion x shift. */
QoffsetY float32 /*!< Quaternion y shift. */
QoffsetZ float32 /*!< Quaternion z shift. */
SrowX [4]float32 /*!< 1st row affine transform. */
SrowY [4]float32 /*!< 2nd row affine transform. */
SrowZ [4]float32 /*!< 3rd row affine transform. */
IntentName [16]byte /*!< 'name' or meaning of data. */
Magic [4]byte /*!< MUST be "ni1\0" or "n+1\0". */
} /**** 348 bytes total ****/
func (header *Nifti1Header) LoadHeader(filepath string) {
reader, err := gzipOpen(filepath)
defer reader.Close()
if err != nil {
fmt.Println(err)
return
}
err = binary.Read(reader, binary.LittleEndian, header)
if err != nil {
fmt.Println(err)
return
}
// fmt.Printf("%+v\n", header)
}
type Nifti1Image struct { /*!< Image storage struct **/
ndim uint32 /*!< last dimension greater than 1 (1..7) */
nx uint32 /*!< dimensions of grid array */
ny uint32 /*!< dimensions of grid array */
nz uint32 /*!< dimensions of grid array */
nt uint32 /*!< dimensions of grid array */
nu uint32 /*!< dimensions of grid array */
nv uint32 /*!< dimensions of grid array */
nw uint32 /*!< dimensions of grid array */
dim [8]uint32 /*!< dim[0]=ndim, dim[1]=nx, etc. */
nvox uint64 /*!< number of voxels = nx*ny*nz*...*nw */
nbyper uint32 /*!< bytes per voxel, matches datatype */
datatype int32 /*!< type of data in voxels: DT_* code */
dx float32 /*!< grid spacings */
dy float32 /*!< grid spacings */
dz float32 /*!< grid spacings */
dt float32 /*!< grid spacings */
du float32 /*!< grid spacings */
dv float32 /*!< grid spacings */
dw float32 /*!< grid spacings */
pixdim [8]float32 /*!< pixdim[1]=dx, etc. */
sclSlope float32 /*!< scaling parameter - slope */
sclInter float32 /*!< scaling parameter - intercept */
calMin float32 /*!< calibration parameter, minimum */
calMax float32 /*!< calibration parameter, maximum */
qformCode int32 /*!< codes for (x,y,z) space meaning */
sformCode int32 /*!< codes for (x,y,z) space meaning */
freqDim int32 /*!< indexes (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 int32 /*!< index for start of slices */
sliceEnd int32 /*!< index for end of slices */
sliceDuration float32 /*!< time between individual slices */
/*! quaternion transform parameters
[when writing a dataset, these are used for qform, NOT qto_xyz] */
quaternB, quaternC, quaternD,
qoffsetX, qoffsetY, qoffsetZ,
qfac float32
qtoXyz [4][4]float32 /*!< qform: transform (i,j,k) to (x,y,z) */
qtoIjk [4][4]float32 /*!< qform: transform (x,y,z) to (i,j,k) */
stoXyz [4][4]float32 /*!< sform: transform (i,j,k) to (x,y,z) */
stoIjk [4][4]float32 /*!< sform: transform (x,y,z) to (i,j,k) */
toffset float32 /*!< 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 (1 file),
2==NIFTI-1 (2 files),
3==NIFTI-ASCII (1 file) */
intentCode int32 /*!< statistic type (or something) */
intentP1 float32 /*!< intent parameters */
intentP2 float32 /*!< intent parameters */
intentP3 float32 /*!< intent parameters */
intentName [16]byte /*!< optional description of intent data */
descrip [80]byte /*!< optional text to describe dataset */
auxFile [24]byte /*!< auxiliary filename */
fname string /*!< header filename (.hdr or .nii) */
iname string /*!< image filename (.img or .nii) */
inameOffset int32 /*!< offset into iname where data starts */
swapsize int32 /*!< swap unit in image data (might be 0) */
byteorder int32 /*!< byte order on disk (MSB_ or LSB_FIRST) */
data []byte /*!< pointer to data: nbyper*nvox bytes */
volumeN int //defined by me, volume vox num
byte2floatF func([]byte) float32 //defined by me, byte2floatF
float2byteF func([]byte, float32)
header Nifti1Header //defined by me, it might be a good idea store the img header in the image structure
numExt int32 /*!< number of extensions in ext_list */
// nifti1_extension *ext_list /*!< array of extension structs (with data) */
// analyze_75_orient_code analyze75_orient /*!< for old analyze files, orient */
}
//filepath ,read data?
//header.VoxOffset, header.Bitpix,header.Dim are very important
func (img *Nifti1Image) LoadImage(filepath string, rdata bool) {
var header Nifti1Header
if filepath == "" {
header = img.header
} else {
header.LoadHeader(filepath)
img.header = header
}
// set dimensions of data array
img.ndim, img.dim[0] = uint32(header.Dim[0]), uint32(header.Dim[0])
img.nx, img.dim[1] = uint32(header.Dim[1]), uint32(header.Dim[1])
img.ny, img.dim[2] = uint32(header.Dim[2]), uint32(header.Dim[2])
img.nz, img.dim[3] = uint32(header.Dim[3]), uint32(header.Dim[3])
img.nt, img.dim[4] = uint32(header.Dim[4]), uint32(header.Dim[4])
img.nu, img.dim[5] = uint32(header.Dim[5]), uint32(header.Dim[5])
img.nv, img.dim[6] = uint32(header.Dim[6]), uint32(header.Dim[6])
img.nw, img.dim[7] = uint32(header.Dim[7]), uint32(header.Dim[7])
img.nvox = 1
for i := int16(1); i <= header.Dim[0]; i++ {
img.nvox *= uint64(header.Dim[i])
}
img.volumeN = int(img.dim[1] * img.dim[2] * img.dim[3])
if header.Bitpix == 0 {
fmt.Println("empty header.Bitpix")
fmt.Println(header)
return
}
//init byte2float32 function
img.nbyper = uint32(header.Bitpix) / 8
// fmt.Println(header)
//setting function to convert float2byte or byte2float
if img.nbyper == 1 {
img.byte2floatF = func(b []byte) float32 {
v := uint8(b[0])
return float32(v)
}
img.float2byteF = func(buff []byte, x float32) {
buff[0] = uint8(x)
}
} else if img.nbyper == 2 {
img.byte2floatF = func(b []byte) float32 {
v := binary.LittleEndian.Uint16(b)
return float32(v)
}
img.float2byteF = func(buff []byte, x float32) {
binary.LittleEndian.PutUint16(buff, uint16(x))
}
} else if img.nbyper == 4 {
img.byte2floatF = func(b []byte) float32 {
v := binary.LittleEndian.Uint32(b)
return math.Float32frombits(v)
}
img.float2byteF = func(buff []byte, x float32) {
v := math.Float32bits(x)
binary.LittleEndian.PutUint32(buff, v)
}
} else if img.nbyper == 8 {
img.byte2floatF = func(b []byte) float32 {
v := binary.LittleEndian.Uint64(b)
return float32(math.Float64frombits(v))
}
img.float2byteF = func(buff []byte, x float32) {
v := math.Float64bits(float64(x))
binary.LittleEndian.PutUint64(buff, v)
}
} else {
fmt.Println("input nbyper:", img.nbyper)
panic("(img *Nifti1Image) byte2float, only support 8 16 32 and 64 bit")
}
// fmt.Println(img.nbyper)
// set the grid spacings
img.dx, img.pixdim[1] = header.Pixdim[1], header.Pixdim[1]
img.dy, img.pixdim[2] = header.Pixdim[2], header.Pixdim[2]
img.dz, img.pixdim[3] = header.Pixdim[3], header.Pixdim[3]
img.dt, img.pixdim[4] = header.Pixdim[4], header.Pixdim[4]
img.du, img.pixdim[5] = header.Pixdim[5], header.Pixdim[5]
img.dv, img.pixdim[6] = header.Pixdim[6], header.Pixdim[6]
img.dw, img.pixdim[7] = header.Pixdim[7], header.Pixdim[7]
if rdata {
var data []byte
var err error
if pkgPath.Ext(filepath) == ".gz" {
reader, err := gzipOpen(filepath)
defer reader.Close()
if err != nil {
fmt.Println("open data error")
return
}
data, err = ioutil.ReadAll(reader)
if err != nil {
fmt.Println("read data error")
return
}
} else {
data, err = ioutil.ReadFile(filepath)
if err != nil {
fmt.Println("open data error")
return
}
}
img.data = data[uint(header.VoxOffset):len(data)]
// fmt.Println(len(img.data), len(data), uint(header.VoxOffset))
}
}
// x,y,z,t,start at zero
func (img *Nifti1Image) GetAt(x, y, z, t uint32) float32 {
tIndex := img.nx * img.ny * img.nz * t
zIndex := img.nx * img.ny * z
yIndex := img.nx * y
xIndex := x
index := uint64(tIndex + zIndex + yIndex + xIndex)
return img.byte2float(img.data[index*uint64(img.nbyper) : (index+1)*uint64(img.nbyper)]) //shift byte
}
func (img *Nifti1Image) GetHeader() Nifti1Header {
return img.header
}
func (img *Nifti1Image) SetNewHeader(newHeader Nifti1Header) {
img.header = newHeader
return
}
func (img *Nifti1Image) SetHeaderDim(dimX, dimY, dimZ, dimT int) {
img.header.SizeofHdr = 348
img.header.Regular = 114
if dimT > 1 {
img.header.Dim[0] = 4
} else {
img.header.Dim[0] = 3
}
img.header.Dim[1] = int16(dimX)
img.header.Dim[2] = int16(dimY)
img.header.Dim[3] = int16(dimZ)
img.header.Bitpix = 32
img.header.Pixdim = [8]float32{-1, 2, 2, 2, 1, 1, 1, 1}
img.header.VoxOffset = 348
img.header.SclSlope = 1
img.header.XyztUnits = 10
img.header.CalMax = 8000
img.header.CalMin = 3000
img.header.QformCode = 4
img.header.SformCode = 4
img.header.QuaternC = 1
img.header.QoffsetX = 90
img.header.QoffsetY = -126
img.header.QoffsetZ = -72
img.header.SrowX = [4]float32{-2, 0, 0, 90}
img.header.SrowY = [4]float32{0, 2, 0, -126}
img.header.SrowZ = [4]float32{0, 0, 2, -72}
img.header.Magic = [4]byte{110, 43, 49, 0}
return
}
func (img *Nifti1Image) SetHeaderDim2(dimX, dimY, dimZ, dimT int) {
if dimT > 1 {
img.header.Dim[0] = 4
} else {
img.header.Dim[0] = 3
}
img.header.Dim[1] = int16(dimX)
img.header.Dim[2] = int16(dimY)
img.header.Dim[3] = int16(dimZ)
img.header.Dim[4] = int16(dimT)
return
}
func (img *Nifti1Image) SetAt(x, y, z, t uint32, elem float32) {
tIndex := img.nx * img.ny * img.nz * t
zIndex := img.nx * img.ny * z
yIndex := img.nx * y
xIndex := x
index := uint64(tIndex + zIndex + yIndex + xIndex)
copy(img.data[index*uint64(img.nbyper):(index+1)*uint64(img.nbyper)], img.float2byte(elem))
}
func (img *Nifti1Image) GetTimeSeries(x, y, z uint32) []float32 {
var timeSeries []float32
timeNum := uint32(len(img.data) / int(img.volumeN*int(img.nbyper)))
var i uint32
for i = 0; i < timeNum; i++ {
timeSeries = append(timeSeries, img.GetAt(x, y, z, i))
}
return timeSeries
}
func (img *Nifti1Image) GetSlice(z, t uint32) [][]float32 {
sliceX := img.nx
sliceY := img.ny
slice := make([][]float32, sliceX)
for i := range slice {
slice[i] = make([]float32, sliceY)
}
var xi uint32
var yi uint32
for xi = 0; xi < sliceX; xi++ {
for yi = 0; yi < sliceY; yi++ {
slice[xi][yi] = img.GetAt(xi, yi, z, t)
// fmt.Println(xi, yi, z, t)
}
}
return slice
}
//Return [x,y,z,t]
func (img *Nifti1Image) GetDims() [4]int {
var dim [4]int
for i, v := range img.dim[1:5] {
dim[i] = int(v)
}
return dim
}
//convert byte to float32,init in LoadImage
func (img *Nifti1Image) byte2float(data []byte) float32 {
v := img.byte2floatF(data)
return v
}
//convert float32 to byte,init in LoadImage
func (img *Nifti1Image) float2byte(data float32) []byte {
buff := make([]byte, img.nbyper)
img.float2byteF(buff, data)
return buff
}
//write nii
func (img *Nifti1Image) Save(filename string) {
padding := int(img.header.VoxOffset) - 348
file, err := os.Create(filename + ".gz")
defer file.Close()
if err != nil {
panic(err)
}
gzipWriter := gzip.NewWriter(file)
defer gzipWriter.Close()
binary.Write(gzipWriter, binary.LittleEndian, img.header)
if padding > 0 {
binary.Write(gzipWriter, binary.LittleEndian, make([]byte, padding)) //padding header
}
binary.Write(gzipWriter, binary.LittleEndian, img.data)
}
func gzipOpen(filepath string) (io.ReadCloser, error) {
f, err := os.Open(filepath)
if err != nil {
return nil, err
}
filepathS := strings.Split(filepath, ".")
if filepathS[len(filepathS)-1] == "gz" {
gzipReader, err := gzip.NewReader(f)
if err != nil {
return nil, err
}
return gzipReader, nil
}
return f, nil
}
//create new image
func NewImg(dimX, dimY, dimZ, dimT int) *Nifti1Image {
var img Nifti1Image
img.header.SizeofHdr = 348
img.header.Regular = 114
if dimT > 1 {
img.header.Dim[0] = 4
} else {
img.header.Dim[0] = 3
}
img.header.Dim[1] = int16(dimX)
img.header.Dim[2] = int16(dimY)
img.header.Dim[3] = int16(dimZ)
img.header.Bitpix = 32
img.header.Pixdim = [8]float32{-1, 2, 2, 2, 1, 1, 1, 1}
img.header.VoxOffset = 348
img.header.SclSlope = 1
img.header.XyztUnits = 10
img.header.CalMax = 8000
img.header.CalMin = 3000
img.header.QformCode = 4
img.header.SformCode = 4
img.header.QuaternC = 1
img.header.QoffsetX = 90
img.header.QoffsetY = -126
img.header.QoffsetZ = -72
img.header.SrowX = [4]float32{-2, 0, 0, 90}
img.header.SrowY = [4]float32{0, 2, 0, -126}
img.header.SrowZ = [4]float32{0, 0, 2, -72}
img.header.Magic = [4]byte{110, 43, 49, 0}
img.LoadImage("", false)
img.data = make([]byte, dimX*dimY*dimZ*dimT*int(img.nbyper))
return &img
}